How to Automate Typesetting of Matrices within Function?

I am toying around with Typst’s scripting functionality, and I was trying to see if it was robust enough to automate the calculation of structure constants for SU(3) in its 3-dimensional representation. However, I am facing a couple issues with rendering and with algebraic numbers. The code I currently have is below.

#let Tr = math.upright("Tr")
#let complex(a, b) = (a, b)
#let complex-add(z1, z2) = (z1.at(0) + z2.at(0), z1.at(1) + z2.at(1))
#let complex-sub(z1, z2) = (z1.at(0) - z2.at(0), z1.at(1) - z2.at(1))
#let complex-mul(z1, z2) = (
  z1.at(0) * z2.at(0) - z1.at(1) * z2.at(1),
  z1.at(0) * z2.at(1) + z1.at(1) * z2.at(0)
)
#let complex-neg(z) = (-z.at(0), -z.at(1))
#let complex-scale(z, s) = (z.at(0) * s, z.at(1) * s)

#let format-complex(z) = {
  let real = z.at(0)
  let imag = z.at(1)
  
  if imag == 0 {
    return str(real)
  } else if real == 0 {
    if imag == 1 { return "i" }
    else if imag == -1 { return "-i" }
    else { return str(imag) + "i" }
  } else {
    let sign = if imag > 0 { "+" } else { "-" }
    let abs-imag = calc.abs(imag)
    if abs-imag == 1 {
      return str(real) + sign + "i"
    } else {
      return str(real) + sign + str(abs-imag) + "i"
    }
  }
}

#let mat-add(A, B) = {
  let rows = A.len()
  let cols = A.at(0).len()
  
  let result = ()
  for i in range(rows) {
    let row = ()
    for j in range(cols) {
      row.push(complex-add(A.at(i).at(j), B.at(i).at(j)))
    }
    result.push(row)
  }
  
  return result
}

#let mat-sub(A, B) = {
  let rows = A.len()
  let cols = A.at(0).len()
  
  let result = ()
  for i in range(rows) {
    let row = ()
    for j in range(cols) {
      row.push(complex-sub(A.at(i).at(j), B.at(i).at(j)))
    }
    result.push(row)
  }
  
  return result
}

#let mat-mul(A, B) = {
  let rows_A = A.len()
  let cols_A = A.at(0).len()
  let cols_B = B.at(0).len()
  
  let result = ()
  for i in range(rows_A) {
    let row = ()
    for j in range(cols_B) {
      let sum = complex(0, 0)
      for k in range(cols_A) {
        sum = complex-add(sum, complex-mul(A.at(i).at(k), B.at(k).at(j)))
      }
      row.push(sum)
    }
    result.push(row)
  }
  
  return result
}

#let mat-commutator(A, B) = {
  let AB = mat-mul(A, B)
  let BA = mat-mul(B, A)
  return mat-sub(AB, BA)
}

#let mat-trace(A) = {
  let sum = complex(0, 0)
  let n = A.len()
  
  for i in range(n) {
    sum = complex-add(sum, A.at(i).at(i))
  }
  
  return sum
}

// Define Gell-Mann matrices
#let t1 = (
  (complex(0, 0), complex(0.5, 0), complex(0, 0)),
  (complex(0.5, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0))
)

#let t2 = (
  (complex(0, 0), complex(0, -0.5), complex(0, 0)),
  (complex(0, 0.5), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0))
)

#let t3 = (
  (complex(0.5, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(-0.5, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0))
)

#let t4 = (
  (complex(0, 0), complex(0, 0), complex(0.5, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(0.5, 0), complex(0, 0), complex(0, 0))
)

#let t5 = (
  (complex(0, 0), complex(0, 0), complex(0, -0.5)),
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0.5), complex(0, 0), complex(0, 0))
)

#let t6 = (
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0.5, 0)),
  (complex(0, 0), complex(0.5, 0), complex(0, 0))
)

#let t7 = (
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, -0.5)),
  (complex(0, 0), complex(0, 0.5), complex(0, 0))
)

#let factor = 1/(2*calc.sqrt(3))
#let t8 = (
  (complex(factor, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(factor, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(-2*factor, 0))
)

#let t = (t1, t2, t3, t4, t5, t6, t7, t8)

#let calculate-f(a, b, c, Cr) = {
  let a_idx = a - 1
  let b_idx = b - 1
  let c_idx = c - 1
  
  let comm = mat-commutator(t.at(a_idx), t.at(b_idx))
  let prod = mat-mul(comm, t.at(c_idx))
  let tr = mat-trace(prod)
  
  // f^(abc) = -i/Cr * Tr([t_a, t_b] t_c)
  let result = complex-mul(complex(0, -1/Cr), tr)
  
  return (
    commutator: comm,
    product: prod,
    trace: tr,
    result: result
  )
}

#let Cr = 1 // Normalization constant

