# core design

This page derives the algebra of `Complex[T]`, states when each
`luna-generic` instance is lawful, and explains the choices behind the
generic type: mutable fields, the unscaled division formula, and the
separation from floating-point analysis.

## Design goal

One generic complex type for every scalar of the ecosystem, whose
instances advertise exactly the algebraic structure the construction has,
and which leaves IEEE special values, branch cuts and transcendental
functions to a backend package.

## Constraints

- `T` ranges over exact rings (`Int`, `BigInt`, exact rationals), the
  floating-point types `Float` and `Double`, and other complex types, so the
  core may use only the operations of the `luna-generic` traits: no order, no
  absolute value, no special values.
- An instance must not claim more structure than the construction has under
  its bound.
- A package cannot add methods or trait instances to a type of another
  package, so whatever needs floating-point semantics has to live in a
  separate package as free functions.

## Main design decisions

- One generic `Complex[T]` with public mutable fields.
- Instances follow the structure of `T`: componentwise operations need only
  the matching operations of `T`, multiplication needs `T : Ring`, inverses
  need `T : Field`.
- The schoolbook product with four multiplications.
- Textbook division through `Inverse::inv` of the norm, because a generic
  field has nothing else.
- No analytic functions in the core; `float_backend` provides them for
  `Double`.

## Mathematical background

### The construction

For a ring $R$, the complex numbers over $R$ are the quotient of the
polynomial ring in a central indeterminate $x$ by the two-sided ideal
generated by $x^2 + 1$:

$$
R[i] = R[x]/(x^2 + 1), \qquad i = x + (x^2 + 1), \qquad i^2 = -1 .
$$

Division by the monic polynomial $x^2 + 1$ leaves a unique remainder of
degree below two, so every element is $a + bi$ with unique $a, b \in R$:
the pair `(re, im)`. Computing with polynomials and reducing $i^2 = -1$
gives the operations of the package:

$$
\begin{aligned}
(a + bi) \pm (c + di) &= (a \pm c) + (b \pm d)i, \\
(a + bi)(c + di) &= ac + (ad + bc)i + bd\,i^2 = (ac - bd) + (ad + bc)i .
\end{aligned}
$$

As a quotient of a ring by a two-sided ideal, $R[i]$ is a ring, with zero
$0 + 0i$ and one $1 + 0i$; $R$ embeds as $a \mapsto a + 0i$. The product
formula keeps the factors of $R$ in their order ($ac$, $ad$, $bc$, $bd$) and
moves only $i$, which commutes with everything because $x$ is central. So
the formula is associative and distributive even when $R$ is not
commutative, which is why `MulMonoid`, `Semiring` and `Ring` need only
`T : Ring`. $R[i]$ is commutative exactly when $R$ is.

The rest of this page assumes $R$ commutative: the norm below lies in $R$
only then. For $z = a + bi$ over a non-commutative $R$,
$z\bar z = (a^2 + b^2) + (ba - ab)i$.

### Conjugation and the norm

$\overline{a + bi} = a - bi$ is the map induced by $x \mapsto -x$, which
fixes the ideal $(x^2 + 1)$, so it is a ring automorphism of order two:

$$
\overline{z + w} = \bar z + \bar w, \qquad \overline{zw} = \bar z\,\bar w,
\qquad \bar{\bar z} = z .
$$

The norm $N(z) = z\bar z = (a + bi)(a - bi) = a^2 + b^2$ lies in $R$ and is
multiplicative, $N(zw) = zw\,\bar z\bar w = N(z)N(w)$.

### Inverses and division

If $N(z)$ is invertible in $R$, then $z \cdot \bar z\,N(z)^{-1} = 1$, so

$$
z^{-1} = \frac{\bar z}{N(z)} = \frac{a}{a^2 + b^2} - \frac{b}{a^2 + b^2}\,i,
\qquad
\frac{a + bi}{c + di} = \frac{(ac + bd) + (bc - ad)i}{c^2 + d^2} .
$$

Conversely, if $z$ is invertible then $N(z)N(z^{-1}) = N(1) = 1$, so $N(z)$
is invertible. Hence $z$ is a unit exactly when $a^2 + b^2$ is.

