fork download
copy 
  1. //The incremental Sieve of Atkin
  2. //adapted from the pseudo-code at the following link
  3. //but with several optimizations and adaptations so as
  4. //to make it run incrementally and functionally
  5. //http://e...content-available-to-author-only...a.org/wiki/Sieve_of_Atkin
  6.  
  7. type cndstate = class val c:uint32 val wi:byte val md12:byte new(cnd,cndwi,mod12) = { c=cnd;wi=cndwi;md12=mod12 } end
  8. type prmsCIS = class val p:uint32 val cont:unit->prmsCIS new(prm,nxtprmf) = { p=prm;cont=nxtprmf } end
  9. type stateCIS<'b> = class val v:uint32 val a:'b val cont:unit->stateCIS<'b> new(curr,aux,cont)= { v=curr;a=aux;cont=cont } end
  10. type allstateCIS<'b> = class val ss:stateCIS<'b> val cont:unit->allstateCIS<'b> new(sbstrm,cont) = { ss=sbstrm;cont=cont } end
  11.  
  12. let primesTFWSA() =
  13.   let count = ref 0u
  14.   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;
  15.                    4uy;8uy;6uy;4uy;6uy;2uy;4uy;6uy;2uy;6uy;6uy;4uy;2uy;4uy;6uy;2uy;6uy;4uy;2uy;4uy;2uy;10uy;2uy;10uy |]
  16.   let rec prmsqrs v sqr = stateCIS(v,sqr,fun() -> let n=v+sqr+sqr in let n=if n<v then 0xFFFFFFFFu else n in prmsqrs n sqr)
  17.   let rec allsqrs (prms:prmsCIS) = let s = prms.p*prms.p in allstateCIS(prmsqrs s s,fun() -> allsqrs (prms.cont()))
  18.   let rec qdrtc v y = stateCIS(v,y,fun() -> let a=(y+1)<<<2 in let a=if a<=0 then (if a<0 then -a else 2) else a
  19.                                             let vn=v+uint32 a in let vn=if vn<v then 0xFFFFFFFFu else vn in qdrtc vn (y+2))
  20.   let rec allqdrtcsX4 x = allstateCIS(qdrtc (((x*x)<<<2)+1u) 1,fun()->allqdrtcsX4 (x+1u))
  21.   let rec allqdrtcsX3 x = allstateCIS(qdrtc (((x*(x+1u))<<<1)-1u) (1 - int x),fun() -> allqdrtcsX3 (x+1u))
  22.   let rec joinT3 (ass:allstateCIS<'b>) = stateCIS<'b>(ass.ss.v,ass.ss.a,fun()->
  23.     let rec (^) (xs:stateCIS<'b>) (ys:stateCIS<'b>) = //union op for CoInductiveStreams
  24.       count := !count + 1u
  25.       match compare xs.v ys.v with
  26.         | 1 -> stateCIS(ys.v,ys.a,fun() -> xs ^ ys.cont())
  27.         | _ -> stateCIS(xs.v,xs.a,fun() -> xs.cont() ^ ys) //<= then keep all the values without combining
  28.     let rec pairs (ass:allstateCIS<'b>) =
  29.       let ys = ass.cont in
  30.       allstateCIS(stateCIS(ass.ss.v,ass.ss.a,fun()->ass.ss.cont()^ys().ss),fun()->pairs (ys().cont()))
  31.     let ys = ass.cont() in let zs = ys.cont() in (ass.ss.cont()^(ys.ss^zs.ss))^joinT3 (pairs (zs.cont())))
  32.   let rec mkprm (cs:cndstate) (sqrs:stateCIS<_>) (qX4:stateCIS<_>) (qX3:stateCIS<_>) tgl =
  33.     let inline advcnd (cs:cndstate) =
  34.       let inline whladv i = if i < 47uy then i + 1uy else 0uy
  35.       let inline modadv m a = let md = m + a in if md >= 12uy then md - 12uy else md
  36.       let a = WHLPTRN.[int cs.wi] in let nc = cs.c+uint32 a
  37.       if nc<cs.c then failwith "Tried to enumerate primes past the numeric range!!!"
  38.       else cndstate(nc,whladv cs.wi,modadv cs.md12 a)
  39.     if cs.c>=sqrs.v then mkprm (if cs.c=sqrs.v then advcnd cs else cs) (sqrs.cont()) qX4 qX3 false //squarefree function
  40.     elif cs.c>qX4.v then mkprm cs sqrs (qX4.cont()) qX3 false
  41.     elif cs.c>qX3.v then mkprm cs sqrs qX4 (qX3.cont()) false
  42.     else match cs.md12 with
  43.            | 7uy -> if cs.c=qX3.v then mkprm cs sqrs qX4 (qX3.cont()) (if qX3.a>0 then not tgl else tgl) //only for a's are positive
  44.                     elif tgl then prmsCIS(cs.c,fun() -> mkprm (advcnd cs) sqrs qX4 qX3 false)
  45.                     else mkprm (advcnd cs) sqrs qX4 qX3 false
  46.            | 11uy -> if cs.c=qX3.v then mkprm cs sqrs qX4 (qX3.cont()) (if qX3.a<0 then not tgl else tgl) //only for a's are negatve
  47.                      elif tgl then prmsCIS(cs.c,fun() -> mkprm (advcnd cs) sqrs qX4 qX3 false)
  48.                      else mkprm (advcnd cs) sqrs qX4 qX3 false
  49.            | _ -> if cs.c=qX4.v then mkprm cs sqrs (qX4.cont()) qX3 (not tgl) //always must be 1uy or 5uy
  50.                   elif tgl then prmsCIS(cs.c,fun() -> mkprm (advcnd cs) sqrs qX4 qX3 false)
  51.                   else mkprm (advcnd cs) sqrs qX4 qX3 false
  52.   let qX4s = joinT3 (allqdrtcsX4 1u) in let qX3s = joinT3 (allqdrtcsX3 1u)
  53.   let rec baseprimes = prmsCIS(11u,fun() -> mkprm (cndstate(13u,1uy,1uy)) initsqrs qX4s qX3s false)
  54.   and initsqrs = joinT3 (allsqrs baseprimes)
  55.   let genseq ps = Seq.unfold (fun (psd:prmsCIS) -> Some(psd.p,psd.cont())) ps //change output n to !count to see # reductions
  56.   seq { count := 0u; yield 2u; yield 3u; yield 5u; yield 7u;
  57.         yield! mkprm (cndstate(11u,0uy,11uy)) initsqrs qX4s qX3s false |> genseq }
  58.  
  59. primesTFWSA() |> Seq.nth 100000 |> printfn "%A"
Success #stdin #stdout 2.44s 13480KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1299721u
https://ideone.com/tldRPw
This code is an implementation of the Sieve of Atkin to an incremental purely functional F# version.
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!