#let display-calculation(a, b, c) = {
  let res = calculate-f(a, b, c, Cr)
  
  [
    == Structure Constant $f^(#a #b #c)$
    
    *Step 1:* Calculate the commutator $[t_(#a), t_(#b)]$
    
    $ [t_(#a), t_(#b)] = t_(#a) t_(#b) - t_(#b) t_(#a) = mat(
        #format-complex(res.commutator.at(0).at(0)), #format-complex(res.commutator.at(0).at(1)), #format-complex(res.commutator.at(0).at(2));
        #format-complex(res.commutator.at(1).at(0)), #format-complex(res.commutator.at(1).at(1)), #format-complex(res.commutator.at(1).at(2));
        #format-complex(res.commutator.at(2).at(0)), #format-complex(res.commutator.at(2).at(1)), #format-complex(res.commutator.at(2).at(2))
      ) $
    
    *Step 2:* Multiply by $t_(#c)$ to get $[t_(#a), t_(#b)]t_(#c)$
    
    $ [t_(#a), t_(#b)]t_(#c) = mat(
        #format-complex(res.product.at(0).at(0)), #format-complex(res.product.at(0).at(1)), #format-complex(res.product.at(0).at(2));
        #format-complex(res.product.at(1).at(0)), #format-complex(res.product.at(1).at(1)), #format-complex(res.product.at(1).at(2));
        #format-complex(res.product.at(2).at(0)), #format-complex(res.product.at(2).at(1)), #format-complex(res.product.at(2).at(2))
      ) $
    
    *Step 3:* Calculate the trace
    
    $ Tr([t_(#a), t_(#b)]t_(#c)) = #format-complex(res.trace) $
    
    *Step 4:* Multiply by $-i/C_r$ to get the structure constant
    
    $ f^(#a #b #c) = frac(-i, C_r) Tr([t_(#a), t_(#b)]t_(#c)) = #format-complex(res.result) $
  ]
}

// Sample calculations for a few structure constants
#display-calculation(1, 2, 3)

The two primary issues I am having is that, when trying to render:

mat(
        #format-complex(res.product.at(0).at(0)), #format-complex(res.product.at(0).at(1)), #format-complex(res.product.at(0).at(2));
        #format-complex(res.product.at(1).at(0)), #format-complex(res.product.at(1).at(1)), #format-complex(res.product.at(1).at(2));
        #format-complex(res.product.at(2).at(0)), #format-complex(res.product.at(2).at(1)), #format-complex(res.product.at(2).at(2))
      ) 

and similarly for the other matrix, it refuses to render as a matrix. What I get is, e.g. (0.5i 0 00 −0.5i 00 0 0), rather than the block form it should have. I cannot figure out what’s going on. I assume the input matrix elements are somehow messing with the rendering.

The second issue I am facing is that I need to use 1/(2*calc.sqrt(3)) for one of the Gell-Mann matrices, but that converts it to a float, making it impossible to render it as 1/sqrt(3) if the output still has the sqrt(3).

Does anyone understand why the matrices aren’t rendering, and if there’s a way to get the sqrt() to still render, even if actual automated calculations require calc.sqrt(3)?

Inserting a space before the semicolon fixes the rendering of the matrix, although I am not sure why the matrix function is behaving like that in the first place. You do not need the space when you create a matrix with (string) elements manually.

Note that the code you shared did not compile because the function Tr() was not defined. For future questions please only share as much code as necessary to reproduce your issue. In this case the function format-complex() and the matrix would have been sufficient.

I don’t really understand your second issue with calc.sqrt(3) and 1/sqrt(3). Could you please create a small (code) example to show what you are trying to achieve?

1 Like

Thanks for pointing out the spacing. Yeah, that is exceptionally odd, but if it works, it works.

And thanks for pointing out the Tr() function issue. I defined that manually via:

#let Tr = math.upright("Tr")

I’ll add that to the original post.

As for what I am trying to get out from the sqrt(3) situation, as an example, if one tries to evaluate:

#display-calculation(1, 2, 8)

one will get on step 2, a diagonal matrix with 0.14433756729740646i as the diagonal element, rather than something like i/(4 sqrt(3)). In short, I am trying to keep the actual sqrt() symbols around, but I don’t know if there’s a way to do that when the actual calculations have to be done with floats. Technically, I would like to get this working with stuff like 0.5 as well, but the sqrt(3) is the really egregious issue.

If you don’t actually need the prefactor 1/sqrt(3) for the computation, could you not just store this separately from the matrix? You would then have an array f that holds the prefactors of the matrices t1 to t8. The matrix computations would then just need integer math and you would include the prefactor(s) again after the computations.

2 Likes

That’s actually a good point. This wouldn’t work for a more general calculation if I needed to be able to square sqrt(3) (assuming the calculation kept the problematic overflow like 3.0000000001), but since f^(a b c) is antisymmetric, that isn’t an issue. I’ll make use of that. Thanks.

You could generalize this further by storing the prefactor in a dictionary where the “value” and the “exponent” are kept separately, e.g.

#let f8 = (value: 3, exponent: decimal("0.5"))

If you need to square this, you would then add the exponents to get decimal("1.0") in which case you would just omit the square root. I know that this is not the most convenient thing to work with, but you can get vary far with integer + decimal math. Even multiplying and reducing fractions is possible without ever touching floating point math.

1 Like

I really like that dictionary idea, and I might have to play around to see if I can get that to work. And while I can see how one might do multiplication of fractions, it’s less obvious how one would do reduction, unless the idea is to do 0 actual multiplication when combining fractions, keeping all the factors separated, and then trying to track down terms with the same base.

The reason is that ; ends a code block, and the function call is done with one (# starts it): #format-complex(...);

with the space, the code block already ends at ) and thus the ; is part of the matrix syntax

1 Like

In case anyone is curious, this is what I eventually got working:

#let Tr = math.upright("Tr")

// ============ SYMBOLIC FRACTION IMPLEMENTATION ============

// Basic symbolic fraction representation
#let sym-frac(num, denom) = {
  // Ensure we always have integer representations
  if type(num) != int { num = calc.round(num) }
  if type(denom) != int { denom = calc.round(denom) }
  
  // Handle special cases
  if denom == 0 { panic("Division by zero") }
  if num == 0 { return (type: "sym-frac", num: 0, denom: 1, sign: 1) }
  
  // Normalize negative signs - always store them separately
  let sign = 1
  if (num < 0 and denom > 0) or (num > 0 and denom < 0) {
    sign = -1 // Negative fraction
  }
  num = calc.abs(num)
  denom = calc.abs(denom)
  
  // Create fraction with normalized sign
  return (
    type: "sym-frac",
    num: num,     // Always positive
    denom: denom, // Always positive
    sign: sign    // -1 if fraction is negative, 1 if positive
  )
}

// Fraction arithmetic operations
#let frac-add(f1, f2) = {
  if f1.type != "sym-frac" or f2.type != "sym-frac" {
    panic("Expected symbolic fractions")
  }
  
  // a/b + c/d = (ad + bc)/bd
  // With signs: (sign_a * a)/b + (sign_c * c)/d
  let sign_a_val = if f1.sign < 0 { -f1.num } else { f1.num }
  let sign_c_val = if f2.sign < 0 { -f2.num } else { f2.num }
  
  let num = sign_a_val * f2.denom + sign_c_val * f1.denom
  let denom = f1.denom * f2.denom
  
  return sym-frac(num, denom)
}

#let frac-sub(f1, f2) = {
  if f1.type != "sym-frac" or f2.type != "sym-frac" {
    panic("Expected symbolic fractions")
  }
  
  // a/b - c/d = (ad - bc)/bd
  // With signs: (sign_a * a)/b - (sign_c * c)/d
  let sign_a_val = if f1.sign < 0 { -f1.num } else { f1.num }
  let sign_c_val = if f2.sign < 0 { -f2.num } else { f2.num }
  
  let num = sign_a_val * f2.denom - sign_c_val * f1.denom
  let denom = f1.denom * f2.denom
  
  return sym-frac(num, denom)
}

#let frac-mul(f1, f2) = {
  if f1.type != "sym-frac" or f2.type != "sym-frac" {
    panic("Expected symbolic fractions")
  }
  
  // (a/b) * (c/d) = (ac)/(bd)
  // Sign: sign_a * sign_c
  let num = f1.num * f2.num
  let denom = f1.denom * f2.denom
  let sign = f1.sign * f2.sign
  
  return (
    type: "sym-frac",
    num: num,
    denom: denom,
    sign: sign
  )
}

#let frac-div(f1, f2) = {
  if f1.type != "sym-frac" or f2.type != "sym-frac" {
    panic("Expected symbolic fractions")
  }
  
  // (a/b) / (c/d) = (ad)/(bc)
  // Sign: sign_a * sign_c
  let num = f1.num * f2.denom
  let denom = f1.denom * f2.num
  let sign = f1.sign * f2.sign
  
  return (
    type: "sym-frac",
    num: num,
    denom: denom,
    sign: sign
  )
}

#let frac-neg(f) = {
  if f.type != "sym-frac" {
    panic("Expected symbolic fraction")
  }
  
  // Negate the sign
  return (
    type: "sym-frac",
    num: f.num,
    denom: f.denom,
    sign: -f.sign
  )
}

#let frac-is-zero(f) = {
  if f.type != "sym-frac" {
    panic("Expected symbolic fraction")
  }
  
  return f.num == 0
}

#let frac-is-one(f) = {
  if f.type != "sym-frac" {
    panic("Expected symbolic fraction")
  }
  
  return f.num == f.denom and f.sign > 0
}

// Format a symbolic fraction
#let format-frac(f) = {
  if f.type != "sym-frac" {
    panic("Expected symbolic fraction")
  }
  
  if f.denom == 1 {
    if f.sign < 0 {
      return [-#f.num]
    } else {
      return [#f.num]
    }
  } else {
    let fraction = math.frac([#f.num], [#f.denom])
    if f.sign < 0 {
      return [-#fraction]
    } else {
      return fraction
    }
  }
}

// ============ ALGEBRAIC NUMBER SYSTEM ============

// Clone a dictionary (since Typst dictionaries don't have a clone method)
#let clone-dict(dict) = {
  let result = (:)
  for (key, value) in dict {
    result.insert(key, value)
  }
  return result
}

// Represent an algebraic number with symbolic fractions for coefficients
#let alg-num(
  real-coef-num: 1, real-coef-denom: 1, real-factors: (:),
  imag-coef-num: 0, imag-coef-denom: 1, imag-factors: (:)
) = {
  return (
    type: "algebraic",
    real: (
      coef: sym-frac(real-coef-num, real-coef-denom),
      factors: clone-dict(real-factors)
    ),
    imag: (
      coef: sym-frac(imag-coef-num, imag-coef-denom),
      factors: clone-dict(imag-factors)
    )
  )
}

// Convenience constructors
#let alg-int(n) = alg-num(real-coef-num: n, real-coef-denom: 1)
#let alg-i = alg-num(real-coef-num: 0, imag-coef-num: 1, imag-coef-denom: 1)

// Create a proper fraction
#let alg-frac(num, denom) = {
  return alg-num(real-coef-num: num, real-coef-denom: denom)
}

// Create a square root
#let alg-sqrt(n) = {
  let factors = (:)
  // Convert n to string if it's an integer
  let base = if type(n) == int { str(n) } else { n }
  factors.insert(base, decimal("0.5"))
  return alg-num(real-factors: factors)
}

// Create a reciprocal square root (1/√n)
#let alg-rsqrt(n) = {
  let factors = (:)
  // Convert n to string if it's an integer
  let base = if type(n) == int { str(n) } else { n }
  factors.insert(base, decimal("-0.5"))
  return alg-num(real-factors: factors)
}

