Gaussian integers

The ring $\mathbb{Z}[i]$ of Gaussian integers is a Euclidean domain with respect to the norm $N(a + bi) = a^2 + b^2$ , and hence a unique factorisation domain. Its units are $\pm 1, \pm i$ , and $1 + i$ is the unique prime above 2, with $2 = -i(1 + i)^2$ . Standard references are [Coh1996, §1.2] for Euclidean rings and [Wei2000] for the binary algorithm.

Division with remainder#

Proposition 1. For $x, y \in \mathbb{Z}[i]$ with $y \ne 0$ , let $q = \lfloor \operatorname{re} (x/y) + 1/2 \rfloor + i \lfloor \operatorname{im}(x/y) + 1/2 \rfloor$ and $r = x - qy$ . Then $N(r) \le N(y)/2$ .

Proof. Each of the real and imaginary parts of $x/y - q$ lies in $(-1/2, 1/2]$ , so $|x/y - q|^2 \le 1/2$ , and $N(r) = N(y)\,|x/y - q|^2$ .

fmpzi_divrem computes this quotient exactly: with $t = 2x\bar y + N(y)(1 + i)$ , it takes the floor quotients of the real and imaginary parts of $t$ by $2N(y)$ , which costs one complex multiplication and two divisions. Ties are rounded upwards in each coordinate, which makes the remainder unique. When $x$ has at least three bits fewer than $y$ (in the larger of its parts), $|x/y| < 1/2$ and the quotient is zero without computation.

fmpzi_divrem_approx only guarantees $N(r) < N(y)$ , which is all the Euclidean algorithm needs. When the bit sizes of $x$ and $y$ differ by less than 45, it computes the quotient from double-precision approximations of $x$ and $y$ (scaled by a common power of 2 if they are large). Each coordinate of the computed quotient then differs from the corresponding coordinate of $x/y$ by at most $1/2 + \varepsilon$ with a tiny $\varepsilon$ , so $|x/y - q|^2 \le 2(1/2 + \varepsilon)^2 < 1$ . Otherwise it calls fmpzi_divrem.

Canonical associates#

Every nonzero $x$ has exactly one associate $i^k x$ with real part positive and imaginary part in $(-\operatorname{re}, \operatorname{re}]$ : the four associates are the rotations of $x$ by multiples of $90°$ , and exactly one of them lies in the sector $-\pi/4 < \arg z \le \pi/4$ . fmpzi_canonical_unit_i_pow returns the exponent $k \in \{0, 1, 2, 3\}$ from two comparisons, and fmpzi_canonicalise_unit applies it. The GCD functions return the canonical associate, so that the GCD is unique; for example $\gcd(2, 1 + i) = 1 + i$ , while $1 - i$ , $-1 + i$ and $-1 - i$ are its other associates.

Greatest common divisors#

Euclidean algorithm#

By Proposition 1 each step of $(x, y) \gets (y, x \bmod y)$ at least halves the norm, so for inputs of $n$ bits there are at most about $2n$ steps, and the total cost is $O(n^2)$ bit operations as for integers. fmpzi_gcd_euclidean uses fmpzi_divrem. fmpzi_gcd_euclidean_improved uses fmpzi_divrem_approx, and as soon as all four parts are below $2^{50}$ in absolute value it finishes in double-precision arithmetic: the approximate quotient may differ from the rounded one, but the remainder $x - qy$ is computed exactly in doubles because all intermediate values stay below $2^{53}$ , and any quotient with $N(r) < N(y)$ gives a correct Euclidean sequence.

The (1+i)-ary algorithm#

Weilert's analogue of the binary GCD [Wei2000] replaces divisions by additions and divisions by $1 + i$ . Since $\mathbb{Z}[i]/(1 + i) \cong \mathbf{F}_2$ , any two elements $x, y$ not divisible by $1 + i$ satisfy $x \equiv y \equiv 1$ , so all four of $x \pm y$ and $x \pm iy$ are divisible by $1 + i$ . fmpzi_gcd_binary first removes the factors $1 + i$ from both inputs, remembering the smaller multiplicity, and then repeatedly replaces the larger of $x, y$ (by norm) by the one of $x \pm y$ , $x \pm iy$ of smallest norm, divided by the largest power of $1 + i$ . The four norms sum to $4(N(x) + N(y))$ , so the smallest is at most $N(x) + N(y) \le 2N(x)$ , and the division by $1 + i$ at least halves it; Weilert shows that the algorithm has quadratic bit complexity, like the binary GCD for integers [Wei2000]. The implementation estimates the norms in double precision and falls back to fmpzi_gcd_euclidean when an input has more than 500 bits.