### When the construction is a field

Let $K$ be a field. $K[x]/(x^2 + 1)$ is a field exactly when $x^2 + 1$ is
irreducible over $K$, that is when $-1$ is not a square in $K$. Directly: a
non-zero $z$ with $N(z) = a^2 + b^2 = 0$ must have $b \ne 0$, and then
$(a/b)^2 = -1$. Conversely, if $s^2 = -1$ then $(s + i)(s - i) = s^2 + 1 =
0$ exhibits zero divisors.

- For $K = \mathbb R$ (`Float`, `Double`), $-1$ is not a square and
  $\mathbb R[i] = \mathbb C$ is a field.
- For $K = \mathbb C$ (`Complex[Double]`), $-1 = i^2$ is a square, so
  `Complex[Complex[Double]]` $\cong \mathbb C[j]/(j^2 + 1) \cong \mathbb C
  \times \mathbb C$ is a ring with zero divisors, such as $(1 + ij)(1 - ij)
  = 1 - i^2 j^2 = 0$.

## Decisions in detail

### Instances follow the structure of `T`

**Problem.** Which `luna-generic` traits may `Complex[T]` implement?

**Choice.** `Zero`, `One`, `AddMonoid` and `AddGroup` need only the
corresponding operations of `T`, because they act componentwise. `MulMonoid`,
`Semiring` and `Ring` need `T : Ring`, because the product uses subtraction
($ac - bd$) and the ring laws of $R[i]$ need those of $R$. `Inverse`,
`MulGroup`, `Field` and the `Div` operator need `T : Field` for $N(z)^{-1}$.
`Conjugate` needs only `Neg`.

The `Field` instance is declared for every `T : Field`. By the previous
section that is lawful when $-1$ is not a square in `T`, which covers the
real scalar types; for `T = Complex[Double]` it is not, and `inv` aborts on
non-zero zero divisors. A trait bound cannot express "−1 is not a square",
so the instance trusts the caller here.[^honesty]

[^honesty]: This is a known gap in the "implement only lawful instances"
rule of Luna Flow. The test suite uses nested complex numbers only on
values with an invertible norm.

### Four multiplications

**Options.** The schoolbook product $(ac - bd) + (ad + bc)i$ with four
multiplications; Gauss's three-multiplication form $k_1 = c(a + b)$, $k_2 =
a(d - c)$, $k_3 = b(c + d)$, product $(k_1 - k_3) + (k_1 + k_2)i$.

**Choice.** Four multiplications. Over floating point the schoolbook form
has a normwise relative error of at most $\sqrt5\,u$,[^bpz] while the
three-multiplication form adds cancellation-prone sums such as $a + b$ and
$d - c$ and is less accurate componentwise; and for generic `T` a multiplication is not
necessarily more expensive than an addition.

[^bpz]: R. Brent, C. Percival and P. Zimmermann, "Error bounds on complex
floating-point multiplication", *Mathematics of Computation* 76 (2007).

