Le tamis d'Eratosthène en F #
https://stackoverflow.com/questions/4629734
-
08-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Je suis intéressé par une mise en œuvre du crible d'Eratosthène purement fonctionnelle F #. Je suis intéressé par une mise en oeuvre du tamis réel, pas la mise en œuvre fonctionnelle naïve ce n'est pas vraiment le tamis , donc pas quelque chose comme ceci:
let rec PseudoSieve list =
match list with
| hd::tl -> hd :: (PseudoSieve <| List.filter (fun x -> x % hd <> 0) tl)
| [] -> []
Le deuxième lien ci-dessus décrit brièvement un algorithme qui nécessiterait l'utilisation d'un multimap, ce qui est disponible en F # pour autant que je sache. La mise en œuvre Haskell donné utilise une carte qui supporte un mode de insertWith, que je ne l'ai pas vu disponible dans la section F # carte fonctionnelle .
Quelqu'un sait-il un moyen de traduire le code de la carte Haskell donnée à F #, ou sait peut-être des méthodes de mise en œuvre alternatives ou des algorithmes tamiser qui sont aussi efficaces et mieux adaptés à une mise en œuvre fonctionnelle ou F #?
La solution
La lecture de cet article, je suis venu avec une idée qui ne nécessite pas multimap. Il gère les clés de carte entrant en collision en déplaçant la clé entrant en collision avant par sa valeur prime encore et encore jusqu'à ce qu'il atteigne une clé qui ne sont pas sur la carte. Ci-dessous primes est une carte avec les clés de la prochaine valeur iterator et les valeurs qui sont des nombres premiers.
let primes =
let rec nextPrime n p primes =
if primes |> Map.containsKey n then
nextPrime (n + p) p primes
else
primes.Add(n, p)
let rec prime n primes =
seq {
if primes |> Map.containsKey n then
let p = primes.Item n
yield! prime (n + 1) (nextPrime (n + p) p (primes.Remove n))
else
yield n
yield! prime (n + 1) (primes.Add(n * n, n))
}
prime 2 Map.empty
Voici l'algorithme de file d'attente prioritaire de cette papier sans l'optimisation carré. J'ai placé les fonctions de file d'attente prioritaire génériques en haut. J'ai utilisé un tuple pour représenter les itérateurs liste paresseux.
let primes() =
// the priority queue functions
let insert = Heap.Insert
let findMin = Heap.Min
let insertDeleteMin = Heap.DeleteInsert
// skips primes 2, 3, 5, 7
let wheelData = [|2L;4L;2L;4L;6L;2L;6L;4L;2L;4L;6L;6L;2L;6L;4L;2L;6L;4L;6L;8L;4L;2L;4L;2L;4L;8L;6L;4L;6L;2L;4L;6L;2L;6L;6L;4L;2L;4L;6L;2L;6L;4L;2L;4L;2L;10L;2L;10L|]
// increments iterator
let wheel (composite, n, prime) =
composite + wheelData.[n % 48] * prime, n + 1, prime
let insertPrime prime n table =
insert (prime * prime, n, prime) table
let rec adjust x (table : Heap) =
let composite, n, prime = findMin table
if composite <= x then
table
|> insertDeleteMin (wheel (composite, n, prime))
|> adjust x
else
table
let rec sieve iterator table =
seq {
let x, n, _ = iterator
let composite, _, _ = findMin table
if composite <= x then
yield! sieve (wheel iterator) (adjust x table)
else
if x = 13L then
yield! [2L; 3L; 5L; 7L; 11L]
yield x
yield! sieve (wheel iterator) (insertPrime x n table)
}
sieve (13L, 1, 1L) (insertPrime 11L 0 (Heap(0L, 0, 0L)))
Voici l'algorithme de file d'attente de priorité à l'optimisation carré. Afin de faciliter l'ajout paresseux nombres premiers à la table de consultation, les décalages de roues devaient être retournés ainsi que les valeurs de choix. Cette version de l'algorithme a O (sqrt (n)) utilisation de la mémoire lorsque la personne est optimisé une O (n).
let rec primes2() : seq<int64 * int> =
// the priority queue functions
let insert = Heap.Insert
let findMin = Heap.Min
let insertDeleteMin = Heap.DeleteInsert
// increments iterator
let wheel (composite, n, prime) =
composite + wheelData.[n % 48] * prime, n + 1, prime
let insertPrime enumerator composite table =
// lazy initialize the enumerator
let enumerator =
if enumerator = null then
let enumerator = primes2().GetEnumerator()
enumerator.MoveNext() |> ignore
// skip primes that are a part of the wheel
while fst enumerator.Current < 11L do
enumerator.MoveNext() |> ignore
enumerator
else
enumerator
let prime = fst enumerator.Current
// Wait to insert primes until their square is less than the tables current min
if prime * prime < composite then
enumerator.MoveNext() |> ignore
let prime, n = enumerator.Current
enumerator, insert (prime * prime, n, prime) table
else
enumerator, table
let rec adjust x table =
let composite, n, prime = findMin table
if composite <= x then
table
|> insertDeleteMin (wheel (composite, n, prime))
|> adjust x
else
table
let rec sieve iterator (enumerator, table) =
seq {
let x, n, _ = iterator
let composite, _, _ = findMin table
if composite <= x then
yield! sieve (wheel iterator) (enumerator, adjust x table)
else
if x = 13L then
yield! [2L, 0; 3L, 0; 5L, 0; 7L, 0; 11L, 0]
yield x, n
yield! sieve (wheel iterator) (insertPrime enumerator composite table)
}
sieve (13L, 1, 1L) (null, insert (11L * 11L, 0, 11L) (Heap(0L, 0, 0L)))
Voici mon programme de test.
type GenericHeap<'T when 'T : comparison>(defaultValue : 'T) =
let mutable capacity = 1
let mutable values = Array.create capacity defaultValue
let mutable size = 0
let swap i n =
let temp = values.[i]
values.[i] <- values.[n]
values.[n] <- temp
let rec rollUp i =
if i > 0 then
let parent = (i - 1) / 2
if values.[i] < values.[parent] then
swap i parent
rollUp parent
let rec rollDown i =
let left, right = 2 * i + 1, 2 * i + 2
if right < size then
if values.[left] < values.[i] then
if values.[left] < values.[right] then
swap left i
rollDown left
else
swap right i
rollDown right
elif values.[right] < values.[i] then
swap right i
rollDown right
elif left < size then
if values.[left] < values.[i] then
swap left i
member this.insert (value : 'T) =
if size = capacity then
capacity <- capacity * 2
let newValues = Array.zeroCreate capacity
for i in 0 .. size - 1 do
newValues.[i] <- values.[i]
values <- newValues
values.[size] <- value
size <- size + 1
rollUp (size - 1)
member this.delete () =
values.[0] <- values.[size]
size <- size - 1
rollDown 0
member this.deleteInsert (value : 'T) =
values.[0] <- value
rollDown 0
member this.min () =
values.[0]
static member Insert (value : 'T) (heap : GenericHeap<'T>) =
heap.insert value
heap
static member DeleteInsert (value : 'T) (heap : GenericHeap<'T>) =
heap.deleteInsert value
heap
static member Min (heap : GenericHeap<'T>) =
heap.min()
type Heap = GenericHeap<int64 * int * int64>
let wheelData = [|2L;4L;2L;4L;6L;2L;6L;4L;2L;4L;6L;6L;2L;6L;4L;2L;6L;4L;6L;8L;4L;2L;4L;2L;4L;8L;6L;4L;6L;2L;4L;6L;2L;6L;6L;4L;2L;4L;6L;2L;6L;4L;2L;4L;2L;10L;2L;10L|]
let primes() =
// the priority queue functions
let insert = Heap.Insert
let findMin = Heap.Min
let insertDeleteMin = Heap.DeleteInsert
// increments iterator
let wheel (composite, n, prime) =
composite + wheelData.[n % 48] * prime, n + 1, prime
let insertPrime prime n table =
insert (prime * prime, n, prime) table
let rec adjust x (table : Heap) =
let composite, n, prime = findMin table
if composite <= x then
table
|> insertDeleteMin (wheel (composite, n, prime))
|> adjust x
else
table
let rec sieve iterator table =
seq {
let x, n, _ = iterator
let composite, _, _ = findMin table
if composite <= x then
yield! sieve (wheel iterator) (adjust x table)
else
if x = 13L then
yield! [2L; 3L; 5L; 7L; 11L]
yield x
yield! sieve (wheel iterator) (insertPrime x n table)
}
sieve (13L, 1, 1L) (insertPrime 11L 0 (Heap(0L, 0, 0L)))
let rec primes2() : seq<int64 * int> =
// the priority queue functions
let insert = Heap.Insert
let findMin = Heap.Min
let insertDeleteMin = Heap.DeleteInsert
// increments iterator
let wheel (composite, n, prime) =
composite + wheelData.[n % 48] * prime, n + 1, prime
let insertPrime enumerator composite table =
// lazy initialize the enumerator
let enumerator =
if enumerator = null then
let enumerator = primes2().GetEnumerator()
enumerator.MoveNext() |> ignore
// skip primes that are a part of the wheel
while fst enumerator.Current < 11L do
enumerator.MoveNext() |> ignore
enumerator
else
enumerator
let prime = fst enumerator.Current
// Wait to insert primes until their square is less than the tables current min
if prime * prime < composite then
enumerator.MoveNext() |> ignore
let prime, n = enumerator.Current
enumerator, insert (prime * prime, n, prime) table
else
enumerator, table
let rec adjust x table =
let composite, n, prime = findMin table
if composite <= x then
table
|> insertDeleteMin (wheel (composite, n, prime))
|> adjust x
else
table
let rec sieve iterator (enumerator, table) =
seq {
let x, n, _ = iterator
let composite, _, _ = findMin table
if composite <= x then
yield! sieve (wheel iterator) (enumerator, adjust x table)
else
if x = 13L then
yield! [2L, 0; 3L, 0; 5L, 0; 7L, 0; 11L, 0]
yield x, n
yield! sieve (wheel iterator) (insertPrime enumerator composite table)
}
sieve (13L, 1, 1L) (null, insert (11L * 11L, 0, 11L) (Heap(0L, 0, 0L)))
let mutable i = 0
let compare a b =
i <- i + 1
if a = b then
true
else
printfn "%A %A %A" a b i
false
Seq.forall2 compare (Seq.take 50000 (primes())) (Seq.take 50000 (primes2() |> Seq.map fst))
|> printfn "%A"
primes2()
|> Seq.map fst
|> Seq.take 10
|> Seq.toArray
|> printfn "%A"
primes2()
|> Seq.map fst
|> Seq.skip 999999
|> Seq.take 10
|> Seq.toArray
|> printfn "%A"
System.Console.ReadLine() |> ignore
Autres conseils
Bien qu'il y ait une réponse donnant un algorithme en utilisant un priorité Queue (PQ) comme dans SkewBinomialHeap , il est peut-être pas le droit PQ pour le travail. Ce que le Crible d'Eratosthène incrémentale (IEOS) exige un PQ qui a d'excellentes performances pour obtenir la valeur minimale et les valeurs réinsérant la plupart du temps un peu plus loin dans la file d'attente, mais n'a pas besoin le nec plus ultra de la performance pour l'ajout de nouvelles valeurs Isoe ajoute que neuf un total des valeurs des nombres premiers jusqu'à la racine carrée de la distance (qui est une petite fraction du nombre de ré-insertion qui se produisent une fois par réduction). Le SkewBinomialHeap PQ ne donne pas vraiment beaucoup plus que d'utiliser la carte intégrée qui utilise un fichier binaire équilibré arbre de recherche - toutes les opérations O (log n) - autres que la pondération change elle des opérations légèrement en faveur des besoins de l'entreprise publique. Cependant, le SkewBinaryHeap nécessite encore beaucoup O (log n) par réduction.
A PQ implémenté comme un Heap plus en particulier en tant que Binary Heap et plus particulièrement en tant que MinHeap à peu près satisfait aux exigences de Isoe avec la performance O (1) à obtenir la minimum et O (log n) performance pour les insertions de re-et ajouter de nouvelles entrées, bien que le rendement est en fait une fraction de O (log n) comme la plupart des insertions de se reproduire dans la partie supérieure de la file d'attente et la plupart des additions de de nouvelles valeurs (qui ne comptent pas car ils sont rares) se produisent à la fin de la file d'attente où ces opérations sont plus efficaces. De plus, le PQ MinHeap peut mettre en œuvre efficacement le minimum supprimer et insérer la fonction dans un (généralement une fraction de) un passage O (log n). Alors, plutôt que pour la carte (qui est mis en œuvre comme AVL arbre ) où il y a une O (log n) opération avec généralement une gamme « log n » complète en raison de la valeur minimum que nous exigeons d'être à la dernière feuille gauche de l'arbre, nous ajoutons en général, et enlever le minimum à la racine et l'insertion sur la moyenne des quelques niveaux bas en un seul passage. Ainsi, le PQ MinHeap peut être utilisé avec seulement une fraction de fonctionnement O (log n) par abattage réduction plutôt que des opérations multiples plus grande fraction de O (log n).
Le PQ MinHeap peut être mis en œuvre avec le code fonctionnel pur (sans « removeMin » mis en œuvre comme Isoe ne l'exige pas, mais il y a une fonction « ajuster » pour une utilisation dans la segmentation), comme suit:
[<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
En utilisant le module ci-dessus, le Isoe peut être écrit avec les optimisations de factorisation des roues et en utilisant efficace Streams Co-inductifs (de la CEI) comme suit:
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 }
Le code ci-dessus calcule les 100.000 premiers nombres premiers dans environ 0,077 secondes, les premiers nombres premiers 1.000.000 en 0,977 secondes, les premiers nombres premiers en 10.000.000 environ 14,33 secondes, et les premiers nombres premiers en 100.000.000 environ 221,87 secondes, le tout sur un i7-2700K ( 3.5GHz) en tant que code 64 bits. Ce code purement fonctionnel est légèrement plus rapide que celle de Dustin Cambell de mutable code basé Dictionnaire avec les optimisations communes ajoutée de factorisation roue, différés ajoutant des nombres premiers de base, et l'utilisation du CID plus efficace est tout ajouté ( et href="http://ideone.com/vjviLS" rel="nofollow noreferrer"> ideone ) mais il est toujours pur code fonctionnel où son utilisation de la classe du dictionnaire est pas . Cependant, pour les gammes plus grandes prime d'environ d'environ deux milliards (environ 100 millions de nombres premiers), le code en utilisant la table de hachage d'un dictionnaire sera plus rapide que le dictionnaire les opérations ne sont pas un facteur O (log n) et ce gain permet de surmonter la complexité de calcul de l'utilisation de tables de hachage Dictionnaire.
Le programme ci-dessus a la particularité en outre que la roue de factorisation est paramétré de telle sorte que, par exemple, on peut utiliser une roue extrêmement importante en mettant WHLPRMS à [| 2u, 3U, 5U, 7U, 11u, 13u, 17u, 19u |] et FSTPRM à 23U pour obtenir un temps de fonctionnement d'environ deux tiers pour les grandes plages à environ 9,34 secondes pour dix millions de nombres premiers, bien qu'il faille noter qu'il faut plusieurs secondes pour calculer le WHLPTRN avant que le programme commence à courir, ce qui est une constante d'importance en tête pas la gamme de choix.
Analyse comparative : Par rapport à cet algorithme est juste un peu plus rapide pur la mise en œuvre de pliage d'arbre incrémental fonctionnelle, car la hauteur utilisée moyenne de l'arbre MinHeap est moins par un facteur de deux à la profondeur de l'arbre plié mais qui est décalé par une perte équivalente de facteur de constante de l'efficacité de la capacité à traverser les niveaux de l'arbre de PQ parce qu'elle est basée sur un tas binaire nécessitant un traitement à la fois de droite et de feuilles de gauche pour chaque niveau de segments de mémoire et une branche ou l'autre manière plutôt qu'une seule comparaison par niveau pour l'arbre pliant avec généralement la branche moins profonde celle prise. Par rapport aux autres PQ et la carte en fonction des algorithmes fonctionnels, des améliorations sont en général d'un facteur constant dans la réduction du nombre d'opérations O (log n) en traversant chaque niveau de la structure d'arbre respective.
Le MinHeap est généralement mis en oeuvre sous forme de tableau mutable tas binaire après un modèle d'arbre généalogique inventé par Michael Eytzinger il y a plus de 400 ans. Je sais que la question a dit qu'il n'y avait pas d'intérêt dans le code mutable non fonctionnel, mais si l'on doit éviter tout sous-code que les utilisations de mutabilité, alors nous ne pouvions pas utiliser de liste ou utilisant la mutabilité LazyList « sous les couvertures » pour des raisons de performance. Alors, imaginez que la version mutable alternative suivante du PQ MinHeap est tel que fourni par une bibliothèque et profiter un autre facteur de plus de deux pour les gammes plus grandes prime de performance:
[<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
Note Geek: J'avais en fait attendre la version mutable d'offrir un rapport de performance bien meilleure amélioration, mais les tourbières dans les insertions de re-en raison de la structure du code imbriqué if-then-else et le comportement aléatoire du premier abattage sélectif ce qui signifie que les valeurs de la prédiction de branchement CPU échoue pour une grande partie des branches qui entraîne de nombreuses années 10 de cycles d'horloge CPU supplémentaires par réduction abattage sélectif reconstruit le cache pré-instruction. chercher
Les seules autres gains de performance de facteur constant sur cet algorithme serait la segmentation et l'utilisation de multi-tâches pour un gain de performance proportionnel au nombre de cœurs de processeur; Cependant, tel qu'il est, c'est l'algorithme fonctionnel le plus rapide SoE pur à ce jour, et même la forme fonctionnelle pure en utilisant la MinHeap fonctionnelle bat des implémentations impératives simplistes telles que Code de Jon Harrop ou crible d'Atkin de Johan Kullbom (qui est en erreurson temps comme il calculé uniquement nombres premiers à 10 millions plutôt que le premier 10 millionième ), mais ces algorithmes serait environ cinq fois plus rapidement si de meilleures optimisations ont été utilisées. Ce rapport d'environ cinq entre le code fonctionnel et impératif sera quelque peu réduite lorsque l'on ajoute le multi-threading de factorisation de roue plus grande que la complexité de calcul des augmentations de code impératives plus rapide que le code fonctionnel et multi-threading aide le code fonctionnel plus lent plus que le un code plus rapide impératif que celui-ci se rapproche de la limite de base du temps nécessaire pour énumérer les nombres premiers trouvés.
Edit_add: Même si l'on pouvait choisir de continuer à utiliser la version fonctionnelle pure de MinHeap, ajoutant efficace segmentation en préparation pour le multi-threading serait légèrement "casser" la « pureté » du code fonctionnel comme suit: 1) la façon la plus efficace de transférer une représentation de nombres premiers composites abattus est une matrice bits emballé la taille du segment, 2) Bien que la taille de la matrice est connue, en utilisant un la compréhension du tableau pour l'initialiser d'une manière fonctionnelle est pas efficace car il utilise « ResizeArray » sous les couvertures qui a besoin de se copier pour tous les ajouts x (je pense que « x » est de huit pour la mise en œuvre actuelle) et en utilisant Array.init doesn « travail t autant de valeurs à des indices particuliers sont ignorés, 3) par conséquent, la meilleure façon de remplir le tableau cueilli-composite est à zeroCreate il de la bonne taille, puis exécuter une fonction d'initialisation qui pourrait écrire à chaque index de tableau mutable plus une fois. Bien que ce ne soit pas strictement « fonctionnelle », il est proche en ce que le tableau est initialisé et jamais modifié à nouveau.
Le code de segmentation ajouté, multi-threading, la circonférence factoriel de roue programmable, et de nombreux réglages de performances est la suivante (à l'exception des nouvelles constantes ajoutés, le code oscillant supplémentaire pour mettre en oeuvre la segmentation et multi-threading est le fond près de la moitié du code à partir de la fonction "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)) }
Notez que les modules MinHeap, à la fois fonctionnelle et basée sur la baie, ont eu une fonction « ajuster » ajoutée pour permettre le réglage de l'état de réforme de la version de chaque fil du PQ au début de chaque nouvelle page de segment. Notez également qu'il est possible d'ajuster le code de sorte que la plupart du calcul est effectué en utilisant 32 plages de bits avec la sortie de la séquence finale uint64 est à peu de frais dans le temps de calcul de sorte qu'actuellement la portée théorique est un peu plus de 100000000000000 (dix élevé à quatorze pouvoir) si l'on veut attendre les trois à quatre mois nécessaires pour calculer cette gamme. Les contrôles de portée numériques ont été supprimés car il est peu probable que quiconque utiliserait cet algorithme pour calculer jusqu'à cette plage et encore moins le dépasser.
Utilisation du MinHeap pur fonctionnelle et 2,3,5,7 roue factorisation, le programme ci-dessus calcule le premier cent mille, un million, dix millions et cent millions de nombres premiers à 0,062, 0,629, 10,53 et 195,62 secondes, respectivement. Utilisation de la MinHeap basée sur la baie des vitesses jusqu'à ce 0,097, 0,276, 3,48 et 51,60 secondes respectivement. En utilisant la molette 2,3,5,7,11,13,17 en changeant WHLPRMS à « [| 2u, 3U, 5U, 7U, 11u, 13u, 17u |] » et FSTPRM à 19U vitesses jusqu'à encore un peu plus de 0,181, 0,308, 2,49, et 36,58 secondes, respectivement (pour améliorer le facteur constant avec une surcharge constant). Ce tweak le plus rapide calcule les nombres premiers 203,280,221 dans la plage de numéros 32 bits dans environ 88,37 secondes. La constante « BFSZ » peut être ajustée avec des compromis entre les temps plus lents pour la version plus petite gammes fois plus rapide pour les gammes plus grandes, avec une valeur de « 1 <<< 14 » recommandé d'être jugé pour les plus grandes plages. Cette constante définit uniquement la taille de tampon minimale, avec le programme de réglage de la taille de la mémoire tampon au-dessus de cette taille automatiquement pour des plages plus grandes tellesque le tampon est si suffisante pour que le plus grand nombre premier de base nécessaire pour la gamme page sera toujours « grève » chaque page au moins une fois; cela signifie que la complexité et les frais généraux d'un « tamis à godets » supplémentaires ne sont pas nécessaires. Cette dernière version entièrement optimisé peut calculer les nombres premiers jusqu'à 10 et 100 milliards environ 256,8 et 3617.4 secondes (un peu plus d'une heure pour les 100 milliards) comme testé à l'aide « primesPQOWSE () |> Seq.takeWhile ((> =) 100000000000UL) |> Seq.fold (fun sp -> s + 1UL) 0UL » pour la sortie. C'est là que les estimations d'environ une demi-journée pour le nombre de nombres premiers à un billion, une semaine pour un maximum de dix billions et environ trois à quatre mois pour jusqu'à cent billions viennent.
Je ne pense pas qu'il est possible de rendre le code fonctionnel ou presque fonctionnelle en utilisant l'algorithme soe supplémentaire pour exécuter beaucoup plus vite que cela. Comme on peut le voir en regardant le code, l'optimisation de l'algorithme incrémental de base a ajouté beaucoup à la complexité du code tel qu'il est probablement un peu plus complexe que équivalente code optimisé basé sur la matrice droite culling avec ce code capable de courir fois environ dix plus vite que cela et sans l'exposant supplémentaire dans la performance qui signifie que ce code fonctionnel supplémentaire a une surcharge de pourcentage supplémentaire de plus en plus.
est cet autre utile que d'un point de vue théorique et intellectuel intéressant? Probablement ce n'est pas. Pour les petites gammes de nombres premiers jusqu'à environ dix millions, le meilleur des années REE fonctionnels supplémentaires de base pas entièrement optimisés sont probablement adéquats et assez simple d'écrire ou ont moins l'utilisation de mémoire RAM que le plus simple impératif SOE. Cependant, ils sont beaucoup plus lents que le code plus impératif à l'aide d'un tableau afin qu'ils « à court de vapeur » pour les gammes ci-dessus que. Bien qu'il ait été démontré ici que le code peut être accéléré par l'optimisation, il est encore 10 années de fois plus lent que une version de tableau pur plus impératif encore a ajouté à la complexité d'être au moins aussi complexe que ce code avec des optimisations équivalentes , et même que le code sous F # sur DotNet est environ quatre fois plus lent que d'utiliser un langage tel que C ++ compilé directement du code natif; si l'on voulait vraiment enquêter sur de grandes gammes de nombres premiers, on utiliserait probablement l'une de ces autres langues et techniques où primesieve peut calculer le nombre de nombres premiers dans la 100000000000000 gamme en moins de quatre heures au lieu des trois mois requis pour ce code. END_EDIT_ADD
Voici une à peu près au maximum optimisé pour carte basée Sieve d'Eratosthène en utilisant des séquences supplémentaires de l'algorithme (et récursive) car il n'y a pas besoin de memoization des valeurs de séquence précédentes (autres que il y a un léger avantage pour la mise en cache des premiers valeurs de base en utilisant Seq.cache), avec les principales optimisations étant qu'il utilise roue factorisation pour la séquence d'entrée et qu'il utilise plusieurs (récursif) flux pour maintenir les nombres premiers de base qui sont inférieures à la racine carrée du dernier numéro étant tamisé, de la manière suivante :
let primesMPWSE =
let whlptrn = [| 2;4;2;4;6;2;6;4;2;4;6;6;2;6;4;2;6;4;6;8;4;2;4;2;
4;8;6;4;6;2;4;6;2;6;6;4;2;4;6;2;6;4;2;4;2;10;2;10 |]
let adv i = if i < 47 then i + 1 else 0
let reinsert oldcmpst mp (prime,pi) =
let cmpst = oldcmpst + whlptrn.[pi] * prime
match Map.tryFind cmpst mp with
| None -> mp |> Map.add cmpst [(prime,adv pi)]
| Some(facts) -> mp |> Map.add cmpst ((prime,adv pi)::facts)
let rec mkprimes (n,i) m ps q =
let nxt = n + whlptrn.[i]
match Map.tryFind n m with
| None -> if n < q then seq { yield (n,i); yield! mkprimes (nxt,adv i) m ps q }
else let (np,npi),nlst = Seq.head ps,ps |> Seq.skip 1
let (nhd,ni),nxtcmpst = Seq.head nlst,n + whlptrn.[npi] * np
mkprimes (nxt,adv i) (Map.add nxtcmpst [(np,adv npi)] m) nlst (nhd * nhd)
| Some(skips) -> let adjmap = skips |> List.fold (reinsert n) (m |> Map.remove n)
mkprimes (nxt,adv i) adjmap ps q
let rec prs = seq {yield (11,0); yield! mkprimes (13,1) Map.empty prs 121 } |> Seq.cache
seq { yield 2; yield 3; yield 5; yield 7; yield! mkprimes (11,0) Map.empty prs 121 |> Seq.map (fun (p,i) -> p) }
Il trouve le 100.000ème nombres premiers jusqu'à 1.299.721 dans environ un 0,445 secondes, mais ne pas être un EoS approprié impératif algorithme n'échelle pas près linéaire avec un nombre croissant de nombres premiers, prend 7.775 secondes pour trouver le premier 1000000 jusqu'à 15.485.867 pour une performance sur cette plage d'environ O (n ^ 1.2) où n est le premier au maximum Trouvées.
Il y a un peu plus l'accord qui pourrait être fait, mais il ne va probablement pas faire beaucoup de différence à un grand pourcentage de meilleures performances comme suit:
-
Comme la bibliothèque de séquence F # est nettement lent, on pourrait utiliser un type d'auto défini qui implémente IEnumerable pour réduire le temps passé dans la séquence interne, mais que les opérations de séquence ne prennent que 20% de la durée totale, même si ceux-ci ont été réduits à zéro du temps le résultat ne serait qu'une réduction de 80% du temps.
-
D'autres formes de stockage de carte pourrait être jugé comme une file d'attente prioritaire comme mentionné par O'Neil ou SkewBinomialHeap tel qu'il est utilisé par @gradbot, mais au moins pour le SkewBinomialHeap, l'amélioration des performances est seulement quelques pour cent . Il semble que dans le choix de différentes implémentations de carte, on est juste commercent meilleure réponse à la recherche et la suppression d'éléments qui sont près du début de la liste contre le temps passé à insérer de nouvelles entrées afin de profiter de ces avantages si le gain net est à peu près un lavage et a encore une performance de O (log n) avec des entrées de plus en plus la carte. L'optimisation ci-dessus en utilisant des flux multiples d'entrées juste à la racine carrée de réduire le nombre d'entrées dans la carte et donc apporter ces améliorations de pas beaucoup d'importance.
Edit_add: Je l'ai fait faire le petit plus de l'optimisation et la performance quelque peu améliorée plus que prévu, probablement en raison de la façon dont l'amélioration de l'élimination de la Seq.skip comme un moyen de faire avancer à travers la séquence de nombres premiers base. Cette optimisation utilise un remplaçant pour la génération de séquence interne comme un tuple de valeur entière et une fonction de maintien servant à l'avance à la valeur suivante dans la séquence, avec la séquence finale F # généré par une opération globale se dérouler. Code se présente comme suit:
type SeqDesc<'a> = SeqDesc of 'a * (unit -> SeqDesc<'a>) //a self referring tuple type
let primesMPWSE =
let whlptrn = [| 2;4;2;4;6;2;6;4;2;4;6;6;2;6;4;2;6;4;6;8;4;2;4;2;
4;8;6;4;6;2;4;6;2;6;6;4;2;4;6;2;6;4;2;4;2;10;2;10 |]
let inline adv i = if i < 47 then i + 1 else 0
let reinsert oldcmpst mp (prime,pi) =
let cmpst = oldcmpst + whlptrn.[pi] * prime
match Map.tryFind cmpst mp with
| None -> mp |> Map.add cmpst [(prime,adv pi)]
| Some(facts) -> mp |> Map.add cmpst ((prime,adv pi)::facts)
let rec mkprimes (n,i) m (SeqDesc((np,npi),nsdf) as psd) q =
let nxt = n + whlptrn.[i]
match Map.tryFind n m with
| None -> if n < q then SeqDesc((n,i),fun() -> mkprimes (nxt,adv i) m psd q)
else let (SeqDesc((nhd,x),ntl) as nsd),nxtcmpst = nsdf(),n + whlptrn.[npi] * np
mkprimes (nxt,adv i) (Map.add nxtcmpst [(np,adv npi)] m) nsd (nhd * nhd)
| Some(skips) -> let adjdmap = skips |> List.fold (reinsert n) (m |> Map.remove n)
mkprimes (nxt,adv i) adjdmap psd q
let rec prs = SeqDesc((11,0),fun() -> mkprimes (13,1) Map.empty prs 121 )
let genseq sd = Seq.unfold (fun (SeqDesc((n,i),tailfunc)) -> Some(n,tailfunc())) sd
seq { yield 2; yield 3; yield 5; yield 7; yield! mkprimes (11,0) Map.empty prs 121 |> genseq }
Les temps requis pour trouver le 100.000ème et les nombres premiers sont 1.000.000ème environ 0,31 et 5,1 secondes, respectivement, donc il y a un gain en pourcentage considérable pour ce petit changement. J'ai essayé ma propre implémentation des interfaces IEnumerable / IEnumerator qui sont la base de séquences, et même si elles sont plus rapides que les versions utilisées par le module de Seq ils font guère de différence plus à cet algorithme, où la plupart du temps est consacré à la carte fonctions. END_EDIT_ADD
Autre que la carte à base de mises en œuvre de EoS supplémentaires, il y a une autre « pure fonctionnelle » mise en œuvre en utilisant l'arbre de pliage qui est dit être un peu plus rapide, mais il a encore un terme O (log n) dans l'arbre de pliage Je soupçonne que ce sera principalement plus rapide (si elle est) en raison de la façon dont l'algorithme est mis en œuvre pour le nombre d'opérations informatiques par rapport à l'aide d'une carte. Si les gens sont intéressés je développerai cette version.
En fin de compte, il faut accepter qu'aucune mise en œuvre fonctionnelle pure des EoS supplémentaires ne viendra jamais proche de la vitesse de traitement brute d'une bonne mise en œuvre impératif pour numeri plus grandesgammes cal. Toutefois, on pourrait arriver à une approche où tout le code est purement fonctionnelle à l'exception du segmentée tamisage de nombres composés sur une plage au moyen d'un réseau (mutable) qui viendrait à proximité de O (n) la performance et dans l'utilisation pratique serait cinquante à une centaine de fois plus rapide que les algorithmes fonctionnels pour les grandes gammes telles que les premiers nombres premiers 200.000.000. Cela a été fait par son blog , mais cela pourrait être accordé plus loin avec très peu de code supplémentaire.
Voici ma tentative d'une traduction raisonnablement fidèle du code Haskell à F #:
#r "FSharp.PowerPack"
module Map =
let insertWith f k v m =
let v = if Map.containsKey k m then f m.[k] v else v
Map.add k v m
let sieve =
let rec sieve' map = function
| LazyList.Nil -> Seq.empty
| LazyList.Cons(x,xs) ->
if Map.containsKey x map then
let facts = map.[x]
let map = Map.remove x map
let reinsert m p = Map.insertWith (@) (x+p) [p] m
sieve' (List.fold reinsert map facts) xs
else
seq {
yield x
yield! sieve' (Map.add (x*x) [x] map) xs
}
fun s -> sieve' Map.empty (LazyList.ofSeq s)
let rec upFrom i =
seq {
yield i
yield! upFrom (i+1)
}
let primes = sieve (upFrom 2)
Prime sur tamis mis en œuvre avec des processeurs de boîte aux lettres:
let (<--) (mb : MailboxProcessor<'a>) (message : 'a) = mb.Post(message)
let (<-->) (mb : MailboxProcessor<'a>) (f : AsyncReplyChannel<'b> -> 'a) = mb.PostAndAsyncReply f
type 'a seqMsg =
| Next of AsyncReplyChannel<'a>
type PrimeSieve() =
let counter(init) =
MailboxProcessor.Start(fun inbox ->
let rec loop n =
async { let! msg = inbox.Receive()
match msg with
| Next(reply) ->
reply.Reply(n)
return! loop(n + 1) }
loop init)
let filter(c : MailboxProcessor<'a seqMsg>, pred) =
MailboxProcessor.Start(fun inbox ->
let rec loop() =
async {
let! msg = inbox.Receive()
match msg with
| Next(reply) ->
let rec filter prime =
if pred prime then async { return prime }
else async {
let! next = c <--> Next
return! filter next }
let! next = c <--> Next
let! prime = filter next
reply.Reply(prime)
return! loop()
}
loop()
)
let processor = MailboxProcessor.Start(fun inbox ->
let rec loop (oldFilter : MailboxProcessor<int seqMsg>) prime =
async {
let! msg = inbox.Receive()
match msg with
| Next(reply) ->
reply.Reply(prime)
let newFilter = filter(oldFilter, (fun x -> x % prime <> 0))
let! newPrime = oldFilter <--> Next
return! loop newFilter newPrime
}
loop (counter(3)) 2)
member this.Next() = processor.PostAndReply( (fun reply -> Next(reply)), timeout = 2000)
static member upto max =
let p = PrimeSieve()
Seq.initInfinite (fun _ -> p.Next())
|> Seq.takeWhile (fun prime -> prime <= max)
|> Seq.toList
Voici mes deux cents, bien que je ne suis pas sûr qu'il répond au critère de l'OP pour être le tamis de truely Eratosthène. Il n'utilise pas la division modulaire et met en œuvre une optimisation du papier cité par l'OP. Il ne fonctionne que pour les listes finies, mais qui me semble être dans l'esprit de la façon dont le tamis a été décrit. En aparté, le document parle de la complexiety en termes de nombre de marques et le nombre de divisions. Il semble que, comme nous devons parcourir une liste liée, que cette ignorant peut-être certains aspects clés des différents algorithmes en termes de performance. En général cependant division modulaire avec des ordinateurs est une opération coûteuse.
open System
let rec sieve list =
let rec helper list2 prime next =
match list2 with
| number::tail ->
if number< next then
number::helper tail prime next
else
if number = next then
helper tail prime (next+prime)
else
helper (number::tail) prime (next+prime)
| []->[]
match list with
| head::tail->
head::sieve (helper tail head (head*head))
| []->[]
let step1=sieve [2..100]
EDIT: Correction d'une erreur dans le code de mon poste d'origine. J'ai essayé le suivi de la logique originale du tamis avec quelques modifications. commencer à savoir avec le premier élément et rayez les multiples de cet élément de l'ensemble. Cet algorithme semble littéralement l'élément suivant qui est un multiple du premier au lieu de faire la division modulaire sur chaque numéro dans l'ensemble. Une optimisation du papier est qu'il commence à chercher des multiples du premier supérieur à p ^ 2.
La partie de la fonction d'aide avec les accords multi-niveaux avec la possibilité que le prochain multiple du premier pourrait déjà être retiré de la liste. Ainsi, par exemple avec le 5 prime, il va essayer de supprimer le numéro 30, mais il ne le trouvera jamais parce qu'il a déjà été retiré par le premier 3. L'espoir que la logique de l'clarifie algorithme.
Pour ce que ça vaut, ce n'est pas un tamis de Erathothenes, mais son très rapide:
let is_prime n =
let maxFactor = int64(sqrt(float n))
let rec loop testPrime tog =
if testPrime > maxFactor then true
elif n % testPrime = 0L then false
else loop (testPrime + tog) (6L - tog)
if n = 2L || n = 3L || n = 5L then true
elif n <= 1L || n % 2L = 0L || n % 3L = 0L || n % 5L = 0L then false
else loop 7L 4L
let primes =
seq {
yield 2L;
yield 3L;
yield 5L;
yield! (7L, 4L) |> Seq.unfold (fun (p, tog) -> Some(p, (p + tog, 6L - tog)))
}
|> Seq.filter is_prime
Il trouve le premier 100.000ème en 1,25 secondes sur ma machine (AMD Phenom II, 3.2GHz quadcore).
Je sais que vous explicitement déclaré que vous étiez intéressé par une mise en œuvre tamis purement fonctionnelle, donc je tenais hors présenter mon tamis jusqu'à présent. Mais en relisant le document que vous avez mentionné, je vois l'algorithme présenté tamis supplémentaire, il est essentiellement la même que la mienne, la seule différence étant les détails de mise en œuvre en utilisant des techniques purement fonctionnelles par rapport à des techniques décidément impératives. Donc, je pense que je au moins la moitié remplissant les conditions nécessaires pour satisfaire votre curiosité. De plus, je dirais que l'utilisation des techniques impératives lorsque des gains de performance peuvent être réalisés, mais cachés par des interfaces fonctionnelles est l'une des techniques les plus puissantes encouragées dans la programmation F #, par opposition à la culture tout pur Haskell. J'ai d'abord publié cette mise en œuvre sur mon Projet Euler F #un un blog mais republier ici avec le dos substitué le code pré-requis et le typage structurel enlevé. primes peut calculer les 100.000 premiers nombres premiers en 0.248 secondes et les premiers 1.000.000 nombres premiers en 4,8 secondes sur mon ordinateur (notez que met en cache primes ses résultats afin que vous aurez besoin de réévaluer chaque fois que vous effectuez une référence).
let inline infiniteRange start skip =
seq {
let n = ref start
while true do
yield n.contents
n.contents <- n.contents + skip
}
///p is "prime", s=p*p, c is "multiplier", m=c*p
type SievePrime<'a> = {mutable c:'a ; p:'a ; mutable m:'a ; s:'a}
///A cached, infinite sequence of primes
let primes =
let primeList = ResizeArray<_>()
primeList.Add({c=3 ; p=3 ; m=9 ; s=9})
//test whether n is composite, if not add it to the primeList and return false
let isComposite n =
let rec loop i =
let sp = primeList.[i]
while sp.m < n do
sp.c <- sp.c+1
sp.m <- sp.c*sp.p
if sp.m = n then true
elif i = (primeList.Count-1) || sp.s > n then
primeList.Add({c=n ; p=n ; m=n*n ; s=n*n})
false
else loop (i+1)
loop 0
seq {
yield 2 ; yield 3
//yield the cached results
for i in 1..primeList.Count-1 do
yield primeList.[i].p
yield! infiniteRange (primeList.[primeList.Count-1].p + 2) 2
|> Seq.filter (isComposite>>not)
}
Voici encore une autre méthode de réalisation de la Sieve d'Eratosthène (SOE) en utilisant uniquement le code incrémental pur F # fonctionnelle. Il est adapté à partir du code Haskell développé comme « Cette idée est due à Dave Bayer, mais il a utilisé une formulation complexe et structure arborescente ternaire équilibrée, l'approfondissement progressif de manière uniforme (formulation simplifiée et biaisée, l'approfondissement de la structure d'arbre binaire droite introduit par Heinrich Apfelmus, encore simplifié par Will Ness) idée de la production en raison de M. Mise en scène O'Neill » comme par le lien suivant: Optimized Arbre de pliage de code en utilisant une roue factoriel Haskell.
Le code suivant a plusieurs optimisations qui le rendent plus approprié pour l'exécution en F #, comme suit:
-
Le code utilise des flux CoInductive au lieu de son LazyList que cet algorithme n'a pas (ou peu) besoin de memoization de LazyList et mes flux CoInductive sont plus efficaces que soit LazyLists (du FSharp.PowerPack) ou la construction dans les séquences. Un autre avantage est que mon code peut être exécuté sur tryFSharp.org et ideone.com sans avoir à copier et coller dans le code source de base Microsoft.FSharp.PowerPack pour le type de LazyList et le module (ainsi que avec mention de copyright)
-
Il a découvert qu'il y est tout à fait beaucoup de frais généraux pour le filtrage de F # sur les paramètres de la fonction de sorte que le précédent type union discriminée plus lisible en utilisant tuples a été sacrifiée au profit de struct par valeur (ou classe fonctionne plus rapidement sur certaines plates-formes) types pour environ un facteur de deux ou plus de vitesse jusqu'à.
-
Les optimisations de Will Ness allant de pliage d'arbre linéaire au pliage bilatéral pliage multivoies et des améliorations à l'aide de la factorisation de roue sont à peu près aussi efficace ratiometrically pour F # comme ils sont pour Haskell, la principale différence entre les deux langues étant que Haskell peut être compilé en code natif et a un compilateur plus optimisé alors que F # a courir plus de frais généraux dans le cadre du système DotNet Framework.
type prmstate = struct val p:uint32 val pi:byte new (prm,pndx) = { p = prm; pi = pndx } end type prmsSeqDesc = struct val v:prmstate val cont:unit->prmsSeqDesc new(ps,np) = { v = ps; cont = np } end type cmpststate = struct val cv:uint32 val ci:byte val cp:uint32 new (strt,ndx,prm) = {cv = strt;ci = ndx;cp = prm} end type cmpstsSeqDesc = struct val v:cmpststate val cont:unit->cmpstsSeqDesc new (cs,nc) = { v = cs; cont = nc } end type allcmpsts = struct val v:cmpstsSeqDesc val cont:unit->allcmpsts new (csd,ncsdf) = { v=csd;cont=ncsdf } end let primesTFWSE = 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; 4uy;8uy;6uy;4uy;6uy;2uy;4uy;6uy;2uy;6uy;6uy;4uy;2uy;4uy;6uy;2uy;6uy;4uy;2uy;4uy;2uy;10uy;2uy;10uy |] let inline whladv i = if i < 47uy then i + 1uy else 0uy let inline advmltpl c ci p = cmpststate (c + uint32 whlptrn.[int ci] * p,whladv ci,p) let rec pmltpls cs = cmpstsSeqDesc(cs,fun() -> pmltpls (advmltpl cs.cv cs.ci cs.cp)) let rec allmltpls (psd:prmsSeqDesc) = allcmpsts(pmltpls (cmpststate(psd.v.p*psd.v.p,psd.v.pi,psd.v.p)),fun() -> allmltpls (psd.cont())) let rec (^) (xs:cmpstsSeqDesc) (ys:cmpstsSeqDesc) = //union op for SeqDesc's match compare xs.v.cv ys.v.cv with | -1 -> cmpstsSeqDesc (xs.v,fun() -> xs.cont() ^ ys) | 0 -> cmpstsSeqDesc (xs.v,fun() -> xs.cont() ^ ys.cont()) | _ -> cmpstsSeqDesc(ys.v,fun() -> xs ^ ys.cont()) //must be greater than let rec pairs (csdsd:allcmpsts) = let ys = csdsd.cont in allcmpsts(cmpstsSeqDesc(csdsd.v.v,fun()->csdsd.v.cont()^ys().v),fun()->pairs (ys().cont())) let rec joinT3 (csdsd:allcmpsts) = cmpstsSeqDesc(csdsd.v.v,fun()-> let ys = csdsd.cont() in let zs = ys.cont() in (csdsd.v.cont()^(ys.v^zs.v))^joinT3 (pairs (zs.cont()))) let rec mkprimes (ps:prmstate) (csd:cmpstsSeqDesc) = let nxt = ps.p + uint32 whlptrn.[int ps.pi] if ps.p >= csd.v.cv then mkprimes (prmstate(nxt,whladv ps.pi)) (csd.cont()) //minus function else prmsSeqDesc(prmstate(ps.p,ps.pi),fun() -> mkprimes (prmstate(nxt,whladv ps.pi)) csd) let rec baseprimes = prmsSeqDesc(prmstate(11u,0uy),fun() -> mkprimes (prmstate(13u,1uy)) initcmpsts) and initcmpsts = joinT3 (allmltpls baseprimes) let genseq sd = Seq.unfold (fun (psd:prmsSeqDesc) -> Some(psd.v.p,psd.cont())) sd seq { yield 2u; yield 3u; yield 5u; yield 7u; yield! mkprimes (prmstate(11u,0uy)) initcmpsts |> genseq } primesLMWSE |> Seq.nth 100000
Edit_add: Ceci a été corrigé en tant que le code d'origine ne gère pas correctement la queue du train et passa la queue du train de paramètre aux paires servent à la fonction joinT3 plutôt que ce qui suit de la queue le flux ZS. Le calendrier ci-dessous a été corrigé en conséquence, avec environ une vitesse de 30% supplémentaire vers le haut. Les codes de lien tryFSharp et ideone ont également été corrigées. END_EDIT_ADD
Le programme ci-dessus fonctionne à environ O (n ^ 1,1) performance avec n le premier maximum calculée soit environ O (n ^ 1,18) lorsque n est le nombre des nombres premiers calculés, et prend environ 2,16 secondes pour calculer le premier million de nombres premiers (environ 0,14 seconde pour les 100.000 premiers nombres premiers) sur un code binaire ordinateur exécutant 64 rapidement en utilisant des types struct plutôt que des classes (il semble que certaines implémentations boîte et Unbox les années struct par valeur dans la formation des fermetures). Je considère que, pour être sur la gamme pratique maximale pour ces purs algorithmes principaux fonctionnels. Ceci est probablement la plus rapide que l'on peut exécuter un algorithme soe fonctionnel pur autre que pour quelques ajustements mineurs de réduire les facteurs constants.
autre que la combinaison segmentation et multi-threading pour partager le calcul entre plusieurs cœurs de processeurs, la plupart des « tweaks » qui pourraient être apportées à cet algorithme sont en augmentation de la circonférence de la roue factorisation pour un gain allant jusqu'à environ 40 % des gains de rendement et mineurs dus à des réglages quant à l'utilisation des structures, des classes, des tuples ou plusieurs paramètres individuels directs dans le passage de l'état entre les fonctions.
EDIT_ADD2: IAvons fait les optimisations ci-dessus avec le résultat que le code est maintenant presque deux fois plus rapide en raison de l'optimisation de la structure avec l'avantage supplémentaire de l'aide en option circonférences roue de factorisation plus grandes pour la réduction plus faible ajoutée. On notera que le dessous évite de code à l'aide de continuations dans la boucle principale de génération de séquence et ne les utilise, le cas échéant pour les nombres premiers flux de base et l'abattage de nombre composite ultérieur ruisseaux dérivée de ces nombres premiers de base. Le nouveau code est le suivant:
type CIS<'T> = struct val v:'T val cont:unit->CIS<'T> new(v,cont) = { v=v;cont=cont } end //Co-Inductive Steam
let primesTFOWSE =
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 inline advcnd c ci = c + uint32 WHLPTRN.[ci]
let inline advmltpl p (c,ci) = (c + uint32 WHLPTRN.[ci] * p,whladv ci)
let rec pmltpls p cs = CIS(cs,fun() -> pmltpls p (advmltpl p cs))
let rec allmltpls k wi (ps:CIS<_>) =
let nxt = advcnd k wi in let nxti = whladv wi
if k < ps.v then allmltpls nxt nxti ps
else CIS(pmltpls ps.v (ps.v*ps.v,wi),fun() -> allmltpls nxt nxti (ps.cont()))
let rec (^) (xs:CIS<uint32*_>) (ys:CIS<uint32*_>) =
match compare (fst xs.v) (fst ys.v) with //union op for composite CIS's (tuple of cmpst and wheel ndx)
| -1 -> CIS(xs.v,fun() -> xs.cont() ^ ys)
| 0 -> CIS(xs.v,fun() -> xs.cont() ^ ys.cont())
| _ -> CIS(ys.v,fun() -> xs ^ ys.cont()) //must be greater than
let rec pairs (xs:CIS<CIS<_>>) =
let ys = xs.cont() in CIS(CIS(xs.v.v,fun()->xs.v.cont()^ys.v),fun()->pairs (ys.cont()))
let rec joinT3 (xs:CIS<CIS<_>>) = CIS(xs.v.v,fun()->
let ys = xs.cont() in let zs = ys.cont() in (xs.v.cont()^(ys.v^zs.v))^joinT3 (pairs (zs.cont())))
let rec mkprm (cnd,cndi,(csd:CIS<uint32*_>)) =
let nxt = advcnd cnd cndi in let nxti = whladv cndi
if cnd >= fst csd.v then mkprm (nxt,nxti,csd.cont()) //minus function
else (cnd,cndi,(nxt,nxti,csd))
let rec pCIS p pi cont = CIS(p,fun()->let (np,npi,ncont)=mkprm cont in pCIS np npi ncont)
let rec baseprimes() = CIS(FSTPRM,fun()->let np,npi = advcnd FSTPRM 0,whladv 0
pCIS np npi (advcnd np npi,whladv npi,initcmpsts))
and initcmpsts = joinT3 (allmltpls FSTPRM 0 (baseprimes()))
let inline genseq sd = Seq.unfold (fun (p,pi,cont) -> Some(p,mkprm cont)) sd
seq { yield! WHLPRMS; yield! mkprm (FSTPRM,0,initcmpsts) |> genseq }
Le code ci-dessus prend environ 0,07, 1,02 et 14,58 secondes pour énumérer le premier cent mille, million, et dix millions de nombres premiers, respectivement, toutes sur la référence de la machine Intel i7-2700K (3.5 GHz) en mode 64 bits. Ce n'est pas beaucoup plus lent que la référence à partir de laquelle la mise en œuvre Haskell ce code a été dérivé, bien qu'il soit un peu plus lent sur et href="http://ideone.com/hgQHVj" rel="nofollow"> ideone en raison d'être en mode 32 bits pour tryfsharp sous Silverlight (environ la moitié à nouveau plus lent) et exécuté sur une machine plus lente sous Mono 2.0 (qui est par nature beaucoup plus lente pour F #) sur ideone alors jusqu'à environ cinq fois plus lente que la machine de référence. Notez que le temps d'exécution rapporté par ideone inclut le temps d'initialisation pour l'aspect intégré des tableaux de table, qui a besoin de temps pour être pris en compte.
Le programme ci-dessus a la particularité en outre que la roue de factorisation est paramétré de telle sorte que, par exemple, on peut utiliser une roue extrêmement importante en mettant WHLPRMS à [| 2u, 3U, 5U, 7U, 11u, 13u, 17u, 19u |] et FSTPRM à 23U pour obtenir un temps de fonctionnement d'environ deux tiers pour les grandes plages à environ 10,02 secondes pour dix millions de nombres premiers, bien qu'il faille noter qu'il faut plusieurs secondes pour calculer le WHLPTRN avant que le programme commence à courir.
Note Geek: Je n'ai pas mis en place un « non-partage fixpoint combinateur pour télescoper la production de nombres premiers en plusieurs étapes » selon la référence du code Haskell, bien que j'ai essayé de le faire, parce qu'on a besoin d'avoir quelque chose comme la liste paresseuse de Haskell pour que cela travail sans fuir dans une boucle infinie et le débordement pile. Bien que mes co-inductifs cours d'eau (CIS de) ont des propriétés de la paresse, ils ne sont pas formellement des listes ou des séquences mises en cache paresseux (ils deviennent des séquences non-cache et peuvent être mises en cache lorsque passé si une fonction telle que la « genseq » que je fournis pour la séquence de sortie final). Je ne voulais pas utiliser la mise en œuvre de PowerPack de LazyList parce qu'il est pas très efficace et il faudrait que je copie que le code source en, qui ne fournit pas tryfsharp et ideone pour les modules importés. En utilisant les séquences intégré (même en cache) est très inefficace lorsque l'on veut utiliser la tête des opérations / arrière qui sont nécessaires à cet algorithme comme le seul moyen d'obtenir la queue d'une séquence est d'utiliser « Seq.skip 1 » qui, utilisations multiples produit une nouvelle séquence en fonction de la séquence d'origine sauté récursive plusieurs fois. Je pourrais mettre en œuvre ma propre classe LazyList efficace basée sur de la CEI, mais il ne semble guère la peine de démontrer un point où les objets récursifs « initcmpsts » et « baseprimes » prennent très peu de code. En outre, le passage d'un LazyList à une fonction pour produire des extensions à cette LazyList qui fonctionnent utilise uniquement des valeurs de près du début de la Lazylist exige que presque toute LazyList est memoized une réduction de l'efficacité de la mémoire: un laissez-passer pour les 10 millions de premiers nombres premiers il faudra une LazyList en mémoire avec près de 180 millions d'éléments. Je pris alors une passe à ce sujet.
Notez que pour les gammes plus larges (10 millions de nombres premiers ou plus), ce code purement fonctionnelle est à peu près la même vitesse que de nombreuses implémentations impératives simplistes du Crible d'Eratosthène ou Atkins, bien que ce soit en raison du manque d'optimisation de ces impératif algorithmes; un moisre mise en œuvre impératif que cela en utilisant des optimisations équivalentes et des réseaux segmentés TAMISAGE seront encore environ dix fois plus rapide que ce que par ma réponse « presque fonctionnelle ».
Notez également que si il est possible de réaliser en utilisant le pliage des arbres tamisant segmenté, il est plus difficile car les algorithmes d'avance cueillant sont enterrés dans les continuations utilisés pour la « union - ^ » opérateur et travailler autour cela signifierait qu'une façon continue avancer gamme abattage doit être utilisé; cela est contrairement à d'autres algorithmes où l'état de la variable d'abattage peut être remise à zéro pour chaque nouvelle page, y compris avec leur gamme réduite, de sorte que si de plus grandes portées de 32 bits sont utilisés, la plage de réforme interne peut toujours être remis à opérer au sein de la 32 -bit même si une gamme gamme de nombres premiers est déterminée à peu de frais 64 bits dans le temps d'exécution par le premier. END_EDIT_ADD2
En fait, j'essayé de faire la même chose, j'ai essayé d'abord la même mise en œuvre naïve comme question, mais il était trop lent. Je me suis alors trouvé cette page YAPES: Problème Sept, Partie 2 , où il a utilisé réel crible d'Eratosthène basé sur Melissa E. O'Neill. Je pris le code à partir de là, juste un peu il a modifié, parce que F # a un peu changé depuis la publication.
let reinsert x table prime =
let comp = x+prime
match Map.tryFind comp table with
| None -> table |> Map.add comp [prime]
| Some(facts) -> table |> Map.add comp (prime::facts)
let rec sieve x table =
seq {
match Map.tryFind x table with
| None ->
yield x
yield! sieve (x+1I) (table |> Map.add (x*x) [x])
| Some(factors) ->
yield! sieve (x+1I) (factors |> List.fold (reinsert x) (table |> Map.remove x)) }
let primes =
sieve 2I Map.empty
primes |> Seq.takeWhile (fun elem -> elem < 2000000I) |> Seq.sum
Cette question demande spécifiquement pour d'autres algorithmes, je fournis la mise en œuvre suivante:
ou sait peut-être des méthodes alternatives de mise en œuvre ou algorithmes tamiser
Non présentation de divers algorithmes d'Eratosthène Sieve (SOE) est vraiment complet sans mentionner Crible d'Atkin (SOA), qui est en fait une variante de SOE en utilisant les solutions à un ensemble d'équations quadratiques pour mettre en œuvre l'abattage composite, ainsi que l'élimination de tous les multiples de carrés des nombres premiers de base (nombre premier inférieur ou égal à la racine carrée du nombre le plus élevé testé pour primalité). En théorie, la SoA est plus efficace que l'entreprise publique en ce qu'il ya des opérations un peu moins sur la plage telle qu'elle devrait avoir avoir environ 20% moins de complexité pour la gamme d'environ 10 à 100 millions, mais en pratique, il est généralement dû plus lents à la les frais généraux de calcul de la complexité de la résolution de plusieurs équations du second degré. Bien que, le très optimisé Daniel J. Bernstein mise en œuvre C est plus rapide que la mise en œuvre SoE contre laquelle il a testé il pour cette gamme particulière de nombres de test , la mise en œuvre SoE contre laquelle il a testé était pas le plus optimal et des versions plus hautement optimisées de SoE droite sont encore plus rapides. Cela semble être le cas ici, même si je reconnais qu'il peut y avoir d'autres optimisations que j'ai manqué.
Depuis O'Neill dans son document sur Soe à l'aide Sieyès sans bornes supplémentaires un but d'abord montrer que le Turner Sieve n'est pas SoE à la fois à l'algorithme et à la performance, elle n'a pas tenu compte d'autres variations de Soe telles que SoA. Faire une recherche rapide de la littérature, je ne trouve aucune application de SoA aux séquences supplémentaires non bornées dont nous parlons ici, donc adapté moi-même comme dans le code suivant.
Tout comme le cas sans bornes pur état de l'environnement peut être considéré comme ayant pour nombres composés d'une séquence non borné de séquences de multiples de nombres premiers, le SoA considère avoir comme nombres premiers potentiels de la séquence sans bornes des séquences non bornées de toutes les expressions de la quadratique équations à l'une des deux variables libres, « x » ou « y » fixé à une valeur de départ et d'une séquence « d'élimination » séparé des séquences de l'ensemble des multiples de nombres premiers de la base, dont la dernière est très semblable au composite des séquences d'élimination de séquences d'état de l'environnement, à l'exception que les séquences avancent plus rapidement par le carré des nombres premiers, plutôt que par un (moindre) multiple des nombres premiers. Je suis efforcé de réduire le nombre de séquences d'équations quadratiques exprimées en reconnaissant que, aux fins d'un tamis supplémentaire, les « 3 * x ^ 2 + y ^ 2 » et « 3 * x ^ 2 - y ^ 2 » séquences sont vraiment la même chose, sauf pour le signe du second terme et l'élimination de toutes les solutions qui ne sont pas bizarre, comme l'application bien 2357 roue factorisation (bien que le SoA a déjà inhérente 235 roue factorisation). Il utilise l'arbre efficace pliage fusion / combinant l'algorithme comme dans l'arbre SoE fusion pour faire face à chaque séquence de séquences, mais avec une simplification que l'opérateur syndical ne se combine pas dans la fusion comme l'algorithme SoA dépend de pouvoir sélectionner l'état premier basé sur la nombre de solutions trouvées quadratique pour une valeur particulière. Le code est plus lent que l'arbre fusion SoE en raison d'environ trois fois le nombre d'opérations aériennes portant sur environ trois fois le nombre de séquences un peu plus complexes, mais il est probable que toute une gamme de tamiser de très grands nombres où il passera SoE en raison de son avantage de performance théorique.
Le code suivant est fidèle à la formulation de la SoA, utilise des types CoInductive cours d'eau plutôt que de LazyList ou de séquences comme memoization est pas nécessaire etla performance est meilleure, aussi ne pas utiliser les syndicats discriminé et évite correspondance de motif pour des raisons de performance:
#nowarn "40"
type cndstate = class val c:uint32 val wi:byte val md12:byte new(cnd,cndwi,mod12) = { c=cnd;wi=cndwi;md12=mod12 } end
type prmsCIS = class val p:uint32 val cont:unit->prmsCIS new(prm,nxtprmf) = { p=prm;cont=nxtprmf } end
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
type allstateCIS<'b> = class val ss:stateCIS<'b> val cont:unit->allstateCIS<'b> new(sbstrm,cont) = { ss=sbstrm;cont=cont } end
let primesTFWSA() =
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;
4uy;8uy;6uy;4uy;6uy;2uy;4uy;6uy;2uy;6uy;6uy;4uy;2uy;4uy;6uy;2uy;6uy;4uy;2uy;4uy;2uy;10uy;2uy;10uy |]
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)
let rec allsqrs (prms:prmsCIS) = let s = prms.p*prms.p in allstateCIS(prmsqrs s s,fun() -> allsqrs (prms.cont()))
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
let vn=v+uint32 a in let vn=if vn<v then 0xFFFFFFFFu else vn in qdrtc vn (y+2))
let rec allqdrtcsX4 x = allstateCIS(qdrtc (((x*x)<<<2)+1u) 1,fun()->allqdrtcsX4 (x+1u))
let rec allqdrtcsX3 x = allstateCIS(qdrtc (((x*(x+1u))<<<1)-1u) (1 - int x),fun() -> allqdrtcsX3 (x+1u))
let rec joinT3 (ass:allstateCIS<'b>) = stateCIS<'b>(ass.ss.v,ass.ss.a,fun()->
let rec (^) (xs:stateCIS<'b>) (ys:stateCIS<'b>) = //union op for CoInductiveStreams
match compare xs.v ys.v with
| 1 -> stateCIS(ys.v,ys.a,fun() -> xs ^ ys.cont())
| _ -> stateCIS(xs.v,xs.a,fun() -> xs.cont() ^ ys) //<= then keep all the values without combining
let rec pairs (ass:allstateCIS<'b>) =
let ys = ass.cont
allstateCIS(stateCIS(ass.ss.v,ass.ss.a,fun()->ass.ss.cont()^ys().ss),fun()->pairs (ys().cont()))
let ys = ass.cont() in let zs = ys.cont() in (ass.ss.cont()^(ys.ss^zs.ss))^joinT3 (pairs (zs.cont())))
let rec mkprm (cs:cndstate) (sqrs:stateCIS<_>) (qX4:stateCIS<_>) (qX3:stateCIS<_>) tgl =
let inline advcnd (cs:cndstate) =
let inline whladv i = if i < 47uy then i + 1uy else 0uy
let inline modadv m a = let md = m + a in if md >= 12uy then md - 12uy else md
let a = WHLPTRN.[int cs.wi] in let nc = cs.c+uint32 a
if nc<cs.c then failwith "Tried to enumerate primes past the numeric range!!!"
else cndstate(nc,whladv cs.wi,modadv cs.md12 a)
if cs.c>=sqrs.v then mkprm (if cs.c=sqrs.v then advcnd cs else cs) (sqrs.cont()) qX4 qX3 false //squarefree function
elif cs.c>qX4.v then mkprm cs sqrs (qX4.cont()) qX3 false
elif cs.c>qX3.v then mkprm cs sqrs qX4 (qX3.cont()) false
else match cs.md12 with
| 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
elif tgl then prmsCIS(cs.c,fun() -> mkprm (advcnd cs) sqrs qX4 qX3 false)
else mkprm (advcnd cs) sqrs qX4 qX3 false
| 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
elif tgl then prmsCIS(cs.c,fun() -> mkprm (advcnd cs) sqrs qX4 qX3 false)
else mkprm (advcnd cs) sqrs qX4 qX3 false
| _ -> if cs.c=qX4.v then mkprm cs sqrs (qX4.cont()) qX3 (not tgl) //always must be 1uy or 5uy
elif tgl then prmsCIS(cs.c,fun() -> mkprm (advcnd cs) sqrs qX4 qX3 false)
else mkprm (advcnd cs) sqrs qX4 qX3 false
let qX4s = joinT3 (allqdrtcsX4 1u) in let qX3s = joinT3 (allqdrtcsX3 1u)
let rec baseprimes = prmsCIS(11u,fun() -> mkprm (cndstate(13u,1uy,1uy)) initsqrs qX4s qX3s false)
and initsqrs = joinT3 (allsqrs baseprimes)
let genseq ps = Seq.unfold (fun (psd:prmsCIS) -> Some(psd.p,psd.cont())) ps
seq { yield 2u; yield 3u; yield 5u; yield 7u;
yield! mkprm (cndstate(11u,0uy,11uy)) initsqrs qX4s qX3s false |> genseq }
Comme indiqué, le code est plus lent que l'arbre pliant roue Optimized SoE telle que publiée dans une autre réponse à une demi-seconde pour les 100.000 premiers nombres premiers, et a à peu près le même O empirique (n ^ 1.2) pour les nombres premiers ont trouvé performance le meilleur des autres pures solutions fonctionnelles. Quelques optimisations supplémentaires qui pourraient être essayées sont que les nombres premiers séquences carrés ne pas utiliser la roue factorisation pour éliminer les 357 multiples des carrés ou même utiliser uniquement les premiers multiples des premiers carrés pour réduire le nombre de valeurs dans les flux de séquence carrés et peut-être d'autres optimisations liées aux courants de séquence d'expression de l'équation quadratique.
Edit_add: Je prend un peu de temps pour se pencher sur les optimisations du SoA et de voir qu'en plus de ce qui précède « quadratfrei » optimisations, ce qui ne sera probablement pas faire beaucoup de différence, que le second degré séquences ont un motif modulo sur chaque 15 éléments qui permettraient un grand nombre des valeurs passées de test composites basculés à une présélection et éliminerait la nécessité de modulo spécifique 12 opérations pour chaque numéro composite. Toutes ces optimisations sont susceptibles d'entraîner une réduction du travail de calcul soumis au pliage des arbres jusqu'à environ 50% pour faire une version légèrement plus optimisée de la course SoA à proximité ou un peu mieux que le meilleur arbre pli fusion SoE. Je ne sais pas quand je pourrais trouver le temps de faire ces quelques jours d'enquête pour déterminer le résultat. END_EDIT_ADD
EDIT_ADD2: En travaillant sur les optimisations ci-dessus qui fait augmenter les performances d'un facteur de deux, je vois pourquoi la performance empirique actuelle avec l'augmentation de n est pas aussi bon que SoE: alors que l'entreprise publique est particulièrement adapté à l'arbre des opérations de pliage en ce que les premières séquences sont plus denses et répéter plus souvent avec des séquences plus tard beaucoup moins denses, les SoA séquences « 4X » sont plus denses pour des séquences plus tard quand ils sont ajoutés et tandis que les séquences « 3X » commencent à moins denses, ils deviennent plus denses que y tend vers zéro, puis obtenir moins dense à nouveau; Cela signifie que les séquences d'appel / retour ne sont pas maintenus à une profondeur minimale que pour SOE mais que cette profondeur augmente au-delà étant proportionnelle à la plage de nombre. Les solutions utilisant le pliage ne sont pas assez comme on peut mettre en oeuvre le pliage gauche pour les séquences que l'augmentation de la densité avec le temps, mais qui laisse encore les parties négatives des séquences « 3X » mal optimisés, comme le fait de briser les séquences « 3X » en positif et parties négatives. La solution la plus simple est susceptible d'enregistrer toutes les séquences à une carte ce qui signifie que le temps d'accès augmentera de quelque chose comme le journal de la racine carrée de la gamme, mais ce sera mieux pour la plage de nombre plus grand que le pliage de l'arborescence. END_EDIT_ADD2
Bien que plus lent, je présente cette solution ici pour montrer comment le code peut être évolué pour exprimer des idées conçues à l'origine du code fonctionnel à impérieusement pur F #. Il fournit des exemples d'utilisation de continuations comme dans CoInductive cours d'eau pour mettre en œuvre la paresse sans utiliser le type Lazy, la mise en œuvre (queue) boucles récursives pour éviter toute exigence de mutabilité, enfiler un accumulateur (TGL) au moyen d'appels récursifs pour obtenir un résultat (le nombre de fois les équations du second degré « frappé » le nombre testé), présentant les solutions d'équations que (séquences paresseux) (ou de flux dans ce cas), etc.
Pour ceux qui voudraient continuer à jouer avec ce code même sans un système de développement basé sur Windows, j'ai posté et href="http://ideone.com/tldRPw" rel="nofollow"> Ideone.com ,même si elle est plus lent sur ces deux plates-formes, à la fois tryfsharp proportionnelle à la vitesse de la machine client local et plus lent en raison de l'exécution sous Silverlight, et Ideone en cours d'exécution sur la CPU du serveur Linux sous Mono-projet 2.0, ce qui est notoirement lent à la fois la mise en œuvre et en particulier aux collections d'ordures.
Je ne pense pas que cette question est complète en considérant que des algorithmes purement fonctionnels que l'état de peau soit une carte ou file d'attente prioritaire dans le cas de quelques réponses ou un arbre de fusion plié dans le cas d'un de mes autres réponses dans que l'un de ceux-ci sont très limitées quant à la facilité d'utilisation pour de grands intervalles de nombres premiers en raison de leur O approximative (n ^ 1,2) rendement ( « ^ » moyen élevé à la puissance de l'endroit où n est le nombre de haut dans la séquence), ainsi que leur temps de calcul pour chaque opération d'abattage. Cela signifie que même pour la plage de nombre de 32 bits, ces algorithmes prendront quelque chose dans la gamme d'au moins plusieurs minutes pour générer les nombres premiers jusqu'à quatre milliards, plus, ce qui est quelque chose qui est très utile.
Il y a eu plusieurs réponses présentant des solutions en utilisant divers degrés de mutabilité, mais ils ont soit pas profité pleinement de leur mutabilité et ont été inefficaces ou ont été seulement les traductions très simplistes de code impératif et fonctionnel laid. Il me semble que le F # tableau mutable est juste une autre forme de cacher l'état mutable l'intérieur d'une structure de données, et qu'un algorithme efficace peut être développé qui a pas d'autre mutabilité utilisé autre que le réseau de tampon mutable utilisé pour l'abattage efficace des nombres composés de segments paginées tampon, le reste du code écrit dans un style fonctionnel pur.
Le code suivant a été développé après avoir vu le code par Jon Harrop , et améliore ces idées comme suit:
-
Code de Jon échoue en termes d'utilisation de la mémoire haute (enregistre tous les nombres premiers générés au lieu de simplement les nombres premiers à la racine carrée du premier candidat le plus élevé, et régénère en permanence tampon des réseaux de plus en plus la taille énorme (égale à la taille du dernier premier Introuvable) quelle que soit la taille du cache du processeur.
-
De plus, son code tel qu'il est présenté n'incluant une séquence de production.
-
En outre, le code tel que présenté ne pas les optimisations au moins ne traitent que des nombres impairs et encore moins sans tenir compte de l'utilisation de l'utilisation factorisation roue.
Si le code de Jon ont été utilisés pour générer la gamme de nombres premiers à la plage de nombres 32 bits de quatre milliards plus, il aurait une exigence de mémoire de giga-octets pour les nombres premiers enregistrés dans la structure de la liste et un autre multi-Gigaoctets pour le tamis tampon, bien qu'il n'y ait aucune raison que ce dernier ne peut pas être d'une taille plus petite fixe. Une fois que le tampon tamis dépasse la taille de la taille du cache du processeur, la performance se détériore rapidement « raclée cache », avec l'augmentation de la perte de performance d'abord la L1, la L2, et enfin les dimensions L3 (si actuelles) sont dépassées.
Ceci est la raison pour laquelle le code de Jon nombres premiers ne calculate jusqu'à environ 25 millions ou même sur ma machine 64 bits avec huit giga-octets de mémoire avant de générer une exception hors de la mémoire et explique aussi pourquoi il y a une plus grande et plus baisse de performance relative que les gammes se supérieur avec environ un O (n ^ 1.4) la performance
avec une gamme de plus en plus et n'est quelque peu sauvé parce qu'il a une telle faible complexité de calcul pour commencer.Les adresses de code suivant toutes ces limitations, en ce qu'elle memoizes seule la base des nombres premiers jusqu'à la racine carrée du nombre maximum dans la plage qui sont calculées selon les besoins (seulement quelques kilo-octets dans le cas du 32 bits numéro de série) et utilise uniquement des tampons de très petites d'environ seize kilo-octets pour chacun du générateur de nombres premiers de base et la page principale segmentés filtre tamis (plus petite que la taille du cache L1 de la plupart CPU) moderne, ainsi que notamment le code de séquence de génération et ( a) étant un peu optimisé pour le tamisage uniquement pour des numéros impairs, ce qui signifie que la mémoire est utilisée de manière plus efficace. En outre, un tableau de bits emballé est utilisé pour futres améliorer l'efficacité de la mémoire; il le coût de calcul est la plupart du temps de compensait en moins des calculs devant être réalisés dans le balayage de la mémoire tampon.
let primesAPF32() =
let rec oddprimes() =
let BUFSZ = 1<<<17 in let buf = Array.zeroCreate (BUFSZ>>>5) in let BUFRNG = uint32 BUFSZ<<<1
let inline testbit i = (buf.[i >>> 5] &&& (1u <<< (i &&& 0x1F))) = 0u
let inline cullbit i = let w = i >>> 5 in buf.[w] <- buf.[w] ||| (1u <<< (i &&& 0x1F))
let inline cullp p s low = let rec cull' i = if i < BUFSZ then cullbit i; cull' (i + int p)
cull' (if s >= low then int((s - low) >>> 1)
else let r = ((low - s) >>> 1) % p in if r = 0u then 0 else int(p - r))
let inline cullpg low = //cull composites from whole buffer page for efficiency
let max = low + BUFRNG - 1u in let max = if max < low then uint32(-1) else max
let sqrtlm = uint32(sqrt(float max)) in let sqrtlmndx = int((sqrtlm - 3u) >>> 1)
if low <= 3u then for i = 0 to sqrtlmndx do if testbit i then let p = uint32(i + i + 3) in cullp p (p * p) 3u
else baseprimes |> Seq.skipWhile (fun p -> //force side effect of culling to limit of buffer
let s = p * p in if p > 0xFFFFu || s > max then false else cullp p s low; true) |> Seq.nth 0 |> ignore
let rec mkpi i low =
if i >= BUFSZ then let nlow = low + BUFRNG in Array.fill buf 0 buf.Length 0u; cullpg nlow; mkpi 0 nlow
else (if testbit i then i,low else mkpi (i + 1) low)
cullpg 3u; Seq.unfold (fun (i,lw) -> //force cull the first buffer page then doit
let ni,nlw = mkpi i lw in let p = nlw + (uint32 ni <<< 1)
if p < lw then None else Some(p,(ni+1,nlw))) (0,3u)
and baseprimes = oddprimes() |> Seq.cache
seq { yield 2u; yield! oddprimes() }
primesAPF32() |> Seq.nth 203280220 |> printfn "%A"
nouveau code calcule les 203,280,221 nombres premiers dans la plage de numéros 32 bits dans environ AJOUTÉE / CORRIGÉ: 25,4 secondes avec des temps en cours d'exécution pour la première 100000, un million, 10 millions, et 100 millions testé 0,01, 0,088, 0,94 et 11,25 secondes respectivement sur un ordinateur de bureau rapide (3,5 GHz i7-2700K @) et peut fonctionner sur et href="http://ideone.com/l8zupF" rel="nofollow"> ideone.com , bien que dans une moindre plage pour ceux-ci en raison des contraintes de temps d'exécution. Il a une performance pire que le code de Jon Harrop pour les petites gammes de quelques milliers de nombres premiers en raison de sa complexité accrue de calcul mais très passe rapidement pour des plages plus grandes en raison de son meilleur algorithme de performance qui compense cette complexité telle qu'il est environ cinq fois plus rapide pour le premier 10 millionième et environ sept fois plus rapide juste avant le code de Jon explose vers le premier 25 millionième.
le temps d'exécution total, plus de la moitié est consacrée à l'énumération de séquence de base et ne serait donc pas aidé dans une large mesure en exécutant des opérations d'abattage Décomposer opérations d'arrière-plan, bien que les optimisations de factorisation des roues en combinaison contribueraient ( bien que plus de calculs, que la complexité serait en arrière-plan) en ce sens qu'ils réduisent le nombre d'opérations de test tampons nécessaires dans l'énumération. D'autres optimisations pourraient être réalisées si l'ordre des séquences n'a pas besoin d'être conservés (comme ne comptant que le nombre de nombres premiers ou en additionnant les nombres premiers) que les séquences peuvent être écrites pour soutenir les interfaces d'énumération parallèle ou le code pourrait être écrit en tant que classe afin que les méthodes membres pourraient faire le calcul sans l'énumération. Ce code pourrait facilement être accordé à l'offre à proximité du même genre de performances que le code C #, mais de façon plus concise exprimé, bien qu'il ne sera jamais atteindre les performances du code natif optimisé de C comme
Je ne suis pas très familier avec multimaps Haskell, mais le F # Power Pack a une classe HashMultiMap, dont résumé XMLDoc est: « tables de hachage, par défaut basée sur F # structurelle « hachage » fonctions et (=) Le tableau peut mapper une seule touche à plusieurs fixations. ». Peut-être que cela pourrait vous aider?
2 * 10 ^ 6 en 1 seconde sur Corei5
let n = 2 * (pown 10 6)
let sieve = Array.append [|0;0|] [|2..n|]
let rec filterPrime p =
seq {for mul in (p*2)..p..n do
yield mul}
|> Seq.iter (fun mul -> sieve.[mul] <- 0)
let nextPrime =
seq {
for i in p+1..n do
if sieve.[i] <> 0 then
yield sieve.[i]
}
|> Seq.tryHead
match nextPrime with
| None -> ()
| Some np -> filterPrime np
filterPrime 2
let primes = sieve |> Seq.filter (fun x -> x <> 0)
