Programming Tips

Red stuff is for editor only...red bits in the code should be deleted in the printed version.

The Magic of Imaginary Factoring

by Stan Wagon

Initialize some options (so that my graphics are as I want them, but my code in the text is a little simpler)

SetOptions[Graphics, DefaultFont->{"Times-Roman", 9},
AspectRatio->Automatic];

Decomplexifying an algorithm for complex numbers yields useful tools for the ordinary integers

Introduction

By using some of Mathematica's built-in algorithms in nonobvious ways, we can often get lots of computing power for very little programming effort. In this column I will show how one of Mathematica's most intriguing algorithms --- the one that creates prime factorizations of Gaussian integers --- can yield powerful methods for ordinary integers.

A Gaussian integer has the form a + b i where a and b are ordinary integers (called rational integers). The theory of primes and factorization in this complex realm proceeds pretty much in the same way as for the rational integers (the units here are ±1 and ±i). The criterion for primality is simple: a + b i is a Gaussian prime (henceforth, G-prime) if and only if either b is nonzero and a2 + b2 is prime, or b = 0 and a is a rational prime congruent to 3 (mod 4). For details of the proofs, see [Wagon 1991, chap. 9].

Examples:

 2 = (1 + i)(1 -- i), so 2 is not a Gaussian prime; neither is 5, which factors as (1 + 2 i)(1 -- 2 i). The prime 3, however, is a G-prime.

 If a prime is congruent to 1 (mod 4) then it has the form a2 + b2 (a very short proof is in [Zagier 1990]) and so it factors into (a + b i)(a -- b i).

 If (c + d i)(e + f i) = (a + b i) then one can show that c2 + d2 divides a2 + b2. This is why the primality of a2 + b2 as a rational integer implies the Gaussian primality of a + b i.

 If p is a rational prime congruent to 3 (mod 4), and p = (c + d i)(e + f i), then c2 + d2 is a proper divisor of p2 so it divides p, contradiction. Thus these ordinary primes are G-primes.

The basic Mathematica commands for dealing with primality over the Gaussian integers are PrimeQ, FactorInteger, and Divisors.

PrimeQ[Range[10], GaussianIntegers -> True]

FactorInteger[466 + 2312 I, GaussianIntegers -> True]

Divisors[4524922, GaussianIntegers -> True]

The preceding example yields 12 divisors, but only one of them corresponds to a 20th-century date: 829 + 1959 I. This is connected to a magic trick involving Stepehen Wolfram's birthday, which will be explained in the next section.

The Underlying Algorithm

The main idea for determining Gaussian factors goes back to the 19th century, when Hermite and Legendre used continued fractions as a way of writing a prime as a sum of two squares. G. Cornacchia, in 1908, showed how this could be done directly from the Euclidean algorithm. We will not give the full details here (see [Wagon 1991, chap. 9]), but point out only that the crux of the method is the ability to write a prime p of the form 4k + 1 as a sum of two squares. That can be done by simply applying the Euclidean algorithm to the pair (p, x) where x is a square root of --1 modulo p. The first two remainders under p are the a and b. Finding the x is not a problem since it cen be easily defined in terms of a quadratic nonresidue, and that can be found by trial and error. We give a brief discussion so the reader can sense the magic of this algorithm.

Compositor: There is a centered dot in next line between 73 and 10.

If p = 73, then, because 272 = 729 = 73·10 -- 1, 27 is a square root of --1 mod 73. Applying the Euclidean algorithm to (73, 27) yields the sequence: 73, 27, 19, 8, 3, 2, 1 and, sure enough, the first two remainders under 73 are 8 and 3, and 64 + 9 = 73. This algorithm is blindingly fast, even for 100-digit primes. As mentioned, this is not the full story, but it is definitely the sine qua non for Gaussian factorization; the omitted details use straightforward ideas of elementary number theory, such as the Chinese remainder theorem. For completeness, we include the following code, which gives a method for getting --1 modulo a prime congruent to 1 (mod 4), and a short routine to get the list of remainders in the Euclidean algorithm.

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

