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.