In section 1.3 of SICP, a method for determining prime numbers is mentioned. The naive prime number determination, known to students around the globe, is if a number n's smallest divisor other than 1 is itself, then n is a prime number.
(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))
;return the smallest divisor of n
(define (smallest-divisor n)
;find a divisor(for 2 to n) of n
(define (find-divisor n test-divisor)
(cond ((< n (sqrt test-divisor)) n)
((= (remainder n test-divisor) 0) test-divisor)
(else (find-divisor n (+ test-divisor 1)))))
(find-divisor n 2))
;check n is prime or not
(define (prime? n)
(= n (smallest-divisor n )))
However, this method is slow, with a complexity of O(n).
SICP then introduces Fermat's Little Theorem:
If a number n is a prime, and for any a<n, a^n≡a(mod n)
Unfortunately, the converse theorem does not hold, but it does not mean it is not useful. We can select a number a smaller than n and perform the aforementioned test. If a^n≡a(mod n) holds, then n is very likely a prime number. The more a's tested against n, the higher the likelihood n is a prime. Thus, it's a probabilistic algorithm; with enough tests, the probability of n not being a prime becomes very small.
This leads us to the fermat test program. Don't worry about the efficiency of calculating a^n, because it's easily solved:
Consider:
a^n = (a^(n/2))^2; when n is even
a^n = a* a^(n-1); when n is odd
Always holds, so with nearly every multiplication, n doubles. If you look at the number line, this is essentially binary search. Thus, the time complexity is reduced. This explains why the pow function provided by programming languages is always faster than writing a*a*a… yourself.
But we shouldn't directly use the built-in pow function. The reason is our goal is not the nth power of a, but to know if the nth power of a modulo m equals a modulo m. Simplified, we only need to verify if this expression holds: a^n mod m = a. If expressed directly as pow(a, n) % m == a, it would be problematic because the value of pow(a, n) can be very large. So, we should use a similar technique:
a^n mod m = ((a^(n / 2))^2) mod m; when n is even
a^n mod m = (a* a^(n-1)) mod m; when n is odd
This way, the calculated value never exceeds m, avoiding computations with excessively large numbers. Hence, we easily have the following program:
(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))
(define (expmod base exps m)
(define (square n) (* n n))
(cond ((= exps 0) 1)
((even? exps)
(remainder (square (expmod base (/ exps 2) m))
m))
(else
(remainder (* base (expmod base (- exps 1) m))
m))))
;fermat test, check n. if (a^n) mod n == a then in all probability n is a prime
(define (fermat-test n)
(define (try-it a)
(= a (expmod a n n)))
;(= a (remainder (expt a n) n)))
(try-it (+ 1 (random (- n 1)))))
;get prime by fermat-test. test n is prime or not, try times time
(define (fermat-prime? n times)
(cond ((= 0 times) #t)
((fermat-test n) (fermat-prime? n (- times 1)))
(else #f)))
However, there's always imperfection.
There exists a type of numbers called Carmichael numbers that can trick the Fermat test. Between 2 and 10000, there are several such numbers:
561, 1105, 1729, 2465, 2821, 6601, 8911
These stubborn numbers are composite, but no matter what a you test them with, they can masquerade as primes in the Fermat test. They are not
few in number and are widely distributed, with a total of 255 Carmichael numbers within 100,000,000. So, if you do not wish to maintain a Carmichael number list, you can strengthen our prime number test method.
Fortunately, in exercise 1.28, Miller and Rabin brought us the Miller-Rabin test method. This method is based on a variation of Fermat's little theorem:
If n is a prime, for any a < n, a^(n-1) ≡ 1(mod n)
To use the MR test, one needs to check for a "non-trivial square root of 1 mod n" during the square step of the Fermat test. If there exists an s = a^2; and s <> 1 and s <> n-1, and s^2 mod m = 1, then n is not a prime.
Matrix67: The Miller-Rabin primality test is also a probabilistic algorithm, and we call composites that pass the Miller-Rabin test with base a as strong pseudoprimes to base a. The first strong pseudoprime to base 2 is 2047. The first strong pseudoprime to bases 2 and 3 is as large as 1,373,653.
So, when you have sufficient test quantity, this method can guarantee safety in industrial applications. Below is my answer to question 1.28:
(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))
(define (mr-expmod base exps m)
(define (square x) (* x x))
(define (check-square-root x n)
(define (check r)
(if (and (not (= x 1)) (not (= x (- n 1))) (= r 1)) 0 r))
(check (remainder (square x) n)))
(cond ((= exps 0) 1)
((even? exps)
(check-square-root(mr-expmod base (/ exps 2) m) m))
(else
(remainder (* base (mr-expmod base (- exps 1) m))
m))))
;Miller-Rabin test, check n. if (a^(n-1)) mod n == 1 then in all probability n is a prime
(define (mr-test n)
(define (try-it a)
(= 1 (mr-expmod a (- n 1) n)))
(try-it (+ 1 (random (- n 1)))))
;get prime by Miller-Rabin-test. test n is prime or not, try times time
(define (mr-prime? n times)
(cond ((= 0 times) #t)
((mr-test n) (mr-prime? n (- times 1)))
(else #f)))
I compiled some data, seeking primes in the intervals [1000, 1100], [10000, 10100], [100000, 100100], and [1000000, 1000100], to see how these methods perform:
Conventional method Conventional method's divisor 2 optimization Miller-Rabin test method
[1000, 1100] 1 ms 1 ms 1 ms
[10000, 10100] 13 ms 8 ms 1 ms
[100000, 100100] 142 ms 78 ms 1 ms
[1000000, 1000100] 1337 ms 830 ms 1 ms
This was run on my Intel Core 1 1.83GHz machine (Ubuntu 8.04 + DrScheme).
I also wrote a C version, with the following results:
Range [2, 10000]:
$ ./prime
Normal: return 0 cost 10ms
Fermat: return 0 cost 24ms
Miller-Rabin: return 0 cost 28ms
Range [200000, 1000000]:
$ ./prime
Normal: return 0 cost 3088ms
Fermat: return 0 cost 824ms
Miller-Rabin: return 0 cost 959ms
The gprof analysis results are mostly as expected, with a portion of the report below:
Each sample counts as 0.01 seconds.
% cumulative self self total
time seconds seconds calls us/call us/call name
63.55 2.58 2.58 800000 3.23 3.23 is_prime(int)
15.76 3.22 0.64 800005 0.80 0.84 fermat_expmod(long, long, long)
10.10 3.63 0.41 800012 0.51 0.88 mr_expmod(long, long, long)
7.27 3.92 0.29 14813786 0.02 0.02 chk_unusual_square_root(long, long)
1.23 3.98 0.05 main
0.86 4.01 0.04 800000 0.04 0.93 mr_test(long, long)
0.74 4.04 0.03 14813658 0.00 0.00 square(long)
0.25 4.05 0.01 800000 0.01 0.94 mr_is_prime(int)
0.25 4.06 0.01 800000 0.01 0.85 fermat_is_prime(int)
0.00 4.06 0.00 800000 0.00 0.84 fermat_test(long, long)
The reason for poor performance with small n values, I guess, is due to the overhead of handling ordinary recursion in C. However, as the data size increases, this overhead becomes relatively insignificant. If someone could convert the above algorithms into iterative form (which I cannot do :P) or tail-recursion form (gcc can optimize tail recursion into iterative form), I expect the performance would be even more pronounced.
