## Pocklington’s N-1 Primality Prover

### September 28, 2012

In previous exercises we have seen two methods for determining if a number is probably prime, the Miller-Rabin and Baillie-Wagstaff methods, and two methods for proving a number prime, trial division and the Lucas N-1 prover, which requires the complete factorization of *n*−1, where *n* is the number to be proved prime. In today’s exercise we examine a method of Pocklington in which a number can be proved prime with only a partial factorization of *n*−1.

We begin by expressing the number to be proved prime as *n* = *f* · *r* + 1, where the complete factorization of *f* is known and exceeds the square root of *n*. Pocklington proved that if, for some *a*, *a*^{n−1} ≡ 1 (mod *n*) and gcd(*a*^{(n−1)/q} − 1, *n*) = 1 for each prime *q*|*f*, then *n* is prime if *f* > √n.

As an example, let’s consider the primality of *n* = 2^{89}−1. A probable-prime test suggests that *n* is prime. Trial division by the primes less than a thousand tells us that *n* = 2 · 3 · 5 · 17 · 23 · 89 · 353 · 397 · 683 · 6194349127121 + 1 = 99924948842910 · 6194349127121 + 1, where the primality of 6194349127121 is indeterminate (it’s composite, but we don’t care). Since *f* = 99924948842910 is greater than the square root of *n*, we will be able to prove the primality of *n* by Pocklington’s theorem.

Consider first the situation *a* = 2; we can choose any *a* randomly from the range 1 < *a* < *n* − 1, but it is easiest to choose *a* from the sequence 2, 3, 4, …. Now 2^{n−1} ≡ 1 (mod *n*), so the first criterion holds, but 2^{(n−1)/2} ≡ 1 (mod *n*), so the gcd is *n*, and the second criterion fails. Next we consider *a* = 3. Now 3^{n−1} ≡ 1 (mod *n*), so the first criterion holds, and the second criterion holds for each *q* ∈ {2, 3, 5, 17, 23, 89, 353, 397, 683}, so *n* = 2^{89} − 1 is prime.

Brillhart, Lehmer and Selfridge later extended Pocklington’s theorem to work with *f* between the cube root and square root of *n*. The idea is to determine the constants *c*_{1} and *c*_{2} such that *n* = *c*_{2} · *f*^{2} + *c*_{1} · *f*^{1} + 1, which is easily done by extracting the “digits” of *n* to the base *f*. Then *n* is composite if *c*_{1}^{2} − 4 *c*_{2} is a perfect square and prime if it is not; Pocklington’s criteria must also hold.

With this extension to Pocklington’s theorem we can prove the primality of 2^{89} − 1 with only the factors less than 500. With *a* = 3, Pocklington’s first criterion holds, and Pocklington’s second criterion holds for each *q* ∈ {2, 3, 5, 17, 23, 89, 353, 397}, and *f* = 146302999770 is greater than the cube root of *n*, so we calculate *n* = 28917 *f*^{2} + 96609474553 *f* + 1, then 96609474553^{2} − 4 · 28917 = 9333390573406754434141 = 12611 · 95783 · 7726832165257 is not a perfect square, so *n* is prime.

Pocklington’s method fails if you can’t find sufficient factors of *n* − 1. We prefer trial division rather than some more powerful method because we have to prove that each of the factors is prime, and trial division does that for us automatically; if you can’t factor *n*−1 by trial division, but you can factor it by Pollard’s rho method or something even more powerful, you could still use Pocklington’s theorem, proving each of the factors prime by using Pocklington’s theorem on each of them recursively. For instance, if you are determined to prove the primality of *n* = 2^{521} − 1, trial division by the primes less than a thousand finds the factors 2, 3, 5, 5, 11, 17, 31, 41, 53, 131, 157 and 521, and Pollard’s rho method quickly finds the factors 1613, 2731, 8191, 42641, 51481, 61681, 409891, and 858001, all of which can be proved prime by Pocklington’s method. Taken together, the product of those factors is greater than the cube root of *n*, so you can prove the primality of 2^{521} − 1 even though you can’t factor *n* − 1 by trial division. Of course, it is still possible that you can’t find enough factors (as an example, try to prove the primality of 2^{607} − 1), in which case you need a more powerful method, but that’s a topic for a different exercise.

