[] module MinHeap = type MinHeapTreeEntry<'T> = class val k:uint32 val v:'T new(k,v) = { k=k;v=v } end type MinHeapTree<'T> = ResizeArray> let empty<'T> = MinHeapTree>() let getMin (pq:MinHeapTree<_>) = if pq.Count > 0 then Some pq.[0] else None let insert k v (pq:MinHeapTree<_>) = if pq.Count = 0 then pq.Add(MinHeapTreeEntry(0xFFFFFFFFu,v)) //add an extra entry so there's always a right max node let mutable nxtlvl = pq.Count in let mutable lvl = nxtlvl <<< 1 //1 past index of value added times 2 pq.Add(pq.[nxtlvl - 1]) //copy bottom entry then do bubble up while less than next level up while ((lvl <- lvl >>> 1); nxtlvl <- nxtlvl >>> 1; nxtlvl <> 0) do let t = pq.[nxtlvl - 1] in if t.k > k then pq.[lvl - 1] <- t else lvl <- lvl <<< 1; nxtlvl <- 0 //causes loop break pq.[lvl - 1] <- MinHeapTreeEntry(k,v); pq let reinsertMinAs k v (pq:MinHeapTree<_>) = //do minify down for value to insert let mutable nxtlvl = 1 in let mutable lvl = nxtlvl in let cnt = pq.Count while (nxtlvl <- nxtlvl <<< 1; nxtlvl < cnt) do let lk = pq.[nxtlvl - 1].k in let rk = pq.[nxtlvl].k in let oldlvl = lvl let k = if k > lk then lvl <- nxtlvl; lk else k in if k > rk then nxtlvl <- nxtlvl + 1; lvl <- nxtlvl if lvl <> oldlvl then pq.[oldlvl - 1] <- pq.[lvl - 1] else nxtlvl <- cnt //causes loop break pq.[lvl - 1] <- MinHeapTreeEntry(k,v); pq let adjust f (pq:MinHeapTree<_>) = //adjust all the contents using the function, then re-heapify if pq <> null then let cnt = pq.Count if cnt > 1 then for i = 0 to cnt - 2 do //change contents using function let e = pq.[i] in let k,v = e.k,e.v in pq.[i] <- MinHeapTreeEntry (f k v) for i = cnt/2 downto 1 do //rebuild by reheapify let kv = pq.[i - 1] in let k = kv.k let mutable nxtlvl = i in let mutable lvl = nxtlvl while (nxtlvl <- nxtlvl <<< 1; nxtlvl < cnt) do let lk = pq.[nxtlvl - 1].k in let rk = pq.[nxtlvl].k in let oldlvl = lvl let k = if k > lk then lvl <- nxtlvl; lk else k in if k > rk then nxtlvl <- nxtlvl + 1; lvl <- nxtlvl if lvl <> oldlvl then pq.[oldlvl - 1] <- pq.[lvl - 1] else nxtlvl <- cnt //causes loop break pq.[lvl - 1] <- kv pq type prmsCIS = class val pg:uint16 val bg:uint16 val pi:int val cont:unit->prmsCIS new(pg,bg,pi,nxtprmf) = { pg=pg;bg=bg;pi=pi;cont=nxtprmf } end type cullstate = struct val p:uint32 val wi:int new(cnd,cndwi) = { p=cnd;wi=cndwi } end let primesPQOWSE() = let WHLPRMS = [| 2u;3u;5u;7u;11u;13u;17u |] in let FSTPRM = 19u in let WHLCRC = int(WHLPRMS |> Seq.fold (*) 1u) let MXSTP = uint64(FSTPRM-1u) in let BFSZ = 1<<<11 in let NUMPRCS = System.Environment.ProcessorCount let WHLLMT = int (WHLPRMS |> Seq.fold (fun o n->o*(n-1u)) 1u) - 1 in let WHLPTRN = Array.zeroCreate (WHLLMT+1) let WHLRNDUP = let gaps (a:int[]) = let rec gap i acc = if a.[i]=0 then gap (i+1) (acc+1) else acc in let b = a |> Array.scan (+) 0 Array.init (WHLCRC>>>1) (fun i-> if a.[i]=0 then 0 else let g=2*gap (i+1) 1 in WHLPTRN.[b.[i]]<-byte g;1) Array.init WHLCRC (fun i->if WHLPRMS |> Seq.forall (fun p->(FSTPRM+uint32(i<<<1))%p<>0u) then 1 else 0) |> gaps |> Array.scan (+) 0 let WHLPOS = WHLPTRN |> Array.map (uint32) |> Array.scan (+) 0u in let advcnd cnd cndi = cnd + uint32 WHLPTRN.[cndi] let MINRNGSTP = if WHLLMT<=31 then uint32(32/(WHLLMT+1)*WHLCRC) else if WHLLMT=47 then uint32 WHLCRC<<<1 else uint32 WHLCRC let MINBFRNG = uint32((BFSZ<<<3)/(WHLLMT+1)*WHLCRC)/MINRNGSTP*MINRNGSTP let MINBFRNG = if MINBFRNG=0u then MINRNGSTP else MINBFRNG let inline whladv i = if i < WHLLMT then i+1 else 0 in let inline culladv c p i = c+uint32 WHLPTRN.[i]*p let rec mkprm (n,wi,pq,(bps:prmsCIS),q,lstp,bgap) = let nxt,nxti = advcnd n wi,whladv wi if n>=q then let p = (uint32 bps.bg<<<16)+uint32 bps.pg let nbps,nxtcmpst,npi = bps.cont(),culladv n p bps.pi,whladv bps.pi let pg = uint32 nbps.pg in let np = p+pg in let sqr = q+pg*((p<<<1)+pg) //only works to p < about 13 million let nbps = prmsCIS(uint16 np,uint16(np>>>16),nbps.pi,nbps.cont) //therefore, algorithm only works to p^2 or about mkprm (nxt,nxti,(MinHeap.insert nxtcmpst (cullstate(p,npi)) pq),nbps,sqr,lstp,(bgap+1us)) //1.7 * 10^14 else match MinHeap.getMin pq with | None -> (uint16(n-uint32 lstp),bgap,wi,(nxt,nxti,pq,bps,q,n,1us)) //fix with q is uint64 | Some kv -> let ca,cs = culladv kv.k kv.v.p kv.v.pi,cullstate(kv.v.p,whladv kv.v.pi) if n>kv.k then mkprm (n,wi,(MinHeap.reinsertMinAs ca cs pq),bps,q,lstp,bgap) elif n=kv.k then mkprm (nxt,nxti,(MinHeap.reinsertMinAs ca cs pq),bps,q,lstp,(bgap+1us)) else (uint16(n-uint32 lstp),bgap,wi,(nxt,nxti,pq,bps,q,n,1us)) let rec pCIS p pg bg pi pq bps q = prmsCIS(pg,bg,pi,fun()-> let (npg,nbg,npi,(nxt,nxti,npq,nbps,nq,nl,ng))=mkprm (p+uint32 WHLPTRN.[pi],whladv pi,pq,bps,q,p,0us) pCIS (p+uint32 npg) npg nbg npi npq nbps nq) let rec baseprimes() = prmsCIS(uint16 FSTPRM,0us,0,fun()-> let np,npi=advcnd FSTPRM 0,whladv 0 pCIS np (uint16 WHLPTRN.[0]) 1us npi MinHeap.empty (baseprimes()) (FSTPRM*FSTPRM)) let prmspg nxt adj pq bp q = //compute next buffer size rounded up to next even wheel circle so at least one each base prime hits the page let rng = max (((uint32(MXSTP+uint64(sqrt (float (MXSTP*(MXSTP+4UL*nxt))))+1UL)>>>1)+MINRNGSTP)/MINRNGSTP*MINRNGSTP) MINBFRNG let nxtp() = async { let rec addprms pqx (bpx:prmsCIS) qx = if qx>=adj then pqx,bpx,qx //add primes to queue for new lower limit else let p = (uint32 bpx.bg<<<16)+uint32 bpx.pg in let nbps = bpx.cont() let pg = uint32 nbps.pg in let np = p+pg in let sqr = qx+pg*((p<<<1)+pg) let nbps = prmsCIS(uint16 np,uint16(np>>>16),nbps.pi,nbps.cont) addprms (MinHeap.insert qx (cullstate(p,bpx.pi)) pqx) nbps sqr let adjcinpg low k (v:cullstate) = //adjust the cull states for the new page low value let p = v.p in let WHLSPN = int64 WHLCRC*int64 p in let db = int64 p*int64 WHLPOS.[v.pi] let db = if k>>1)%(WHLSPN>>>1))+int64 p-1L)/int64 p)] in let x = int64 WHLPOS.[r]*int64 p let r = if r>WHLLMT then 0 else r in let x = if x>>5) let rec nxtp' wi cnt = let _,nbg,_,ncnt = mkprm cnt in let nwi = wi + int nbg if nwi < bfbtsz then nbuf.[nwi>>>5] <- nbuf.[nwi>>>5] ||| (1u<<<(nwi&&&0x1F)); nxtp' nwi ncnt else let _,_,pq,bp,q,_,_ = ncnt in nbuf,pq,bp,q //results incl buf and cont parms for next page let npq,nbp,nq = addprms pq bp q return nxtp' 0 (0u,0,MinHeap.adjust (adjcinpg adj) npq,nbp,nq-adj,0u,0us) } rng,nxtp() |> Async.StartAsTask let nxtpg nxt (cont:(_*System.Threading.Tasks.Task<_>)[]) = //(len,pq,bp,q) = let adj = (cont |> Seq.fold (fun s (r,_) -> s+r) 0u) let _,tsk = cont.[0] in let _,pq,bp,q = tsk.Result let ncont = Array.init (NUMPRCS+1) (fun i -> if i (Array.empty,MinHeap.empty,baseprimes(),FSTPRM*FSTPRM-FSTPRM)) let nxtcond n = prmspg (uint64 n) (n-FSTPRM) MinHeap.empty (baseprimes()) (FSTPRM*FSTPRM-FSTPRM) let initcont = Seq.unfold (fun (n,((r,_)as v))->Some(v,(n+r,nxtcond (n+r)))) (FSTPRM,initcond) |> Seq.take (NUMPRCS+1) |> Seq.toArray let rec nxtprm (c,ci,i,buf:uint32[],cont) = let rec nxtprm' c ci i = let nc = c + uint64 WHLPTRN.[ci] in let nci = whladv ci in let ni = i + 1 in let nw = ni>>>5 if nw >= buf.Length then let (npg,ncont)=nxtpg nc cont in nxtprm (c,ci,-1,npg,ncont) elif (buf.[nw] &&& (1u <<< (ni &&& 0x1F))) = 0u then nxtprm' nc nci ni else nc,nci,ni,buf,cont nxtprm' c ci i seq { yield! WHLPRMS |> Seq.map (uint64); yield! Seq.unfold (fun ((c,_,_,_,_) as cont)->Some(c,nxtprm cont)) (nxtprm (uint64 FSTPRM-uint64 WHLPTRN.[WHLLMT],WHLLMT,-1,Array.empty,initcont)) } primesPQOWSE() |> Dump |> Ignore