// Check if a decimal is approximately an integer
#let is-integer(dec) = {
  let int_value = calc.round(dec)
  if type(dec) == decimal {
    return calc.abs(dec - decimal(str(int_value))) < decimal("0.00001")
  } else {
    return calc.abs(dec - int_value) < 0.00001
  }
}

// Simplify one part (real or imaginary) of an algebraic number
#let simplify-part(part) = {
  let coef = part.coef
  let factors = clone-dict(part.factors)
  
  // Handle zero coefficient
  if frac-is-zero(coef) {
    return (coef: sym-frac(0, 1), factors: (:))
  }
  
  // Simplify factors with zero or integer exponents
  let keys = factors.keys()
  for base in keys {
    let exp = factors.at(base)
    
    // Convert to decimal if not already
    if type(exp) != decimal {
      exp = decimal(str(exp))
      factors.at(base) = exp
    }
    
    // Remove factors with zero exponent
    if exp == decimal("0") {
      factors.remove(base)
    }
  }
  
  return (coef: coef, factors: factors)
}

// Simplify an algebraic number
#let simplify-alg(num) = {
  if num.type != "algebraic" {
    return alg-int(num)
  }
  
  // Simplify both parts
  let real-part = simplify-part(num.real)
  let imag-part = simplify-part(num.imag)
  
  // Create the simplified number
  return (
    type: "algebraic",
    real: real-part,
    imag: imag-part
  )
}