Your task is to write a program that proves the primality of an input number *n*, or shows that it is composite, using the method of Pocklington with extensions from Brillhart, Lehmer and Selfridge. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

## Essay: Programming With Prime Numbers

### September 25, 2012

We introduce today a new feature of Programming Praxis: essays. An essay isn’t an exercise; it doesn’t present a task for you to perform, but instead provides extended discussion and complete source code. Though essays may be based on exercises, the idea is that an essay can delve more deeply into a topic than is possible in the limited space and time of an exercise.

Our first essay is about a topic dear to my heart: programming with prime numbers. The essay presents two versions of the Sieve of Eratosthenes, uses trial division to identify primes and factor composites, and provides the Miller-Rabin pseudoprime checker and Pollard rho factoring algorithm. Source code is given in five languages: C, Haskell, Java, Python, and Scheme.

Please read the essay, and feel free to comment on it below; comments on the essay itself are closed. Let me know if you see any errors. And feel free to link to the essay on the internet if you know of places where it is appropriate.

## Autumn Equinox

### September 21, 2012

Today is the autumn equinox, when the day and night are of equal lengths at the equator.

Your task is to write a program that calculates the length of the day and the night to show that they are of equal length. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

## ABC Conjecture

### September 18, 2012

The ABC Conjecture has recently been in the news on math blogs because of the claim that it has been proved by Shinichi Mochizuki. Though the proof is being taken seriously, due to Mochizuki’s reputation, it is five hundred pages long, and confirmation will take several months. The conjecture states that given two positive integers *a* and *b* and their sum *c* with no common factors, the product of the distinct prime factors of *abc* is rarely much smaller than *c*.

The *radical* of a number *n* is the product of the distinct prime factors of *n*; for instance, rad(18) = 6 because 18 = 2 × 3 × 3 and, eliminating the duplicate occurrence of 3, 2 × 3 = 6. The *quality* of an (*a*,*b*,*c*) triple is given by *q*(*a*,*b*,*c*) = log(*c*) / log(rad(*abc*)). The precise statement of the ABC conjecture is that for every ε > 0 there are only finitely many triples (*a*,*b*,*c*) with *a* and *b* coprime positive integers and *a* + *b* = *c* such that rad(*abc*)^{1+ε} < *c*, or equivalently, such that *q*(*a*,*b*,*c*) > 1.

Your task is to write the functions rad and q and find all the triples with *c* < 1000 for which *q*(*a*,*b*,*c*) > 1. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

## Tribonacci Numbers

### September 14, 2012

You will recall that fibonacci numbers are formed by a sequence starting with 0 and 1 where each succeeding number is the sum of the two preceeding numbers; that is, F[n] = F[n-1] + F[n-2] with F[0] = 0 and F[1] = 1. We studied fibonacci numbers in a previous exercise.

Tribonacci numbers are like fibonacci numbers except that the starting sequence is 0, 0 and 1 and each succeeding number is the sum of the three preceeding numbers; that is, T[n] = T[n-1] + T[n-2] + T[n-3] with T[-1] = 0, T[0] = 0 and T[1] = 1. The powering matrix for tribonacci numbers, used similarly to the powering matrix for fibonacci numbers, is:

The first ten terms of the tribonacci sequence, ignoring the two leading zeros, are 1, 1, 2, 4, 7, 13, 24, 44, 81 and 149 (A000073).

Your task is to write two functions that calculate the first *n* terms of the tribonacci sequence by iteration and the *n*th term by matrix powering; you should also calculate the tribonacci constant, which is the limit of the ratio between successive tribonacci numbers as *n* tends to infinity. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

## The Sum Of The First Billion Primes

### September 11, 2012

The first ten prime numbers are 2, 3, 5, 7, 11, 13, 17, 19, 23, and 29. Their sum is 129.

Your task is to write a program that calculates the sum of the first billion primes. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

## The First Two Programs

### September 7, 2012

In their book *The C Programming Language*, Brian Kernigan and Dennis Ritchie say that the first program every programmer should write is a program that writes “Hello, world!” to the console. Then they give the second program that produces a fahrenheit/celsius temperature conversion table, with fahrenheit temperatures every 20 degrees from 0 to 300 and the corresponding celsius temperature to the nearest degree, each pair written on its own line with the two temperatures separated by tabs.

Your task is to write the first two programs. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

## Fountain Codes

### September 4, 2012

