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

@@ -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])
+			// a is real but b is not: complex division still requires the
-			c[0].Quo(a[0], b[0])
+			// 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() {

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