// Multiply a part (real or imaginary) by another part
#let multiply-parts(a, b) = {
  // If either coefficient is zero, result is zero
  if frac-is-zero(a.coef) or frac-is-zero(b.coef) {
    return (coef: sym-frac(0, 1), factors: (:))
  }
  
  // Multiply coefficients
  let coef = frac-mul(a.coef, b.coef)
  
  // Combine factors by adding exponents
  let result_factors = clone-dict(a.factors)
  for (base, exp) in b.factors {
    if base in result_factors {
      result_factors.at(base) = result_factors.at(base) + exp
    } else {
      result_factors.insert(base, exp)
    }
  }
  
  return (coef: coef, factors: result_factors)
}

// Multiply two algebraic numbers
#let alg-mul(a, b) = {
  if a.type != "algebraic" {
    a = alg-int(a)
  }
  if b.type != "algebraic" {
    b = alg-int(b)
  }
  
  // (a + bi)(c + di) = (ac - bd) + (ad + bc)i
  
  // Real part: ac - bd
  let real_ac = multiply-parts(a.real, b.real)
  let real_bd = multiply-parts(a.imag, b.imag)
  
  // Imaginary part: ad + bc
  let imag_ad = multiply-parts(a.real, b.imag)
  let imag_bc = multiply-parts(a.imag, b.real)
  
  // For the real part, need to subtract bd from ac
  let real_result
  if frac-is-zero(real_bd.coef) {
    // If bd is zero, real part is just ac
    real_result = real_ac
  } else if frac-is-zero(real_ac.coef) {
    // If ac is zero, real part is -bd
    real_result = (
      coef: frac-neg(real_bd.coef),
      factors: real_bd.factors
    )
  } else {
    // Both non-zero, but we can only combine if they have identical factors
    let factors_match = true
    if real_ac.factors.keys().sorted() == real_bd.factors.keys().sorted() {
      for (base, exp_ac) in real_ac.factors {
        let exp_bd = real_bd.factors.at(base, default: none)
        if exp_bd == none or exp_ac != exp_bd {
          factors_match = false
          break
        }
      }
    } else {
      factors_match = false
    }
    
    if factors_match {
      real_result = (
        coef: frac-sub(real_ac.coef, real_bd.coef),
        factors: real_ac.factors
      )
    } else {
      // If we can't combine, just use ac for now
      real_result = real_ac
    }
  }
  
  // For the imaginary part, need to add ad and bc
  let imag_result
  if frac-is-zero(imag_ad.coef) {
    // If ad is zero, imaginary part is just bc
    imag_result = imag_bc
  } else if frac-is-zero(imag_bc.coef) {
    // If bc is zero, imaginary part is just ad
    imag_result = imag_ad
  } else {
    // Both non-zero, but we can only combine if they have identical factors
    let factors_match = true
    if imag_ad.factors.keys().sorted() == imag_bc.factors.keys().sorted() {
      for (base, exp_ad) in imag_ad.factors {
        let exp_bc = imag_bc.factors.at(base, default: none)
        if exp_bc == none or exp_ad != exp_bc {
          factors_match = false
          break
        }
      }
    } else {
      factors_match = false
    }
    
    if factors_match {
      imag_result = (
        coef: frac-add(imag_ad.coef, imag_bc.coef),
        factors: imag_ad.factors
      )
    } else {
      // If we can't combine, just use ad for now
      imag_result = imag_ad
    }
  }

  return simplify-alg((
    type: "algebraic",
    real: real_result,
    imag: imag_result
  ))
}

