Imaginary Numbers in Cryptography

§ 1

Why Do Imaginary Numbers Appear in Cryptography?

At first glance, cryptography seems to live entirely in the world of integers — modular arithmetic, prime factorization, discrete logarithms. But imaginary and complex numbers weave through it at every level: as the underlying algebraic structure of the fastest algorithms, as the language of elliptic curve theory, and as the foundation of post-quantum cryptography.

There are six distinct ways imaginary numbers enter cryptography, each contributing something different:

Gaussian Integers ℤ[i]

Integers extended with i. A unique factorization domain with "Gaussian primes." Basis for Gaussian-integer variants of RSA and DH.

√−1 in Finite Fields

For prime p ≡ 1 (mod 4), −1 has a square root in ℤ_p. Enables the GLV scalar decomposition trick that speeds up ECC by ~50%.

Roots of Unity & NTT

The Number Theoretic Transform replaces complex roots of unity (e^2πi/n) with modular ones. Heart of Kyber and Dilithium performance.

Cyclotomic Polynomials

Φ_n(x) whose roots are primitive nth roots of unity. Ring-LWE lives in ℤ[x]/Φ_n(x) — the algebraic structure that makes lattice crypto efficient.

Complex Multiplication

Elliptic curves over imaginary quadratic fields ℚ(√−d). Used to construct curves with specific group orders for pairing-friendly cryptography.

Bilinear Pairings

Weil and Tate pairings map to roots of unity in field extensions. Enable IBE, BLS signatures, and zero-knowledge proofs.

§ 2

Gaussian Integers ℤ[i]

The Gaussian integers are complex numbers a + bi where both a and b are ordinary integers. They extend the familiar integer number line into a 2D lattice on the complex plane.

Definition — Gaussian Integers

ℤ[i] = { a + bi : a, b ∈ ℤ, i² = −1 }
Addition: (a+bi) + (c+di) = (a+c) + (b+d)i
Multiplication: (a+bi)(c+di) = (ac−bd) + (ad+bc)i
Norm: N(a+bi) = a² + b² (always a non-negative integer)

Gaussian Primes

A Gaussian integer π is a Gaussian prime if it cannot be factored into smaller Gaussian integers. The rules for which ordinary primes stay prime in ℤ[i] are determined by whether −1 is a quadratic residue mod p:

Gaussian Prime Classification p = 2 → 2 = −i(1+i)² (ramifies; 1+i is a Gaussian prime, N(1+i)=2) p ≡ 1 (mod 4) → p splits: p = π·π̄ into two conjugate Gaussian primes e.g., 5 = (2+i)(2−i), 13 = (3+2i)(3−2i) p ≡ 3 (mod 4) → p stays prime in ℤ[i] (inert) e.g., 3, 7, 11, 19, 23 are all Gaussian primes

Theorem — Fermat's Two-Square

An odd prime p is the sum of two squares p = a² + b² if and only if p ≡ 1 (mod 4). This is exactly when p factors in ℤ[i] as p = (a+bi)(a−bi). Example: 29 = 25+4 = 5²+2² = (5+2i)(5−2i).

Gaussian integers Gaussian primes Split primes (p≡1 mod 4)

Division with Remainder in ℤ[i]

ℤ[i] is a Euclidean domain: you can always divide with remainder, enabling a Gaussian GCD algorithm. Division uses the complex norm as the "size" measure:

Gaussian Division Algorithm // Divide α by β: find q, r such that α = q·β + r, N(r) < N(β) q = round( α/β ) // round each of real and imaginary parts to nearest integer r = α − q·β // Example: (11+3i) ÷ (2+5i) α/β = (11+3i)(2−5i)/29 = (37−49i)/29 ≈ 1.28 − 1.69i → round to 1−2i q = 1−2i, r = (11+3i) − (1−2i)(2+5i) = (11+3i)−(12+i) = −1+2i N(r) = 1+4 = 5 < N(2+5i) = 29 ✓

GNU Octave Gaussian integer arithmetic, GCD, Gaussian prime test

%% Gaussian Integer Arithmetic in GNU Octave
%% Represent a+bi as [a, b]; all operations return [real_part, imag_part]

function r = gn(z)
  r = z(1)^2 + z(2)^2;       % Gaussian norm N(a+bi) = a²+b²
