From 351ba05ff413a9d5d9383641349ecf31ed0e8c48 Mon Sep 17 00:00:00 2001 From: DualFroz Date: Wed, 2 Sep 2026 05:04:22 +0200 Subject: [PATCH] fix(bignum): ComplexMultiplier.Quo wrong result and panic for real / non-real The a-real / b-non-real branch of Quo divided the real numerator by each component of b independently, which is not complex division: (1+0i)/(1+i) returned 1+1i instead of 0.5-0.5i, (1+0i)/i returned +Inf+1i, and (0+0i)/i divided 0/0 and panicked. Use the full a*conj(b)/|b|^2 formula (the same one the general branch uses), which reduces correctly when a is real. Closes #577 --- utils/bignum/complex.go | 20 ++++++++++++++++-- utils/bignum/complex_test.go | 39 ++++++++++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 2 deletions(-) create mode 100644 utils/bignum/complex_test.go diff --git a/utils/bignum/complex.go b/utils/bignum/complex.go index 6e74b8b..bbe5d4f 100644 --- a/utils/bignum/complex.go +++ b/utils/bignum/complex.go @@ -197,8 +197,24 @@ func (cEval *ComplexMultiplier) Quo(a, b, c *Complex) { c[0].Quo(a[0], b[0]) c[1].SetFloat64(0) } else { - c[1].Quo(a[0], b[1]) - c[0].Quo(a[0], b[0]) + // a is real but b is not: complex division still requires the + // full formula a/b = a*conj(b) / |b|^2. Dividing the real + // numerator by each component of b independently is incorrect + // (and divides 0/0 when a is zero). + cEval.tmp0.Mul(a[0], b[0]) + cEval.tmp1.Mul(a[1], b[1]) + cEval.tmp2.Mul(a[1], b[0]) + cEval.tmp3.Mul(a[0], b[1]) + + cEval.tmp0.Add(cEval.tmp0, cEval.tmp1) + cEval.tmp1.Sub(cEval.tmp2, cEval.tmp3) + + cEval.tmp2.Mul(b[0], b[0]) + cEval.tmp3.Mul(b[1], b[1]) + cEval.tmp2.Add(cEval.tmp2, cEval.tmp3) + + c[0].Quo(cEval.tmp0, cEval.tmp2) + c[1].Quo(cEval.tmp1, cEval.tmp2) } } else { if b.IsReal() { diff --git a/utils/bignum/complex_test.go b/utils/bignum/complex_test.go new file mode 100644 index 0000000..d57a165 --- /dev/null +++ b/utils/bignum/complex_test.go @@ -0,0 +1,39 @@ +package bignum + +import ( + "math/big" + "testing" + + "github.com/stretchr/testify/require" +) + +func TestComplexMultiplierQuo(t *testing.T) { + const prec = 128 + + newComplex := func(re, im float64) *Complex { + return &Complex{ + new(big.Float).SetPrec(prec).SetFloat64(re), + new(big.Float).SetPrec(prec).SetFloat64(im), + } + } + f := func(v float64) *big.Float { return new(big.Float).SetPrec(prec).SetFloat64(v) } + + cm := NewComplexMultiplier() + for _, tc := range []struct { + name string + a, b *Complex + wantRe, wantIm float64 + }{ + // Real numerator, non-real denominator: the previously broken branch. + {"1/(1+i)", newComplex(1, 0), newComplex(1, 1), 0.5, -0.5}, + {"1/i", newComplex(1, 0), newComplex(0, 1), 0, -1}, + {"0/i", newComplex(0, 0), newComplex(0, 1), 0, 0}, // must not panic + // General complex division (regression guard for the other branch). + {"(2+3i)/(1+i)", newComplex(2, 3), newComplex(1, 1), 2.5, 0.5}, + } { + c := newComplex(0, 0) + cm.Quo(tc.a, tc.b, c) + require.Zerof(t, c[0].Cmp(f(tc.wantRe)), "%s real: got %v, want %v", tc.name, c[0], tc.wantRe) + require.Zerof(t, c[1].Cmp(f(tc.wantIm)), "%s imag: got %v, want %v", tc.name, c[1], tc.wantIm) + } +}