// Calculate the complex conjugate of an algebraic number
#let alg-conjugate(a) = {
  if a.type != "algebraic" {
    a = alg-int(a)
  }
  
  // Negate the imaginary part
  return (
    type: "algebraic",
    real: a.real,
    imag: (
      coef: frac-neg(a.imag.coef),
      factors: clone-dict(a.imag.factors)
    )
  )
}

// Divide two algebraic numbers using a/(b+ci) = a(b-ci)/(b²+c²)
#let alg-div(a, b) = {
  if a.type != "algebraic" {
    a = alg-int(a)
  }
  if b.type != "algebraic" {
    b = alg-int(b)
  }
  
  // Special case: denominator is pure real
  if frac-is-zero(b.imag.coef) {
    // Create copies of a's factors
    let real_factors = clone-dict(a.real.factors)
    let imag_factors = clone-dict(a.imag.factors)
    
    // For each factor in the denominator, add it with negated exponent
    for (base, exp) in b.real.factors {
      if base in real_factors {
        real_factors.at(base) = real_factors.at(base) - exp
      } else {
        real_factors.insert(base, -exp) // Add with negative exponent
      }
      
      if base in imag_factors {
        imag_factors.at(base) = imag_factors.at(base) - exp
      } else {
        imag_factors.insert(base, -exp) // Add with negative exponent
      }
    }
    
    return (
      type: "algebraic",
      real: (
        coef: frac-div(a.real.coef, b.real.coef),
        factors: real_factors  // Now contains correct combined factors
      ),
      imag: (
        coef: frac-div(a.imag.coef, b.real.coef),
        factors: imag_factors  // Now contains correct combined factors
      )
    )
  }
  
  // For complex denominators, multiply by conjugate
  // a/(b+ci) = a(b-ci)/((b+ci)(b-ci)) = a(b-ci)/(b²+c²)
  let b_conj = alg-conjugate(b)
  
  // Multiply numerator by conjugate of denominator
  let num = alg-mul(a, b_conj)
  
  // Calculate denominator: (b+ci)(b-ci) = b² + c²
  let denom = alg-mul(b, b_conj)
  
  // Now denom should be real-only
  if not frac-is-zero(denom.imag.coef) {
    // This should not happen mathematically
    return simplify-alg(num)
  }
  
  // Divide by the real denominator
  return alg-div(num, denom)
}