Nonresidue[p_?OddQ] := Module[{n = 1},
	While[JacobiSymbol[Prime[n], p] == 1, n++]; 
	Prime[n]]
	
EuclideanRemainders[a_, b_] := EuclideanRemainders[a, b, {a}]

EuclideanRemainders[a_, b_?Positive, rems_List] :=
	 EuclideanRemainders[b, Mod[a, b], Append[rems, b]]

EuclideanRemainders[a_, 0, rems_]  := rems

EuclideanRemainders[73, 29]

We can use much larger primes and the solution to the sum-of-two-squares problem is still pretty much instantaneous. We show only the first two remainders under p.

p = Prime[10^9]  (* happens to be 1 modulo 4 *)

x = SqrtNegOne[p]

Take[Select[EuclideanRemainders[p, x], #^2 < p &], 2]

Apply[Plus, %^2]

Getting Sums of Two Squares

The intricacies of the preceding section can be ignored by the user who just wishes to make efficient use of these algorithms. For since Gaussian factorization is built into Mathematica, we can turn the algorithmic idea on its head and use Gaussian factorization to write any integer as a sum of two squares. If an integer n, not necessarily prime, equals a2 + b2, then b = 
(a + b i)(a -- bi), which means that a + bi appears among the Gaussian divisors of n. The following code selects the Gaussian divisors that work, and returns their real and imaginary parts, sorted and with duplicates removed.

SumTwoSquares[n_] := Union[
  Sort[{Re[#], Im[#]}] & /@ 
    Select[Divisors[n, GaussianIntegers->True], 
      # * Conjugate[#] == n &]]

SumTwoSquares[4524922]

829^2 + 1959^2

In fact, we can extend this code to get all representations of an integer as a sum of d squares. Our method is not very efficient, however. To get all representations as 3 squares we just repeatedly subtract squares z2 from n and get all the representations of n ­ z2 as a sum of two squares. We use the same idea, via recursion, to deal with larger d. This is hopeless for large values of n and d, but we can get, in only a few seconds, all the representations of, say, 1000 as a sum of 4 squares. We cache intermediate results for efficiency.

SumSquares[n_, 2] := SumSquares[n, 2] = Union[
  Map[Sort[{Re[#], Im[#]}] &, Select[
    Divisors[n, GaussianIntegers->True], 
      # * Conjugate[#] == n &]]]

SumSquares[n_, d_] := (SumSquares[n, d] = Union[Map[Sort, Join[
   Append[#, 0] & /@ SumSquares[n, d - 1],
   If[IntegerQ[Sqrt[n]], {Append[Array[0&, d-1], Sqrt[n]]}, {}],
   Flatten[Table[Append[#, z]& /@ SumSquares[n - z^2, d-1], 
                  {z, Sqrt[n - 1]}], 1]] ]])

SumSquares[1000, 4]

I find it remarkable that so much can be accomplished with so little code! Here is a cute magic trick using these ideas [Wagon 1994]. The magician asks a member of the audience to take his day and month of birth, make a 4-digit number out of them (July 9, 1951, becomes 0709), and square that number. Then square the year of birth and add these two squares, presenting the sum to the magician. The magician, with the help of Mathematica, comes back with the birth date in almost no time. The idea is to select the appropriate possibility (there are a handful of cases for which the answer is not unique!) from all the 2-square representations. Here is some code that outputs the date automatically. Potentially good complex numbers are selected by a pure function that checks real and imaginary parts. Then Complex[] :> StringForm[] is used to replace the answer, a complex number here, by a string representation of the date.

BirthDate[n_Integer] := Select[
  Divisors[n, GaussianIntegers->True], 
        1900 < Im[#] < 1985 && 
        Mod[Re[#], 100] <= 31 && 
        # * Conjugate[#] == n &] /. 
  Complex[monthday_, year_] :> 
      StringForm["`` ``, ``", 
    	   {January, February, March, April, May, June,
            July, August, September, October, 
            November, December}[[Quotient[monthday, 100]]], 
        Mod[monthday, 100], year] /. {ans_} :> ans

Stephen Wolfram replied to my query with the 7-digit number 4524922. This number has 12 Gaussian divisors, but Mathematica filters these down to the correct birthdate in an instant.

BirthDate[4524922]

Perhaps the reader can develop a version of this trick that uses small enough numbers that the 3- or 4-squares algorithm can be used. We mention in passing that these and analogous routines for the number of representations of an integer as a sum of 2 or 3 or more squares are useful in investigations related to the Madelung constant of a chemical lattice. See [Crandall 1994, chap. 2]; a discussion of the Madelung constant will appear in my next column in this journal.

Visualizing the Gaussian Primes

An image of the G-primes is quite pleasing visually and has been used in some cloth weavings manufactured in the 1950s. In fact, I have used this pattern to tile my shower walls. To see the pattern, we simply make a table of the G-primes a + b i where a and b are nonengative. Since b + a i is prime iff a + b i is, we then use a transformation rule to turn a single G-prime into two Rectangle objects. This picks up the rest of the first-quadrant primes.

GPrimes[n_] := Select[
  Flatten[Table[{a, b}, {a, 0, n}, {b, 0, a}], 1],
    PrimeQ[#[[1]] + I #[[2]], GaussianIntegers->True]&];

GPrimes[8]

Compositor: The last line of code in following can be omitted
from published version. (except  "];"  of course)

Show[Graphics[
  GPrimes[50] /. {a_Integer, b_Integer} :> 
    {Rectangle[{a, b}, {a+1, b+1}], 
     Rectangle[{b, a}, {b+1, a+1}]}],
  Axes->Automatic,
  Ticks->Table[{i + 0.5, i}, {2}, {i, 10, 50, 10}]];

Figure 1. The Gaussian primes in the first quadrant.

The Moat Problem

Images of the G-primes bring to the fore the infamous moat problem. Everyone knows that a finite person cannot walk along the real axis from 2 to infinity stepping only on the primes. For if the walker can take steps of size K or less, then he is bound to run into a ³moat² of composites, such as L! + 2, L! + 3, . . . , L! + L where L  K + 2. But is such a prime walk to infinity possible in the G-primes? That is the infamous unsolved problem posed by Basil Gordon in 1962. Using extensions of ideas presented here, one can show that step-size 26 is not sufficient to get to infinity [Gethner, Wagon, and Wick]; Jordan and Rabung [Jordan and Rabung 1970] had settled step-size 10 in 1970. Let's use Mathematica to show that such a walk doesn't exist for steps of size 2 or less (i.e., to display a 2-moat). To whet your appetite, have a look at the moat images in Figure 2.

Figure 2: Individual and GraphicsArray

Code for figure 2

GPrimesInDisk[r_] := Select[Flatten[
  Table[{a, b}, {a, 1, r}, {b, 0, Min[a, Sqrt[r^2 - a^2]]}], 1],
    PrimeQ[#[[1]] + I #[[2]], GaussianIntegers->True]&]

poss[{1, 1}, 1] := {{2, 1}, {1, 2}};

poss[pt_, 1] := Map[pt + # &, {{1, 1}, {1, -1}, {-1, 1}, {-1, -1}}];

poss[{1, 1}, 2] := {{2, 1}, {1, 2}, {-1, 1}, {1, -1}};

poss[pt_, 2] := Map[pt + # &, {{1, 1}, {1, -1}, {-1, 1}, {-1, -1},
                            {2, 0}, {0, 2}, {0, -2}, {-2, 0}}];
poss[{1, 1}, 3] := {{2, 1}, {1, 2}, {-1, 1}, {1, -1}};

poss[pt_, 3] := Map[pt + # &, {{1, 1}, {1, -1}, {-1, 1}, {-1, -1},
                  {2, 0}, {0, 2}, {0, -2}, {-2, 0},
                  {2, 2}, {2, -2}, {-2, 2}, {-2, -2}}];

neighbors[Gp_] := Sort[{Gp, #}] & /@ Select[poss[Gp, case], 
           PrimeQ[#[[1]] + I #[[2]], GaussianIntegers->True] &]
           
symmetrize[{a_?NumberQ, b_}] := {{a, b}, {-a, b}, {a, -b}, 
 {-a, -b}, {b, a}, {-b, a}, {b, -a}, {-b, -a}};
 
symmetrize[{p_, q_}] := Transpose[{symmetrize[p], symmetrize[q]}];
   
symmetrize[pts_] := Flatten[symmetrize /@ pts, 1] /; Depth[pts] == 3;

bounds = {15, 50, 100};
case = 1;
  (* The three cases are denoted 1, 2, and 3 *)

roads = Union[Flatten[Map[neighbors, GPrimesInDisk[bounds[[case]]]], 1]];

Clear[mark];
mark[_] := False; 
mark[{1, 1}] = True;

fix[pairs_] := Scan[(If[mark[#[[1]]], mark[#[[2]]] = True,
    If[mark[#[[2]]], mark[#[[1]]] = True]]) &, pairs]
    
    
unmarkedPairs = roads;
FixedPointList[(
  fix[unmarkedPairs]; 
  Length[unmarkedPairs = 
      Select[unmarkedPairs, !mark[#[[1]]] || !mark[#[[2]]] &]])&, 
   Length[unmarkedPairs]]

firstOctant = Select[roads, mark[#[[1]]] &];

roadNetwork = Flatten[symmetrize /@ firstOctant, 1];

unreachablePrimes = symmetrize[Select[
  Flatten[roads, 1], 
    !mark[#] && #[[1]]^2 + #[[2]]^2 <= bounds[[case]]^2 &]];

g[case]=Show[Graphics[{
    Circle[{0, 0}, bounds[[case]]], 
    {AbsolutePointSize[1], Point /@ unreachablePrimes},
    {AbsoluteThickness[0.4], Line /@ roadNetwork}}],
    AspectRatio->Automatic];

Show[GraphicsArray[{g[1], g[2], g[3]}]];

Figure 2.   These iamges show how far one can walk from the origin using steps of size 2, 2, or 8, respectively, and stepping only on Gaussian primes.

We can break down the problem into four steps.

1.  Identify the G-primes in a given disk (first octant only, because of symmetry).
2.  For each identified G-prime, find all  its 2-neighbors: the G-primes at distance two or less.
3.  Form the road network consisting of all lines between G-primes and their 2-neighbors.
4.  Determine the connected component of 1 + i in this network.

Step 1 is easy. We just select the G-primes from all possibilities.

GPrimesInDisk[r_] := Select[Flatten[
  Table[{a, b}, {a, 1, r}, {b, 0, Min[a, Sqrt[r^2 - a^2]]}], 1],
    PrimeQ[#[[1]] + I #[[2]], GaussianIntegers->True]&]

GPrimesInDisk[8]

Step 2 is easy too. Because of parity considerations, there are essentially two possibilities for a 2-neighbor: add (2, 0) or (1, 1). But signs and order yield 8 possibilities, and {1, 1} is a special case. The function poss forms the set of possibilities.

poss[{1, 1}] := {{2, 1}, {1, 2}, {-1, 1}, {1, -1}};
poss[pt_] := Map[pt + # &, {{1, 1}, {1, -1}, {-1, 1}, {-1, -1},
                            {2, 0}, {0, 2}, {0, -2}, {-2, 0}}];

Step 3 is easy too: Just agglomerate all the neighbors of all the G-primes in a disk of radius 50. We use Sort so that Union, which we use in roads, will eliminate duplicates. As a reality check, we use Line to show the result.

neighbors[Gp_] := Sort[{Gp, #}] & /@ Select[poss[Gp], 
           PrimeQ[#[[1]] + I #[[2]], GaussianIntegers->True] &]
           
roads = Union[Flatten[Map[neighbors, GPrimesInDisk[50]], 1]];

Compositor: In code below, erase the Ticks option from printed version.

Show[Graphics[{Map[Line, roads]}],
   Axes->Automatic, Ticks->Table[Range[10, 50, 10], {2}]];

Figure 3. The road network in the Gaussian primes using steps of size 2 or less.


Now, step 4 is the tricky part. We do it by first defining mark to be False on everything and True on {1, 1}. Then we can scan (via fix[]) the pairs that make up the roads, marking a road-end as True if the other end is so marked. When this process ceases making a change, we will have the G-primes in the connected component marked as True. To decide when this process is finished, we call on FixedPointList (note that FixedPoint would suffice) using a function that returns the number of unmarked pairs.

Clear[mark];
mark[_] := False;
mark[{1, 1}] = True;

fix[pairs_] := Scan[Which[
    mark[#[[1]]], mark[#[[2]]] = True,
    mark[#[[2]]], mark[#[[1]]] = True] &, pairs]

unmarkedPairs = roads;
FixedPointList[(
  fix[unmarkedPairs]; 
  Length[unmarkedPairs = 
      Select[unmarkedPairs, !mark[#[[1]]] || !mark[#[[2]]] &]])&, 
   Length[unmarkedPairs]]

This shows that 111 of our original 266 pairs are unmarked. To see the result in its symmetric form we define symmetrize to apply the 8 transformations, and then symmetrize the road network consisting of reachable G-primes and use just dots for the unreachable ones. The following code generates the image of the 2-moat shown in Figure 2.

symmetrize[{a_?NumberQ, b_}] := {{a, b}, {-a, b}, {a, -b}, 
 {-a, -b}, {b, a}, {-b, a}, {b, -a}, {-b, -a}};
 
symmetrize[{p_, q_}] := Transpose[{symmetrize[p], symmetrize[q]}];
   
symmetrize[pts_] := 
   Flatten[symmetrize /@ pts, 1] /; Depth[pts] == 3;

firstOctant = Select[roads, mark[#[[1]]] &];

roadNetwork = Flatten[symmetrize /@ firstOctant, 1];

unreachablePrimes = symmetrize[Select[
  Flatten[roads, 1], !mark[#] && #[[1]]^2 + #[[2]]^2 <= 50^2 &]];

Show[Graphics[{Circle[{0, 0}, 50], 
    {AbsolutePointSize[1], Point /@ unreachablePrimes},
    {AbsoluteThickness[0.4], Line /@ roadNetwork}}],
 AspectRatio -> Automatic];

References

Richard E. Crandall, Projects in Scientific Computation,TELOS/Springer-Verlag,
                         New York, 1994.

Ellen Gethner, Stan Wagon, and Brian Wick, A stroll through the Gaussian integers, forthcoming.

J. H. Jordan and J. R. Rabung, A conjecture of Paul Erdös, Mathematics of 
                         Computation 24 (1970) 221-223.

Stan Wagon, Mathematica in Action, W. H. Freeman, New York, 1991.

Stan Wagon, A Mathematica'l magic trick, College Mathematics Journal 25 (1994) 
                         325­326.
Don Zagier, A one-sentence proof that every prime p  1 (mod 4) is a sum of two 
                          squares, American Mathematical Monthly 97 (1990) 144.

About the author (use same copy as in 4-4)

Some color images of 2-moats and 3-moats for possible cover use. Editor: If you want code, please ask.

A color 2-moat

A color 3-moat

A poster-quality image with the first-quadrant Gaussian primes in the
background.