fork download
copy 
  1. //This work improves on the ideas of Jon Harrop as of
  2. //http://f...content-available-to-author-only...t.ca/2010/02/sieve-of-eratosthenes.html
  3. //which fails in terms of high memory use (saves all generated primes instead of
  4. //just the primes to the square root of the highest candidate prime, and continuously
  5. //regenerating buffer arrays of ever increasing huge size irrespective of the CPU cache sizes)
  6. //as well not including a generating sequence and not having the optimizations to
  7. //only deal with odd numbers let alone using wheel factorization
  8.  
  9. (*let primesJH() =
  10.   let primes =
  11.     let a = ResizeArray[2]
  12.     let grow() =
  13.       let p0 = a.[a.Count-1]+1
  14.       let b = Array.create p0 true
  15.       for di in a do
  16.         let rec loop i =
  17.           if i<b.Length then
  18.             b.[i] <- false
  19.             loop(i+di)
  20.         let i0 = p0/di*di
  21.         loop(if i0<p0 then i0+di-p0 else i0-p0)
  22.       for i=0 to b.Length-1 do
  23.         if b.[i] then a.Add(p0+i)
  24.     fun n ->
  25.       while n >= a.Count do
  26.         grow()
  27.       a.[n]
  28.   Seq.initInfinite id |> Seq.map primes
  29.  
  30. primesJH() |> Seq.nth 10000000  *)
  31.  
  32. //#time
  33.  
  34. let primesAPF() =
  35.   let rec oddprimes() =
  36.     let BUFSZ = 1<<<17 in let buf = Array.zeroCreate (BUFSZ>>>5) in let BUFRNG = uint32 BUFSZ<<<1
  37.     let inline testbit i = (buf.[i >>> 5] &&& (1u <<< (i &&& 0x1F))) = 0u
  38.     let inline cullbit i = let w = i >>> 5 in buf.[w] <- buf.[w] ||| (1u <<< (i &&& 0x1F))
  39.     let inline cullp p s low = let rec cull' i = if i < BUFSZ then cullbit i; cull' (i + int p)
  40.                                cull' (if s >= low then int((s - low) >>> 1)
  41.                                       else let r = ((low - s) >>> 1) % p in if r = 0u then 0 else int(p - r))
  42.     let inline cullpg low = //cull composites from whole buffer page for efficiency
  43.       let max = low + BUFRNG - 1u in let max = if max < low then uint32(-1) else max
  44.       let sqrtlm = uint32(sqrt(float max)) in let sqrtlmndx = int((sqrtlm - 3u) >>> 1)
  45.       if low <= 3u then for i = 0 to sqrtlmndx do if testbit i then let p = uint32(i + i + 3) in cullp p (p * p) 3u
  46.       else baseprimes |> Seq.skipWhile (fun p -> //force side effect of culling to limit of buffer
  47.           let s = p * p in if p > 0xFFFFu || s > max then false else cullp p s low; true) |> Seq.nth 0 |> ignore
  48.     let rec mkpi i low =
  49.       if i >= BUFSZ then let nlow = low + BUFRNG in Array.fill buf 0 buf.Length 0u; cullpg nlow; mkpi 0 nlow
  50.       else (if testbit i then i,low else mkpi (i + 1) low)
  51.     cullpg 3u; Seq.unfold (fun (i,lw) -> //force cull the first buffer page then doit
  52.         let ni,nlw = mkpi i lw in let p = nlw + (uint32 ni <<< 1)
  53.         if p < lw then None else Some(p,(ni+1,nlw))) (0,3u)
  54.   and baseprimes = oddprimes() |> Seq.cache
  55.   seq { yield 2u; yield! oddprimes() }
  56.  
  57. primesAPF() |> Seq.nth 10000000 |> printfn "%A"
Success #stdin #stdout 4.35s 12552KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
179424691u
https://ideone.com/l8zupF
Segmented functional Sieve of Eratosthenes not using mutable state other than the mutable buffer array contents; sieving odds in a packed bit array
language:
F# (mono 4.1)
created:
10 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!