[Edu-sig] Miller-Rabin (first draft, slow)

Tim Peters tim.one@home.com
Thu, 14 Dec 2000 22:38:34 -0500


[posted & mailed]

[Kirby Urner, with a supernaturally slow Miller-Rabin program]

Haven't had time to look at this in detail, but two things you need to pick
up if you're going to sling bigints in Python:

1)     1L << N
   computes the same result as
       2L ** N
   (for int N >= 0) but shifting is much faster than exponentiating.

2)     pow(i, j, k)
   computes the same result as
       (i**j) % k
   but can be unboundedly faster.

#2 is crucial.  To see why, compute, e.g., 2L**1000000 % 3 by hand,
computing 2L**1000000 first <wink>.  Then do it smarter:  "OK,
(2L**1000000)%3 == (2L**2 % 3)**500000 % 3 == 1L**500000 % 3 == 1L % 3 ==
1L".  That's the kind of savings 3-argument pow can buy you (it does mods
repeatedly under the covers, so that it never actually needs to do
arithmetic with numbers larger than k**2; and the latter can be unboundedly
smaller than i**j, so 3-argument pow can be unboundedly faster).

>     # if any of the primes is a factor, we're done
>     if 0 in map(lambda i,n=n,b=bases:
>                 primes.gcd(n,b[i])==1, range(len(bases))):
>        print "Definitely not prime"
>        return

*Try* writing a nice, readable Python loop -- it might grow on you <wink>.
Note that since you're using primes as bases, gcd(n,b[i]) == 1 is equivalent
to the simpler n % b[i] != 0.  So:

    for base in bases:
        if n % base == 0:
            print "Definitely not prime"
            return

from-the-convoluted-to-the-obvious-ly y'rs  - tim