finding primes

S.D.Mechveliani mechvel@math.botik.ru
Thu, 21 Dec 2000 09:35:54 +0300


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. 


Shlomi Fish <shlomif@vipe.technion.ac.il>  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 <Colin.Runciman@cs.york.ac.uk>  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?
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"

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.


------------------
Sergey Mechveliani
mechvel@botik.ru