This section describes two methods for checking the primality of an
integer $n$, one with order of growth
$\Theta(\sqrt{n})$, and a
probabilistic algorithm with order of growth
$\Theta(\log n)$. The exercises at the end of
this section suggest programming projects based on these algorithms.
Searching for divisors
Since ancient times, mathematicians have been fascinated by problems
concerning prime numbers, and many people have worked on the problem
of determining ways to test if numbers are prime. One way
to test if a number is prime is to find the number's divisors. The
following program finds the smallest integral divisor (greater than 1)
of a given number $n$. It does this in a
straightforward way, by testing $n$ for
divisibility by successive integers starting with 2.
Original
JavaScript
(define (smallest-divisor n)
(find-divisor n 2))
(define (find-divisor n test-divisor)
(cond ((> (square test-divisor) n) n)
((divides? test-divisor n) test-divisor)
(else (find-divisor n (+ test-divisor 1)))))
(define (divides? a b)
(= (remainder b a) 0))
function smallest_divisor(n) {
return find_divisor(n, 2);
}
function find_divisor(n, test_divisor) {
return square(test_divisor) > n
? n
: divides(test_divisor, n)
? test_divisor
: find_divisor(n, test_divisor + 1);
}
function divides(a, b) {
return b % a === 0;
}
We can test whether a number is prime as follows:
$n$ is prime if and only if
$n$ is its own smallest divisor.
Original
JavaScript
(define (prime? n)
(= n (smallest-divisor n)))
function is_prime(n) {
return n === smallest_divisor(n);
}
The end test for
find-divisorfind_divisor
is based on the fact that if $n$ is not prime it
must have a divisor less than or equal to
$\sqrt{n}$.[1]
This means that the algorithm need only test divisors between 1 and
$\sqrt{n}$. Consequently, the number of steps
required to identify $n$ as prime will have order
of growth $\Theta(\sqrt{n})$.
The Fermat test
The $\Theta(\log n)$ primality test is based on
a result from number theory known as
Fermat's Little
Theorem.[2]
Fermat's Little Theorem:
If $n$ is a prime number and
$a$ is any positive integer less than
$n$, then $a$ raised
to the $n$th power is congruent to
$a$ modulo $n$.
(Two numbers are said to be
congruent modulo
$n$ if they both have the same remainder when
divided by $n$. The remainder of a number
$a$ when divided by
$n$ is also referred to as the
remainder of $a$ modulo
$n$, or simply as $a$
modulo $n$.)
If $n$ is not prime, then, in general, most of
the numbers $a < n$ will not satisfy the above
relation. This leads to the following algorithm for testing primality:
Given a number $n$, pick a
random number $a < n$ and compute the
remainder of $a^n$ modulo
$n$. If the result is not equal to
$a$, then $n$ is
certainly not prime. If it is $a$, then chances
are good that $n$ is prime. Now pick another
random number $a$ and test it with the same
method. If it also satisfies the equation, then we can be even more
confident that $n$ is prime. By trying more and
more values of $a$, we can increase our
confidence in the result. This algorithm is known as the Fermat test.
To implement the Fermat test, we need a
procedurefunction
that computes the
exponential of a number modulo another number:
Original
JavaScript
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder
(square (expmod base (/ exp 2) m))
m))
(else
(remainder
(* base (expmod base (- exp 1) m))
m))))
function expmod(base, exp, m) {
return exp === 0
? 1
: is_even(exp)
? square(expmod(base, exp / 2, m)) % m
: (base * expmod(base, exp - 1, m)) % m;
}
This is very similar to the
fast-exptfast_exptprocedurefunction
of section 1.2.4. It uses successive
squaring, so that the number of steps grows logarithmically with the
exponent.[3]
The Fermat test is performed by choosing at random a number
$a$ between 1 and
$n-1$ inclusive and checking whether the remainder
modulo $n$ of the
$n$th power of $a$ is
equal to $a$. The random number
$a$ is chosen using the
procedure
random,
primitive function
math_random,
which we assume is included as a primitive in Scheme.Random
returns a nonnegative integer less than its integer input. Hence, to obtain
a random number between 1 and $n-1$, we call
random with an input of
$n-1$ and add 1 to the result:
which returns a nonnegative number less than 1. Hence, to obtain
a random number between 1 and $n-1$, we multiply
the return value of math_random by
$n-1$, round down the result with the
primitive function
math_floor,
and add 1:
Original
JavaScript
(define (fermat-test n)
(define (try-it a)
(= (expmod a n n) a))
(try-it (+ 1 (random (- n 1)))))
function fermat_test(n) {
function try_it(a) {
return expmod(a, n, n) === a;
}
return try_it(1 + math_floor(math_random() * (n - 1)));
}
The following
procedurefunction
runs the test a given number of times, as specified by a parameter. Its
value is true if the test succeeds every time, and false otherwise.
Original
JavaScript
(define (fast-prime? n times)
(cond ((= times 0) true)
((fermat-test n)
(fast-prime? n (- times 1)))
(else false)))
function fast_is_prime(n, times) {
return times === 0
? true
: fermat_test(n)
? fast_is_prime(n, times - 1)
: false;
}
Probabilistic methods
The Fermat test differs in character from most familiar algorithms, in which
one computes an answer that is guaranteed to be correct. Here, the answer
obtained is only probably correct. More precisely, if
$n$ ever fails the Fermat test, we can be certain
that $n$ is not prime. But the fact that
$n$ passes the test, while an extremely strong
indication, is still not a guarantee that $n$ is
prime. What we would like to say is that for any number
$n$, if we perform the test enough times and find
that $n$ always passes the test, then the
probability of error in our primality test can be made as small as we like.
Unfortunately, this assertion is not quite correct. There do exist numbers
that fool the Fermat test: numbers $n$ that are
not prime and yet have the property that $a^n$ is
congruent to $a$ modulo
$n$ for all integers
$a < n$. Such numbers are extremely rare, so
the Fermat test is quite reliable in practice.[4]
There are variations of the Fermat test that cannot be fooled. In these
tests, as with the Fermat method, one tests the primality of an integer
$n$ by choosing a random integer
$a < n$ and checking some condition that
depends upon $n$ and
$a$. (See
exercise 1.28 for an example of such a test.)
On the other hand, in contrast to the Fermat test, one can prove that, for
any $n$, the condition does not hold for most of
the integers $a < n$ unless
$n$ is prime. Thus, if
$n$ passes the test for some random choice
of $a$, the chances are better than even
that $n$ is prime. If
$n$ passes the test for two random choices of
$a$, the chances are better than 3 out of 4 that
$n$ is prime. By running the test with more and
more randomly chosen values of $a$ we can make
the probability of error as small as we like.
The existence of tests for which one can prove that the chance of error
becomes arbitrarily small has sparked interest in algorithms of this type,
which have come to be known as probabilistic algorithms. There is
a great deal of research activity in this area, and probabilistic algorithms
have been fruitfully applied to many fields.[5]
Exercise 1.21
Use the
smallest-divisorsmallest_divisorprocedurefunction
to find the smallest divisor of each of the following numbers: 199, 1999,
19999.
Exercise 1.22
Most Lisp implementations include a primitive called
runtime
that returns an integer that specifies the amount of time the system has
been running (measured, for example, in microseconds). The following
timed-prime-test
procedure,
Assume a primitive function
get_time of no arguments
that returns the number of milliseconds that have passed since 00:00:00 UTC
on Thursday, 1 January, 1970.[6]
The following
timed_prime_test
function,
when called with an integer $n$, prints
$n$ and checks to see if
$n$ is prime. If $n$
is prime, the
procedurefunction
prints three asterisks[7] followed by the amount of
time used in performing the test.
Using this
procedure,function,
write a
procedurefunctionsearch-for-primessearch_for_primes
that checks the primality of consecutive odd integers in a specified range.
Use your
procedure
function
to find the three smallest primes larger than 1000; larger than 10,000;
larger than 100,000; larger than 1,000,000. Note the time needed to test
each prime. Since the testing algorithm has order of growth of
$\Theta(\sqrt{n})$, you should expect that testing
for primes around 10,000 should take about
$\sqrt{10}$ times as long as testing for primes
around 1000. Do your timing data bear this out? How well do the data for
100,000 and 1,000,000 support the $\sqrt{n}$
prediction? Is your result compatible with the notion that programs on
your machine run in time proportional to the number of steps required for
the computation?
Original
JavaScript
function search_for_primes(start, times) {
return times === 0
? true
: start > 2 && start % 2 === 0
? search_for_primes(start + 1, times)
// if we get true, it's a prime
: timed_prime_test(start)
? search_for_primes(start + 2, times - 1)
: search_for_primes(start + 2, times);
}
The timing data quite clearly supports the $\sqrt{n}$ for prediction that is sufficiently big, such as 100,000 and 1,000,000.
Exercise 1.23
The
smallest-divisorsmallest_divisorprocedurefunction
shown at the start of this section does lots of needless testing: After it
checks to see if the number is divisible by 2 there is no point in checking
to see if it is divisible by any larger even numbers. This suggests that
the values used for
test-divisortest_divisor
should not be 2, 3, 4, 5, 6, … but rather 2, 3, 5, 7, 9,
…. To implement this change,
define a proceduredeclare a functionnext that returns 3 if its input is equal to 2
and otherwise returns its input plus 2. Modify the
smallest-divisorsmallest_divisorprocedurefunction
to use
(next test-divisor)next(test_divisor)
instead of
(+ test-divisor 1).test_divisor + 1.
With
timed-prime-testtimed_prime_test
incorporating this modified version of
smallest-divisor,smallest_divisor,
run the test for each of the 12 primes found in
exercise 1.22.
Since this modification halves the number of test steps, you should expect
it to run about twice as fast. Is this expectation confirmed? If not, what
is the observed ratio of the speeds of the two algorithms, and how do you
explain the fact that it is different from 2?
Original
JavaScript
function next(input) {
return input === 2
? 3
: input + 2;
}
function find_divisor(n, test_divisor) {
return square(test_divisor) > n
? n
: divides(test_divisor, n)
? test_divisor
: find_divisor(n, next(test_divisor));
}
The ratio of the speeds of the two algorithms is not exactly 2, but this might be due to hardware / network issues.
It is about 1.5 times faster compared to previous solution.
Exercise 1.24
Modify the
timed-prime-testtimed_prime_testprocedurefunction
of exercise 1.22 to use
fast-prime?fast_is_prime
(the Fermat method), and test each of the 12 primes you found in that
exercise. Since the Fermat test has
$\Theta(\log n)$ growth, how would you expect
the time to test primes near 1,000,000 to compare with the time needed to
test primes near 1000? Do your data bear this out? Can you explain any
discrepancy you find?
Original
JavaScript
function timed_prime_test(n) {
display(n);
return start_prime_test(n, get_time());
}
function start_prime_test(n, start_time) {
return fast_is_prime(n, math_floor(math_log(n)))
? report_prime(get_time() - start_time)
: true;
}
function report_prime(elapsed_time) {
display(" *** ");
display(elapsed_time);
}
The time to test primes near 1,000,000 using fast_is_prime
is about 4 ms, 4 times the time needed to test primes near 1,000. This is faster compared to 8 ms
that we achieved if we use the is_prime. However, despite being 4 times slower,
this fact cannot lead us to believe that it has a greater growth than $\Theta(\log n)$,
as it should be tested with greater numbers to gain a more accurate understanding of the growth of the function.
Exercise 1.25
Alyssa P. Hacker complains that we went to a lot of extra work in writing
expmod. After all, she says, since we already
know how to compute exponentials, we could have simply written
Original
JavaScript
(define (expmod base exp m)
(remainder (fast-expt base exp) m))
function expmod(base, exp, m) {
return fast_expt(base, exp) % m;
}
Is she correct?
Would this
procedurefunction
serve as well for our fast prime tester? Explain.
Alyssa's suggestion is correct at first sight: her
expmod function computes
$\textit{base}^{\textit{exp}}$ and
then finds its remainder modulo $m$, as
required in the Fermat test.
However, for large bases, Alyssa's method will quickly bump into
limitations because JavaScript uses 64 bits to represent numbers,
following the double-precision floating point standard. When the
numbers become so large that they cannot be represented precisely
any longer in this standard, the results become unreliable. Even
worse, the method might exceed the largest number that can be
represented in this standard, and the computation leads to an
error.
For small bases, however, Alyssa's method may be even faster than
the original expmod function,
because it will carry out only one single remainder operation.
Exercise 1.26
Louis Reasoner is having great difficulty doing
exercise 1.24.
His
fast-prime?fast_is_prime
test seems to run more slowly than his
prime?is_prime
test. Louis calls his friend Eva Lu Ator over to help. When they examine
Louis's code, they find that he has rewritten the
expmodprocedurefunction
to use an explicit multiplication, rather than calling
square:
Original
JavaScript
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder
(* (expmod base (/ exp 2) m)
(expmod base (/ exp 2) m))
m))
(else
(remainder
(* base (expmod base (- exp 1) m))
m))))
function expmod(base, exp, m) {
return exp === 0
? 1
: is_even(exp)
? ( expmod(base, exp / 2, m)
* expmod(base, exp / 2, m)) % m
: (base * expmod(base, exp - 1, m)) % m;
}
I don't see what difference that could make,
says Louis. I do. says Eva. By writing the
procedurefunction
like that, you have transformed the
$\Theta(\log n)$ process into a
$\Theta(n)$ process. Explain.
Eva is correct: by evaluating the expression:
(expmod(base, exp / 2, m) * expmod(base, exp / 2, m)) % m
the expression
expmod(base, exp / 2, m)
is evaluated twice at each step in the computation when the exponent is
even, eliminating the benefit of the fast exponentiation
algorithm—which halves the exponent when the exponent is
even—therefore eliminating the feature of the algorithm that makes
it faster.
Exercise 1.27
Demonstrate that the
Carmichael numbers listed in
footnote 4 really do fool the Fermat
test. That is, write a
procedurefunction
that takes an integer $n$ and tests whether
$a^n$ is congruent to
$a$ modulo $n$ for
every $a < n$, and try your
procedurefunction
on the given Carmichael numbers.
Original
JavaScript
function carmichael(n) {
function expmod(base, exp, m) {
return exp === 0
? 1
: is_even(exp)
? square(expmod(base, exp / 2, m)) % m
: (base * expmod(base, exp - 1, m)) % m;
}
function fermat_test(n, a) {
return expmod(a, n, n) === a;
}
function iter(n, i) {
return i === n
? true
: fermat_test(n, i)
? iter(n, i + 1)
: false;
}
return iter(n, 2);
}
Exercise 1.28
One variant of the Fermat test that cannot be fooled is called the
Miller–Rabin test (Miller 1976;
Rabin 1980). This starts from
an alternate form of
Fermat's Little Theorem, which states that if
$n$ is a prime number and
$a$ is any positive integer less than
$n$, then $a$ raised
to the $(n-1)$st power is congruent to 1
modulo $n$. To test the primality of a
number $n$ by the Miller–Rabin test, we pick a
random number $a < n$ and raise
$a$ to the $(n-1)$st
power modulo $n$ using the
expmodprocedure.function.
However, whenever we perform the squaring step in
expmod, we check to see if we have discovered a
nontrivial square root of 1
modulo $n$,
that is, a number not equal to 1 or $n-1$ whose
square is equal to 1 modulo $n$. It is
possible to prove that if such a nontrivial square root of 1 exists, then
$n$ is not prime. It is also possible to prove
that if $n$ is an odd number that is not prime,
then, for at least half the numbers $a < n$,
computing $a^{n-1}$ in this way will reveal a
nontrivial square root of 1 modulo $n$.
(This is why the Miller–Rabin test cannot be fooled.) Modify the
expmodprocedurefunction
to signal if it discovers a nontrivial square root of 1, and use this to
implement the Miller–Rabin test with a
procedurefunction
analogous to
fermat-test.fermat_test.
Check your
procedurefunction
by testing various known primes and non-primes. Hint: One convenient way to
make expmod signal is to have it return 0.
Original
JavaScript
function random(n) {
return math_floor(math_random() * n);
}
function miller_rabin_test(n) {
function expmod(base, exp, m) {
return exp === 0
? 1
: is_even(exp)
? square(trivial_test(expmod(base,
exp / 2,
m),
m))
% m
: (base * expmod(base, exp - 1, m))
% m;
}
function trivial_test(r, m) {
return r === 1 || r === m - 1
? r
: square(r) % m === 1
? 0
: r;
}
function try_it(a) {
return expmod(a, n - 1, n) === 1;
}
return try_it(1 + random(n - 1));
}
function do_miller_rabin_test(n, times) {
return times === 0
? true
: miller_rabin_test(n)
? do_miller_rabin_test(n, times - 1)
: false;
}
[1]
If
$d$ is a divisor of
$n$, then so is $n/d$.
But $d$ and $n/d$
cannot both be greater than $\sqrt{n}$.
[2]
Pierre
de Fermat (1601–1665) is considered to be
the founder of modern
number theory. He obtained many important number-theoretic results,
but he usually announced just the results, without providing his proofs.
Fermat's Little Theorem was stated in a letter he wrote in 1640.
The first published proof was given by
Euler in 1736 (and an
earlier, identical proof was discovered in the unpublished manuscripts
of
Leibniz). The most famous of Fermat's results—known as
Fermat's Last Theorem—was jotted down in 1637 in his copy of
the book Arithmetic (by the third-century Greek mathematician
Diophantus) with the remark I have discovered a truly remarkable
proof, but this margin is too small to contain it. Finding a proof
of Fermat's Last Theorem became one of the most famous challenges in
number theory. A complete
solution was finally given in 1995 by
Andrew Wiles of Princeton
University.
[3]
The reduction steps in the cases where the exponent
$e$ is greater than 1 are based on the fact that,
for any integers $x$,
$y$, and $m$, we can
find the remainder of $x$ times
$y$ modulo $m$ by
computing separately the remainders of $x$ modulo
$m$ and $y$ modulo
$m$, multiplying these, and then taking the
remainder of the result modulo $m$. For
instance, in the case where $e$ is even, we
compute the remainder of $b^{e/2}$ modulo
$m$, square this, and take the remainder modulo
$m$. This technique is useful because it means
we can perform our computation without ever having to deal with numbers much
larger than $m$. (Compare
exercise 1.25.)
[4]
Numbers that fool the
Fermat test are called
Carmichael numbers, and little is known
about them other than that they are extremely rare. There are 255
Carmichael numbers below 100,000,000. The smallest few are 561, 1105,
1729, 2465, 2821, and 6601. In testing primality of very large
numbers chosen at random, the chance of stumbling upon a value that
fools the Fermat test is less than the chance that
cosmic radiation will cause the computer to make an error in carrying out a
correct algorithm. Considering an algorithm to be inadequate
for the first reason but not for the second illustrates the difference
between
mathematics and engineering.
[5]
One of the most
striking applications of
probabilistic prime testing has been to the field of
cryptography.
Although it is now computationally infeasible to factor an arbitrary 200-digit
number, the primality of such a number can be checked in a few seconds with the
Fermat test.
Although it is computationally infeasible to factor an arbitrary 300-digit
number as of this writing (2021), the primality of such a number can be checked
in a few seconds with the Fermat test.
This fact forms the basis of a technique for constructing
unbreakable codes suggested by
Rivest,
Shamir, and
Adleman (1977). The resulting
RSA algorithm has become a widely used technique for enhancing the
security of electronic communications. Because of this and related
developments, the study of
prime numbers, once considered the epitome of a topic in pure
mathematics to be studied only for its own sake, now turns out to have
important practical applications to cryptography, electronic funds transfer,
and information retrieval.
[6]
This date is called the
UNIX epoch and is part of the specification of functions
that deal with time in the
UNIX$^{\textrm{TM}}$ operating
system.
[7]
The primitive
function display returns its
argument, but also prints it. Here
"***"
is a
string, a sequence of characters that we pass as argument
to the display function.
Section 2.3.1 introduces strings more
thoroughly.