Pythagorean Triples
October 26, 2012
Today’s exercise feels like a Project Euler problem:
A pythagorean triple consists of three positive integers a, b and c with a < b < c such that a2 + b2 = c2. For example, the three numbers 3, 4 and 5 form a pythagorean triple because 32 + 42 = 9 + 16 = 25 = 52.
The perimeter of a pythagorean triple is the sum of the three numbers that make up the pythagorean triple. For example, the perimeter of the 3, 4, 5 pythagorean triple is 3 + 4 + 5 = 12. There are 17 pythagorean triples with perimeter not exceeding 100. Ordered by ascending perimeter, they are: 3, 4, 5; 6, 8, 10; 5, 12, 13; 9, 12, 15; 8, 15, 17; 12, 16, 20; 7, 24, 25; 15, 20, 25; 10, 24, 26; 20, 21, 29; 18, 24, 30; 16, 30, 34; 21, 28, 35; 12, 35, 37; 15, 36, 39; 9, 40, 41; and 24, 32, 40.
How many pythagorean triples exist with perimeter not exceeding one million?
Your task is to write a program to compute the count of pythagorean triples with perimeter not exceeding one million. 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.
Universal Hash Function
October 23, 2012
Recently, when writing a program that shall remain anonymous, I needed a hash table with keys that could be either strings or integers. That may sound weird to you — it felt weird enough to me that I rearranged things so that all the keys had the same type. But the experience got me thinking about a universal hash function that could be used with keys of any type.
Your task is to write a hash function, suitable for your normal programming environment, that can take a value of any type and return a thirty-two bit integer suitable for use in a hash table. 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.
Prime Partitions
October 19, 2012
Today’s exercise is my penance for a wrong answer at /r/math.
A partition of a number is a way to add numbers to equal the target number; for instance, 1+1+2+3+5 is a partition of 11. We studied partitions in two previous exercises. If all the numbers used in the summation are prime, it is known as a prime partition; for instance, 2+2+2+2+3 = 2+3+3 = 2+2+2+5 = 3+3+5 = 2+2+7 = 11 are the 6 prime partitions for 11. The number of prime partitions of a number is a function known by the Greek letter kappa in number theory, so κ(11)=6. You can see the sequence of prime partitions at A000607.
The usual computation of the number of prime partitions is done in two parts. The first part is a function to compute the sum of the prime factors of a number; for instance, 42=2·3·7, so sopf(42)=12. Note that sopf(12)=5 because the sum is over only the distinct factors of a number, so even though 12=2·2·3, only one 2 is included in the calculation of the sum. You can see the sequence of the sums of the distinct primes dividing n at A008472.
Given the function to compute the sum of the prime factors of a number, the number of prime partitions can be calculated using the following formula:
Your task is to write a function that computes the number of prime partitions of a number; use it to calculate κ(1000). 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.
Version Control
October 16, 2012
We looked at the make programmer’s tool in a previous exercise, and today we look at another programmer’s tool, a version control system. Version control systems allow programmers to keep multiple versions of text files in an economical manner, and there are many of them available: git, subversion, rcs, and the grand-daddy of all of them, sccs.
Our version control system provides three commands. Put takes a single filename argument, stores the file in a history file, and removes the original file from the file system; the history file has the same name as the original file with an added suffix .hist. Get takes a filename argument plus an optional version argument and writes to the file system the requested version of the file, defaulting to the current version of the file if no version number is given; versions are numbered 0 for the current version, 1 for the previous version, 2 for the version before that, and so on. Dir takes a single filename argument and shows a directory listing all available versions of a file. Unlike the bigger version control systems, there are no branches, only a single trunk of history.
The key to the version control system is the history file, which stores the most recent version of the file in its entirety, preceeded by a header line. Then follows a succession of deltas, one for each historical version, each preceeded by its own header line. Here’s a sample history file primes.ss.hist:
@@@ phil Mon Oct 15 18:52:22 CDT 2012 Factorization by trial division
(define (primes n)
(let ((bits (make-vector (+ n 1) #t)))
(let loop ((p 2) (ps (list)))
(cond ((< n p) (reverse ps))
((vector-ref bits p)
(do ((i p (+ i p))) ((< n i))
(vector-set! bits i #f))
(loop (+ p 1) (cons p ps)))
(else (loop (+ p 1) ps))))))(define (factors n)
(let loop ((n n) (f 2) (fs (list)))
(cond ((< n (* f f))
(reverse (cons n fs)))
((zero? (modulo n f)
(loop (/ n f) f (cons f fs)))
(else (loop n (+ f 1) fs)))))
@@@ phil Mon Oct 15 18:50:32 CDT 2012 Fix typo s/vectr/vector
10,17d
@@@ phil Mon Oct 15 18:49:57 CDT 2012 Sieve of Eratosthenes
7c
(vectr-set! bits i #f))
.
The header line that preceeds each delta, including the current file version, starts with three asterisks. Then comes the name of the user that created the delta, and the date and time when it was created. The rest of the line is a comment written by the user when the delta is added to the history file.
Deltas are computed by the Unix diff command with the -e option (we wrote our own version in a previous exercise) comparing the old and new versions and applied by the Unix ed command by applying edit commands to the newest version to recreate older versions. We also use the Unix wc command to count the lines in a delta.
Your task is to write a version control system as described above. 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.
Birthday Paradox
October 12, 2012
In the study of probability, the Birthday Paradox states that in a group of 23 people, there is a 50% chance that two will have the same birthday; in a group of 57 people, the odds rise to 99%.
Your task is to simulate the birthday paradox over many trials and verify the odds. 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.
Two Word Games
October 9, 2012
It’s been a while since we played word games. We have two today:
1) Find all the words in a dictionary that contain exactly five vowels (a, e, i, o and u) in ascending order.
2) Find all the words in a dictionary of length at least six that contain letters in strictly ascending alphabetical order.
These games are easy to play using regular expressions, so you should solve them without regular expressions, using only simple string manipulations. If your system doesn’t provide a dictionary, you can find one at http://icon.shef.ac.uk/Moby/mwords.html.
Your task is to play the two games. 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.
AKS Primality Prover, Part 2
October 5, 2012
In the previous exercise we wrote some auxiliary functions needed for the implementation of the AKS primality prover. Today we will implement the AKS algorithm:
AKS Primality Prover: Given an integer n > 1, determine if it is prime or composite.
1. If n = ab for integers a > 0 and b > 1, output COMPOSITE.
2. Find the smallest r such that ordr(n) > (log2 n)2.
3. If 1 < gcd(a, n) < n for some a ≤ r, output COMPOSITE.
4. If n ≤ r, output PRIME.
5. For each a from 1 to floor √φ(r) · log2 n), if (x + a)n &neq; xn + a (mod xr − 1, n), output COMPOSITE.
6. Output PRIME.
Here ordr(n) is the multiplicative order of n modulo r, both logarithms are to the base 2, φ(r) is Euler’s totient function, and the polynomial modular exponentiation is done as in the previous exercise.
Your task is to write a program to prove primality using the AKS algorithm. 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.
AKS Primality Prover, Part 1
October 2, 2012
In the previous exercise, we looked at Pocklington’s algorithm for proving the primality of a number n; the algorithm depends on at least a partial factorization of n−1. In today’s exercise we look at an algorithm of Manindra Agrawal, Neeraj Kayal, and Nitin Saxena of the Indian Institute of Technology in Kanpur, which does not rely on n to be of any special form, is purely deterministic, of provable polynomial time, not conditional on the Riemann Hypothesis, and unfortunately very slow.
The entire mathematical community was astonished in 2002 when Agrawal, Kayal and Saxena published their algorithm. In the first place, Kayal and Saxena were undergraduate students at the time; Agrawal was their professor. And in the second place, their proof is simple, not relying on any complicated math; the general reaction among mathematicians was “How could we have missed that? What else have we missed?” The three won numerous awards, and their algorithm has been thoroughly studied and analyzed.
We will look at the AKS algorithm in the next exercise. Today, we write two functions that will be needed to implement the AKS algorithm.
The first function computes the multiplicative order of n modulo r, written as ordr(n), which is the smallest number k such that nk ≡ 1 (mod r). The algorithm is straight forward: starting from k = 2, increment k until you find the answer. The multiplicative order exists only when n and r are co-prime.
The second function computes the modular exponentiation of a polynomial; it’s not hard to write, but requires some explanation. To compute the modular exponentiation of a polynomial is simple; it uses the same square-and-multiply algorithm that is used for the modular exponentiation of an integer, which is part of our Standard Prelude and is described in Algorithm 3.C of the essay Programming with Prime Numbers. All we need is the modular multiplication of two polynomials, which we will explain by example.
Consider the problem of squaring polynomial x3 + 4x2 + 12x + 3 modulo (x5 − 1, 17). Polynomial multiplication is exactly the same as grade-school multiplication, except there are no carries, so the process looks like this:
1 4 12 3
1 4 12 3
--- --- --- ---
3 12 36 9
12 48 144 36
4 16 48 12
1 4 12 3
--- --- --- --- --- --- ---
1 8 40 102 168 72 9
Thus, x3 + 4x2 + 12x + 3 squared is x6 + 8x5 + 40x4 + 102x3 + 168x2 + 72x + 9. To reduce the result modulo x5 − 1 we divide using the grade-school long-division algorithm and take the remainder, which gives x + 8 with a remainder of 40x4 + 102x3 + 168x2 + 73x + 17:
1 8
+ --- --- --- --- --- --- ---
1 0 0 0 0 -1 | 1 8 40 102 168 72 9
1 0 0 0 0 -1
--- --- --- --- --- ---
8 40 102 168 73 9
8 0 0 0 0 -8
--- --- --- --- --- ---
40 102 168 73 17
We can confirm the calculation by multiplying and adding the remainder:
1 0 0 0 0 -1
1 8
--- --- --- --- --- ---
8 0 0 0 0 -8
1 0 0 0 0 -1
--- --- --- --- --- --- ---
1 8 0 0 0 -1 -8
40 102 168 73 17
--- --- --- --- --- --- ---
1 8 40 102 168 72 9
Then we simply reduce each of the coefficients of the remainder modulo 17, giving the result 6x4 + 0x3 + 15x2 + 5x + 0. The whole computation can be organized as shown below. Note how the division and reduction modulo x5 − 1 is accomplished, eliminating the high-order coefficients and adding them back to the low-order coefficients; we are relying here on the simple form of the polynomial modulus, and would have to do more work if it was more complicated:
1 4 12 3 multiplicand
1 4 12 3 multiplier
--- --- --- ---
3 12 36 9 3 * 1 4 12 3 * x^0
12 48 144 36 12 * 1 4 12 3 * x^1
4 16 48 12 4 * 1 4 12 3 * x^2
1 4 12 3 1 * 1 4 12 3 * x^3
--- --- --- --- --- --- ---
1 8 40 102 168 72 9 1 4 12 3 * 1 4 12 3
-1 -8 1 8 reduce modulo x^5 - 1
--- --- --- --- --- --- ---
40 102 168 73 17 reduce modulo 17
6 0 15 5 0 final result
Your task is to write functions that compute the multiplicative order of n modulo r and the modular exponentiation of polynomials. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution in the comments below.
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, an−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 = 289−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 2n−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 3n−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 = 289 − 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 c1 and c2 such that n = c2 · f2 + c1 · f1 + 1, which is easily done by extracting the “digits” of n to the base f. Then n is composite if c12 − 4 c2 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 289 − 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 f2 + 96609474553 f + 1, then 966094745532 − 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 = 2521 − 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 2521 − 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 2607 − 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.
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.
Programming Praxis 