Componentwise, $\mathrm{fl}(ac - bd) = ac(1 + \theta_2) - bd(1 + \theta_2')$
with $|\theta_2| \le \gamma_2 = 2u/(1 - 2u)$, so

$$
|\mathrm{fl}(ac - bd) - (ac - bd)| \le \gamma_2\,(|ac| + |bd|),
$$

which is small relative to the result unless $ac$ and $bd$ cancel.

### Textbook division in the generic core

**Problem.** Division needs $N(w)^{-1}$; for floating point, $c^2 + d^2$
overflows or underflows long before $c + di$ does.

**Options.** Scaled algorithms (Smith's, or rescaling by a power of two)
need comparisons and absolute values that a generic field does not have;
the textbook formula needs only field operations.

**Choice.** The generic `Div` and `Inverse` use the textbook formula with
`Inverse::inv` of `T`. It is exact for exact fields. For `Double` it is
correct while $c^2 + d^2$ stays in range, and `Inverse` of `Double` aborts
on zero, including an underflowed norm. Scaled, special-value-aware
division lives in [`float_backend`](float_backend.md#smiths-division).

### Mutable fields

**Problem.** Complex values are often updated in loops (accumulators,
recurrences), where allocating a new value per step is wasteful.

**Choice.** `Complex[T]` is a `pub(all)` struct with `mut re` and `mut im`
plus `set`, `set_re` and `set_im`. Every operation returns a new value and
never mutates its inputs, so code that does not call the setters can treat
complex numbers as values. Callers that do mutate must remember that the
struct is shared by reference.

### No analytic functions in the core

$\sqrt z$, $\log z$ or $\sin z$ need an order, absolute values,
transcendental functions of the real scalar, a choice of branch, and IEEE
special values. None of these exist for a generic field. The root package
therefore stops at algebra, and the [float_backend](float_backend.md)
package supplies analysis for `Double`.

### Text and debug forms

`Show` prints `re + imi` using `T`'s text (`1 + -2i` for $1 - 2i$), which
is the canonical text of the type and the reason `to_string` stays a
promoted method. `Debug` prints the record form for tests and diagnostics.
MoonBit 0.10 promotion of the other trait methods is explicit in
[`src/extends.mbt`](../../../src/extends.mbt).

## Correctness and invariants

- Ring laws hold for `T : Ring` because $R[i]$ is a quotient ring; the test
  suite checks associativity, identities and $\overline{zw} = \bar z\bar w$
  on integers.
- `z * z.conjugate()` equals $N(z) + 0i$ exactly for exact `T`.
- `z * w / w == z` and `z.inv() * z == 1` hold for exact fields when
  $N(w)$, $N(z)$ are units, and up to rounding for `Double` within range.
- No operation mutates its arguments; only the three setters do.

### Floating-point instances

For `T = Double` every law holds only up to rounding, and the unscaled norm
adds hard limits. Let $\Omega \approx 1.80 \times 10^{308}$ be the largest
finite `Double`, $\lambda = 2^{-1022} \approx 2.2 \times 10^{-308}$ the
smallest normal one and $\mu = 2^{-1074} \approx 4.9 \times 10^{-324}$ the
smallest subnormal one. The norm $n = \mathrm{fl}(c^2 + d^2)$ of $w = c + di$
behaves as follows:

$$
\begin{aligned}
n = \infty &\iff c^2 + d^2 \gtrsim \Omega
&&\Longleftarrow\ |w| > \sqrt{\Omega} \approx 1.34 \times 10^{154}, \\
n < \lambda &\iff |w| \lesssim \sqrt{\lambda} \approx 1.49 \times 10^{-154}, \\
\mathrm{fl}(1/n) = \infty &\iff n < 1/\Omega
&&\Longleftarrow\ |w| \lesssim \Omega^{-1/2} \approx 7.46 \times 10^{-155}, \\
n = 0 &\iff c^2 \le \mu/2 \text{ and } d^2 \le \mu/2
&&\Longleftarrow\ \max(|c|, |d|) \lesssim \sqrt{\mu/2} \approx 1.5 \times 10^{-162}.
\end{aligned}
$$

The last line uses that a square below half the smallest subnormal rounds to
zero. In the first case $n^{-1} = 0$ and `inv` returns zero parts, so
$z z^{-1} = 0 \ne 1$; in the third the parts are $\pm\infty$ or NaN; in the
fourth `Inverse::inv` of `Double` aborts. The field law $z z^{-1} = 1$ for
every $z \ne 0$ therefore fails outside roughly
$10^{-154} \le |z| \le 10^{154}$, which is the honest domain of the `Field`
instance for `Complex[Double]`. [`float_backend`](float_backend.md#smiths-division)
divides without forming $c^2 + d^2$.

## Alternatives rejected

- **A `Double`-only complex type.** It would duplicate the type for `Float`,
  integers (Gaussian integers) and exact rationals.
- **Immutable fields.** Cleaner value semantics, but a breaking change for
  existing callers that update `re` and `im` in place.
- **Scaled division in the core.** It needs order and magnitude operations
  that generic fields lack.

## Boundaries

- No analytic or transcendental functions, no special-value handling, no
  branch cuts; see [float_backend](float_backend.md).
- No polar representation in the core.
- No check that `T` makes $R[i]$ a field; the `Field` instance trusts it.
- No overflow-safe division for floating-point `T`.