SICP Exercise 1.24: Faster Prime Time

Modify the timed-prime-test procedure of exercise 1.22 to use fast-prime? (the Fermat method), and test each of the 12 primes you found in that exercise. Since the Fermat test has Θ(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?

Okay, for this one I've gone back to the original procedure definitions of exercise 1.22. I'll run timed-prime-test directly instead of searching for primes.

Now the implementation of fast-prime? executes the fermat-test the number of times indicated by the second operand. This is useful for us as we've already seen that we need to test the numbers a large number of times in order for us to be able to get useful timing data. After a wee bit of experimentation I opted to use 250,000 times. Similar to exercise 1.23, the exercise here says to "Modify the timed-prime-test procedure of exercise 1.22 to use fast-prime?", but it's actually start-prime-test that was making the call to prime? in the original implementation, and so is the procedure we actually modify. Here's my updated start-prime-test procedure:

(define (start-prime-test n start-time)
  (if (fast-prime? n 250000)
      (report-prime (- (runtime) start-time))
      (display "")))

And here's what happens when I run the tests:

> (timed-prime-test 1009)
1009 *** 1588.0
> (timed-prime-test 1013)
1013 *** 1663.0
> (timed-prime-test 1019)
1019 *** 1689.0
> (timed-prime-test 10007)
10007 *** 2064.0
> (timed-prime-test 10009)
10009 *** 1974.0
> (timed-prime-test 10037)
10037 *** 2048.0
> (timed-prime-test 100003)
100003 *** 6108.0
> (timed-prime-test 100019)
100019 *** 5984.0
> (timed-prime-test 100043)
100043 *** 5986.0
> (timed-prime-test 1000003)
1000003 *** 9727.0
> (timed-prime-test 1000033)
1000033 *** 8330.0
> (timed-prime-test 1000037)
1000037 *** 8772.0

Now the Fermat test has Θ(log n) growth so we can work out the expected ratios between the times taken for the first three primes larger than 1,000, and for the subsequent orders of magnitude. Also, we can take the average times for each of these tuples of primes and work out the actual ratios:

Start	1st Prime Time (ms)	2nd Prime Time (ms)	3rd Prime Time (ms)	Average Time (ms)	Expected Ratio	Actual Ratio
1,000	1588.0	1663.0	1689.0	1646.7	-	-
10,000	2064.0	1974.0	2048.0	2028.7	1.333	1.232
100,000	6108.0	5984.0	5986.0	6026.0	1.667	3.659
1,000,000	9727.0	8330.0	8772.0	8943.0	2.0	5.431

For the first three primes above 10,000 the ratio seems to be about right. But for 100,000 and 1,000,000 the ratios are much higher than expected. Why would this be?

A hint at an answer is perhaps given in the footnote to section 1.2.3:

"...if we count process steps as "machine operations" we are making the assumption that the number of machine operations needed to perform, say, a multiplication is independent of the size of the numbers to be multiplied, which is false if the numbers are sufficiently large."

I've noted before that I'm using DrRacket for most of my SICP work. It so happens that I'm using a 32-bit version of DrRacket. The user manual indicates that it uses 30-bits plus one sign bit to represent fixnums, and if a value is too large to represent as a fixnum then it will be represented as a scheme_bignum_type - an arbitrarily large integer value. Now 2³⁰ = 1073741824, so provided the calculation of expmod remains within this range DrRacket will be able to do all the calculations using fixnums and so should be nice and efficient. However, if we look at the procedure expmod we see that it contains a line where it squares the results of a recursive call to expmod. If the result returned by this recursive call to expmod exceeds √2³⁰ = 2¹⁵ = 32,768 then DrRacket will have to represent its square as a scheme_bignum_type, which is likely to have inferior performance to fixnums for calculations. So when can this happen?

Well, the result of this recursive call will be in the range 0..n, where n is the number we are testing the primality of. So it can't happen for the first three primes larger than either 1,000 or 10,000, but it can happen for the first three primes larger than 100,000 and 1,000,000. Whether it does or not will depend upon the random number chosen for the test. However, as we are running the fermat test 250,000 times with random numbers in the range 0..n, chances are we're going to hit on quite a few that result in a recursive call to expmod returning such a value. As soon as it does that, the interpreter has to use the less efficient representation and so the performance degrades.