// Format an algebraic number for display
#let format-alg(num) = {
  if num.type != "algebraic" {
    return str(num)
  }
  
  num = simplify-alg(num)
  
  // Special case for zero
  if frac-is-zero(num.real.coef) and frac-is-zero(num.imag.coef) {
    return "0"
  }
  
  // Function to format a single part (real or imaginary)
  let format-part(part) = {
    let coef = part.coef
    let factors = part.factors
    
    // If coefficient is zero, return empty content
    if frac-is-zero(coef) {
      return []
    }
    
    // Separate factors into positive and negative exponents
    let pos_factors = (:)
    let neg_factors = (:)
    for (base, exp) in factors {
      if exp > decimal("0") {
        pos_factors.insert(base, exp)
      } else if exp < decimal("0") {
        neg_factors.insert(base, -exp) // Store positive value in neg_factors
      }
    }
    
    // Get the sign for the entire expression
    let sign_content = []
    if coef.sign < 0 {
      sign_content = math.minus
    }
    
    // Format the numerator (coefficient and positive exponent factors)
    let num_content = []
    
    // Add coefficient absolute value if not 1 or if there are no positive factors
    if coef.num != coef.denom or pos_factors.keys().len() == 0 {
      if coef.num == coef.denom {
        // Coefficient is 1, don't show it unless there are no factors
        if pos_factors.keys().len() == 0 {
          num_content = [1]
        }
      } else {
        // Show the coefficient (always positive now)
        num_content = [#coef.num]
      }
    }
    
    // Add positive exponent factors
    let sorted_pos_factors = pos_factors.keys().sorted()
    for base in sorted_pos_factors {
      let exp = pos_factors.at(base)
      if exp == decimal("0.5") {
        num_content = [#num_content #math.sqrt(base)]
      } else if exp == decimal("1") {
        num_content = [#num_content #base]
      } else {
        num_content = [#num_content #base^(#exp)]
      }
    }
    
    // Check if we need a fraction
    // If we have no negative factors and denominator is 1, no fraction needed
    if neg_factors.keys().len() == 0 and coef.denom == 1 {
      return [#sign_content #num_content]
    }
    
    // Format the denominator (coefficient denominator and negative exponent factors)
    let denom_content = []
    
    // Add coefficient denominator if not 1
    if coef.denom != 1 {
      denom_content = [#coef.denom]
    }
    
    // Add negative exponent factors
    let sorted_neg_factors = neg_factors.keys().sorted()
    for base in sorted_neg_factors {
      let exp = neg_factors.at(base)
      if exp == decimal("0.5") {
        denom_content = [#denom_content #math.sqrt(base)]
      } else if exp == decimal("1") {
        denom_content = [#denom_content #base]
      } else {
        denom_content = [#denom_content #base^(#exp)]
      }
    }
    
    // Return the fraction with sign outside
    let fraction = math.frac([#num_content], [#denom_content])
    return [#sign_content #fraction]
  }
  
  // Format real and imaginary parts
  let real_result = format-part(num.real)
  let imag_result = format-part(num.imag)
  
  // If imaginary part has content, add 'i'
  if not frac-is-zero(num.imag.coef) {
    imag_result = [#imag_result i]
  }
  
  // Combine real and imaginary parts
  if frac-is-zero(num.real.coef) {
    return imag_result
  } else if frac-is-zero(num.imag.coef) {
    return real_result
  } else {
    // Need + sign between real and imaginary parts (but only if imag is positive)
    if num.imag.coef.sign > 0 {
      return [#real_result + #imag_result]
    } else {
      // The minus sign is already included in imag_result
      return [#real_result #imag_result]
    }
  }
}

// ============ COMPLEX NUMBER IMPLEMENTATION (FOR MATRICES) ============

// Basic complex number operations (used for matrices)
#let complex(a, b) = (a, b)
#let complex-add(z1, z2) = (z1.at(0) + z2.at(0), z1.at(1) + z2.at(1))
#let complex-sub(z1, z2) = (z1.at(0) - z2.at(0), z1.at(1) - z2.at(1))
#let complex-mul(z1, z2) = (
  z1.at(0) * z2.at(0) - z1.at(1) * z2.at(1),
  z1.at(0) * z2.at(1) + z1.at(1) * z2.at(0)
)
#let complex-scale(z, s) = (z.at(0) * s, z.at(1) * s)

// Convert complex to algebraic
#let complex-to-alg(z) = {
  let real = z.at(0)
  let imag = z.at(1)
  
  return alg-num(
    real-coef-num: real, real-coef-denom: 1,
    imag-coef-num: imag, imag-coef-denom: 1
  )
}

// Format a complex number for matrix display
#let format-complex(z) = {
  let real = z.at(0)
  let imag = z.at(1)
  
  if imag == 0 {
    return str(real)
  } else if real == 0 {
    if imag == 1 { return "i" }
    else if imag == -1 { return "-i" }
    else { return str(imag) + "i" }
  } else {
    let sign = if imag > 0 { "+" } else { "-" }
    let abs-imag = calc.abs(imag)
    if abs-imag == 1 {
      return str(real) + sign + "i"
    } else {
      return str(real) + sign + str(abs-imag) + "i"
    }
  }
}

// ============ MATRIX OPERATIONS ============

#let mat-add(A, B) = {
  let rows = A.len()
  let cols = A.at(0).len()
  
  let result = ()
  for i in range(rows) {
    let row = ()
    for j in range(cols) {
      row.push(complex-add(A.at(i).at(j), B.at(i).at(j)))
    }
    result.push(row)
  }
  
  return result
}

#let mat-sub(A, B) = {
  let rows = A.len()
  let cols = A.at(0).len()
  
  let result = ()
  for i in range(rows) {
    let row = ()
    for j in range(cols) {
      row.push(complex-sub(A.at(i).at(j), B.at(i).at(j)))
    }
    result.push(row)
  }
  
  return result
}

#let mat-mul(A, B) = {
  let rows_A = A.len()
  let cols_A = A.at(0).len()
  let cols_B = B.at(0).len()
  
  let result = ()
  for i in range(rows_A) {
    let row = ()
    for j in range(cols_B) {
      let sum = complex(0, 0)
      for k in range(cols_A) {
        sum = complex-add(sum, complex-mul(A.at(i).at(k), B.at(k).at(j)))
      }
      row.push(sum)
    }
    result.push(row)
  }
  
  return result
}

#let mat-commutator(A, B) = {
  let AB = mat-mul(A, B)
  let BA = mat-mul(B, A)
  return mat-sub(AB, BA)
}

#let mat-trace(A) = {
  let sum = complex(0, 0)
  let n = A.len()
  
  for i in range(n) {
    sum = complex-add(sum, A.at(i).at(i))
  }
  
  return sum
}

// ============ GELL-MANN MATRICES ============

// Define Gell-Mann matrices without prefactors
#let t1_base = (
  (complex(0, 0), complex(1, 0), complex(0, 0)),
  (complex(1, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0))
)

#let t2_base = (
  (complex(0, 0), complex(0, -1), complex(0, 0)),
  (complex(0, 1), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0))
)

#let t3_base = (
  (complex(1, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(-1, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0))
)

#let t4_base = (
  (complex(0, 0), complex(0, 0), complex(1, 0)),
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(1, 0), complex(0, 0), complex(0, 0))
)

#let t5_base = (
  (complex(0, 0), complex(0, 0), complex(0, -1)),
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 1), complex(0, 0), complex(0, 0))
)

#let t6_base = (
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(1, 0)),
  (complex(0, 0), complex(1, 0), complex(0, 0))
)

#let t7_base = (
  (complex(0, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(0, -1)),
  (complex(0, 0), complex(0, 1), complex(0, 0))
)

#let t8_base = (
  (complex(1, 0), complex(0, 0), complex(0, 0)),
  (complex(0, 0), complex(1, 0), complex(0, 0)),
  (complex(0, 0), complex(0, 0), complex(-2, 0))
)