end
function r = gmul(u, v)
  r = [u(1)*v(1)-u(2)*v(2),  u(1)*v(2)+u(2)*v(1)];
end
function [q, r] = gdivrem(alpha, beta)
  % α = q·β + r,  N(r) < N(β)
  nb = gn(beta);
  qr = (alpha(1)*beta(1) + alpha(2)*beta(2)) / nb;
  qi = (alpha(2)*beta(1) - alpha(1)*beta(2)) / nb;
  q  = [round(qr), round(qi)];
  r  = alpha - gmul(q, beta);
end
function g = ggcd(a, b)
  while any(b ~= 0)
    [~, r] = gdivrem(a, b);
    a = b; b = r;
  end
  g = a;
end
function tf = is_gaussian_prime(z)
  % A Gaussian integer z is prime if N(z) is prime, or
  % z = unit * p where p is an ordinary prime ≡ 3 (mod 4)
  n   = gn(z);
  tf  = isprime(n) || (z(2)==0 && isprime(abs(z(1))) && mod(abs(z(1)),4)==3) || ...
         (z(1)==0 && isprime(abs(z(2))) && mod(abs(z(2)),4)==3);
end

% Factoring 5 in Z[i]: 5 = (2+i)(2-i)
a = [2,1]; b = [2,-1];
prod = gmul(a,b);
printf('(2+i)(2-i) = %d+%di  [expect 5+0i]\n', prod(1), prod(2));
printf('2+i Gaussian prime? %d\n', is_gaussian_prime([2,1]));
printf('3   Gaussian prime? %d  (3≡3 mod 4)\n', is_gaussian_prime([3,0]));
printf('5   Gaussian prime? %d  (splits)\n',   is_gaussian_prime([5,0]));

% Gaussian GCD: gcd(11+3i, 1+8i)
g = ggcd([11,3], [1,8]);
printf('gcd(11+3i, 1+8i) = %d+%di\n', g(1), g(2));

% Print Gaussian primes with norm ≤ 25
printf('\nGaussian primes a+bi with 0≤a,b≤5:\n');
for a = 0:5; for b = 0:5
  if (a+b>0) && is_gaussian_prime([a,b])
    printf('  %d+%di (N=%d)\n', a, b, a^2+b^2);
  end
end; end

§ 3

The Square Root of −1 in Finite Fields

In the reals, √−1 doesn't exist. But in a finite field 𝔽_p, it exists exactly when −1 is a quadratic residue — i.e., when p ≡ 1 (mod 4). This is not philosophical: it's a concrete, computable number modulo p.

Existence Criterion √−1 exists in 𝔽_p ⟺ p ≡ 1 (mod 4) // By Euler's criterion: (−1)^((p-1)/2) ≡ 1 (mod p) ⟺ (p-1)/2 is even ⟺ p≡1(mod 4) // Computing √−1: use the Tonelli-Shanks algorithm, or // the shortcut when p ≡ 1 (mod 4): Find generator g of 𝔽_p*, then i = g^(p−1)/4 mod p // Example: p = 17 ≡ 1 (mod 4); generator g = 3 i = 3^(17−1)/4 = 3⁴ = 81 ≡ 13 (mod 17) Verify: 13² = 169 = 9·17 + 16 ≡ −1 (mod 17) ✓

Wilson's Theorem Shortcut

When p ≡ 1 (mod 4), a formula for √−1 mod p is i ≡ ((p−1)/2)! (mod p). This is Wilson's theorem in disguise: (p−1)! ≡ −1 (mod p), which splits nicely when (p−1)/2 is even.

GNU Octave Compute √−1 mod p; verify quadratic residuosity

%% Finding sqrt(-1) mod p for primes p ≡ 1 (mod 4)

function r = mod_pow(base, exp, m)
  r = 1; base = mod(base, m);
  while exp > 0
    if bitand(exp,1); r=mod(r*base,m); end
    exp=bitshift(exp,-1); base=mod(base^2,m);
  end
end

function g = find_generator(p)
  for g = 2:p-1
    if mod_pow(g, (p-1)/2, p) ~= 1; return; end
  end
end

