?

Log in

No account? Create an account
In which Ben is evil to a poor little Haskell interpreter. - Adventures in Engineering — LiveJournal
The wanderings of a modern ronin.

Ben Cantrick
  Date: 2007-02-14 22:37
  Subject:   In which Ben is evil to a poor little Haskell interpreter.
Public
  Mood:der uber-nerd
  Music:Type O: Gravitational Constant G = 6.67 x 10^8 Cm^3/Gm Sec^2
  Tags:  haskell, prime numbers, programming



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)

Post A Comment | 7 Comments | | Link






A Queen Of The Shadow Realm
  User: mistress_wench
  Date: 2007-02-15 07:40 (UTC)
  Subject:   (no subject)
*swoons*
oooh I KNEW there was a reason I put you on the flist ;)**


**sarcasm not intended
Reply | Thread | Link



Trevor Stone: mathnet - to cogitate and to solve
  User: flwyd
  Date: 2007-02-15 15:44 (UTC)
  Subject:   (no subject)
Keyword:mathnet - to cogitate and to solve
I assume you want &lt; there instead of &gt;, i.e. 2 <= K < N. Otherwise, only negative numbers, 0, and 1 can be prime.

I've never read a Haskell tutorial (though I've installed an Eclipse plugin in a fit of "find every language Eclipse can develop in"), but I was able to read all of that. That's more than can be said for LISP.
Reply | Thread | Link



Ben Cantrick
  User: mackys
  Date: 2007-02-15 19:04 (UTC)
  Subject:   (no subject)
Ah, duh! Thanks for the correction! Fixed in the post now.
Reply | Parent | Thread | Link



cgibbard
  User: cgibbard
  Date: 2007-02-15 16:06 (UTC)
  Subject:   Hugs isn't built for speed either.
If you want your Haskell programs to run quickly, grab a copy of GHC (get a binary release, you can't compile it without GHC). Hugs doesn't even really try to do things quickly. I also don't know why you're getting illegal operation errors. Those are most likely bugs in Hugs.

From the screenshot, it looks like you're running windows, so you'd get that from http://haskell.org/ghc/download_ghc_66.html#windows

Compiled with optimisations using GHC 6.6, your program ran in 3.504s on 502560410469881 and 33.078s on 21167397384723757 on my 2.4 GHz P4.

Join us at #haskell on Freenode if you haven't yet! Newbie questions are quite welcome there. :)
Reply | Thread | Link



Ben Cantrick
  User: mackys
  Date: 2007-02-15 19:08 (UTC)
  Subject:   Re: Hugs isn't built for speed either.
If you want your Haskell programs to run quickly, grab a copy of GHC (get a binary release, you can't compile it without GHC).

I'll try that, thanks!

Hugs doesn't even really try to do things quickly. I also don't know why you're getting illegal operation errors. Those are most likely bugs in Hugs.

I have another theory - it doesn't optimize tail recursion, and so it's allocating a stack frame for every function call.

Join us at #haskell on Freenode if you haven't yet! Newbie questions are quite welcome there. :)

Cool. I was thinking there must be an IRC channel out there where I could ask how to do integer square root... ended up having to google for it instead.
Reply | Parent | Thread | Link



  User: (Anonymous)
  Date: 2007-02-16 02:22 (UTC)
  Subject:   Re: Hugs isn't built for speed either.
Oh, Hugs certainly would optimise tail recursion. It's something which is required for Haskell implementations. However, due to lazy evaluation, you can still end up with tail recursions that while implemented as tight loops, build huge unevaluated values whose evaluation creates a stack overflow. It's possible that's happening. I'll have a look at it at some point.
Reply | Parent | Thread | Link



Ben Cantrick
  User: mackys
  Date: 2007-02-16 03:56 (UTC)
  Subject:   Re: Hugs isn't built for speed either.
Do you want me to enter a bug in the Hugs bug tracker?
Reply | Parent | Thread | Link



browse
May 2015