// Create the algebraic prefactors with symbolic fractions
#let half = alg-frac(1, 2)
#let sqrt3 = alg-sqrt(3)
#let half_sqrt3 = alg-div(half, sqrt3)

// Define prefactors for each matrix
#let prefactors = (
  half,           // t1: 1/2
  half,           // t2: 1/2
  half,           // t3: 1/2
  half,           // t4: 1/2
  half,           // t5: 1/2
  half,           // t6: 1/2
  half,           // t7: 1/2
  half_sqrt3      // t8: 1/(2√3)
)

// Scaled versions of the matrices for display
#let t_bases = (t1_base, t2_base, t3_base, t4_base, t5_base, t6_base, t7_base, t8_base)

// ============ STRUCTURE CONSTANT CALCULATION ============

// Debug function to examine the structure of an algebraic number
#let debug-alg(num) = {
  if num.type != "algebraic" {
    return "not an algebraic number: " + str(num)
  }
  
  let real_coef = "{ num: " + str(num.real.coef.num) + ", denom: " + str(num.real.coef.denom) + " }"
  let imag_coef = "{ num: " + str(num.imag.coef.num) + ", denom: " + str(num.imag.coef.denom) + " }"
  
  let result = "algebraic: {\n"
  result += "  real: { coef: " + real_coef + ", factors: " + str(num.real.factors) + " }\n"
  result += "  imag: { coef: " + imag_coef + ", factors: " + str(num.imag.factors) + " }\n"
  result += "}"
  
  return result
}

// Create the structure constants calculation function
#let calculate-f(a, b, c, Cr) = {
  let a_idx = a - 1
  let b_idx = b - 1
  let c_idx = c - 1
  
  // Calculate using base matrices (without prefactors)
  let comm = mat-commutator(t_bases.at(a_idx), t_bases.at(b_idx))
  let prod = mat-mul(comm, t_bases.at(c_idx))
  let tr = mat-trace(prod)
  
  // Get prefactors for each matrix
  let pa = prefactors.at(a_idx)
  let pb = prefactors.at(b_idx)
  let pc = prefactors.at(c_idx)
  
  // Combined prefactor (pa * pb * pc)
  let combined_prefactor = alg-mul(pa, alg-mul(pb, pc))
  
  // Convert trace to algebraic number
  let trace_alg = complex-to-alg(tr)
  
  // Structure constant formula: f^abc = -i/Cr * Tr([t_a, t_b] t_c) * combined_prefactor
  let i_factor = alg-div(alg-mul(alg-int(-1), alg-i), alg-int(Cr))
  
  // Calculate each intermediate result for debugging
  let trace_with_prefactor = alg-mul(trace_alg, combined_prefactor)
  
  // Final result with everything multiplied
  let final_result = alg-mul(i_factor, trace_with_prefactor)
  
  return (
    commutator: comm,
    product: prod,
    trace: tr,
    trace_alg: trace_alg,
    prefactor_a: pa,
    prefactor_b: pb,
    prefactor_c: pc,
    combined_prefactor: combined_prefactor,
    i_factor: i_factor,
    trace_with_prefactor: trace_with_prefactor,
    final_result: final_result
  )
}

#let Cr = 1 // Normalization constant

// ============ DISPLAY CALCULATIONS ============