function i = sqrt_neg1_mod_p(p)
  if mod(p, 4) ~= 1
    error('sqrt(-1) does not exist mod %d (p ≡ 3 mod 4)', p); 
  end
  g = find_generator(p);
  i = mod_pow(g, (p-1)/4, p);
end

% Test several primes
test_primes = [5, 13, 17, 29, 41, 61, 101, 10007];
printf('%-8s %-6s %-10s %-8s\n', 'p', 'p mod 4', 'sqrt(-1)', 'verify');
printf('%s\n', repmat('-',36,1));
for p = test_primes
  if mod(p,4)==1
    s   = sqrt_neg1_mod_p(p);
    chk = mod(s^2, p);   % should be p-1
    printf('%-8d %-6d %-10d %-8d\n', p, mod(p,4), s, chk==p-1);
  else
    printf('%-8d %-6d %-10s %-8s\n', p, mod(p,4), 'N/A', '—');
  end
end

%% The secp256k1 (Bitcoin) field prime: p = 2^256 - 2^32 - 977
%% p ≡ 3 (mod 4), so sqrt(-1) does NOT exist in 𝔽_p directly.
%% But the cube root of 1 does, enabling the GLV endomorphism (see §7).
printf('\nsecp256k1 field prime mod 4 = %d (no sqrt(-1) in F_p)\n', ...
       mod(3, 4));  % symbolic only; actual prime too large for floating point

§ 4

Gaussian Integer RSA

Since ℤ[i] is a unique factorization domain with a notion of "Gaussian primes," it is possible to construct an RSA-like scheme over Gaussian integers. The security rests on the hardness of factoring a Gaussian composite — knowing N(π₁·π₂) but not the factors.

Gaussian RSA Key Generation // Choose two Gaussian primes π₁, π₂ with large norms n = π₁ · π₂ // Gaussian composite φ(n) = (N(π₁)−1)(N(π₂)−1) // analogue of Euler totient e = choose with gcd(e, φ(n)) = 1 d = e⁻¹ mod φ(n) // Encryption: for Gaussian integer message M = a+bi C = M^e mod n // Gaussian modular exponentiation M = C^d mod n

GNU Octave Gaussian integer modular exponentiation and toy RSA

%% Gaussian Integer RSA (toy example)
%% Gaussian modular reduction: (a+bi) mod n means reduce norm

function r = gmod(alpha, n)
  % Reduce Gaussian integer alpha mod ordinary integer n
  r = mod(alpha, n);
end
function r = gmul_mod(u, v, n)
  r = mod([u(1)*v(1)-u(2)*v(2),  u(1)*v(2)+u(2)*v(1)], n);
end
function r = gpow_mod(base, exp, n)
  % Square-and-multiply for Gaussian integers
  r = [1, 0];
  base = mod(base, n);
  while exp > 0
    if bitand(exp, 1)
      r = gmul_mod(r, base, n);
    end
    base = gmul_mod(base, base, n);
    exp  = bitshift(exp, -1);
  end
end

% Key setup: Gaussian primes π₁ = 2+3i (N=13), π₂ = 4+i (N=17)
% n = π₁·π₂;  use ordinary n = N(π₁)·N(π₂) = 13·17 = 221 as modulus
N1 = 13; N2 = 17;    % norms of the two Gaussian primes
n  = N1 * N2;          % = 221
phi= (N1-1)*(N2-1);   % = 192
e  = 5;                 % gcd(5,192) = 1 ✓
% d = e^{-1} mod phi
function inv = modinv(a,m)
  [g,x,~]=gcd(mod(a,m),m); inv=mod(x,m);
end
d  = modinv(e, phi);
printf('Gaussian RSA: n=%d, e=%d, d=%d, phi=%d\n', n, e, d, phi);

% Encrypt Gaussian message M = 5+7i
M   = [5, 7];
C   = gpow_mod(M, e, n);
Dec = gpow_mod(C, d, n);
printf('M = %d+%di\n', M(1), M(2));
printf('C = M^e mod n = %d+%di\n', C(1), C(2));
printf('Decrypted = C^d mod n = %d+%di  OK=%d\n', Dec(1), Dec(2), isequal(Dec,M));

§ 5

Cyclotomic Polynomials & Roots of Unity

The n-th roots of unity are the n complex solutions to zⁿ = 1. The primitive n-th roots of unity satisfy no lower-degree equation; they are the roots of the cyclotomic polynomial Φ_n(x).

