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) + } +}