(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "NeXT Mathematica Notebook Front End Version 2.2"; NeXTStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L3, e8, 24, "Times"; ; fontset = subtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; ; fontset = subsubtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 14, "Times"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 14, "Times"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 12, "Times"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 10, "Times"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 14, "Times"; ; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 10, "Courier"; ; fontset = name, inactive, noPageBreakInGroup, nowordwrap, nohscroll, preserveAspect, M7, italic, B65535, 10, "Times"; ; fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = leftheader, 12; fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, 12; fontset = leftfooter, 12; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12; fontset = completions, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12; fontset = special1, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12; fontset = special2, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, 12; fontset = special3, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, right, M7, 12; fontset = special4, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12; fontset = special5, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12; paletteColors = 128; currentKernel; ] :[font = title; inactive; preserveAspect] Chapter Nine: Imaginary Primes and Prime Imaginaries :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.1 The Complex Euclidean Algorithm :[font = input; preserveAspect] ComplexGCD[z_, u_] := Block[{q = z/u}, Print[z]; ComplexGCD[u, z - u (Round[Re[q]] + I Round[Im[q]])]] ComplexGCD[z_, 0_] := z :[font = input; Cclosed; preserveAspect; startGroup] ComplexGCD[73, 27+I] :[font = print; inactive; preserveAspect] 73 27 + I :[font = output; output; inactive; preserveAspect; endGroup] -8 - 3*I ;[o] -8 - 3 I :[font = input; Cclosed; preserveAspect; startGroup] ComplexGCD[-1350180 + 3266646 I, -4416 + 1040 I] :[font = print; inactive; preserveAspect] -1350180 + 3266646 I -4416 + 1040 I 780 - 1882 I 128 + 718 I -400 - 190 I -82 + 128 I -16 + 56 I 6 + 32 I :[font = output; output; inactive; preserveAspect; endGroup] 4 - 14*I ;[o] 4 - 14 I :[font = input; Cclosed; preserveAspect; startGroup] ComplexGCD[Product[a + b I, {a, 1, 8}, {b, 1, 7}], Product[a + b I, {a, 12, 16}, {b, 12, 15}]] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] 270015132441600000 + 4044308486400000*I ;[o] 270015132441600000 + 4044308486400000 I :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.2 Quadratic Residues :[font = input; Cclosed; preserveAspect; startGroup] (* PowerMod must be made Listable, see CRT section in Chapter 8 *) PowerMod[Range[1, 18], 2, 19] :[font = output; output; inactive; preserveAspect; endGroup] {1, 4, 9, 16, 6, 17, 11, 7, 5, 5, 7, 11, 17, 6, 16, 9, 4, 1} ;[o] {1, 4, 9, 16, 6, 17, 11, 7, 5, 5, 7, 11, 17, 6, 16, 9, 4, 1} :[font = input; Cclosed; preserveAspect; startGroup] SetAttributes[JacobiSymbol, Listable] JacobiSymbol[Range[18], 19] :[font = output; output; inactive; preserveAspect; endGroup] {1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, -1, 1, 1, -1} ;[o] {1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, -1, 1, 1, -1} :[font = input; Cclosed; preserveAspect; startGroup] Flatten[Position[%, 1]] :[font = output; output; inactive; preserveAspect; endGroup] {1, 4, 5, 6, 7, 9, 11, 16, 17} ;[o] {1, 4, 5, 6, 7, 9, 11, 16, 17} :[font = input; Cclosed; preserveAspect; startGroup] Select[Prime[Range[2, 50]], JacobiSymbol[-1, #] == 1 &] :[font = output; output; inactive; preserveAspect; endGroup] {5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137, 149, 157, 173, 181, 193, 197, 229} ;[o] {5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137, 149, 157, 173, 181, 193, 197, 229} :[font = input; initialization; preserveAspect] *) Nonresidue[p_?OddQ] := Block[{n = 0}, While[JacobiSymbol[Prime[++n], p] == 1]; Prime[n]] Attributes[Nonresidue] = Listable (* :[font = input; Cclosed; preserveAspect; startGroup] Nonresidue[Prime[Range[100000, 100030]]] :[font = output; output; inactive; preserveAspect; endGroup] {2, 7, 5, 2, 11, 2, 5, 2, 2, 3, 19, 2, 2, 2, 5, 2, 2, 3, 2, 3, 2, 2, 2, 2, 13, 2, 3, 3, 3, 5, 7} ;[o] {2, 7, 5, 2, 11, 2, 5, 2, 2, 3, 19, 2, 2, 2, 5, 2, 2, 3, 2, 3, 2, 2, 2, 2, 13, 2, 3, 3, 3, 5, 7} :[font = input; Cclosed; preserveAspect; startGroup] JacobiSymbol[Prime[Range[10]], Prime[100010]] :[font = output; output; inactive; preserveAspect; endGroup] {1, 1, 1, 1, 1, 1, 1, -1, -1, -1} ;[o] {1, 1, 1, 1, 1, 1, 1, -1, -1, -1} :[font = input; initialization; preserveAspect] *) SqrtNegOne[p_] := PowerMod[Nonresidue[p], (p-1)/4, p] /; Mod[p, 4]==1 SqrtNegOne[2] := 1 (* :[font = input; Cclosed; preserveAspect; startGroup] SqrtNegOne[73] :[font = output; output; inactive; preserveAspect; endGroup] 27 ;[o] 27 :[font = input; Cclosed; preserveAspect; startGroup] primes = Select[Prime[Range[100]], Mod[#, 4] == 1 &] PowerMod[SqrtNegOne /@ primes, 2, primes] - primes :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1} ;[o] {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1} :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.3 Sums of Two Squares via Complex GCDs :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] ComplexGCD with Print deleted :[font = input; preserveAspect; endGroup] ComplexGCD[z_, u_] := Block[{q = z/u}, ComplexGCD[u, z - u (Round[Re[q]] + I Round[Im[q]])]] ComplexGCD[z_, 0_] := z :[font = input; preserveAspect] Sum2SquaresViaComplexGCD[p_] := Block[{g = ComplexGCD[p, SqrtNegOne[p] + I]}, Abs[{Re[g], Im[g]}]] /; Mod[p, 4] == 1 || p == 2 Attributes[Sum2SquaresViaComplexGCD] = Listable :[font = input; Cclosed; preserveAspect; startGroup] Sum2SquaresViaComplexGCD[73] :[font = output; output; inactive; preserveAspect; endGroup] {8, 3} ;[o] {8, 3} :[font = input; Cclosed; preserveAspect; startGroup] Sum2SquaresViaComplexGCD[848654483879497562821] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {28440994650, 6305894639} ;[o] {28440994650, 6305894639} :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.4 Gaussian Primes :[font = input; preserveAspect] 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] :[font = input; preserveAspect] GaussianPrimes[10]; :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] version 2.0 Notes :[font = text; inactive; preserveAspect; endGroup] 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. ;[s] 5:0,0;82,1;106,2;115,3;128,4;292,-1; 5:1,13,10,Times,0,14,0,0,0;1,11,9,Courier,1,14,0,0,0;1,13,10,Times,0,14,0,0,0;1,11,9,Courier,1,14,0,0,0;1,13,10,Times,0,14,0,0,0; :[font = input; preserveAspect] (* 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] :[font = input; Cclosed; preserveAspect; startGroup] GaussianFactors[1488311206705990 + 8526860574095268 I] :[font = output; output; inactive; preserveAspect; endGroup] {{1 + I, 2}, {6 + I, 11}, {6 - I, 3}, {40 + 17*I, 1}, {994 - 335*I, 1}} ;[o] {{1 + I, 2}, {6 + I, 11}, {6 - I, 3}, {40 + 17 I, 1}, {994 - 335 I, 1}} :[font = input; Cclosed; preserveAspect; startGroup] GaussianFactors[21!] :[font = output; output; inactive; preserveAspect; endGroup] {{1 + I, 36}, {3, 9}, {2 + I, 4}, {2 - I, 4}, {7, 3}, {11, 1}, {3 + 2*I, 1}, {3 - 2*I, 1}, {4 + I, 1}, {4 - I, 1}, {19, 1}} ;[o] {{1 + I, 36}, {3, 9}, {2 + I, 4}, {2 - I, 4}, {7, 3}, {11, 1}, {3 + 2 I, 1}, {3 - 2 I, 1}, {4 + I, 1}, {4 - I, 1}, {19, 1}} :[font = input; Cclosed; preserveAspect; startGroup] GaussianFactors[270015132441600000 + 4044308486400000 I] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {{1 + I, 22}, {3, 4}, {2 + I, 5}, {2 - I, 5}, {7, 1}, {3 + 2*I, 2}, {3 - 2*I, 1}, {4 + I, 1}, {4 - I, 2}, {6 + I, 1}, {4 + 5*I, 1}, {4 - 5*I, 1}, {8 + 3*I, 1}, {8 + 7*I, 1}} ;[o] {{1 + I, 22}, {3, 4}, {2 + I, 5}, {2 - I, 5}, {7, 1}, {3 + 2 I, 2}, {3 - 2 I, 1}, {4 + I, 1}, {4 - I, 2}, {6 + I, 1}, {4 + 5 I, 1}, {4 - 5 I, 1}, {8 + 3 I, 1}, {8 + 7 I, 1}} :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.5 Sums of Two Squares via Real GCDs :[font = input; preserveAspect] (* 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 :[font = input; Cclosed; preserveAspect; startGroup] Sum2Squares[848654483879497562821] :[font = output; output; inactive; preserveAspect; endGroup] 28440994650 ;[o] 28440994650 :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] Bonus: New Code :[font = text; inactive; preserveAspect] The following version of MiddleRemainders is shorter, more direct, and 20% faster. Recursion is tempting, but often not best! :[font = input; preserveAspect; endGroup; endGroup] Attributes[MiddleRemainder1] = Listable MiddleRemainder1[x_, y_] := Block[{xx = x, yy = y}, While[xx^2 > x, {xx, yy} = {yy, Mod[xx,yy]}]; xx] :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.6 Square Roots Modulo an Integer :[font = input; preserveAspect] SqrtModPrime[a_, 2] := Mod[a, 2] SqrtModPrime[a_, p_] := {} /; JacobiSymbol[a, p] == -1 SqrtModPrime[a_, p_] := 0 /; Mod[a, p] == 0 :[font = input; preserveAspect] 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 :[font = input; preserveAspect] SqrtModPrime[a_, p_] := Min[PowerMod[a, (p+1)/4, p]*{-1,1} + {p,0}] /; Mod[p, 4] == 3 :[font = input; Cclosed; preserveAspect; startGroup] Mod[PowerMod[49, 7, 3329]*PowerMod[3, 702, 3329], 3329] :[font = output; output; inactive; preserveAspect; endGroup] 3322 ;[o] 3322 :[font = input; preserveAspect] 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 :[font = input; Cclosed; preserveAspect; startGroup] testprimes = Prime[Range[10000, 10020]] SqrtModPrime[49, testprimes] :[font = output; output; inactive; preserveAspect; endGroup] {7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7} ;[o] {7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7} :[font = input; Cclosed; preserveAspect; startGroup] SqrtModPrime[48, testprimes] :[font = output; output; inactive; preserveAspect; endGroup] {{}, {}, 30153, 10443, 36940, {}, {}, {}, {}, {}, 43550, {}, {}, 10977, 28462, 561, {}, 7126, {}, {}, 39920} ;[o] {{}, {}, 30153, 10443, 36940, {}, {}, {}, {}, {}, 43550, {}, {}, 10977, 28462, 561, {}, 7126, {}, {}, 39920} :[font = input; Cclosed; preserveAspect; startGroup] Mod[%^2, testprimes] :[font = output; output; inactive; preserveAspect; endGroup] {{}, {}, 48, 48, 48, {}, {}, {}, {}, {}, 48, {}, {}, 48, 48, 48, {}, 48, {}, {}, 48} ;[o] {{}, {}, 48, 48, 48, {}, {}, {}, {}, {}, 48, {}, {}, 48, 48, 48, {}, 48, {}, {}, 48} :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] SqrtModPrimeShanks from Appendix (needed for next routine) :[font = input; preserveAspect; endGroup] 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 :[font = input; preserveAspect] 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 :[font = input; preserveAspect] bruteforce[a_, n_] := Select[Range[0, n-1], Mod[#^2, n] == Mod[a, n] &] Attributes[bruteforce] = Listable :[font = input; Cclosed; preserveAspect; startGroup] a = 16; p = 2; r = 5 {SqrtModPrimePower[a, p, r], bruteforce[a, p^r]} :[font = output; output; inactive; preserveAspect; endGroup] {{4, 12, 20, 28}, {4, 12, 20, 28}} ;[o] {{4, 12, 20, 28}, {4, 12, 20, 28}} :[font = input; Cclosed; preserveAspect; startGroup] a = Range[32]; p = 2; r = 5 SqrtModPrimePower[a, p, r] == bruteforce[a, p^r] :[font = output; output; inactive; preserveAspect; endGroup] True ;[o] True :[font = input; preserveAspect] 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]] :[font = input; Cclosed; preserveAspect; startGroup] {SqrtMod[1, 3*5*7], bruteforce[1, 3*5*7]} :[font = output; output; inactive; preserveAspect; endGroup] {{1, 29, 34, 41, 64, 71, 76, 104}, {1, 29, 34, 41, 64, 71, 76, 104}} ;[o] {{1, 29, 34, 41, 64, 71, 76, 104}, {1, 29, 34, 41, 64, 71, 76, 104}} :[font = input; Cclosed; preserveAspect; startGroup] SqrtMod[4, 234535662466] :[font = output; output; inactive; preserveAspect; endGroup] {2, 20429237894, 24102849916, 30319719362, 38739333488, 44532087808, 50748957254, 55983196106, 59168571380, 74851807172, 76412433998, 83271421298, 95281045064, 100515283916, 103700659190, 113591140658, 120944521808, 130835003276, 134020378550, 139254617402, 151264241168, 158123228468, 159683855294, 175367091086, 178552466360, 183786705212, 190003574658, 195796328978, 204215943104, 210432812550, 214106424572, 234535662464} ;[o] {2, 20429237894, 24102849916, 30319719362, 38739333488, 44532087808, 50748957254, 55983196106, 59168571380, 74851807172, 76412433998, 83271421298, 95281045064, 100515283916, 103700659190, 113591140658, 120944521808, 130835003276, 134020378550, 139254617402, 151264241168, 158123228468, 159683855294, 175367091086, 178552466360, 183786705212, 190003574658, 195796328978, 204215943104, 210432812550, 214106424572, 234535662464} :[font = input; Cclosed; preserveAspect; startGroup] Mod[%^2, 234535662466] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4} ;[o] {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4} :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.7 A General Solution to the Sum-of-Two-Squares Problem :[font = input; preserveAspect] 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 :[font = input; Cclosed; preserveAspect; startGroup] PrimitiveReps[Prime[Range[10000, 10005]]] :[font = output; output; inactive; preserveAspect; endGroup] {{{323, 20}}, {}, {}, {{269, 180}}, {{322, 33}}, {}} ;[o] {{{323, 20}}, {}, {}, {{269, 180}}, {{322, 33}}, {}} :[font = input; Cclosed; preserveAspect; startGroup] SqrtsSquareDivisors[2 3^2 5^2 23] :[font = output; output; inactive; preserveAspect; endGroup] {1, 5, 3, 15} ;[o] {1, 5, 3, 15} :[font = input; Cclosed; preserveAspect; startGroup] n = 2454^2 + 4 GeneralSum2Squares[n] :[font = output; output; inactive; preserveAspect; endGroup] {{1838, 1626}, {1962, 1474}, {2266, 942}, {2322, 794}, {2334, 758}, {2378, 606}, {2446, 198}, {2454, 2}} ;[o] {{1838, 1626}, {1962, 1474}, {2266, 942}, {2322, 794}, {2334, 758}, {2378, 606}, {2446, 198}, {2454, 2}} :[font = input; Cclosed; preserveAspect; startGroup] Map[Plus@@# &, %^2] GeneralSum2Squares[10^10] :[font = output; output; inactive; preserveAspect; endGroup] {{80000, 60000}, {84320, 53760}, {93600, 35200}, {96000, 28000}, {99712, 7584}, {100000, 0}} ;[o] {{80000, 60000}, {84320, 53760}, {93600, 35200}, {96000, 28000}, {99712, 7584}, {100000, 0}} :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] PrimitiveReps from Appendix :[font = input; preserveAspect; endGroup] 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 :[font = input; Cclosed; preserveAspect; startGroup] PrimitiveReps[13904, 1, 7] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] {{31, 43}, {73, 35}, {101, 23}, {109, 17}} ;[o] {{31, 43}, {73, 35}, {101, 23}, {109, 17}} :[font = section; inactive; Cclosed; preserveAspect; startGroup] BONUS on Lagrange's Theorem on 4 squares :[font = text; inactive; Cclosed; preserveAspect; startGroup] 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. :[font = input; preserveAspect; endGroup] (* 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 :[font = input; Cclosed; preserveAspect; startGroup] Sum4Squares[10^10+17] :[font = output; output; inactive; preserveAspect; endGroup] {1864, 7648, 99681, 1316} ;[o] {1864, 7648, 99681, 1316} :[font = input; preserveAspect; endGroup] Plus @@ (%^2) :[font = section; inactive; Cclosed; preserveAspect; startGroup] 9.8 Eisenstein Primes :[font = input; preserveAspect] Representation[p_] := ({u = MiddleRemainder[p, SqrtModPrimeShanks[-3, p]], Sqrt[(p - u^2)/3]}) /; Mod[p, 3] == 1 Attributes[Representation] = Listable :[font = input; Cclosed; preserveAspect; startGroup] Representation[{7, 13, 19, 31, 37, 43, 1299709}] :[font = output; output; inactive; preserveAspect; endGroup] {{2, 1}, {1, 2}, {4, 1}, {2, 3}, {5, 2}, {4, 3}, {791, 474}} ;[o] {{2, 1}, {1, 2}, {4, 1}, {2, 3}, {5, 2}, {4, 3}, {791, 474}} :[font = input; Cclosed; preserveAspect; startGroup] Map[#[[1]]^2 + 3 #[[2]]^2 &, %] :[font = output; output; inactive; preserveAspect; endGroup] {7, 13, 19, 31, 37, 43, 1299709} ;[o] {7, 13, 19, 31, 37, 43, 1299709} :[font = input; preserveAspect] 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]) :[font = input; preserveAspect; endGroup] EisensteinPrimes[10]; ^*)