Ben Cantrick (mackys) wrote,
Ben Cantrick
mackys

  • Mood:
  • Music:

In which Ben is evil to a poor little Haskell interpreter.




So I decided I wanted to dork around with Haskell, just for grins. I downloaded WinHugs and started looking at a tutorial.

I got as far as the recursive definition of prime (a function to test if a number is prime or not) before my evil side manifested. You can see the code for this function in the lower-right of the screen shot. It probably isn't obvious to any non-Haskeller, but the algorithm the tutorial gave was:

"A number N is non-prime if there is another number K such that: 2 < K < N, and N is evenly divisible by K." In other words, the approach was to take your number, N, and try and integer divide it by every number from 2 to N-1 and see if any of them divided it evenly!

Now I quite understand this was only for tutorial purposes, and nobody would actually implement a prime test function this way. Still, it was screechingly obvious to me that this algorithm was no AKS primality test. Heck, it wasn't even a Sieve of Eratosthenes! I couldn't resist making the obvious optimization of only checking from 2 up to (N/2). (I know you only strictly need to check up to sqrt(N), but I don't know how to do integer square roots in Haskell yet.) Even so I was pretty sure this thing was going to be dog-slow. Unsure of exactly HOW slow it was going to be on my "mere" 850 MHz Pentium 3, and also unsure of what format Haskell uses internally for integers, I decided to start off with a few random tests.

So I started typing in large random numbers, just for kicks. About the third time I tried it, WinHugs crashed and burned with an "Illegal Operation" error. So much for the theoretical inherent superiority of functional languages!! I have no idea why it crashed. I could have overflowed an integer somewhere and scribbled stack, or a million other things. My laughter at this point was loud enough to annoy the neighbors.

Throughly amused, I restarted WinHugs and found the excellent List Of Small Primes page, and started up again, timing the runs by counting cursor flashes, which I assume are 1 per second. You can see me first testing a few primes that test the boundaries around 16 bits, and then seeing that all seemed well, I kicked it up a notch...

So, want to know the punchline? It took it 78 seconds to tell me that 1,299,827 (the 100,000th prime) was actually prime. Interestingly enough, because of the way the algorithm works - attempting to check for even divisiblity by 2 as the very first test - it takes no perceptible time at all to tell you that 1,299,828 (or any other even number) is non-prime. If you look at the not(factors 2 x) line of the code, you'll see that this is basically a NON prime-number tester. In short, it's doing some math to every number less than N, in order to prove that N has a certain property. Not exactly the optimal way of doing things.

Still, it is pretty cool that even a rank beginner who just installed the interpreter five minutes ago can do a prime number tester in 5 lines of Haskell. (Not that this would be a particularly hard algorithm to implement in C, natch.)


There are a couple more obvious optimizations I should make to the code I have before I toss it out and implement a better algorithm. First, as mentioned before, only check numbers from 2 to sqrt(N). Secondly, only test the odd numbers in that range, since any even number is by definition also divisible by 2 and we already checked that.

If I get ambitious tonite, I may figure out how to pass a list to a function, and implement the Sieve of Eratosthenes, or maybe even the Sieve of Atkin.


Edit: Holy crap, man! I did the two obvious optimizations, and now it figures out that 1299827 is prime before my finger comes off the Enter key! Same with 15485863! There's a barely perceptible delay, maybe 150 ms, to confirm the 15 millionth prime, 275604541. At this point I'm suspecting it's broken and I need to find a huge composite number to test it with...

Ahahahaha!! Tried prime 70610771126426561 (= 256,203,221 * 275,604,541 ; the 14 millionth + 1 prime and the 15 millionth prime) and it illegal op crashed instantly!

Illegal op's after about 3 minutes on 502560410469881 (= 15,485,867 * 32,452,843 ; the 2nd millionth prime and the third millionth - 1 prime). Running out of memory?

Takes about half a second to confirm non-primality on 62615533 = 7907 * 7919. (These are the 9,999th and 10,000th primes.)

Four seconds for 10967535067 (= 104723 * 104729, the 9,999th and 10,000th primes).

51 seconds for 1689542430967 (= 1299821 * 1299827, the 100,007th and 100,008th primes).

I quit and restarted the Hugs window, then set it running on 502560410469881 again. It did it in 9:48. I think that's about the outer edge of what most people would be willing to sit there waiting for. So 15 digit composite numbers are about the limit. Very nice...

One last test, just for grins. 17 digit prime 21167397384723757. Illegal Op after just a few seconds less than an hour. I suspect some kind of garbage collection or memory issue when the last iteration of hasFactors returns and the call stack is being unwound.

FYI, code is now:

isqrt :: Integral a => a -> a
isqrt = floor . sqrt . fromIntegral

hasFactors sqrtn n m | m > sqrtn = False
                     | n `mod` m == 0 = True
                     | True = hasFactors sqrtn n (m+2)

prime :: Integral x => x -> Bool
prime x | x < 4 = True
        | x `mod` 2 == 0 = False
        | True = not(hasFactors (isqrt x) x 3)

Tags: haskell, prime numbers, programming
Subscribe
  • Post a new comment

    Error

    default userpic

    Your reply will be screened

    Your IP address will be recorded 

    When you submit the form an invisible reCAPTCHA check will be performed.
    You must follow the Privacy Policy and Google Terms of use.
  • 7 comments