Re: Elementary group theory: Proof of Fermat-Maas primality-test (was: correcting *** ...)

From: "*** T. Winter" <***.Win...@xxxxxx>
The program I am using for primality tests initially does trial
division by primes upto 100,000; next it performs Fermat for a
single number; ...

I got to thinking about that, whether it's really cost effective to
do trial division up to such a high divisor. So I did some timing
tests here on the Unix shell system. I measured repeat loops of (1)
modular multiplication of two bignums in the same range of 10**15
upward to 10**100, (2) modular exponentation with all three bignums
in the same range as before, and (3) trial division of bignum in
that same range divided by single-word integer. Here are the three
test numbers I used, only the power of ten varied between different
bignum-size timings. The high-order digits shown correspond to
bignums in the ballpark of 1E50 to 4E50 (except single-float
doesn't allow such large exponents, as I discovered later when I
tried to take logarithms to calculate density of primes, so you'd
need to specify these as 1D50 and 4D50 respectively).

(setq bigval1 (isqrt (* 2 (expt 10 200)))) ;Starting value, changed by each
; calculation in *mod or exptmod loop
;141421356237309504880168872420969807856967187537694...
(setq bigval2 (isqrt (* 2 (expt 10 201)))) ;Modulus (trial prime), constant
;447213595499957939281834733746255247088123671922305...
(setq bigval3 (isqrt (* bigval1 bigval2))) ;Exponent (not used for trial division)
;2514866859365870816635531009317715676551936728029772...

internal-time-units-per-second ;100
That's what it says, the ratio between compute units given by
get-internal-run-time and actual clock time in seconds.

Here's the timing loop for trial division:
(progn (gc) (setq t0 (get-internal-run-time))
(let ((val1 bigval1) (val3 bigval3))
(dotimes (ix 50000) (setq val1 (mod val3 (+ 1 ix))))
(setq t1 (get-internal-run-time)) (- t1 t0)))
I did the *mod and exptmod timings earlier, and stripped the code
to make the trial division timing loop, sigh, instead of making a
backup copy. But I think the earlier versions said:
... (*mod val1 val1 val2)...
... (exptmod val1 val3 val2)...
also different number instead of 50000 for number of iterations,
and binding for val2, otherwise identical code for all three tests.

Don't complain about details. This is just to get a ballpark idea
how long each of the three kinds of operations typically takes
relatively to each other.

So here are the timing results:

;LoopRep op 10**? timings in hundredths of second time per operation
;30000 *mod E15: 197 197 200+-2 201 202 ;(/ 2.00 30000)=6.67e-5
;30000 *mod E25: 206 207 209+-2 209 215 ;(/ 2.09 30000)=6.97e-5
;30000 *mod E50: 236 237 238+-2 240 244 ;(/ 2.38 30000)=7.93e-5
;30000 *mod E100: 291 293 294+-2 295 304 ;(/ 2.94 30000)=9.80e-5
Notice that bignum multiplication and modulus (remainder) doesn't
vary much with size of numbers. I don't know why. Maybe this
version of CMU Common Lisp uses FFT algorithms??

;LoopRep op 10**? timings in hundredths of second time per operation
;1000 exptmod E15: 535 536 536+-4 545 553 ;(/ 5.36 1000)=0.00536
;1000 exptmod E25: 891 911 914+-3 916 919 ;(/ 9.14 1000)=0.00914
;1000 exptmod E50: 2066 2082 2084+-10 2103 2108 ;(/ 20.84 1000)=0.02084
(It was taking a half minute realtime, so I switched to fewer times
through the loop. Loop overhead seems to have increased the per-op
cost by about 1 or 2 percent, nothing to worry about when the timing
varies between runs by that much anyway.)
;500 exptmod E50: 1035 1048 1057+-5 1058 1065 ;(/ 10.57 500)=0.02114
;500 exptmod E100: 2438 2443 2459+-21 2486 2488 ;(/ 24.59 500)=0.04918
(Notice time per operation is slightly worse than linear in number
of digits, due to 1.5 * number of digits *mod-ops needed for each
one exptmod-op, and near-constant time for single *mod-op.)

;LoopRep op 10**? timings in hundredths of second time per operation
;50000 /int E15: 201 203 204+-1 205 206 ;(/ 2.04 50000)=4.08e-5
;50000 /int E25: 202 202 204+-1 204 206 ;(/ 2.04 50000)=4.08e-5
;50000 /int E50: 210 210 211+-2 213 220 ;(/ 2.11 50000)=4.22e-5
;50000 /int E100: 206 208 209+-1 210 210 ;(/ 2.09 50000)=4.18e-5
(This surprized me: A single division to get remainder, of bignum
by fixnum, takes almost as long as *mod which is a bignum*bignum
multiplication followed by a bignum/bignum division to get
remainder. Maybe the loop overhead dominates over bignum/fixnum
division. But this would happen in a real program too, so I'm
sticking with these timing figures for now.)


Now we ignore the *mod timings, and compare timings of exptmod with /int:
Number of bignum/int that can be done for cost of just one exptmod(bignums):
E15: (/ 0.00536 4.08e-5)= 131
E25: (/ 0.00914 4.08e-5)= 224
E50: (/ 0.021 4.22e-5)= 498
E100: (/ 0.04918 4.18e-5)=1176

But note the expectation of success (learning number is not prime and
thereby avoiding all further work) with Fermat's exptmod, is the
density of composites within all odd integers in the range of
interest, which follows the Prime Number Theorem:
Interval Number primes Density primes in odds Density composites in odds
[1E15, 3E15] 5.679537e+13 0.05679537 0.94320464
[1E25, 3E25] 3.4347323e+23 0.034347318 0.9656527
[1D50, 3D50] 1.7273...d+48 0.01727337333499136d0 0.9827266266650087d0
[1D100, 3D100] 8.6613...d+97 0.008661368757353492d0 0.9913386312426465d0
But as you see, that density is essentially equal to 1.0.

Whereas expectation of success is 1/n when using n as trial divisor.
(Assuming trial divisors only by primes, a bit worse if using all odd
numbers as trial divisors to avoid the overhead of keeping a list of
lots of primes to test.)

So for trial divisors smaller than numbers above (131 ... 1176), it's
more cost-effective to do another trial division and hope to get
lucky, whereas for trial divisors larger than those numbers above,
it's better to stop trial division and perform Fermat's test
immediately for better benefit/cost ratio. As a simple
easy-to-remember rule of thumb, stop trial divisors at one thousand
for testing large primes (100 decimal digits) typical in RSA
cryptography, maybe two thousand if going for 200-digit primes.

Note: I failed to take into account density of Carmichael numbers
and other false primes where Fermat test fails to catch the
composite. This would lower the effectiveness of the Fermat test a
little and thereby increase the cutoff point for trial divisors.
But the density of Carmichael numbers is extremely low, much lower
than density of prime numbers which are already less than one
percent in the 1E100 region, and other false primes have only small
chance of randomBase**(composite-1) =modComposite= 1, so this isn't
a big enough effect to worry about, right?
.