Fountain codes provide an effective method of file transmission over lossy channels. Invented by Michael Luby in 1998, but not published unti 2002, they are used to send bit torrents and streaming media, transmit map and traffic data to in-car receivers, and probably for the recent “brain transplant” to the Curiosity rover on Mars (I’m speculating, but wouldn’t you expect to lose a few packets between here and Mars?). The idea is to send many chunks of a file, larger than the original file, like a fountain spraying droplets of water, and for the receiver to collect drops in a bucket, which can be reassembled into the original file without collecting all the drops sprayed by the fountain.

The encoding process is simple. Pick a random number *d* between 1 and *n*, the number of blocks in the file. Pick *d* blocks from the file at random (no duplicates) and xor them together. Send the xor of the combined blocks, along with the list of blocks used to make the xor. The distribution of *d* is important; it is not just a uniform random number. We choose *d* with probability *d* / (1/2)*n*(*n*+1), so if the file has 20 blocks, the denominator is 210 (the sum of the integers from 1 to 20), the degree *d* is 1 with probability 20/210, 2 with probability 19/210, and so on, until it is 20 with probability 1/210.

The decoding process isn’t much harder. Start with an array of blocks of length *n*, with an indication that each block is not currently decoded. Then each packet is analyzed as it is received. If the packet has a single block, it can be added immediately to the array of blocks; if not, any known blocks in the input list can be xor’ed into the packet, reducing its degree. A packet that is reduced to a single block is added to the output, others are added to a hold list. Any time a new block is discovered, the entire hold list is scanned for any further reductions in any of its packets.

Let’s illustrate with a fountain that sends blocks of a single byte from the “Programming Praxis” input string. We begin with an empty array and hold list and fetch the first packet from the fountain (22 15 9). That packet has two blocks, so we put it in the hold and go on. The next packet is (80 0), so we put a “P” in block 0 and go on. The next packet is (55 12 10), which we add to the hold. Next we get a packet (14 10 8), which doesn’t help immediately, but it shares a block with (55 12 10), which will be useful later; we add it to the hold, which now contains (22 15 9), (55 12 10) and (14 10 8). The next packet is (103 10). We can immediately add a “g” at block 10, then xor 55 with 103 to get the “P” at block 12 and xor 14 with 103 to get the “i” at block 8; our string is now “P…….i.g.P…..” and our hold is now (22 15 9). We next receive packets (19 5 1), (88 11 15), (45 1 9 14 12), (105 16) and (114 1), so we can add an “i” at block 16, an “r” at block 1, and, by virtue of the (19 5 1) packet in the hold, an “a” at block 5; now the string is “Pr…a..i.g.P…..” and the hold contains (22 15 9), (88 11 15), and (45 1 9 14 12). The next packet is (97 14), which gives us another block, the “a” of “Praxis,” and also lets us reduce the packet (45 1 9 14 12) in the hold to (76 1 9 12). The next packet is (99 16 4 10 1 7), which we can reduce to (118 16 4 7) by applying blocks 10 and 1. Next we receive (120 15), which gives us the “x” in block 15 directly, and also gives us the “n” in block 9 and the space character in block 11 by applying block 15 to the (22 15 9) and (88 11 15) items in the hold; we now have “Pr…a..ing P..x..” and packets (76 1 9 12) and (118 16 4 7) in the hold, with nine blocks left to be decoded. And so the process continues. Eventually the entire string is revealed, on average after *O*(*n* log_{e} *n*) packets have been received, which is about 52 for our 18-block string.

In practice, it is common to do things a little bit differently. Certainly it is normal to use blocks larger than a single byte. Instead of sending the list of blocks, it is normal to send the seed to a simple random number generator; if both ends of the transmission use the same random number generator, both can create the list of blocks from the same seed. And it is common to send a checksum with each packet; if the packet fails the checksum, it is simply ignored, as if it had never been received.

It is also common to use a different degree-distribution than the one we used here. It turns out that a bimodal degree-distribution with lots of 1-, 2- or 3-block packets, combined with an n-block packet and a few n-minus-two-or-three block packets, with very few packets in between, makes the whole process go much quicker, an average *O*(*n*) with just a small overhead of 5% or 10% over the size of the file being transmitted. Use of this degree-distribution is called a *raptor code*.

Your task is to write a program that creates a fountain and decoder. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.