#let display-calculation(a, b, c) = {
  let res = calculate-f(a, b, c, Cr)
  
  [
    == Structure Constant $f^(#a #b #c)$
    
    *Step 1:* Calculate the commutator $[t_(#a), t_(#b)]$
    
    $ [t_(#a), t_(#b)] = #format-alg(res.prefactor_a) dot.c #format-alg(res.prefactor_b) (mat(delim: "[",
      #format-complex(t_bases.at(a - 1).at(0).at(0)), #format-complex(t_bases.at(a - 1).at(0).at(1)), #format-complex(t_bases.at(a - 1).at(0).at(2)) ;
      #format-complex(t_bases.at(a - 1).at(1).at(0)), #format-complex(t_bases.at(a - 1).at(1).at(1)), #format-complex(t_bases.at(a - 1).at(1).at(2)) ;
      #format-complex(t_bases.at(a - 1).at(2).at(0)), #format-complex(t_bases.at(a - 1).at(2).at(1)), #format-complex(t_bases.at(a - 1).at(2).at(2))
    ) mat(delim: "[",
      #format-complex(t_bases.at(b - 1).at(0).at(0)), #format-complex(t_bases.at(b - 1).at(0).at(1)), #format-complex(t_bases.at(b - 1).at(0).at(2)) ;
      #format-complex(t_bases.at(b - 1).at(1).at(0)), #format-complex(t_bases.at(b - 1).at(1).at(1)), #format-complex(t_bases.at(b - 1).at(1).at(2)) ;
      #format-complex(t_bases.at(b - 1).at(2).at(0)), #format-complex(t_bases.at(b - 1).at(2).at(1)), #format-complex(t_bases.at(b - 1).at(2).at(2))
    ) - mat(delim: "[",
      #format-complex(t_bases.at(b - 1).at(0).at(0)), #format-complex(t_bases.at(b - 1).at(0).at(1)), #format-complex(t_bases.at(b - 1).at(0).at(2)) ;
      #format-complex(t_bases.at(b - 1).at(1).at(0)), #format-complex(t_bases.at(b - 1).at(1).at(1)), #format-complex(t_bases.at(b - 1).at(1).at(2)) ;
      #format-complex(t_bases.at(b - 1).at(2).at(0)), #format-complex(t_bases.at(b - 1).at(2).at(1)), #format-complex(t_bases.at(b - 1).at(2).at(2))
    ) mat(delim: "[",
      #format-complex(t_bases.at(a - 1).at(0).at(0)), #format-complex(t_bases.at(a - 1).at(0).at(1)), #format-complex(t_bases.at(a - 1).at(0).at(2)) ;
      #format-complex(t_bases.at(a - 1).at(1).at(0)), #format-complex(t_bases.at(a - 1).at(1).at(1)), #format-complex(t_bases.at(a - 1).at(1).at(2)) ;
      #format-complex(t_bases.at(a - 1).at(2).at(0)), #format-complex(t_bases.at(a - 1).at(2).at(1)), #format-complex(t_bases.at(a - 1).at(2).at(2))
    )) = #format-alg(alg-mul(res.prefactor_a, res.prefactor_b)) mat(delim: "[",
      #format-complex(res.commutator.at(0).at(0)), #format-complex(res.commutator.at(0).at(1)), #format-complex(res.commutator.at(0).at(2)) ;
      #format-complex(res.commutator.at(1).at(0)), #format-complex(res.commutator.at(1).at(1)), #format-complex(res.commutator.at(1).at(2)) ;
      #format-complex(res.commutator.at(2).at(0)), #format-complex(res.commutator.at(2).at(1)), #format-complex(res.commutator.at(2).at(2))
    ) $
    
    *Step 2:* Multiply by $t_(#c)$ to get $[t_(#a), t_(#b)]t_(#c)$
    
    $ [t_(#a), t_(#b)]t_(#c) = #format-alg(res.prefactor_c) dot.c #format-alg(alg-mul(res.prefactor_a, res.prefactor_b)) mat(delim: "[",
      #format-complex(res.commutator.at(0).at(0)), #format-complex(res.commutator.at(0).at(1)), #format-complex(res.commutator.at(0).at(2)) ;
      #format-complex(res.commutator.at(1).at(0)), #format-complex(res.commutator.at(1).at(1)), #format-complex(res.commutator.at(1).at(2)) ;
      #format-complex(res.commutator.at(2).at(0)), #format-complex(res.commutator.at(2).at(1)), #format-complex(res.commutator.at(2).at(2))
    ) mat(delim: "[",
      #format-complex(t_bases.at(c - 1).at(0).at(0)), #format-complex(t_bases.at(c - 1).at(0).at(1)), #format-complex(t_bases.at(c - 1).at(0).at(2)) ;
      #format-complex(t_bases.at(c - 1).at(1).at(0)), #format-complex(t_bases.at(c - 1).at(1).at(1)), #format-complex(t_bases.at(c - 1).at(1).at(2)) ;
      #format-complex(t_bases.at(c - 1).at(2).at(0)), #format-complex(t_bases.at(c - 1).at(2).at(1)), #format-complex(t_bases.at(c - 1).at(2).at(2))
    ) = #format-alg(res.combined_prefactor) mat(delim: "[",
      #format-complex(res.product.at(0).at(0)), #format-complex(res.product.at(0).at(1)), #format-complex(res.product.at(0).at(2)) ;
      #format-complex(res.product.at(1).at(0)), #format-complex(res.product.at(1).at(1)), #format-complex(res.product.at(1).at(2)) ;
      #format-complex(res.product.at(2).at(0)), #format-complex(res.product.at(2).at(1)), #format-complex(res.product.at(2).at(2))
    ) $
    
    *Step 3:* Calculate the trace
    
    $ Tr([t_(#a), t_(#b)]t_(#c)) = #format-alg(res.combined_prefactor) dot.c Tr(mat(delim: "[",
      #format-complex(res.product.at(0).at(0)), #format-complex(res.product.at(0).at(1)), #format-complex(res.product.at(0).at(2)) ;
      #format-complex(res.product.at(1).at(0)), #format-complex(res.product.at(1).at(1)), #format-complex(res.product.at(1).at(2)) ;
      #format-complex(res.product.at(2).at(0)), #format-complex(res.product.at(2).at(1)), #format-complex(res.product.at(2).at(2))
    )) = #format-alg(res.combined_prefactor) dot.c #format-complex(res.trace) = #format-alg(res.trace_with_prefactor) $
    
    *Step 4:* Multiply by $-i/C_r$ to get the structure constant
    
    $ f^(#a #b #c) = -frac(i, C_r) dot.c #format-alg(res.trace_with_prefactor) = #format-alg(res.final_result) $
  ]
}

// Calculate and display all structure constants where a != b != c
#let display-all-distinct-constants() = [
  = All SU(3) Structure Constants with Distinct Indices

  #for a in range(1, 9) {
    for b in range(1, 9) {
      for c in range(1, 9) {
        // Skip if any two indices are equal
        if a != b and b != c and a != c {
          display-calculation(a, b, c)
          v(1em) // Add vertical space between calculations
        }
      }
    }
  }
]

// Add this to the end of your document
#display-all-distinct-constants()

This works for any case where a != b != c, but as soon as one equality holds, some weird symbolic failures appear, and I cannot be bothered at this point getting something more robust working.

At this point, beyond the symbolic issues, the one problem I have is that reducing fractions using the Euclidean algorithm is off the table, since I wasn’t sure how to handle fractions with radicals involved.