Cyclotomic Polynomial Definition Φₙ(x) = ∏ (x − ω^k) over primitive n-th roots ω = e^2πi/n Φ₁(x) = x − 1 Φ₂(x) = x + 1 Φ₄(x) = x² + 1 ← roots are ±i (MOST IMPORTANT for crypto) Φ₈(x) = x⁴ + 1 Φ₁₂(x) = x⁴ − x² + 1 Φₚ(x) = x^p−1 + x^p−2 + ⋯ + x + 1 for prime p xⁿ − 1 = ∏ Φd(x) over all divisors d of n

Why Φ₂ₙ(x) = xⁿ + 1 Matters

Ring-LWE and Kyber/Dilithium use the ring R = ℤ_q[x] / (xⁿ + 1) where n is a power of 2. This is the ring of integers in the cyclotomic number field ℚ(ζ_2n). The polynomial xⁿ + 1 = Φ_2n(x) is irreducible over ℚ, making R a valid field extension.

Ring-LWE Ring Structure Rq = ℤ_q[x] / (xⁿ + 1) n = 256 in Kyber, 512/1024 in Dilithium // Elements are polynomials of degree < n with coefficients in ℤ_q // Multiplication is polynomial mult, then reduce mod (xⁿ + 1) // Key: xⁿ ≡ −1 (mod xⁿ+1), so xⁿ⁺¹ ≡ −x, etc. This −1 is exactly the √−1 connection: the cyclotomic ring wraps complex roots into algebra

GNU Octave Compute cyclotomic polynomials; polynomial ring multiplication mod xⁿ+1

%% Cyclotomic Polynomials and Ring Arithmetic mod x^n + 1

function phi = cyclotomic(n)
  % Compute Φ_n(x) using the Möbius inversion formula
  % Φ_n = (x^n - 1) / ∏ Φ_d(x) for d | n, d < n
  phi_table = {}; phi_table{1} = [-1, 1];  % Φ₁ = x-1
  for k = 2:n
    % x^k - 1
    num = zeros(1, k+1); num(1) = -1; num(end) = 1;
    den = [1];
    for d = 1:k-1
      if mod(k,d)==0
        den = conv(den, phi_table{d});
      end
    end
    [phi_table{k}, ~] = deconv(num, den);
    phi_table{k} = round(phi_table{k});
  end
  phi = phi_table{n};
end

% Display first several cyclotomic polynomials
for n = [1,2,3,4,5,6,8,10,12]
  c = cyclotomic(n);
  printf('Φ_%2d(x) = ', n);
  for k = length(c):-1:1
    deg=k-1; co=c(length(c)-deg);
    if co~=0
      if deg==0; printf('%+d',co);
      elseif deg==1; printf('%+dx',co);
      else printf('%+dx^%d',co,deg); end
    end
  end; printf('\n');
end

%% Polynomial multiplication in ℤ_q[x]/(x^n + 1) — Ring-LWE style
function r = ring_mul(a, b, n, q)
  % Multiply polynomials a, b in Z_q[x]/(x^n + 1)
  r = zeros(1, n);
  for i = 0:n-1; for j = 0:n-1
    deg = i+j;
    sign= 1 - 2*(floor(deg/n));   % +1 or -1: xⁿ ≡ -1
    r(mod(deg,n)+1) = r(mod(deg,n)+1) + sign*a(i+1)*b(j+1);
  end; end
  r = mod(r, q);
end

% Demo in Z_17[x]/(x^4+1): multiply [1,2,3,4] * [5,6,7,8]
a = [1,2,3,4]; b = [5,6,7,8]; q = 17; n = 4;
printf('\nRing mult in Z_17[x]/(x^4+1):\n');
printf('a*b = '); disp(ring_mul(a, b, n, q));

§ 6

Number Theoretic Transform (NTT)

The NTT is the Discrete Fourier Transform over a finite field. Instead of complex roots of unity e^2πi/n, it uses modular roots of unity — integers ω in ℤ_q such that ωⁿ ≡ 1 (mod q). It powers the O(n log n) polynomial multiplication at the heart of Kyber and Dilithium.

