2013/01/23

Fun with primes

Every year, on new year's holidays I tend to play with writing algorithms that aren't directly related to work. (In Japanese tradition, we write calligraphy in the beginning of each year, called 書き初め (first writing). Mine is a bit like that---first coding.)

Last year I played with procedural modeling using Markov-chain Monte Carlo (WiLiKi:Gauche:MetropolisProceduralModeling). Initially I planned to pick some graphics-related papers again. However, @cddddr wrote a nice entry of generating infinite prime sequence taking advantage of Gauche's generators and I thought that'd be something nice to have in Gauche, so I decided to add math.prime module and to throw in some utilities related to prime numbers.

First thing the module offer is an infinite lazy sequence of prime numbers *primes*. I optimize the original code by @cddddr a bit, notably to avoid using continuations; currently it can generate first 10^7 primes in 14 seconds on 2.4GHz Core2 machine.

gosh> (time (list-ref *primes* #e1e7))
;(time (list-ref *primes* 10000000))
; real  14.297
; user  14.060
; sys    0.240
179424691

Primality tests are also suitable to be in the math.prime module, right? I only knew Miller-Rabin probabilistic primality test, so I googled a bit to hunt around other options. I found out a few interesting things.

• Miller-Rabin tests doesn't need to be probabilistic. Up to certain numbers it is known to give definite answer with the selected "witness" values. Wikipedia lists up to 341550071728321---for a number below that we can say definitely it's a prime or not, by only using up to first 7 primes (2 to 17) as witnesses. If I increase the number of witnesses this bound can be higher, though I found out that the following BPSW test becomes faster than 8 or more rounds of Miller-Rabin tests.
• Baillie-PSW primality test is, like Miller-Rabin, a test that may tell a pseudoprime as a prime, but it is verified that no pseudoprime exist under 2^64, according to the page (the author of the page says independent check is needed, though). So up to that range it can be used as a deterministic test.
• Then there's AKS primality test, fastest known determinisic primality test.

I implemented Miller-Rabin and BPSW tests but haven't done AKS test; running out of time.

What else would be nice to have? How about prime factorization? This, again, is very interesting topic to search for various algorithms. I found This answer in stackoverflow and used it as a starting point. I ended up writing Monte-Carlo factorization (Brent's modification of Pollard's Rho alrogithm), as well as the naive "divide by every prime" approach.

There's something fascinating in prime factorization. I found myself blankly typing random numbers to factorize, again and again...

gosh> (mc-factorize 435728572856238754559)
(104579 1018217 4091958413)
gosh> (mc-factorize 984597945623658375937)
(79 12463265134476688303)

If the number goes around 20digits in decimal or so, sometimes Monte-Carlo factorization takes very long. Factoring 19740274219868223167 took several minutes.

gosh> (mc-factorize 19740274219868223167)
(2971215073 6643838879)

I ran mc-factorize on 7-th fermat number (2^128+1) for more than a day but it couldn't find factors.

I thought I'd stop here, but when I tried WolframAlpha it factorized those big numbers without sweat, in seconds. Ok, I'd proceed to Elliptic curve factorization... some day...

The module also provides Jacobi symbol (used in BPSW test), and also Euler's totient function. Following Wikipedia links reveals more and more interesting functions (Mobius function? What's a cool name!) but it'd be endless process so I had to stop, for now.

I have no prospect to use this module in any practical work so far---if you need serious calculation in these fields, there already exist more general and optimized packages like PARI/GP, so I think math.prime module has little practical importance. Nonetheless, it was fun to make those things run---and who knows someday it might come in handy.

Tags: math.prime, Prime