fork(1) 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 prmstate = struct val p:uint32 val pi:byte new (prm,pndx) = { p = prm; pi = pndx } end
  29. type prmsSeqDesc = struct val v:prmstate val cont:unit->prmsSeqDesc new(ps,np) = { v = ps; cont = np } end
  30. type cmpststate = struct val cv:uint32 val ci:byte val cp:uint32 new (strt,ndx,prm) = {cv = strt;ci = ndx;cp = prm} end
  31. type cmpstsSeqDesc = struct val v:cmpststate val cont:unit->cmpstsSeqDesc new (cs,nc) = { v = cs; cont = nc } end
  32. type allcmpsts = struct val v:cmpstsSeqDesc val cont:unit->allcmpsts new (csd,ncsdf) = { v=csd;cont=ncsdf } end
  33.  
  34. let primesTFOWSE =
  35.   let whlptrn = [| 2uy;4uy;2uy;4uy;6uy;2uy;6uy;4uy;2uy;4uy;6uy;6uy;2uy;6uy;4uy;2uy;6uy;4uy;6uy;8uy;4uy;2uy;4uy;2uy;
  36.                    4uy;8uy;6uy;4uy;6uy;2uy;4uy;6uy;2uy;6uy;6uy;4uy;2uy;4uy;6uy;2uy;6uy;4uy;2uy;4uy;2uy;10uy;2uy;10uy |]
  37.   let inline whladv i = if i < 47uy then i + 1uy else 0uy
  38.   let inline advmltpl c ci p = cmpststate (c + uint32 whlptrn.[int ci] * p,whladv ci,p)
  39.   let rec pmltpls cs = cmpstsSeqDesc(cs,fun() -> pmltpls (advmltpl cs.cv cs.ci cs.cp))
  40.   let rec allmltpls (psd:prmsSeqDesc) =
  41.     allcmpsts(pmltpls (cmpststate(psd.v.p*psd.v.p,psd.v.pi,psd.v.p)),fun() -> allmltpls (psd.cont()))
  42.   let rec (^) (xs:cmpstsSeqDesc) (ys:cmpstsSeqDesc) = //union op for SeqDesc's
  43.     match compare xs.v.cv ys.v.cv with
  44.       | -1 -> cmpstsSeqDesc (xs.v,fun() -> xs.cont() ^ ys)
  45.       | 0 -> cmpstsSeqDesc (xs.v,fun() -> xs.cont() ^ ys.cont())
  46.       | _ -> cmpstsSeqDesc(ys.v,fun() -> xs ^ ys.cont()) //must be greater than
  47.   let rec pairs (csdsd:allcmpsts) =
  48.     let ys = csdsd.cont in
  49.     allcmpsts(cmpstsSeqDesc(csdsd.v.v,fun()->csdsd.v.cont()^ys().v),fun()->pairs (ys().cont()))
  50.   let rec joinT3 (csdsd:allcmpsts) = cmpstsSeqDesc(csdsd.v.v,fun()->
  51.     let ys = csdsd.cont() in let zs = ys.cont() in (csdsd.v.cont()^(ys.v^zs.v))^joinT3 (pairs (zs.cont())))
  52.   let rec mkprimes (ps:prmstate) (csd:cmpstsSeqDesc) =
  53.     let nxt = ps.p + uint32 whlptrn.[int ps.pi]
  54.     if ps.p >= csd.v.cv then mkprimes (prmstate(nxt,whladv ps.pi)) (csd.cont()) //minus function
  55.     else prmsSeqDesc(prmstate(ps.p,ps.pi),fun() -> mkprimes (prmstate(nxt,whladv ps.pi)) csd)
  56.   let rec baseprimes = prmsSeqDesc(prmstate(11u,0uy),fun() -> mkprimes (prmstate(13u,1uy)) initcmpsts)
  57.   and initcmpsts = joinT3 (allmltpls baseprimes)
  58.   let genseq sd = Seq.unfold (fun (psd:prmsSeqDesc) -> Some(psd.v.p,psd.cont())) sd
  59.   seq { yield 2u; yield 3u; yield 5u; yield 7u; yield! mkprimes (prmstate(11u,0uy)) initcmpsts |> genseq }
  60.  
  61. primesTFOWSE |> Seq.nth 100000 |> printfn "%A"
  62.  
Success #stdin #stdout 0.5s 12760KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1299721u
https://ideone.com/RwC66u
Optimized Tree Folding Wheel based SoE, corrected
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!