NTT vs FFT FFT: X̂[k] = Σ_j=0^n-1 x[j] · e^−2πijk/n NTT: X̂[k] = Σ_j=0^n-1 x[j] · ω^jk (mod q) // Requirement: q must be a prime with q ≡ 1 (mod n) // so that a primitive n-th root of unity ω exists in Z_q Kyber-768: q = 3329, n = 256, ω = 17 (17^((3329-1)/256) = 17^13 mod 3329) Dilithium: q = 8380417, n = 256, ω = 1753

Polynomial mult O(n²)

→

NTT (pointwise in freq domain)

→

O(n log n) !

The key insight is identical to the FFT: a polynomial product in the "time domain" (coefficient space) becomes a pointwise product in the "frequency domain" (NTT domain), and the NTT can be computed in O(n log n) via the Cooley-Tukey butterfly structure — with all arithmetic mod q instead of over ℂ.

Interactive NTT Demo

Input polynomial coefficients (comma-separated, length must be power of 2) Modulus q (prime, q ≡ 1 mod n recommended)

Press "Run NTT" to see the transform

GNU Octave Full NTT + INTT + fast polynomial multiplication (Kyber parameters)

%% Number Theoretic Transform — the imaginary-number-free FFT used in Kyber

function A = ntt(a, q, omega)
  % Cooley-Tukey NTT (iterative, bit-reversal permutation)
  n = length(a);
  A = a(bitrevorder(1:n));  % bit-reversal permutation
  len = 2;
  while len <= n
    w = mod_pow(omega, n/len, q);
    for i = 1:len:n
      wn = 1;
      for j = 0:len/2-1
        u = A(i+j);
        t = mod(wn * A(i+j+len/2), q);
        A(i+j)          = mod(u+t, q);
        A(i+j+len/2) = mod(u-t+q, q);
        wn = mod(wn*w, q);
      end
    end
    len *= 2;
  end
end

function a = intt(A, q, omega)
  % Inverse NTT: use ω^{-1} and multiply by n^{-1}
  n     = length(A);
  omega_inv = mod_pow(omega, q-2, q);  % ω^{-1} = ω^{q-2} mod q (Fermat)
  n_inv     = mod_pow(n, q-2, q);
  a     = ntt(A, q, omega_inv);
  a     = mod(a * n_inv, q);
end

% NTT-based polynomial multiplication in Z_q (q=17, n=8, ω=2)
% Check: 2^8 = 256 ≡ 1 (mod 17) ← 2 is a primitive 8th root of unity mod 17
q = 17; n = 8; omega = 2;
printf('Verify ω^n mod q: 2^8 mod 17 = %d  (must be 1)\n', mod_pow(omega,n,q));

a = [1,2,3,4,0,0,0,0];  % polynomial 1+2x+3x²+4x³
b = [5,6,7,8,0,0,0,0];  % polynomial 5+6x+7x²+8x³

% Method 1: direct convolution (O(n²))
direct = mod(conv(a,b), q);

% Method 2: NTT-based (O(n log n))
A   = ntt(a, q, omega);
B   = ntt(b, q, omega);
C   = mod(A .* B, q);           % pointwise multiply!
ntt_result = intt(C, q, omega);

printf('Direct:  '); disp(direct(1:8));
printf('NTT:     '); disp(ntt_result);
printf('Match: %d\n', isequal(direct(1:8), ntt_result));

§ 7

Complex Multiplication on Elliptic Curves

An elliptic curve E over a field K has an endomorphism ring End(E) — the set of algebraic maps from E to itself. For most curves End(E) ≅ ℤ. But special curves have End(E) ≅ an order in an imaginary quadratic field K = ℚ(√−d). These are called CM curves.

Imaginary Quadratic Fields Used in Crypto ℚ(√−1) → ℤ[i] endomorphism ring // secp256k1-style curves, j-invariant 1728 ℚ(√−2) → ℤ[√−2] ℚ(√−3) → ℤ[ζ₃] (Eisenstein integers) // j-invariant 0, BLS12-381 curve ℚ(√−7) → order in ℤ[(1+√−7)/2] // CM theory gives control over the group order #E(𝔽_p): #E(𝔽_p) = p + 1 − t, where t = π + π̄ (trace of Frobenius) π · π̄ = p and π is a Gaussian/quadratic integer in ℚ(√−d)

Why This Matters: Pairing-Friendly Curves

