Re: finding primes

21 Dec 2000

      On Thu, 21 Dec 2000, S.D.Mechveliani wrote:
...
Hello,
On generating prime numbers, people wrote
...
|     import Array
|     primes n = [ i | i <- [2 ..n], not (primesMap ! i)] where
| 	primesMap   = accumArray (||) False (2,n) multList
| [..]
I think that it is good for functional programming to avoid arrays, 
as possible.
Yes. So I realized. It's a pity though, because sometimes arrays are very 
convenient.
...
Shlomi Fish   writes
...
[..]
primes :: Int -> [Int]
primes how_much = sieve [2..how_much] where
        sieve (p:x) = 
            p : (if p <= mybound
                then sieve (remove (p*p) x)
                else x) where
            remove what (a:as) | what > how_much = (a:as)
                               | a < what = a:(remove what as)
                               | a == what = (remove (what+step) as)
                               | a > what = a:(remove (what+step) as)
            remove what [] = []
            step = (if (p == 2) then p else (2*p)) 
        sieve [] = []
        mybound = ceiling(sqrt(fromIntegral how_much))
I optimized it quite a bit, but the concept remained the same.
Anyway, this code can scale very well to 100000 and beyond. But it's not
exactly the same algorithm.
[..]
C.Runciman   gives a paper reference 
about this.
Aslo may I ask what do you mean by "can scale to 100000",
the value of a prime or its position No in the list?
I mean it can easily generate the list of primes from 2 up to 100000. Not
100000 separate primes.
...
Anyway, here are my attempts:
-----------------------------------
1.
  primes1 = s [(2::Int)..] :: [Int]
                     where  s (p:ns) = p: (s (filter (notm p) ns))
                            notm p n = (mod n p) /= 0
main = putStr $ shows (primes1!!9000) "\n"
This strikes me as the sieve version which I gave as the first example.
It's not very efficient because it uses modulo which is a costy operation.
...
After compiling by  GHC-4.08 -O2
it yields           93187  in  115 sec   on  Intel-586, 160 MHz.
-----------------------------------
...
2.
The DoCon program written in Haskell (again, no arrays) gives
about 10 sec (for Integer values). 
Also it finds, for example, first 5 primes after 10^9 as follows:
take 5 $ filter isPrime [(10^9 ::Integer) ..] 
  -->
  [1000000007,1000000009,1000000021,1000000033,1000000087]
This takes  0.05 sec.
But DoCon uses a particular  isPrime  test method:
  Pomerance C., Selfridge J.L., Wagstaff S.S.:
  The pseudoprimes to 25*10^9.
  Math.Comput., 1980, v.36, No.151, pp.1003-1026.
After 25*10^9 it becomes again, expensive.
Where can I find this DoCon program?

Regards,

	Shlomi Fish
...
------------------
Sergey Mechveliani
mechvel@botik.ru
----------------------------------------------------------------------
Shlomi Fish        shlomif@vipe.technion.ac.il 
Home Page:         http://t2.technion.ac.il/~shlomif/
Home E-mail:       shlomif@techie.com

The prefix "God Said" has the extraordinary logical property of 
converting any statement that follows it into a true one.