Chapter Nine:

Imaginary Primes and Prime Imaginaries

9.1 The Complex Euclidean Algorithm

ComplexGCD[z_, u_] := Block[{q = z/u}, Print[z];
          ComplexGCD[u, z - u (Round[Re[q]] + I Round[Im[q]])]]
ComplexGCD[z_, 0_] := z

ComplexGCD[73, 27+I]

ComplexGCD[-1350180 + 3266646 I, -4416 + 1040 I]

ComplexGCD[Product[a + b I, {a, 1, 8},   {b,  1,  7}],
           Product[a + b I, {a, 12, 16}, {b, 12, 15}]]

9.2 Quadratic Residues

(* PowerMod must be made Listable, see CRT section in Chapter 8 *)

PowerMod[Range[1, 18], 2, 19]

SetAttributes[JacobiSymbol, Listable]
JacobiSymbol[Range[18], 19]

Flatten[Position[%, 1]]

Select[Prime[Range[2, 50]], JacobiSymbol[-1, #] == 1 &]

Nonresidue[p_?OddQ] := Block[{n = 0},
	While[JacobiSymbol[Prime[++n], p] == 1]; Prime[n]]
Attributes[Nonresidue] = Listable

Nonresidue[Prime[Range[100000, 100030]]]

JacobiSymbol[Prime[Range[10]], Prime[100010]]

SqrtNegOne[p_] := 
		PowerMod[Nonresidue[p], (p-1)/4, p]  /; Mod[p, 4]==1
SqrtNegOne[2] := 1

SqrtNegOne[73]

primes = Select[Prime[Range[100]], Mod[#, 4] == 1 &]
PowerMod[SqrtNegOne /@ primes, 2, primes] - primes

9.3 Sums of Two Squares via Complex GCDs

ComplexGCD with Print deleted

ComplexGCD[z_, u_] := Block[{q = z/u},
          ComplexGCD[u, z - u (Round[Re[q]] + I Round[Im[q]])]]
ComplexGCD[z_, 0_] := z

Sum2SquaresViaComplexGCD[p_] :=
	Block[{g = ComplexGCD[p, SqrtNegOne[p] + I]},
	Abs[{Re[g], Im[g]}]]           /; Mod[p, 4] == 1 || p == 2
Attributes[Sum2SquaresViaComplexGCD] = Listable

Sum2SquaresViaComplexGCD[73]

Sum2SquaresViaComplexGCD[848654483879497562821]

9.4 Gaussian Primes

boxes[q_] := Map[Rectangle[#, #+1] &,
              { q, Reverse[q],   q*{1,-1}, Reverse[q*{1,-1}],
               -q, Reverse[-q], -q*{1,-1}, Reverse[-q*{1,-1}]}]

units = {GrayLevel[.5], boxes[{1,0}]}

GaussianPrimes[n_] := Show[Graphics[{units,
 Map[boxes,
  {{1,1}} ~Join~
  Map[{#,0} &, Select[Range[3, n, 4], PrimeQ]] ~Join~
  Sum2SquaresViaComplexGCD[Select[Range[5,n^2,4], PrimeQ]]]}],
AspectRatio->1]

GaussianPrimes[10];

version 2.0 Notes

the factoring of Gaussian integers has been built into version 2.0 via an option, GaussianIntegers -> True, to the FactorInteger command. Note that this means that the representation of primes as a sum of two squares is essentially built-in: just factor the prime over the Gaussian integers.

(* Uses Sum2SquaresViaComplexGCD, SqrtNegOne, Nonresidue, and   		ComplexGCD *)

divides[u_, z_] := IntegerQ[Re[z/u]] && IntegerQ[Im[z/u]]

GaussianFactors[z_] := Block[
	{norm = Abs[z]^2, plist, factors = {}, w, s, a, b},
	plist = FactorInteger[norm];
	If[EvenQ[norm], AppendTo[factors, {1 + I, plist[[1,2]]}];
		            plist = Rest[plist]];
	Scan[If[Mod[#[[1]], 4] == 3, 
	  AppendTo[factors, {#[[1]], #[[2]]/2}],
		{a, b} = Sum2SquaresViaComplexGCD[#[[1]]];
		w =  a + b I;
		s = 0;
		While[divides[w, z], s++; w *= a + b I];
		If[s > 0, AppendTo[factors, {a + b I, s}]];
		If[s < #[[2]], AppendTo[factors,
								{a - b I, #[[2]] - s}]]] &,
	plist];
factors]

GaussianFactors[1488311206705990 + 8526860574095268 I]

GaussianFactors[21!]

GaussianFactors[270015132441600000 + 4044308486400000 I]

9.5 Sums of Two Squares via Real GCDs

(* Uses SqrtNegOne and Nonresidue *)

Attributes[MiddleRemainder] = Listable
Attributes[Sum2Squares] = Listable

MiddleRemainder[p_, m_, x_] := x /; x^2 < p
MiddleRemainder[p_, m_, x_] :=
           MiddleRemainder[p, x, Mod[m, x]] /; x^2 >= p
MiddleRemainder[p_, x_]     := MiddleRemainder[p, p, x]

Sum2Squares[p_] :=
     MiddleRemainder[p, SqrtNegOne[p]] /; Mod[p, 4] == 1

Sum2Squares[848654483879497562821]

Bonus: New Code

The following version of MiddleRemainders is shorter, more direct, and 20% faster. Recursion is tempting, but often not best!

Attributes[MiddleRemainder1] = Listable

MiddleRemainder1[x_, y_] := Block[{xx = x, yy = y},
 While[xx^2 > x, {xx, yy} = {yy, Mod[xx,yy]}];
 xx]

9.6 Square Roots Modulo an Integer

SqrtModPrime[a_, 2]  := Mod[a, 2]
SqrtModPrime[a_, p_] := {}  /;  JacobiSymbol[a, p] == -1
SqrtModPrime[a_, p_] := 0   /;  Mod[a, p] == 0

SqrtModPrime::usage = "SqrtModPrime[a, p] assumes p is prime and returns the least nonnegative square root of a mod p if such exists, {} otherwise."
Attributes[SqrtModPrime] = Listable

SqrtModPrime[a_, p_] :=
  Min[PowerMod[a, (p+1)/4, p]*{-1,1} + {p,0}] /; Mod[p, 4] == 3

Mod[PowerMod[49, 7, 3329]*PowerMod[3, 702, 3329], 3329]

SqrtModPrime[a_, p_] := (k = (p - 1)/4; s = 0;
  h = Nonresidue[p];
  While[EvenQ[k], k /= 2; s++];
  e2 = p-1;
  Scan[(e2 /= 2; 
     If[Mod[# * PowerMod[h, e2, p], p] != 1, e2 += (p-1)/2]) &,
     Reverse[NestList[Mod[# #, p]&, PowerMod[a, k, p], s]]];
  Min[Mod[PowerMod[a, (k+1)/2, p]*PowerMod[h, e2/2, p], p] * 
	{-1,1} + {p,0}]) 						/; Mod[p, 4] == 1

testprimes = Prime[Range[10000, 10020]]
SqrtModPrime[49, testprimes]

SqrtModPrime[48, testprimes]

Mod[%^2, testprimes]

SqrtModPrimeShanks from Appendix (needed for next routine)

SqrtModPrimeShanks::usage = "SqrtModPrimeShanks[a, p] assumes p is prime and returns, via Shanks's algorithm, the least nonnegative square root of a mod p if such exists, {} otherwise. It requires the routine Nonresidue"

Attributes[SqrtModPrimeShanks] = Listable
SetAttributes[PowerMod, Listable] 

SqrtModPrimeShanks[a_, 2]  := Mod[a, 2]

SqrtModPrimeShanks[a_, p_] := {}             /; JacobiSymbol[a, p] == -1

SqrtModPrimeShanks[a_, p_] := 0                        /; Mod[a, p] == 0

SqrtModPrimeShanks[a_, p_] :=
  Min[PowerMod[a, (p+1)/4, p]*{-1, 1} + {p, 0}]        /; Mod[p, 4] == 3

SqrtModPrimeShanks[a_, p_] :=
  Block[{k = (p - 1)/2, i, L, n, r, s = 1, c, b},
	While[EvenQ[k], k /= 2; s++];
	{c, r, n} = PowerMod[{Nonresidue[p], a, a}, {k, (k+1)/2,k}, p];
	While[n != 1, 
		L = s; b = n;
		Do[
		If[b == 1, b = c; s = i - 1,  b = Mod[b b, p]], {i, L}];
		{c, r, n} = Mod[b {b, r, b n}, p]];
	Min[r, p - r]]                                   /; Mod[p, 4] == 1

Clear[SqrtModPrimePower]


SqrtModPrimePower::usage = "SqrtModPrimePower[a, p, r] returns the set of all mod-p^r square roots of a, in increasing order; it uses SqrtModPrimeShanks, which in turn uses Nonresidue."

Attributes[SqrtModPrimePower] = Listable

(* Exponent = 1 case reduces to previous routine *)
SqrtModPrimePower[a_, p_, 1] :=
Block[{root = SqrtModPrimeShanks[a, p]},	
If[!NumberQ[root], {}, Union[Mod[{p, 0} + {-1, 1} root, p]]]]

(* Non-coprime case, any prime, reduces to coprime case *)
SqrtModPrimePower[a_, p_, r_] := Block[{s = 1, j},
 If[Mod[a, p^r] == 0, Return[Range[0, p^r-1, p^Ceiling[r/2]]]];
 While[Mod[a, p^s] == 0, s++]; s--;
 If[OddQ[s], {},
    Flatten@Table[
    (SqrtModPrimePower[a/p^s, p, r-s] + (j-1) p^(r-s)) p^(s/2),
		       {j, p^(s/2)}]]]      /; Mod[a, p] == 0 && r > 1

(* Coprime case, odd prime *)
SqrtModPrimePower[a_, p_, r_] :=
 Block[{j, root = SqrtModPrimeShanks[a, p]},
 If[!NumberQ[root], Return[{}], 
 Do[root = Mod[root - PowerMod[2 root, -1, p]*(root^2-a), p^j],
								{j, 2, r}];
 Return[Sort[{-1,1}*root + {p^r, 0}]]]] /;
			p != 2 && Mod[a, p] != 0 && r > 1

(* Coprime case, prime is 2, no solution *)
SqrtModPrimePower[a_, 2, r_] := {}   /;
(r==2 && Mod[a, 4] == 3) || (3<=r && OddQ[a] && Mod[a, 8] != 1)

(* Coprime case modulo 4 *)
SqrtModPrimePower[a_, 2, 2] := {1, 3} /; Mod[a, 4] == 1

(* Coprime case, prime is 2, exponent is at least 3, solutions ÊÊÊÊexist *)
SqrtModPrimePower[a_, 2, r_] :=
 Block[{root = Nest[Mod[(#^3 + (2-a)#)/2, 2^r]&, 1, r-3]},
 Union[
    {root, 2^r-root}, Mod[{root, 2^r-root} + 2^(r-1), 2^r]]] /;
  										OddQ[a] && 3 <= r

bruteforce[a_, n_] := Select[Range[0, n-1], 
				Mod[#^2, n] == Mod[a, n] &]
Attributes[bruteforce] = Listable

a = 16; p = 2; r = 5
{SqrtModPrimePower[a, p, r], bruteforce[a, p^r]}

a = Range[32]; p = 2; r = 5
SqrtModPrimePower[a, p, r] == bruteforce[a, p^r]

SqrtMod::usage = "SqrtMod[a, n] returns the set of all mod-n square roots of a, in increasing order; it uses CRT, SqrtModPrimePower, SqrtModPrimeShanks, and Nonresidue.";
Attributes[SqrtMod] = Listable

SqrtMod[a_, n_] := Block[{factors = Transpose@FactorInteger[n]},
  rootlist = SqrtModPrimePower[a, factors[[1]], factors[[2]]];
  If[MemberQ[rootlist, {}], Return[{}]];
  PrependTo[rootlist, 
  	CRT[List[##], factors[[1]]^factors[[2]]]&];
  Union@Flatten[Outer @@ rootlist, Length[rootlist] - 2]]

{SqrtMod[1, 3*5*7], bruteforce[1, 3*5*7]}

SqrtMod[4, 234535662466]

Mod[%^2, 234535662466]

9.7 A General Solution to the Sum-of-Two-Squares Problem

PrimitiveReps::usage = "PrimitiveReps[n] returns a list of all pairs of nonnegative integers {a, b} such that a >= b,
n = a^2 + b^2, and GCD[a, b] = 1."

SqrtsSquareDivisors::usage = "SqrtsSquareDivisors[n] returns a list of all positive integers a whose square divides n"

GeneralSum2Squares::usage = "GeneralSum2Squares[n] returns a list of all pairs of nonnegative integers {a, b} such that
a >= b and n = a^2 + b^2"

Attributes[PrimitiveReps] = Listable
Attributes[GeneralSum2Squares] = Listable

PrimitiveReps[n_] := Map[{#, Sqrt[n-#^2]}&,
MiddleRemainder[n, Select[SqrtMod[-1,n], # <= n/2 &]]] /; n > 1

PrimitiveReps[1] = {{1,0}};					

SqrtsSquareDivisors[n_] :=
  Block[{factorlist = Transpose@FactorInteger[n]},
    Flatten [Outer @@ Prepend[
	 factorlist[[1]]^(Range[0, factorlist[[2]]/2]), Times]]]

GeneralSum2Squares[1] = {{1, 0}}
GeneralSum2Squares[0] = {{0, 0}}

GeneralSum2Squares[n_] :=
 Block[{temp = SqrtsSquareDivisors[n]},
 Sort[Complement[
  Flatten[PrimitiveReps[n/temp^2] temp, 1], {{}}]]] /; n > 1

PrimitiveReps[Prime[Range[10000, 10005]]]

SqrtsSquareDivisors[2 3^2 5^2 23]

n = 2454^2 + 4
GeneralSum2Squares[n]

Map[Plus@@# &, %^2]

GeneralSum2Squares[10^10]

PrimitiveReps from Appendix

PrimitiveReps::usage = "PrimitiveReps[n, f, g] returns the set of representations of n as f a^2 + g b^2, where gcd(a, b) = 1. It is assumed that gcd(n, f g) = 1 and n > f + g. The routines SqrtMod, Nonresidue, and CRT must be loaded."

MiddleRemainder::usage = "MiddleRemainder[n, x, f, g] returns the first remainder r that satisfies f r^2 < n, when the Euclidean algorithm is applied to n, x."

SetAttributes[{MiddleRemainder, PrimitiveReps}, Listable]

MiddleRemainder[n_, m_, x_, f_, g_] := x /; f x^2 < n
MiddleRemainder[n_, m_, x_, f_, g_] :=
           MiddleRemainder[n, x, Mod[m, x], f, g] /; f x^2 >= n
MiddleRemainder[n_, x_, f_, g_] := MiddleRemainder[n, n, x, f, g]

PrimitiveReps[n_, f_, g_] := Block[{temp},
   Complement[
      Map[If[IntegerQ[temp = Sqrt[(n-f #^2)/g]], {#, temp}, {}] &,
         MiddleRemainder[n,
          Select[SqrtMod[-g PowerMod[f, -1, n], n], # <= n/2 &], f, g]],
  {{}}]]                                                    /; n > f + g

PrimitiveReps[13904, 1, 7]

BONUS on Lagrange's Theorem on 4 squares

Not in book, but here is my implementation of a routine from the paper [RS] that represents any integer as a sum of 4 squares. It could probably be improved by first checking if the input is itself a square or sum of two squares.

(* Needs Sum2Squares from chapter 9 *)

Sum4Squares[n_] := Block[{x=0,y=0,aa},
	While[!PrimeQ[n-x^2-y^2],
 x = Random[Integer, {0,Floor@Sqrt[n]//N}];
 y = Random[Integer, {0,Floor@Sqrt[n-x^2]}//N]];
 Sort[{x,y,aa=Sum2Squares[n-x^2-y^2],
             Sqrt[n-x^2-y^2-aa^2]},
                Mod[#1,2] < Mod[#2,2]&]] /;
          EvenQ[n] && OddQ[n/2]

Sum2Squares[2] = 1

Sum4Squares[n_?OddQ] :=
Block[{ans = Sum4Squares[2n]},
 Abs[{ans[[1]]-ans[[2]], ans[[1]]+ans[[2]],
      ans[[3]]+ans[[4]], ans[[3]]-ans[[4]]}/2]]

Sum4Squares[n_?EvenQ] :=  Block[{d=0,ans},
 While[EvenQ[n/2^d], d++];
 If[OddQ[d], 2^((d-1)/2) * Sum4Squares[2 n/2^d],
 	ans = Sum4Squares[2n];
 	Abs[{ans[[1]]-ans[[2]], ans[[1]]+ans[[2]],
      ans[[3]]+ans[[4]], ans[[3]]-ans[[4]]}/2]]]

Attributes[Sum4Squares]=Listable

Sum4Squares[10^10+17]

Plus @@ (%^2)

9.8 Eisenstein Primes

Representation[p_] :=
 ({u = MiddleRemainder[p, SqrtModPrimeShanks[-3, p]],
   Sqrt[(p - u^2)/3]}) 					/; Mod[p, 3] == 1
Attributes[Representation] = Listable

Representation[{7, 13, 19, 31, 37, 43, 1299709}]

Map[#[[1]]^2 + 3 #[[2]]^2 &, %]

EisensteinPrimes::usage = "EisensteinPrimes[n] uses hexagons to generate an image of the Eisenstein primes with norm less than n^2."

SetAttributes[{coords, factor, associates, hexagon}, Listable]

omega = N[Exp[2 I Pi /3]];
units = {1, -omega^2, omega, -1, omega^2, -omega};
hexagon[0] = N[Exp[I Range[Pi/6, 2 Pi, Pi/3]] / Sqrt[3]]
hexagon[z_] := z + hexagon[0];
coords[z_] := {Re[z], Im[z]}; 
associates[z_] := coords[hexagon[units * z]]

factor[p_] :=
 ({u, v} = Representation[p]; u + v + 2 v omega^{1,2}) 

EisensteinPrimes[n_] := (
 realprimes = Flatten[
	associates[Append[Select[Range[5, n, 6], PrimeQ], 2]], 1];
 nonrealprimes = Flatten[associates[factor[
	Select[Range[7, n^2, 6], PrimeQ]]], 2];
 Show[Graphics[{Polygon /@ associates[1 - omega],
                Polygon /@ realprimes,
                Polygon /@ nonrealprimes,
                GrayLevel[.5],
                Polygon /@ associates[1]}],	AspectRatio->1])

EisensteinPrimes[10];