Pairing-based cryptography (BLS signatures, zkSNARKs) requires curves with very specific group orders and embedding degrees. The CM method constructs these curves by choosing a desired order, finding the CM discriminant D, computing the Hilbert class polynomial H_D(x), and setting j(E) to be a root of H_D(x) mod p.

CM Method — Constructing Curves with Prescribed Order 1. Choose order m = p+1−t with t²−4p = D·y² (D < 0 square-free) 2. Compute Hilbert class polynomial H_D(x) over ℤ 3. Find root j₀ of H_D(x) mod p → j₀ = j-invariant of desired curve 4. Recover a, b from j₀: a = 3j₀/(j₀−1728), b = 2j₀/(j₀−1728) 5. Check #E(𝔽_p) = m; twist if needed

GNU Octave Frobenius trace, group order, and j-invariant for CM curves

%% CM Curves: Counting Points with Complex Multiplication
%% For curves y² = x³ + ax + b over F_p, we use Schoof's baby approach

function count = count_points_naive(a, b, p)
  % Count affine points on y² = x³ + ax + b mod p, plus point at infinity
  count = 1;   % include point at infinity
  for x = 0:p-1
    rhs = mod(x^3 + a*x + b, p);
    if rhs == 0
      count += 1;
    elseif mod_pow(rhs, (p-1)/2, p) == 1   % QR: Euler's criterion
      count += 2;
    end
  end
end

function j = j_invariant(a, b, p)
  % j(E) = 1728 · 4a³ / (4a³ + 27b²) mod p
  num = mod(1728 * 4 * mod_pow(a, 3, p), p);
  den = mod(4*mod_pow(a,3,p) + 27*mod_pow(b,2,p), p);
  j   = mod(num * mod_pow(den, p-2, p), p);
end

p = 127;   % small prime for demo
printf('Curve  a    b     #E    t=p+1-#E  j-invariant\n');
printf('%s\n', repmat('-',52,1));
for [a,b] = [[1;0], [0;1], [-1;0], [3;5], [0;7]]
  ord = count_points_naive(a, b, p);
  t   = p + 1 - ord;
  j   = j_invariant(mod(a,p), mod(b,p), p);
  printf('y²=x³%+dx%+d  %-5d %+d   j=%d\n', a,b,ord,t,j);
end
%% j=0:    curve with CM by ℤ[ω₃] (cube roots of unity, BLS-type)
%% j=1728: curve with CM by ℤ[i]  (square roots of -1)

§ 8

GLV Endomorphism: Using √−1 to Speed Up ECC

The Gallant-Lambert-Vanstone (GLV) method exploits a curve's endomorphism to split an n-bit scalar multiplication into two n/2-bit multiplications, roughly halving computation time. For secp256k1 (Bitcoin/Ethereum), this endomorphism comes from a complex cube root of unity in the base field.

secp256k1 GLV Endomorphism // Curve: y² = x³ + 7 over 𝔽_p (j = 0, CM by ℤ[ω₃]) φ: (x, y) ↦ (β·x, y) where β³ ≡ 1 (mod p), β ≠ 1 // β is a non-trivial cube root of unity in 𝔽_p (not √-1, but related) // Since φ is an endomorphism: φ(P) = λ·P for some λ with λ² + λ + 1 ≡ 0 (mod n) // For any scalar k, decompose: k = k₁ + k₂·λ (both k₁, k₂ ≈ n/2 bits) k·P = k₁·P + k₂·φ(P) = k₁·P + k₂·(β·x, y) // Use Simultaneous Double-and-Add: halves the number of doublings

The √−1 Connection for Other Curves

For curves where p ≡ 1 (mod 4) and End(E) ≅ ℤ[i], the endomorphism is literally (x, y) ↦ (−x, i·y) where i = √−1 mod p. This gives φ(P) = i·P (the imaginary unit acts as a scalar!).

GNU Octave GLV scalar decomposition on a CM curve with ℤ[i] endomorphism

%% GLV Scalar Decomposition using the ℤ[i] endomorphism
%% Curve: y² = x³ + x  over F_p (j=1728, CM by ℤ[i])
%% Endomorphism: φ(x,y) = (-x, i·y) where i² ≡ -1 (mod p)

