fork download
copy 
  1. //Following is the Haskell code from which this algorithm was adapted:
  2. //from http://w...content-available-to-author-only...l.org/haskellwiki/Prime_numbers#Tree_merging_with_Wheel
  3. //As a reference, here is the fastest pure functional sieve in Haskell produced to date:
  4. //
  5. //primesTMWE = 2:3:5:7: gapsW 11 wheel (joinT3 $ rollW 11 wheel primes')
  6. //  where
  7. //    primes' = 11: gapsW 13 (tail wheel) (joinT3 $ rollW 11 wheel primes')
  8. // 
  9. //gapsW k ws@(w:t) cs@(c:u) | k==c  = gapsW (k+w) t u 
  10. //                          | True  = k : gapsW (k+w) t cs  
  11. //rollW k ws@(w:t) ps@(p:u) | k==p  = scanl (\c d->c+p*d) (p*p) ws 
  12. //                                      : rollW (k+w) t u 
  13. //                          | True  = rollW (k+w) t ps    
  14. //joinT3 ((x:xs): ~(ys:zs:t)) = x : union xs (union ys zs) //NOTE:  '~' means deferred evaluation    
  15. //                                   `union` joinT3 (pairs t)  
  16. //wheel = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:
  17. //        4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel
  18. //
  19. //union (x:xs) (y:ys) = case (compare x y) of 
  20. //           LT -> x : union  xs  (y:ys)
  21. //           EQ -> x : union  xs     ys 
  22. //           GT -> y : union (x:xs)  ys
  23. //union  xs     []    = xs
  24. //union  []     ys    = ys
  25. //
  26. //pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t
  27.  
  28. type CIS<'T> = struct val v:'T val cont:unit->CIS<'T> new(v,cont) = { v=v;cont=cont } end //Co-Inductive Steam
  29. let primesTFOWSE =
  30.   let WHLPRMS = [| 2u;3u;5u;7u |] in let FSTPRM = 11u in let WHLCRC = int (WHLPRMS |> Seq.fold (*) 1u) >>> 1
  31.   let WHLLMT = int (WHLPRMS |> Seq.fold (fun o n->o*(n-1u)) 1u) - 1
  32.   let WHLPTRN =
  33.     let wp = Array.zeroCreate (WHLLMT+1)
  34.     let gaps (a:int[]) = let rec gap i acc = if a.[i]=0 then gap (i+1) (acc+1uy) else acc
  35.                          {0..WHLCRC-1} |> Seq.fold (fun s i->
  36.                            let ns = if a.[i]<>0 then wp.[s]<-2uy*gap (i+1) 1uy;(s+1) else s in ns) 0 |> ignore
  37.     Array.init (WHLCRC+1) (fun i->if WHLPRMS |> Seq.forall (fun p->(FSTPRM+uint32(i<<<1))%p<>0u)
  38.                                   then 1 else 0) |> gaps;wp
  39.   let inline whladv i = if i < WHLLMT then i+1 else 0 in let inline advcnd c ci = c + uint32 WHLPTRN.[ci]
  40.   let inline advmltpl p (c,ci) = (c + uint32 WHLPTRN.[ci] * p,whladv ci)
  41.   let rec pmltpls p cs = CIS(cs,fun() -> pmltpls p (advmltpl p cs))
  42.   let rec allmltpls k wi (ps:CIS<_>) =
  43.     let nxt = advcnd k wi in let nxti = whladv wi
  44.     if k < ps.v then allmltpls nxt nxti ps
  45.     else CIS(pmltpls ps.v (ps.v*ps.v,wi),fun() -> allmltpls nxt nxti (ps.cont()))
  46.   let rec (^) (xs:CIS<uint32*_>) (ys:CIS<uint32*_>) = 
  47.     match compare (fst xs.v) (fst ys.v) with //union op for composite CIS's (tuple of cmpst and wheel ndx)
  48.       | -1 -> CIS(xs.v,fun() -> xs.cont() ^ ys)
  49.       | 0 -> CIS(xs.v,fun() -> xs.cont() ^ ys.cont())
  50.       | _ -> CIS(ys.v,fun() -> xs ^ ys.cont()) //must be greater than
  51.   let rec pairs (xs:CIS<CIS<_>>) =
  52.     let ys = xs.cont() in CIS(CIS(xs.v.v,fun()->xs.v.cont()^ys.v),fun()->pairs (ys.cont()))
  53.   let rec joinT3 (xs:CIS<CIS<_>>) = CIS(xs.v.v,fun()->
  54.     let ys = xs.cont() in let zs = ys.cont() in (xs.v.cont()^(ys.v^zs.v))^joinT3 (pairs (zs.cont())))
  55.   let rec mkprm (cnd,cndi,(csd:CIS<uint32*_>)) =
  56.     let nxt = advcnd cnd cndi in let nxti = whladv cndi
  57.     if cnd >= fst csd.v then mkprm (nxt,nxti,csd.cont()) //minus function
  58.     else (cnd,cndi,(nxt,nxti,csd))
  59.   let rec pCIS p pi cont = CIS(p,fun()->let (np,npi,ncont)=mkprm cont in pCIS np npi ncont)
  60.   let rec baseprimes() = CIS(FSTPRM,fun()->let np,npi = advcnd FSTPRM 0,whladv 0
  61.                                            pCIS np npi (advcnd np npi,whladv npi,initcmpsts))
  62.   and initcmpsts = joinT3 (allmltpls FSTPRM 0 (baseprimes()))
  63.   let inline genseq sd = Seq.unfold (fun (p,pi,cont) -> Some(p,mkprm cont)) sd
  64.   seq { yield! WHLPRMS; yield! mkprm (FSTPRM,0,initcmpsts) |> genseq }
  65.  
  66. primesTFOWSE |> Seq.nth 1000000 |> printfn "%A"
Success #stdin #stdout 5.38s 13168KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
15485867u
https://ideone.com/hgQHVj
Optimized Tree Folding Wheel based SoE, further optimized
language:
F# (mono 4.1)
created:
11 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!