The lattice algorithm#

The ideal $(x, y) = g\mathbb{Z}[i]$ generated by $x$ and $y$ is a lattice in $\mathbb{Z}^2$ , generated over $\mathbb{Z}$ by $x, ix, y, iy$ (as vectors of real and imaginary parts). Its Hermite normal form is $\begin{pmatrix} c & b \\ 0 & a \end{pmatrix}$ with $ac = N(g)$ , and $g$ is, up to a unit, a shortest nonzero vector of the lattice: every element is $gw$ with $N(w) \ge 1$ . fmpzi_gcd_shortest, due to D. Schultz, computes the Hermite form by integer GCDs of the coordinates and then the shortest vector of the two-dimensional lattice in the $\ell^\infty$ norm from the continued fraction expansion of an interval around $a/b$ , computed by the quasi-linear ball continued fraction routine _fmpq_ball_get_cfrac of fmpq. A special case handles purely real or purely imaginary $g$ , where a vector $g(1 + i)$ has the same $\ell^\infty$ norm as $g$ .

Choice of algorithm#

fmpzi_gcd returns the canonical associate of the other input if one is zero. If all parts are below $2^{50}$ in absolute value it uses the double-precision Euclidean algorithm directly. If either input has both parts of one word, or either input has at most 30000 bits, it uses fmpzi_gcd_euclidean_improved. Otherwise it uses fmpzi_gcd_shortest, first reducing the larger input modulo the smaller with one approximate division if their sizes differ by more than a factor 1.5.

The prime 1 + i#

For $x = a + bi \ne 0$ , let $s = \min(\operatorname{ord}_2 a, \operatorname{ord}_2 b)$ (with $\operatorname{ord}_2 0 = \infty$ ). Then $x = 2^s x'$ with $x'$ not divisible by 2, and $x'$ is divisible by $1 + i$ , exactly once, if and only if both parts of $x'$ are odd. Since $2^s = (1 + i)^{2s} i^{-s}$ , the multiplicity of $1 + i$ in $x$ is $2s$ or $2s + 1$ . fmpzi_remove_one_plus_i computes it from the 2-adic valuations of the parts, divides by $2^s$ with shifts, multiplies by $i^{s}$ , and divides by $1 + i$ once more using $(a + bi)/(1 + i) = ((a + b) + (b - a)i)/2$ if needed.

Primality#

The Gaussian primes are, up to units, $1 + i$ , the rational primes $p \equiv 3 \pmod 4$ , and the elements $\pi$ with $N(\pi) = p$ a rational prime $p \equiv 1 \pmod 4$ [Coh1996, §1.2]. Hence $a + bi$ is prime if and only if either $ab = 0$ and the nonzero part is $\pm p$ with $p$ prime and $p \equiv 3 \pmod 4$ , or $ab \ne 0$ and $a^2 + b^2$ is prime. fmpzi_is_prime applies this criterion with fmpz_is_prime, which proves primality, and fmpzi_is_probabprime with fmpz_is_probabprime. FLINT does not provide factorisation in $\mathbb{Z}[i]$ , which reduces to factoring the norm and, for each prime $p \equiv 1 \pmod 4$ dividing it, finding $\pi = \gcd(p, t + i)$ with $t^2 \equiv -1 \pmod p$ .

Use in FLINT#

Division: fmpzi_divrem, fmpzi_divrem_approx, fmpzi_divexact.
Normalisation: fmpzi_canonical_unit_i_pow, fmpzi_canonicalise_unit.
GCD: fmpzi_gcd_euclidean, fmpzi_gcd_euclidean_improved, fmpzi_gcd_binary, fmpzi_gcd_shortest, fmpzi_gcd; the last is the GCD of the generic ring gr_ctx_init_fmpzi.
fmpzi_remove_one_plus_i, fmpzi_is_prime, fmpzi_is_probabprime.
Multiplication: fmpzi_mul and fmpzi_sqr use the complex multiplication flint_mpn_mul_complex when all parts are multiword; polynomial products use _fmpzi_poly_mulmid.