% Toy parameters: p=9001 (prime ≡ 1 mod 4 so sqrt(-1) exists)
p   = 9001;
a_c = 1; b_c = 0;   % y² = x³ + x
printf('p mod 4 = %d\n', mod(p, 4));  % must be 1

% Compute sqrt(-1) mod p
i_mod_p = sqrt_neg1_mod_p(p);
printf('sqrt(-1) mod %d = %d  (verify: %d)\n', p, i_mod_p, ...
       mod(i_mod_p^2, p));  % should be p-1

% ℤ[i] endomorphism: φ(x,y) = (-x mod p, i_mod_p * y mod p)
function Q = phi_endo(P, i_val, p)
  Q = [mod(-P(1), p),  mod(i_val * P(2), p)];
end

% GLV decomposition: given k, find k1, k2 s.t. k ≡ k1 + k2*λ (mod n)
% where λ is the eigenvalue of φ (λ ≡ i_mod_p for ℤ[i] CM)
function [k1,k2] = glv_decompose(k, lambda, n)
  % Simple decomposition (exact for ℤ[i] case, round-based)
  k2 = round(k * mod_inverse(lambda,n) / n);
  k1 = k - k2 * lambda;
end

% Find a point on y² = x³ + x mod 9001
G = [];
for x = 1:p
  rhs = mod(x^3 + x, p);
  y   = mod_pow(rhs, floor((p+1)/4), p);
  if mod(y^2,p)==rhs && y>0; G=[x y]; break; end
end
phiG = phi_endo(G, i_mod_p, p);
printf('G       = (%d, %d)\n', G(1), G(2));
printf('φ(G)    = (%d, %d)  [= (-x, i·y)]\n', phiG(1), phiG(2));
printf('φ(G) lies on curve: y²-x³-x ≡ %d (mod %d)\n', ...
       mod(phiG(2)^2 - phiG(1)^3 - phiG(1), p), p);  % should be 0

§ 9

Bilinear Pairings & Roots of Unity

A bilinear pairing maps two groups of elliptic curve points to a group of roots of unity in a field extension. This is where complex roots of unity re-emerge, now living in a finite field extension 𝔽_p^k.

Weil / Tate Pairing e: G₁ × G₂ → μ_r ⊂ 𝔽_p^k* μ_r = { ζ ∈ 𝔽_p^k* : ζ^r = 1 } // r-th roots of unity! Bilinearity: e(aP, bQ) = e(P, Q)^ab Non-degeneracy: e(P, Q) ≠ 1 for P,Q ≠ ∞ // k = embedding degree: smallest k s.t. r | (p^k − 1) // BLS12-381: k=12, r≈255-bit prime, p≈381-bit prime

Applications Enabled by Pairings

BLS Signatures

Sign: σ = x·H(m). Verify: check e(σ, G) = e(H(m), xG). Supports signature aggregation — combine n signatures into one. Used in Ethereum 2.0.

Identity-Based Encryption

Encrypt to an email address without a pre-distributed key. Uses pairings to bind identity strings to public keys. Boneh-Franklin IBE (2001) was the first practical IBE scheme.

zkSNARKs (Groth16)

Zero-knowledge proofs with constant-size proofs (192 bytes). The pairing check e(A, B) = e(C, δ)·e(D, γ) is the final verification step. Powers Zcash and many ZK rollups.

Attribute-Based Encryption

Encrypt data with a policy ("CEO AND CTO"). Decryption works only if attributes satisfy the policy. Pairings make the policy matching algebraically verifiable.

Miller's Algorithm — Computing the Weil Pairing // Core loop: builds a function f_{r,P} on the curve f = 1; T = P for bit in binary(r) from MSB-1 to 0: f ← f² · l_T,T(Q) / v_2T(Q) // line function at doubling T ← 2T if bit = 1: f ← f · l_T,P(Q) / v_T+P(Q) // line function at addition T ← T+P e(P,Q) = f^{(p^k − 1)/r} (final exponentiation to μ_r)

§ 10

Ring-LWE: Imaginary Quadratic Fields Meet Lattices

Ring-LWE (Ring Learning With Errors) is the algebraic number theory problem underlying Kyber and Dilithium. The ring R = ℤ[x]/(xⁿ+1) is the ring of integers of the cyclotomic field ℚ(ζ_2n), and its properties — rooted in the arithmetic of complex roots of unity — provide both the security and efficiency of post-quantum lattice schemes.

