Хотя один ответ давал алгоритм, использующий Priority Queue (PQ) , как в SkewBinomialHeap , это, возможно, не правильный PQ дляработа.Для инкрементного сита эратосфена (iEoS) требуется PQ, обладающий превосходной производительностью для получения минимального значения и повторной вставки значений, в основном, чуть дальше по очереди, но не требует максимальной производительности для добавления новых значений, поскольку iSoE добавляет только как новыеоценивает общее количество простых чисел вплоть до квадратного корня из диапазона (что составляет крошечную долю от числа повторных вставок, которые происходят один раз за сокращение).SkewBinomialHeap PQ на самом деле не дает намного больше, чем использование встроенной карты, которая использует сбалансированное двоичное дерево поиска - все операции O (log n) - кроме того, что она слегка меняет вес операций вв пользу требований SoE.Однако SkewBinaryHeap по-прежнему требует много операций O (log n) для каждого сокращения.
PQ, реализованный как Heap , более конкретно, как Binary Heap и даже большев частности, MinHeap в значительной степени удовлетворяет требованиям iSoE с производительностью O (1) для получения минимальной производительности и производительности O (log n) для повторных вставок и добавления новых записей, хотя производительность на самом деле составляет долю O (log n), так как большинствоповторные вставки происходят в верхней части очереди, и большинство добавлений новых значений (которые не имеют значения, поскольку они встречаются редко) происходят в конце очереди, где эти операции наиболее эффективны.Кроме того, MinHeap PQ может эффективно реализовать минимум удаления и функцию вставки за один (как правило, часть) один O (log n) проход.Затем вместо карты (которая реализована в виде AVL-дерева ), где есть одна операция O (log n) с полным диапазоном 'log n' из-за минимального значения, которое нам требуется прина крайнем левом последнем листе дерева мы обычно добавляем и удаляем минимум в корне и вставляем в среднем несколько уровней вниз за один проход.Таким образом, MinHeap PQ может использоваться только с одной долей операции O (log n) на сокращение выборки, а не с несколькими операциями с большей долей O (log n).
MinHeap PQ может быть реализован с помощью чистого функционального кода (без «removeMin», реализованного, поскольку iSoE не требует этого, но есть функция «Adjust» для использования в сегментации), следующим образом:
[<RequireQualifiedAccess>]
module MinHeap =
type MinHeapTreeEntry<'T> = class val k:uint32 val v:'T new(k,v) = { k=k;v=v } end
[<CompilationRepresentation(CompilationRepresentationFlags.UseNullAsTrueValue)>]
[<NoEquality; NoComparison>]
type MinHeapTree<'T> =
| HeapEmpty
| HeapOne of MinHeapTreeEntry<'T>
| HeapNode of MinHeapTreeEntry<'T> * MinHeapTree<'T> * MinHeapTree<'T> * uint32
let empty = HeapEmpty
let getMin pq = match pq with | HeapOne(kv) | HeapNode(kv,_,_,_) -> Some kv | _ -> None
let insert k v pq =
let kv = MinHeapTreeEntry(k,v)
let rec insert' kv msk pq =
match pq with
| HeapEmpty -> HeapOne kv
| HeapOne kv2 -> if k < kv2.k then HeapNode(kv,pq,HeapEmpty,2u)
else let nn = HeapOne kv in HeapNode(kv2,nn,HeapEmpty,2u)
| HeapNode(kv2,l,r,cnt) ->
let nc = cnt + 1u
let nmsk = if msk <> 0u then msk <<< 1
else let s = int32 (System.Math.Log (float nc) / System.Math.Log(2.0))
(nc <<< (32 - s)) ||| 1u //never ever zero again with the or'ed 1
if k <= kv2.k then if (nmsk &&& 0x80000000u) = 0u then HeapNode(kv,insert' kv2 nmsk l,r,nc)
else HeapNode(kv,l,insert' kv2 nmsk r,nc)
else if (nmsk &&& 0x80000000u) = 0u then HeapNode(kv2,insert' kv nmsk l,r,nc)
else HeapNode(kv2,l,insert' kv nmsk r,nc)
insert' kv 0u pq
let private reheapify kv k pq =
let rec reheapify' pq =
match pq with
| HeapEmpty -> HeapEmpty //should never be taken
| HeapOne kvn -> HeapOne kv
| HeapNode(kvn,l,r,cnt) ->
match r with
| HeapOne kvr when k > kvr.k ->
match l with //never HeapEmpty
| HeapOne kvl when k > kvl.k -> //both qualify, choose least
if kvl.k > kvr.k then HeapNode(kvr,l,HeapOne kv,cnt)
else HeapNode(kvl,HeapOne kv,r,cnt)
| HeapNode(kvl,_,_,_) when k > kvl.k -> //both qualify, choose least
if kvl.k > kvr.k then HeapNode(kvr,l,HeapOne kv,cnt)
else HeapNode(kvl,reheapify' l,r,cnt)
| _ -> HeapNode(kvr,l,HeapOne kv,cnt) //only right qualifies
| HeapNode(kvr,_,_,_) when k > kvr.k -> //need adjusting for left leaf or else left leaf
match l with //never HeapEmpty or HeapOne
| HeapNode(kvl,_,_,_) when k > kvl.k -> //both qualify, choose least
if kvl.k > kvr.k then HeapNode(kvr,l,reheapify' r,cnt)
else HeapNode(kvl,reheapify' l,r,cnt)
| _ -> HeapNode(kvr,l,reheapify' r,cnt) //only right qualifies
| _ -> match l with //r could be HeapEmpty but l never HeapEmpty
| HeapOne(kvl) when k > kvl.k -> HeapNode(kvl,HeapOne kv,r,cnt)
| HeapNode(kvl,_,_,_) when k > kvl.k -> HeapNode(kvl,reheapify' l,r,cnt)
| _ -> HeapNode(kv,l,r,cnt) //just replace the contents of pq node with sub leaves the same
reheapify' pq
let reinsertMinAs k v pq =
let kv = MinHeapTreeEntry(k,v)
reheapify kv k pq
let adjust f (pq:MinHeapTree<_>) = //adjust all the contents using the function, then rebuild by reheapify
let rec adjust' pq =
match pq with
| HeapEmpty -> pq
| HeapOne kv -> HeapOne(MinHeapTreeEntry(f kv.k kv.v))
| HeapNode (kv,l,r,cnt) -> let nkv = MinHeapTreeEntry(f kv.k kv.v)
reheapify nkv nkv.k (HeapNode(kv,adjust' l,adjust' r,cnt))
adjust' pq
Используя вышеупомянутый модуль, iSoE можно записатьс оптимизацией факторизации колеса и использованием эффективных коиндуктивных потоков (CIS) следующим образом:
type CIS<'T> = class val v:'T val cont:unit->CIS<'T> new(v,cont) = { v=v;cont=cont } end
type cullstate = struct val p:uint32 val wi:int new(cnd,cndwi) = { p=cnd;wi=cndwi } end
let primesPQWSE() =
let WHLPRMS = [| 2u;3u;5u;7u |] in let FSTPRM = 11u in let WHLCRC = int (WHLPRMS |> Seq.fold (*) 1u) >>> 1
let WHLLMT = int (WHLPRMS |> Seq.fold (fun o n->o*(n-1u)) 1u) - 1
let WHLPTRN =
let wp = Array.zeroCreate (WHLLMT+1)
let gaps (a:int[]) = let rec gap i acc = if a.[i]=0 then gap (i+1) (acc+1uy) else acc
{0..WHLCRC-1} |> Seq.fold (fun s i->
let ns = if a.[i]<>0 then wp.[s]<-2uy*gap (i+1) 1uy;(s+1) else s in ns) 0 |> ignore
Array.init (WHLCRC+1) (fun i->if WHLPRMS |> Seq.forall (fun p->(FSTPRM+uint32(i<<<1))%p<>0u)
then 1 else 0) |> gaps;wp
let inline whladv i = if i < WHLLMT then i + 1 else 0 in let advcnd c i = c + uint32 WHLPTRN.[i]
let inline culladv c p i = let n = c + uint32 WHLPTRN.[i] * p in if n < c then 0xFFFFFFFFu else n
let rec mkprm (n,wi,pq,(bps:CIS<_>),q) =
let nxt = advcnd n wi in let nxti = whladv wi
if nxt < n then (0u,0,(0xFFFFFFFFu,0,MinHeap.empty,bps,q))
elif n>=q then let bp,bpi = bps.v in let nc,nci = culladv n bp bpi,whladv bpi
let nsd = bps.cont() in let np,_ = nsd.v in let sqr = if np>65535u then 0xFFFFFFFFu else np*np
mkprm (nxt,nxti,(MinHeap.insert nc (cullstate(bp,nci)) pq),nsd,sqr)
else match MinHeap.getMin pq with | None -> (n,wi,(nxt,nxti,pq,bps,q))
| Some kv -> let ca,cs = culladv kv.k kv.v.p kv.v.wi,cullstate(kv.v.p,whladv kv.v.wi)
if n>kv.k then mkprm (n,wi,(MinHeap.reinsertMinAs ca cs pq),bps,q)
elif n=kv.k then mkprm (nxt,nxti,(MinHeap.reinsertMinAs ca cs pq),bps,q)
else (n,wi,(nxt,nxti,pq,bps,q))
let rec pCID p pi pq bps q = CIS((p,pi),fun()->let (np,npi,(nxt,nxti,npq,nbps,nq))=mkprm (advcnd p pi,whladv pi,pq,bps,q)
pCID np npi npq nbps nq)
let rec baseprimes() = CIS((FSTPRM,0),fun()->let np=FSTPRM+uint32 WHLPTRN.[0]
pCID np (whladv 0) MinHeap.empty (baseprimes()) (FSTPRM*FSTPRM))
let genseq sd = Seq.unfold (fun (p,pi,pcc) ->if p=0u then None else Some(p,mkprm pcc)) sd
seq { yield! WHLPRMS; yield! mkprm (FSTPRM,0,MinHeap.empty,baseprimes(),(FSTPRM*FSTPRM)) |> genseq }
Приведенный выше код вычисляет первые 100 000 простых чисел примерно за 0,077 секунды, первые 1 000 000 простых чисел за 0,977 секунды, первое10 000 000 простых чисел за 14,33 секунды и первые 100 000 000 простых чисел за 221,87 секунды, все на i7-2700K (3,5 ГГц) в виде 64-битного кода.Этот чисто функциональный код немного быстрее, чем код изменчивого кода на основе словаря Дастина Камбелла с добавленной общей оптимизацией факторизации колес, отложенным добавлением базовых простых чисел и использованием всех более эффективных CID (tryfsharp и ideone ) , но все еще является чисто функциональным кодом, в котором он использует класс Dictionary не .Однако для больших диапазонов простых чисел, составляющих около двух миллиардов (около 100 миллионов простых), код, использующий словарь на основе хеш-таблицы, будет работать быстрее, поскольку операции словаря не имеют коэффициента O (log n), и этот выигрыш преодолевает вычислительныйсложность использования словарных хеш-таблиц.
Вышеуказанная программа имеет еще одну особенность, заключающуюся в том, что колесо факторизации параметризовано так, что, например, можно использовать очень большое колесо, установив WHLPRMS в [| 2u; 3u; 5u; 7u; 11u; 13u; 17u; 19u |] и FSTPRM до 23u, чтобы получить время выполнения около двух третей для больших диапазонов при 9,34 секунды для десяти миллионов простых чисел, хотя следует учесть, что для этого требуется несколько секунд вычислите WHLPTRN до запуска программы, что является постоянным расходом независимо от простого диапазона.
Сравнительный анализ : По сравнению с чисто функциональной реализацией пошагового сворачивания дерева этот алгоритм немного быстрее, потому что средняя используемая высота дерева MinHeap меньше в два раза, чем глубина сложенное дерево, но это компенсируется эквивалентной постоянной потерей коэффициента эффективности в способности обходить уровни дерева PQ, поскольку оно основано на двоичной куче, требующей обработки как правого, так и левого листьев для каждого уровня кучи и ветви, так или иначе чем одно сравнение на уровень для свертывания дерева с обычно менее глубокой ветвью, чем взятая. По сравнению с другими функциональными алгоритмами, основанными на PQ и карте, улучшения, как правило, являются постоянным фактором сокращения числа операций O (log n) при обходе каждого уровня соответствующих древовидных структур.
MinHeap обычно реализуется как изменяемый массив двоичная куча после модели на основе генеалогического дерева, изобретенной Майклом Эйтцингером более 400 лет назад. Я знаю, что вопрос сказал, что не интересен нефункциональный изменяемый код, но если нужно избегать всех вложенных кодов, которые используют изменчивость, то мы не могли бы использовать списки или LazyList, которые используют изменчивость «под прикрытием» из соображений производительности. Итак, представьте, что следующая альтернативная изменяемая версия MinHeap PQ предоставлена библиотекой и имеет еще один коэффициент более двух для больших простых диапазонов производительности:
[<RequireQualifiedAccess>]
module MinHeap =
type MinHeapTreeEntry<'T> = class val k:uint32 val v:'T new(k,v) = { k=k;v=v } end
type MinHeapTree<'T> = ResizeArray<MinHeapTreeEntry<'T>>
let empty<'T> = MinHeapTree<MinHeapTreeEntry<'T>>()
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
Замечание Geek: на самом деле я ожидал, что изменяемая версия предложит намного лучшее улучшенное соотношение производительности, но она затормозится при повторных вставках из-за вложенной структуры кода if-then-else и случайного поведения основной выборки значения означают, что предсказание ветвления ЦП не выполняется для большой части ветвей, что приводит к множеству дополнительных 10-ти тактовых циклов ЦП на единицу сокращения для восстановления кеша предварительной выборки команд.
Единственным другим приростом производительности с постоянным коэффициентом по этому алгоритму будет сегментирование и использование многозадачности для увеличения производительности, пропорционального количеству ядер ЦП; однако, в его нынешнем виде, это самый быстрый чисто функциональный алгоритм SoE на сегодняшний день, и даже чистая функциональная форма с использованием функционала MinHeap превосходит упрощенные императивные реализации, такие как код Джона Харропа или сито Йохана Кульбома Аткин (что является ошибкой в его времени, так как он вычислил только простых чисел до 10 миллионов, а не 10 миллионовных простых ), но эти алгоритмы были бы примерно в пять раз быстрее, если бы использовались лучшие оптимизации. Это соотношение около пяти между функциональным и императивным кодом будет несколько уменьшено, когда мы добавим многопоточность с большей факторизацией колеса, так как вычислительная сложность императивного кода увеличивается быстрее, чем функциональный код, а многопоточность помогает более медленному функциональному коду больше, чем более быстрый императивный код, поскольку последний приближается к базовому пределу времени, необходимого для перечисления через найденные простые числа.
EDIT_ADD: Несмотря на то, что можно было бы продолжать использовать чисто функциональную версию MinHeap, добавление эффективной сегментации при подготовке к многопоточности немного «сломало бы» чистоту «функционального кода следующим образом: 1) Наиболее эффективный способ передачи представления составных отбракованных чисел - это массив упакованных битов размером сегмента, 2) хотя размер массива известен с использованием понимания массива функционально инициализировать его неэффективно, так как он использует «ResizeArray» под крышками, который должен копировать себя для каждого добавления x (я думаю, «x» - восемь для текущей реализации), а использование Array.init не работать так, как пропускаются многие значения по определенным индексам, 3) Следовательно, самый простой способ заполнить составной массив - это обнулить его, создать правильный размер и затем запустить функцию инициализации, которая может записывать в каждый индекс изменяемого массива не более одного раза , Хотя это не является строго «функциональным», оно близко к тому, что массив инициализируется, а затем никогда не изменяется снова.
Код с добавленной сегментацией, многопоточностью, программируемой факториальной окружностью колес и многими настройками производительности выглядит следующим образом (кроме некоторых добавленных новых констант, дополнительный настроенный код для реализации сегментации и многопоточности составляет примерно половину кода, начинающегося с функции "prmspg"):
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.wi,cullstate(kv.v.p,whladv kv.v.wi)
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.wi]
let db = if k<low then let nk = int64(low-k)+db in nk-((nk/WHLSPN)*WHLSPN)
else let nk = int64(k-low) in if db<nk then db+WHLSPN-nk else db-nk
let r = WHLRNDUP.[int((((db>>>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<db then x+WHLSPN-db else x-db in uint32 x,cullstate(p,r)
let bfbtsz = int rng/WHLCRC*(WHLLMT+1) in let nbuf = Array.zeroCreate (bfbtsz>>>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<NUMPRCS then cont.[i+1]
else prmspg (nxt+uint64 adj) adj pq bp q)
let _,tsk = ncont.[0] in let nbuf,_,_,_ = tsk.Result in nbuf,ncont
//init cond buf[0], no queue, frst bp sqr offset
let initcond = 0u,System.Threading.Tasks.Task.Factory.StartNew (fun()->
(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)) }
Обратите внимание, что в модулях MinHeap, как функциональных, так и основанных на массивах, была добавлена функция «Adjust», позволяющая регулировать состояние отбраковки версии PQ каждого потока в начале каждой новой страницы сегмента. Также обратите внимание, что было возможно настроить код так, чтобы большая часть вычислений выполнялась с использованием 32-битных диапазонов с конечной последовательностью, выводимой как uint64, с небольшими затратами времени вычислений, так что в настоящее время теоретический диапазон составляет что-то более 100 триллионов (десять повышено до четырнадцатая степень), если кто-то готов подождать около трех-четырех месяцев, необходимых для вычисления этого диапазона. Проверки числового диапазона были удалены, поскольку маловероятно, чтобы кто-либо использовал этот алгоритм для вычисления этого диапазона, не говоря уже о его пропуске.
Используя чистый функционал MinHeap и факторизацию колес 2,3,5,7, вышеуказанная программа вычисляет первые сто тысяч, миллион, десять миллионов и сто миллионов простых чисел за 0,062, 0,629, 10,53 и 195,62 секунды, соответственно. Использование MinHeap на основе массива ускоряет это до 0,097, 0,276, 3,48 и 51,60 секунд соответственно. Использование колеса 2,3,5,7,11,13,17 путем изменения WHLPRMS на «[| 2u; 3u; 5u; 7u; 11u; 13u; 17u |]» и FSTPRM до 19u ускоряет эту скорость еще больше до 0,181, 0,308, 2,49 и 36,58 секунд соответственно (для постоянного улучшения коэффициента при постоянных издержках). Эта самая быстрая настройка вычисляет 203 280 221 простых чисел в диапазоне 32-разрядных чисел за 88,37 секунды. Константа «BFSZ» может быть скорректирована с помощью компромиссов между более медленным временем для версии с меньшими диапазонами и более быстрым временем для больших диапазонов, при этом значение «1 <<< 14» рекомендуется использовать для больших диапазонов. Эта константа устанавливает только минимальный размер буфера, при этом программа автоматически настраивает размер буфера выше этого размера для больших диапазонов, чтобы буфер был достаточным для того, чтобы наибольшее базовое простое число, требуемое для диапазона страниц, всегда "ударяло" по каждой странице хотя бы один раз ; это означает, что сложность и дополнительные затраты на дополнительное «сито ведра» не требуются. Эта последняя полностью оптимизированная версия может вычислять простые числа до 10 и 100 миллиардов примерно за 256,8 и 3617,4 секунды (чуть более часа для 100 миллиардов), как было протестировано с использованием "primesPQOWSE () |> Seq.takeWhile ((> =) 100000000000UL) |> Seq.fold (fun sp -> s + 1UL) 0UL "для вывода. Вот откуда берутся примерно полдня для подсчета простых чисел до триллиона, неделя для десяти триллионов и от трех до четырех месяцев для ста триллионов.
Я не думаю, что возможно сделать функциональный или почти функциональный код, используя инкрементный алгоритм SoE, чтобы работать намного быстрее, чем этот.Как видно из кода, оптимизация базового алгоритма приращения значительно увеличила сложность кода, так что он, вероятно, немного сложнее, чем эквивалентно оптимизированный код, основанный на отбраковке прямого массива, причем этот код способен работать примерно в десять раз быстрее, чемэто и без дополнительного показателя в производительности, означающего, что этот функциональный инкрементный код имеет постоянно увеличивающиеся дополнительные процентные издержки.
Так полезно ли это, кроме как с интересной теоретической и интеллектуальной точки зрения?Вероятно, это не так.Для меньших диапазонов простых чисел, примерно до десяти миллионов, лучшие из базовых не полностью оптимизированных инкрементальных функциональных SoE, вероятно, являются адекватными и довольно простыми для записи или имеют меньше использования оперативной памяти, чем самые простые императивные SoE.Однако они намного медленнее, чем более императивный код, использующий массив, поэтому они «выдыхаются» для диапазонов выше этого.Несмотря на то, что здесь было продемонстрировано, что код может быть ускорен за счет оптимизации, он все еще в 10 раз медленнее, чем более императивная версия на основе чистого массива, но добавляет сложности, по крайней мере, такой же сложной, как этот код с эквивалентными оптимизациями.и даже этот код под F # в DotNet примерно в четыре раза медленнее, чем при использовании языка, такого как C ++, скомпилированного непосредственно в нативный код;если бы кто-то действительно хотел исследовать большие диапазоны простых чисел, он, вероятно, использовал бы один из тех других языков и методов, где primesieve может вычислить число простых чисел в диапазоне сто триллионов менее чем за четыре часа вместо примерно трехмесяцев, необходимых для этого кода. END_EDIT_ADD