home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / h / hugs / hs / mersenne < prev    next >
Encoding:
Text File  |  1995-02-14  |  3.1 KB  |  72 lines
  1. -----------------------------------------------------------------------------
  2. -- Mersenne.hs                                                 Mark P. Jones
  3. --                                                          February 7, 1995
  4. --
  5. -- Here is a Hugs program to calculate the 30th Mersenne prime using the
  6. -- builtin bignum arithmetic.
  7. --
  8. -- For those who don't know, a Mersenne prime is a prime number of the form:
  9. --
  10. --                               n
  11. --                              2  - 1
  12. --
  13. -- The first few Mersenne primes are for:
  14. --   n = 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203,
  15. --       2281, 3217, 4253, 4423, ...
  16. --
  17. -- The thirtieth Mersenne prime occurs for n = 216,091.  In decimal
  18. -- notation, this number has 65050 digits.
  19. --
  20. --
  21. -- A little story about me and this number:
  22. --
  23. -- As I recall, this fact was discovered nearly ten years ago.  I
  24. -- wrote an Intel 8080 assembly language program to calculate this
  25. -- number.  Running on a Z80A based machine, it used a 32K array --
  26. -- more than half of the total memory available -- with each byte
  27. -- containing two binary coded decimal digits.   The array was
  28. -- initialized to contain the number 1 and a loop was used to double
  29. -- the value in the array a total of 216091 times, before the final 1
  30. -- was subtracted.  Using the timings for individual Z80A
  31. -- instructions, I estimated the running time for the program and,
  32. -- when it finished on Thursday April 17, 1986, after running for a
  33. -- little under 18 hours, I was delighted that my predictions were
  34. -- within 10 seconds of the actual running time.  Of course, now I
  35. -- understand a little more about error bounds and tolerances, I realize
  36. -- that this was more by luck than judgement, but at the time, I was
  37. -- delighted!  I don't remember if I knew the O(log n) algorithm for
  38. -- exponentials at the time, but it wouldn't have been easy to apply
  39. -- with the limited amount of memory at my disposal back then.  (Of
  40. -- course, it wouldn't have been O(log n) in practice either because
  41. -- the individual multiplications can hardly be considered O(1)!)
  42. --
  43. -- Now I can run this program, written in Hugs (or to be accurate,
  44. -- written using calls to Hugs primitive functions), on the machine
  45. -- on my desk while I'm editing files and reading mail in other
  46. -- windows, and it still finishes in under 7 minutes.  Of course,
  47. -- it did use 6M of heap (though not all at the same time), but
  48. -- who's counting?  :-)
  49.  
  50. p         :: Integer
  51. p          = 2 ^ 216091 - 1
  52.  
  53. digitsInP :: Integer
  54. digitsInP  = genericLength (show p)
  55.  
  56. -- Here are the smaller Mersenne primes listed above:
  57.  
  58. smallMPindices :: [Int]
  59. smallMPindices  = [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127,
  60.                    521, 607, 1279, 2203, 2281, 3217, 4253, 4423 ]
  61.  
  62. smallMP  :: [Integer]
  63. smallMP   = [ 2 ^ n - 1 | n <- smallMPindices ]
  64.  
  65. -- Does an incremental algorithm buy us anything?  Not much, it would seem!
  66.  
  67. smallMP' :: [Integer]
  68. smallMP'  = map (subtract 1) (scanl (\x i -> x * 2^i) (2^n) ns)
  69.             where (n:ns) = zipWith (-) smallMPindices (0:smallMPindices)
  70.  
  71. -----------------------------------------------------------------------------
  72.