Ring-LWE Problem R = ℤ[x]/(xⁿ+1), R_q = R/(q) = ℤ_q[x]/(xⁿ+1) Secret: s ← R_q (small coefficients) Sample: (a, b = a·s + e) where a ← R_q uniform, e small // Ring-LWE Hard Problem: distinguish (a, as+e) from (a, uniform) // Why xⁿ+1? It's irreducible over ℤ → field extension // Its roots are ζ, ζ³, ζ⁵, ... (primitive 2n-th roots of unity) // These complex roots guarantee algebraic structure needed for proofs

The Canonical Embedding

The canonical embedding maps a ring element f ∈ R to n complex numbers by evaluating at each of the primitive 2n-th roots of unity. This embedding is the bridge between the abstract ring and complex analysis — it is what makes the geometric notion of "small" (short vectors in ℂⁿ) correspond to "small coefficients" in the polynomial ring.

Canonical Embedding σ: R → ℂⁿ σ(f) = ( f(ζ), f(ζ³), f(ζ⁵), …, f(ζ²ⁿ⁻¹) ) where ζ = e^iπ/n (primitive 2n-th root of unity) ‖σ(f)‖² = n · Σ fᵢ² // connects ℂⁿ norm to polynomial coefficients The security of Ring-LWE ultimately reduces to the hardness of finding short vectors in lattices defined by these complex evaluations.

GNU Octave Canonical embedding, Ring-LWE key exchange, visualizing complex roots

%% Ring-LWE Key Exchange: canonical embedding and Kyber-style demo

function emb = canonical_embedding(f, n)
  % Evaluate f at primitive 2n-th roots of unity ζ^(2k+1), k=0..n-1
  % ζ = exp(i*π/n)
  emb = zeros(1, n);
  for k = 0:n-1
    zeta  = exp(1i * pi * (2*k+1) / n);  % ζ^(2k+1)
    val   = 0;
    for j = 0:n-1
      val += f(j+1) * zeta^j;
    end
    emb(k+1) = val;
  end
end

% Show that xⁿ+1 roots are on the unit circle in C
n = 8;
printf('Roots of x^%d + 1 (primitive 2n-th roots of unity):\n', n);
for k = 0:n-1
  zeta = exp(1i * pi * (2*k+1) / n);
  printf('  k=%d: angle=%.2f°  ζ^%d + 1 = %.6f\n', ...
         k, (2*k+1)*180/n, n, abs(zeta^n + 1));  % should be ~0
end

% Ring-LWE toy key exchange in Z_q[x]/(x^n+1)
q = 97; n = 8;
% Alice's key pair
s_a = randi([-1,1], 1, n);          % small secret
a   = randi([0, q-1], 1, n);        % public random polynomial
e_a = randi([-1,1], 1, n);          % small error
b_a = ring_mul(a, s_a, n, q);
b_a = mod(b_a + e_a, q);             % public key: b_a = a*s_a + e_a
% Bob's key pair
s_b = randi([-1,1], 1, n);
e_b = randi([-1,1], 1, n);
b_b = mod(ring_mul(a, s_b, n, q) + e_b, q);
% Shared secrets (approximately equal due to small errors)
K_a = ring_mul(b_b, s_a, n, q);    % Alice computes b_b * s_a
K_b = ring_mul(b_a, s_b, n, q);    % Bob computes b_a * s_b
diff= max(abs(double(K_a) - double(K_b)));
printf('\nRing-LWE key exchange:\n');
printf('K_a = '); disp(K_a);
printf('K_b = '); disp(K_b);
printf('Max coefficient difference: %d  (small ≈ error magnitude)\n', diff);
% In real Kyber, reconciliation mechanism ensures exact agreement

The Full Circle

We started with i = √−1 — an abstract concept with "no real meaning." We end with it as the organizing principle of the most important family of post-quantum algorithms. The roots of unity e^2πi/n define the algebraic structure of ℚ(ζ_n). The rings ℤ[x]/(Φ_n(x)) inherited from those roots are precisely the ones that make lattice cryptography simultaneously secure and fast. The imaginary unit is not an abstraction bolted onto cryptography — it is woven into its foundations.