Le tamis D'Eratosthène en F#

bien qu'il y ait eu une réponse donnant un algorithme utilisant une file D'attente prioritaire (PQ) comme dans un SkewBinomialHeap , ce N'est peut-être pas le bon PQ pour le travail. Ce que le tamis incrémentiel D'Eratosthène (iEoS) exige est un PQ qui a d'excellentes performances pour obtenir la valeur minimale et réinsérer des valeurs la plupart du temps un peu plus loin dans la file d'attente, mais n'a pas besoin de l'ultime dans la performance pour ajouter de nouveaux valeurs comme iSoE ajoute seulement comme nouvelles valeurs un total des nombres premiers jusqu'à la racine carrée de la gamme (qui est une infime fraction du nombre de Ré-insertions qui se produisent une fois par réduction). Le SkewBinomialHeap PQ ne donne pas vraiment beaucoup plus que l'utilisation de la carte intégrée qui utilise un arbre de recherche binaire équilibré - toutes les opérations O(log n) - autre que cela change la pondération des opérations légèrement en faveur des exigences de la SoE. Cependant, le SkewBinaryHeap nécessite encore de nombreuses opérations O (log n) par réduction.

UN PQ mis en œuvre comme un en Tas", 151960920" en plus particulier comme un Tas Binaire et encore plus particulièrement en tant que MinHeap assez bien satisfait à l'iSoE des exigences actuelles en O(1) rendement à obtenir le minimum et le O(log n) de la performance pour la ré-insertion et d'ajouter de nouvelles entrées, bien que la performance est en fait une fraction de O(log n), puisque la plupart de la re-insertions se produire à proximité de la haut de la file d'attente et la plupart des ajouts de nouvelles valeurs (qui n'ont pas d'importance car ils sont peu fréquents) se produisent près de la fin de la file d'attente où ces opérations sont les plus efficaces. De plus, le MinHeap PQ peut mettre en œuvre efficacement la fonction delete minimum et insert en un seul passage (généralement une fraction de) un o(log n). Ensuite, plutôt que pour la carte (qui est implémentée comme un "arbre AVL ) où il y a une opération O (log n) Avec généralement une gamme complète "log n" due à la valeur minimum que nous exigeons étant à l'extrémité gauche de l'arbre, nous ajoutons et enlevons généralement le minimum à la racine et insérons en moyenne quelques niveaux vers le bas dans un passage. Ainsi, le PQ de MinHeap peut être utilisé avec seulement une fraction D'opération O(log n) par réduction d'abattage plutôt qu'avec plusieurs opérations de plus grande fraction d'opération O(log n).

le PQ MinHeap peut être implémenté avec du code fonctionnel pur (sans" removeMin " implémenté comme l'iSoE ne le fait pas)

[<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, l'iSoE peut être écrit avec les optimisations de factorisation de roue et en utilisant des flux Co-inductifs efficaces (CIS's) comme suit:

type CIS<'T> = class val v:'T val cont:unit->CIS<'T> new(v,cont) = { v=v;cont=cont } 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.pi,cullstate(kv.v.p,whladv kv.v.pi)
                                                   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 en environ 0.077 secondes, les 1.000.000 premiers nombres premiers en 0.977 secondes, les 10.000.000 premiers nombres premiers en environ 14.33 secondes, et les 100 millions de premiers nombres premiers en environ 221,87 secondes, tous sur un i7-2700K (3,5 GHz) comme code 64 bits. Ce code purement fonctionnel est légèrement plus rapide que celui de Dustin Cambell mutable Dictionary basé code avec les optimisations communes ajoutées de factorisation de la roue, l'ajout différé des amorces de base, et l'utilisation de la plus efficace CID tous ajoutés ( tryfsharp et ideone ) mais est encore le code fonctionnel pur où son utilisation de la classe du dictionnaire n'est pas . Cependant, pour des gammes de prime plus grandes d'environ deux milliards (environ 100 millions de nombres premiers), le code utilisant le dictionnaire de table de hachage sera plus rapide que les opérations de dictionnaire n'ont pas un facteur O(log n) et ce gain surmonte la complexité de calcul de L'utilisation des tables de hachage de dictionnaire.

le programme ci-dessus a la caractéristique supplémentaire que la roue de factorisation est paramétrée de sorte que, par exemple, on peut utiliser une roue extrêmement grande en réglant WHLPRMS à [|2u;3u;5u;7u;11u;13u;17u;19u/] et FSTPRM à 23u pour obtenir un temps d'exécution d'environ deux tiers pour les grandes portées à 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 ne commence à courir, ce qui est une constante au-dessus, peu importe la gamme première.

analyse Comparative : par rapport à l'arbre purement fonctionnel incrémentiel mise en œuvre de pliage, cet algorithme est juste un peu plus rapide parce que la hauteur moyenne utilisée de L'arbre de Menheap est moins d'un facteur de deux que la profondeur de l'arbre plié, mais qui est compensée par un facteur constant équivalent perte d'efficacité dans la capacité de traverser les niveaux d'arbre PQ en raison d'elle étant basée sur un tas binaire exigeant le traitement des feuilles droite et gauche pour chaque niveau de tas et une branche dans un sens ou dans l'autre plutôt qu'une seule comparaison par niveau pour le pliage de l'arbre avec généralement la branche moins profonde la prise. Par rapport à d'autres algorithmes fonctionnels basés sur PQ et Map, les améliorations sont généralement par un facteur constant dans la réduction du nombre D'opérations O(log n) dans la traversée de chaque niveau des structures d'arbre respectives.

La MinHeap est généralement , mis en œuvre en tant que mutable tableau tas binaire après l'arbre généalogique de base du modèle inventé par Michael Eytzinger il y a plus de 400 ans. Je sais que la question disait qu'il n'y avait aucun intérêt pour le code mutable non fonctionnel, mais si l'on doit éviter tous les sous-codes qui utilisent la mutabilité, alors nous ne pourrions pas utiliser ceux de list ou de LazyList qui utilisent la mutabilité "sous les couvertures" pour des raisons de performance. Imaginez donc que la version mutable alternative suivante du PQ de MinHeap soit fournie par une bibliothèque et bénéficie d'un autre facteur de plus de deux pour les plus grandes gammes de premier ordre dans la 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

Geek note: je m'attendais en fait à ce que la version mutable offre un bien meilleur rapport de performance, mais elle s'embourbe dans les ré-insertions en raison de la structure de code imbriquée si-then-else et du comportement aléatoire des valeurs premières de cull, ce qui signifie que la prédiction de branche CPU échoue pour une grande proportion des branches résultant en de nombreux 10's supplémentaires de cycles D'horloge CPU par réduction de cull pour reconstruire le cache d'instruction pre-fetch.

le seul autre les gains de performance à facteur constant sur cet algorithme seraient la segmentation et l'utilisation de multitâches pour un gain de performance proportionnel au nombre de cœurs CPU; cependant, en l'état actuel, il s'agit de l'algorithme fonctionnel SOE le plus rapide à ce jour, et même la forme fonctionnelle pure en utilisant le MinHeap fonctionnel Bat mises en œuvre impératives simplistes telles que code de Jon Harrop ou Sieve de Johan Kullbom d'Atkin (qui est une erreur dans son timing comme il a seulement calculé les primes à 10 millions plutôt que le 10 millionième prime ), mais ces algorithmes seraient environ cinq fois plus rapide si de meilleures optimisations ont été utilisées. Ce rapport d'environ cinq entre le code fonctionnel et le code impératif sera quelque peu réduit lorsque nous ajouterons le multi-filetage de la plus grande factorisation de roue que la complexité computationnelle du code impératif augmente plus rapidement que le code fonctionnel et multi-filetage aide le code fonctionnel plus lent plus que le code impératif plus rapide car ce dernier se rapproche de la limite de base du temps requis pour énumérer à travers les nombres premiers trouvés.

EDIT_ADD: même si l'on pouvait choisir de continuer à utiliser la version fonctionnelle pure de MinHeap, en ajoutant efficace segmentation en préparation pour multi-filetage serait légèrement "briser" la "pureté" du code fonctionnel comme suit: 1) la manière la plus efficace de transférer une représentation de nombres premiers composés est un tableau à emboîtement de la taille du segment, 2) alors que la taille du tableau est connue, utiliser une compréhension de tableau pour l'initialiser d'une manière fonctionnelle n'est pas efficace car il utilise "ResizeArray" sous les couvertures qui doit se Copier pour chaque x additions (je pense que 'x' est huit pour l'implémentation actuelle) et en utilisant le tableau.init ne fonctionne pas car de nombreuses valeurs à des index particuliers sont sautées, 3) par conséquent, la façon la plus facile de remplir le le tableau composite sélectionné doit le zérocréater de la taille correcte et ensuite exécuter une fonction d'initialisation qui pourrait écrire à chaque indice de tableau mutable pas plus d'une fois. Bien que ce ne soit pas strictement "fonctionnel", il est proche en ce que le tableau est initialisé et puis jamais modifié à nouveau.

le code avec segmentation ajoutée, multi-threading, circonférence factorielle de roue programmable, et de nombreux code supplémentaire accordé pour mettre en œuvre la segmentation et multi-threading est la moitié inférieure environ 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.pi,cullstate(kv.v.p,whladv kv.v.pi)
                        if n>kv.k then mkprm (n,wi,(MinHeap.reinsertMinAs ca cs pq),bps,q,lstp,bgap)
                        elif n=kv.k then mkprm (nxt,nxti,(MinHeap.reinsertMinAs ca cs pq),bps,q,lstp,(bgap+1us))
                        else (uint16(n-uint32 lstp),bgap,wi,(nxt,nxti,pq,bps,q,n,1us))
  let rec pCIS p pg bg pi pq bps q = prmsCIS(pg,bg,pi,fun()->
    let (npg,nbg,npi,(nxt,nxti,npq,nbps,nq,nl,ng))=mkprm (p+uint32 WHLPTRN.[pi],whladv pi,pq,bps,q,p,0us)
    pCIS (p+uint32 npg) npg nbg npi npq nbps nq)
  let rec baseprimes() = prmsCIS(uint16 FSTPRM,0us,0,fun()->
                           let np,npi=advcnd FSTPRM 0,whladv 0
                           pCIS np (uint16 WHLPTRN.[0]) 1us npi MinHeap.empty (baseprimes()) (FSTPRM*FSTPRM))
  let prmspg nxt adj pq bp q =
    //compute next buffer size rounded up to next even wheel circle so at least one each base prime hits the page
    let rng = max (((uint32(MXSTP+uint64(sqrt (float (MXSTP*(MXSTP+4UL*nxt))))+1UL)>>>1)+MINRNGSTP)/MINRNGSTP*MINRNGSTP) MINBFRNG
    let nxtp() = async {
      let rec addprms pqx (bpx:prmsCIS) qx = 
        if qx>=adj then pqx,bpx,qx //add primes to queue for new lower limit
        else let p = (uint32 bpx.bg<<<16)+uint32 bpx.pg in let nbps = bpx.cont()
             let pg = uint32 nbps.pg in let np = p+pg in let sqr = qx+pg*((p<<<1)+pg)
             let nbps = prmsCIS(uint16 np,uint16(np>>>16),nbps.pi,nbps.cont)
             addprms (MinHeap.insert qx (cullstate(p,bpx.pi)) pqx) nbps sqr
      let adjcinpg low k (v:cullstate) = //adjust the cull states for the new page low value
        let p = v.p in let WHLSPN = int64 WHLCRC*int64 p in let db = int64 p*int64 WHLPOS.[v.pi]
        let db = if k<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 fonctionnels et basés sur des tableaux, ont une fonction" adjust " ajoutée pour permettre le réglage de l'état de cull de la version du PQ de chaque thread au début de chaque nouvelle page de segment. Notez également qu'il était possible de modifier le code afin que le calcul est fait en utilisant des gammes de 32 bits avec la séquence finale de production comme uint64 à peu de coût en temps de calcul de sorte que actuellement la gamme théorique est quelque chose plus de 100 trillions (dix soulevées à la puissance de quatorze) si l'on était prêt à attendre les environ trois à quatre mois nécessaires pour calculer cette gamme. Les vérifications numériques des fourchettes ont été supprimées, car il est peu probable que quelqu'un utilise cet algorithme pour calculer jusqu'à cette fourchette, et encore moins au-delà.

utilisant la fonction pure MinHeap et 2,3,5 et 7 de la roue de la factorisation, le programme ci-dessus calcule les premiers cent mille, un million, dix millions de dollars, et une centaine de millions de nombres premiers dans 0.062, 0.629, 10.53, et 195.62 secondes, respectivement. En utilisant le MinHeap basé sur un réseau, on accélère ceci jusqu'à 0.097, 0.276, 3.48 et 51.60 secondes, respectivement. En utilisant la roue 2,3,5,7,11,13,17 en changeant WHLPRMS à "[|2u;3u;5u;7u;11u;13u;17u|] " et fstprm à 19u vitesses qui montent encore un peu plus à 0,181, 0,308, 2,49, et 36,58 secondes, (pour une amélioration constante du facteur avec un dépassement constant). Ce tweak le plus rapide calcule les 203.280.221 nombres premiers dans la gamme de 32 bits en environ 88.37 secondes. La constante" BFSZ "peut être ajustée avec des compromis entre des temps plus lents pour des gammes plus petites version des temps plus rapides pour des gammes plus grandes, avec une valeur de" 1<<14 " recommandé d'être essayé pour les gammes plus grandes. Cette constante ne définit que la taille minimale du tampon, le programme ajustant la taille du tampon au-dessus de cette taille automatiquement pour des portées plus grandes de sorte que le tampon est suffisant de sorte que la prime de base la plus grande requise pour la gamme de page sera toujours "rayer" chaque page au moins une fois; cela signifie que la complexité et les frais généraux d'un "tamis à seau" supplémentaire ne sont pas nécessaires. Cette dernière version entièrement optimisée peut calculer les primes jusqu'à 10 et 100 milliards en environ 256,8 et 3617,4 secondes (un peu plus d'une heure pour les 100 milliards) comme testé en utilisant "primesPQOWSE() |> Seq.takeWhile (( > =) 1000000000UL) / > Seq.fold (fun S p - > s + 1UL) 0UL" for output. C'est de là que proviennent les estimations d'environ une demi-journée pour le compte de nombres premiers à un trillion, une semaine pour un maximum de dix trillions et environ trois à quatre mois pour un maximum de cent trillions.

Je ne pense pas qu'il soit possible de faire du code fonctionnel ou presque fonctionnel en utilisant l'algorithme incrémental SoE pour courir 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é grandement à la complexité du code, de sorte qu'il est probablement un peu plus complexe que le code optimisé de façon équivalente basé sur l'abattage en réseau rectiligne avec ce code capable de courir environ dix fois plus vite que cela et sans l'exposant supplémentaire dans la performance, ce qui signifie que ce code incrémental fonctionnel a un pourcentage supplémentaire sans cesse croissant de frais généraux.

est-ce donc utile autrement que d'un point de vue théorique et intellectuel intéressant? Probablement, il ne l'est pas. Pour de plus petites gammes de nombres premiers jusqu'à environ dix millions, les meilleurs de l'incrémental functional SOE fonctionnelle de base qui n'est pas entièrement optimisée sont probablement adéquates et assez simples à écrire ou ont moins D'utilisation de mémoire vive que le plus simple impératif SoE. Cependant, ils sont beaucoup plus lents que le code plus impératif en utilisant un tableau de sorte qu'ils "s'épuisent de vapeur" pour les gammes au-dessus de cela. Alors qu'il a été démontré ici que le code peut être accéléré par l'optimisation, c'est encore 10 fois plus lent qu'un plus impératif pur la version basée sur des tableaux a cependant ajouté à la complexité d'être au moins aussi complexe que ce code avec des optimisations équivalentes, et même ce code sous F# sur DotNet est environ quatre fois plus lent que l'utilisation d'un langage tel que C++ compilé directement au code natif; si on voulait vraiment étudier de grandes gammes de nombres premiers, on utiliserait probablement l'un de ces autres langages et techniques où primessive peut calculer le nombre de nombres premiers dans la gamme de centaines de trillions quatre heures au lieu des trois mois requis pour ce code. END_EDIT_ADD

voici ma tentative d'une traduction raisonnablement fidèle du code de Haskell en 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)

Premier tamis mis en œuvre avec la boîte aux lettres de processeurs:

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 un tamis cartographique incrémental (et récursif) d'Eratosthènes utilisant des séquences car il n'y a pas besoin de mémoïser les valeurs de séquence précédentes (sauf qu'il y a un léger avantage à cacher les valeurs de base des nombres premiers en utilisant Seq.cache), avec les principales optimisations étant qu'il utilise la factorisation de la roue pour la séquence d'entrée et qu'il utilise plusieurs flux (récursifs) pour maintenir les nombres premiers de base qui sont inférieurs à racine carrée du dernier numéro étant tamisé, comme suit:

  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 les 100.000 E primes jusqu'à 1.299.721 dans environ une 0.445 seconde, mais n'étant pas un algorithme EOS impératif approprié il ne se met pas à l'échelle près linéairement avec un nombre accru de primes, prend 7.775 secondes pour trouver le 1.000.000 prime jusqu'à 15.485.867 pour une performance sur cette gamme d'environ O(N^1.2) où n est le prime maximum trouvé.

il y a un peu plus d'accordage que cela pourrait être fait, mais cela ne fera probablement pas une grande différence quant à un pourcentage élevé de meilleurs résultats comme suit:

1. comme la bibliothèque de séquence F# est très lente, on pourrait utiliser un type défini qui met en œuvre IEnumerable pour réduire le temps passé dans la séquence interne, mais comme les opérations de séquence ne prennent qu'environ 20% du temps total, même si ceux-ci étaient réduits à temps zéro, le résultat ne serait qu'une réduction à 80% du temps.

2. D'autres formes de stockage de cartes pourraient être essayées telles qu'une file d'attente prioritaire comme mentionné par O'Neil ou le SkewBinomialHeap comme utilisé par @gradbot, mais au moins pour le SkewBinomialHeap, l'amélioration de la performance est seulement quelques pour cent. Il semble que dans le choix de différentes implémentations de carte, on échange Juste une meilleure réponse dans la recherche et la suppression des éléments qui sont près du début de la liste contre le temps passé à insérer de nouveaux entrées afin de profiter de ces avantages de sorte que le gain net est à peu près un lavage et a toujours une performance O(log n) avec des entrées croissantes dans la carte. L'optimisation ci-dessus en utilisant plusieurs flux d'entrées juste à la racine carrée réduire le nombre d'entrées dans la carte et donc faire ces améliorations de peu d'importance.

EDIT_ADD: j'ai fait le peu supplémentaire d'optimisation et la performance a amélioré un peu plus que prévu, probablement en raison de la meilleure façon d'éliminer les Seq.sauter comme un moyen d'avancer à travers la séquence des nombres premiers de base. Cette optimisation utilise un remplacement pour la génération de la séquence interne comme un tuple de valeur entière et une fonction de continuation utilisée pour avancer à la valeur suivante dans la séquence, avec la séquence finale F# générée par une opération de déploiement global. Le Code est le suivant:

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 nécessaires pour trouver le 100,000 ème et 1.000.000 TH premiers sont d'environ 0,31 et 5,1 secondes, respectivement, il ya donc 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 des séquences, et bien qu'elles soient plus rapides que les versions utilisées par le module Seq, elles ne font guère plus de différence à cet algorithme où la plupart du temps est passé dans les fonctions de mappage. END_EDIT_ADD

autres par rapport aux implémentations EOS incrémentielles basées sur la carte, il y a une autre implémentation "purement fonctionnelle" utilisant le pliage D'arbre qui est dit être légèrement plus rapide, mais comme il a encore un terme O(log n) dans le pliage d'arbre, je soupçonne qu'il sera principalement plus rapide (si c'est le cas) en raison de la façon dont l'algorithme est mis en œuvre quant au nombre d'opérations informatiques par rapport à l'utilisation d'une carte. Si les gens sont intéressés, je développerai aussi cette version.

en fin de compte, il faut accepter qu'aucun pur la mise en œuvre fonctionnelle de L'EoS incrémental sera jamais proche de la vitesse de traitement brute d'une bonne mise en œuvre impérative pour des gammes numériques plus grandes. Cependant, on pourrait proposer une approche où tout le code est purement fonctionnel, sauf pour le tamisage segmenté des nombres composites sur une plage en utilisant un tableau (mutable) qui serait proche de la performance O(n) et en pratique serait de cinquante à cent fois plus rapide que les algorithmes fonctionnels pour de grandes gammes telles que les premiers 200 millions de nombres premiers. Cela a été fait par @Jon Harrop dans son blog , mais cela pourrait être accordé plus loin avec très peu de code supplémentaire.

Voici mes deux cents, bien que je ne suis pas sûr qu'il répond au critère de L'OP pour être vraiment le tamis d'Eratosthène. Il n'utilise pas la division modulaire et met en œuvre une optimisation à partir du papier cité par L'OP. Il ne fonctionne que pour les listes finies, mais cela me semble être dans l'esprit de la façon dont le tamis a été décrit à l'origine. De plus, le document traite de la complexité en termes de nombre de marques et de divisions. Il semble que, comme nous devons traverser un liked list, que cela peut-être en ignorant certains aspects clés des divers algorithmes en termes de performance. En général, bien que la 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 message original. J'ai essayé de suivre la logique originale du tamis avec quelques modifications. À savoir commencer par le premier élément et rayez les multiples de cet élément de l'ensemble. Cet algorithme cherche littéralement l'article suivant qui est un multiple du premier au lieu de faire modulaire division de chaque nombre dans l'ensemble. Une optimisation du papier est qu'il commence à chercher des multiples du prime supérieur à p^2.

la partie de la fonction helper à plusieurs niveaux traite de la possibilité que le prochain multiple du premier puisse déjà être retiré de la liste. Donc par exemple avec le prime 5, il va essayer de supprimer le nombre 30, mais il ne le trouvera jamais parce qu'il a été déjà supprimé par le prime 3. Espérons que précise l'algorithme de logique.

je sais que vous avez explicitement déclaré que vous étiez intéressé par une implémentation de tamis purement fonctionnelle, donc j'ai attendu jusqu'à présent pour présenter mon tamis. Mais en relisant l'article que vous avez mentionné, je vois que l'algorithme de crible incrémentiel présenté est essentiellement le même que le mien, la seule différence étant les détails de mise en œuvre de l'utilisation de techniques purement fonctionnelles par rapport à des techniques décidément impératives. Donc je pense que je suis au moins à moitié qualifié pour satisfaire ta curiosité. En Outre, Je serait d'avis que l'utilisation de techniques impératives lorsque des gains de performance significatifs 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 Haskell pure. J'ai d'abord publié cette mise en œuvre sur mon projet Euler pour F#un blog mais re-publier ici avec le code pré-requis substitué retour dans et structural typing supprimé. primes peut calculer les 100 000 premiers primes en 0,248 secondes et les 1,000,000 premiers primes en 4,8 secondes sur mon ordinateur (notez que primes cache ses résultats de sorte que vous aurez besoin de le réévaluer chaque fois que vous effectuez un benchmark).

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 pour accomplir le tamis incrémentiel D'Eratosthènes (SoE) en utilisant seulement le code F fonctionnel pur. Il est adapté du code Haskell développé comme "cette idée est due à Dave Bayer, bien qu'il ait utilisé une formulation complexe et la structure équilibrée de l'arbre ternaire, s'approfondissant progressivement de manière uniforme (formulation simplifiée et une asymétrique, l'approfondissement à la bonne structure de l'arbre binaire introduit par Heinrich Apfelmus, encore simplifié par Will Ness). Mise en scène idée de production due à M. O'Neill" selon le lien suivant: code de pliage optimisé de L'arbre à l'aide d'une roue factorielle en Haskell .

le code suivant a plusieurs optimisations qui le rendent plus approprié pour l'exécution en F#, comme suit:

1. le code utilise des flux coinductifs au lieu de Lazylist's comme cet algorithme n'a pas (ou peu) besoin de la mémoization de LazyList et mes flux coinductifs sont plus efficace que les Lazylistes (du FSharp.PowerPack) ou les séquences intégrées. Un autre avantage est que mon code peut être exécuté sur tryFSharp.org et ideone.com sans avoir à copier-coller dans le Microsoft.FSharp.PowerPack Base de code source pour le LazyList type et le module (avec l'avis de droit d'auteur)

2. il a été découvert qu'Il ya beaucoup de frais généraux pour le modèle de F# l'appariement sur les paramètres de fonction de sorte que le précédent plus lisible type d'union discriminé en utilisant tuples a été sacrifié en faveur des types de structure de sous-valeur (ou classe comme des courses plus rapides sur certaines plates-formes) pour environ un facteur de deux ou plus d'accélération.

3. les optimisations de Will Ness allant du repliement linéaire de l'arbre au repliement bilatéral en passant par le repliement multidirectionnel et les améliorations utilisant la factorisation de la roue sont à peu près aussi efficaces rationométriquement pour F# qu'elles le sont pour Haskell, avec la principale différence entre les deux langues étant que Haskell peut être compilé en code natif et a un compilateur plus fortement optimisé tandis que F# a plus de fonctionnement supérieur sous le système de cadre DotNet.

   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é car le code original ne gérait pas correctement la queue du flux et passait la queue du paramètre stream à la fonction de paires joinT3 fonctionne plutôt que la queue suivant le flux zs. Le timing ci-dessous a également été corrigé en conséquence, avec environ 30% d'accélération supplémentaire. Les codes de liaison tryFSharp et ideone ont également été corrigés. END_EDIT_ADD

le programme ci-dessus fonctionne à environ O (N^1.1) la performance avec n la prime maximale calculée ou environ O (N^1.18) quand n est le nombre de primes calculé, et prend environ 2.16 secondes pour calculer le premier million des nombres premiers (environ 0,14 seconde pour les 100 000 premiers nombres premiers) sur un ordinateur rapide exécutant du code 64 bits en utilisant des types de struct plutôt que des classes (il semble que certaines implémentations box et unbox les structures de by-value dans la formation de fermetures). Je considère que c'est à propos de la portée maximale pratique pour n'importe lequel de ces algorithmes purs premiers fonctionnels. Il s'agit probablement de la plus rapide que l'on puisse exécuter un algorithme pure fonctionnelle SoE autre que pour quelques retouches mineures pour réduire les facteurs constants.

mis à part la combinaison de la segmentation et du multi-threading pour partager le calcul entre plusieurs cœurs CPU, la plupart des "retouches" qui pourraient être apportées à cet algorithme sont dans l'augmentation de la circonférence de la factorisation de la roue pour un gain de performance jusqu'à environ 40% et des gains mineurs dus aux retouches quant à l'utilisation de structures, classes, tuples, ou des paramètres individuels plus directs dans le passage de l'état entre les fonctions.

EDIT_ADD2: j'ai fait les optimisations ci-dessus avec le résultat que le code est maintenant presque deux fois plus rapide en raison des optimisations de la structure avec le bonus ajouté d'utiliser optionnellement des circonférences de factorisation de roues plus grandes pour la réduction plus petite ajoutée. Notez que le code ci-dessous évite d'utiliser des continuations dans la boucle de génération de séquence principale et ne les utilise que si nécessaire pour les flux de base des nombres premiers et les flux de cull subséquents composés dérivés de ces nombres premiers de base. Nouvelle le 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 les premiers cent mille, million, et dix millions de nombres premiers, respectivement, tous sur la machine de référence Intel i7-2700K (3.5 GHz) en mode 64 bits. Ce n'est pas beaucoup plus lent que l'implémentation Haskell de référence à partir de laquelle ce code a été dérivé, bien qu'il soit légèrement plus lent sur tryfsharp et ideone en raison d'être en mode 32 bits pour tryfsharp sous Silverlight (environ la moitié de nouveau plus lent) et courant sur une machine plus lente sous Mono 2.0 (qui est intrinsèquement beaucoup plus lent pour F#) sur ideone est donc jusqu'à cinq fois plus lent que la machine de référence. Notez que le temps d'exécution rapporté par ideone inclut le temps d'initialisation pour les tableaux de tables de recherche intégrées, dont le temps doit être pris en compte.

le programme ci-dessus a la caractéristique supplémentaire que la roue de factorisation est paramétré de façon à ce que, par exemple, on puisse utiliser une roue extrêmement grande en réglant WHLPRMS à [| 2u;3u;5u;7u;11u;13u;17u;19u |] et FSTPRM à 23u pour obtenir un temps d'exécution d'environ deux tiers pour les grandes portées à 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 ne commence à fonctionner.

Geek note: je n'ai pas mis en œuvre un "non-sharing fixpoint combinator for telescoping multi-étages production de primes" comme selon le code de référence 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 fonctionne sans s'enfuir dans une boucle infinie et débordement de la pile. Bien que mes flux Co-inductifs (CIS) aient certaines propriétés de paresse, ils ne sont pas formellement des listes paresseuses ou des séquences mises en cache (elles deviennent des séquences non-mises en cache et pourraient être mises en cache lorsqu'elles sont passées ainsi une fonction telle que celle de "genseq" que je fournis pour la séquence de sortie finale). Je ne voulais pas utilisez L'implémentation PowerPack de LazyList car elle n'est pas très efficace et exigerait que je copie ce code source dans tryfsharp et ideone, qui ne fournissent pas de modules importés. L'utilisation des séquences intégrées (même mises en cache) est très inefficace lorsque l'on veut utiliser des opérations tête/queue comme celles requises pour cet algorithme car la seule façon d'obtenir la queue d'une séquence est d'utiliser "Seq.skip 1" qui, lors d'utilisations multiples, produit une nouvelle séquence basée sur la séquence originale de façon récursive sauté de nombreuses fois. Je pourrais mettre en œuvre ma propre classe LazyList efficace basée sur les CIS, mais il semble à peine la peine de démontrer un point où les objets récursifs "initcmpsts" et "baseprimes" prennent très peu de code. De plus, le passage D'un Lazyliste à une fonction pour produire des extensions à ce Lazyliste qui n'utilise des valeurs que depuis le début du Lazyliste exige que presque tout le Lazyliste est memoized pour une réduction de l'efficacité de la mémoire: une passe pour les 10 premiers des millions de nombres premiers nécessiteront un Lazyliste en mémoire avec près de 180 millions d'éléments. J'ai donc pris un pass sur ce.

noter que pour des gammes plus grandes (10 millions de nombres premiers ou plus), ce code purement fonctionnel est à peu près la même vitesse que de nombreuses implémentations impératives simplistes du tamis D'Eratosthène ou D'Atkins, bien que cela soit dû au manque d'optimisation de ces algorithmes impératifs; une implémentation plus impérative que cela en utilisant des optimisations équivalentes et les réseaux de tamisage segmentés seront encore environ dix fois plus rapides que ceci selon ma réponse "presque fonctionnelle".

notent également que s'il est possible de mettre en œuvre le tamisage segmenté en utilisant le repliement de l'arbre, il est plus difficile puisque les algorithmes d'avance de l'abattage sont enfouis à l'intérieur des continuations utilisées pour l'opérateur "union - ^" et que travailler autour de cela signifierait qu'une gamme d'abattage en progression continue doit être utilisée; ceci est à la différence d'autres algorithmes où l'état de la variable de cull peut être réinitialisée pour chaque nouvelle page, y compris avoir leur portée réduite, de sorte que si des portées plus grandes que 32 bits sont utilisés, la portée interne de culled peut encore être réinitialisée pour fonctionner dans la gamme de 32 bits, même si une gamme de 64 bits de primes est déterminée à peu de coût dans le temps d'exécution par prime. END_EDIT_ADD2

pour ce que sa valeur, ce n'est pas un tamis D'Erathothènes, mais son très fast:

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 100000 ème prime en 1.25 secondes sur ma machine (AMD Phenom II, 3.2 GHZ quadcore).

en fait j'ai essayé de faire la même chose, j'ai d'abord essayé la même implémentation naïve qu'en question, mais c'était trop lent. J'ai alors trouvé cette page YAPES: Problem Seven, Part 2 , où il a utilisé le tamis réel D'Eratosthène basé sur Melissa E. O'Neill. J'ai pris le code de là, juste un peu modifié, parce que F # a changé un peu 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

comme cette question demande spécifiquement d'autres algorithmes, je fournis l'implémentation suivante:

ou connaît peut-être d'autres méthodes de mise en œuvre ou des algorithmes de tamisage

aucune soumission de tamis divers d'algorithmes D'Eratosthène (SoE) n'est vraiment complète sans mentionner le tamis D'Atkin (SoA), qui est en fait une variation de SoE en utilisant les solutions à un ensemble de équations quadratiques pour mettre en œuvre l'abattage composite ainsi que l'élimination de tous les multiples des carrés des nombres premiers de base (nombres premiers inférieurs ou égaux à la racine carrée du nombre le plus élevé testé pour la primalité). En théorie, la SoA est plus efficace que la SoE en ce qu'il y a un peu moins d'opérations dans la gamme, de sorte qu'elle aurait dû avoir environ 20% moins de complexité pour la gamme d'environ 10 à 100 millions, mais en pratique, elle est généralement plus lente en raison des frais généraux de la complexité de résoudre plusieurs équations quadratiques. Bien que L'implémentation C de Daniel J. Bernstein soit plus rapide que L'implémentation SoE avec laquelle il l'a testé pour cette gamme particulière de numéros de test , l'implémentation SoE avec laquelle il a testé n'était pas la plus optimale et les versions plus optimisées de SOE sont encore plus rapides. Cela semble être le cas ici, même si je reconnais qu'il peut y être d'autres optimisations que j'ai manqué.

comme O'Neill, dans son article sur la SoE, utilise des tamis incrémentiels non limités pour montrer principalement que le tamis Turner n'est pas une SoE tant au niveau de l'algorithme qu'au niveau de la performance, elle n'a pas tenu compte de nombreuses autres variations de la SoE telles que la SoA. En faisant une recherche rapide de la littérature, Je ne trouve aucune application de SoA aux séquences incrémentielles illimitées dont nous discutons ici, donc je l'ai adapté moi-même comme dans le code suivant.

tout comme le cas PUR sans limite de SoE peut être considéré comme ayant comme des nombres composés une séquence sans limite de séquences de nombres premiers multiples, le SoA considère comme des nombres premiers potentiels la séquence sans limite des séquences sans limite de toutes les expressions des équations quadratiques avec une des deux variables libres, 'x' ou 'y' fixée à une valeur de départ et avec une séquence "élimination" séparée des séquences de tous les multiples des nombres premiers de base, qui la dernière est très similaire aux séquences d'élimination composées de séquences pour SoE sauf que les séquences progressent plus rapidement par le carré des nombres premiers plutôt que par un (plus petit) multiple des nombres premiers. J'ai essayé de réduire le nombre d'équations quadratiques séquences exprimées en reconnaissant que, pour les fins d'un différentiel, tamis, la "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 étranges, ainsi que l'application 2357 factorisation de roue (bien que L'aos a déjà implicite 235 factorisation de roue). Il utilise l'algorithme efficace de fusion/combinaison de pliage d'arbre comme dans la fusion d'arbre SoE pour traiter chaque séquence de séquences, mais avec une simplification que l'opérateur de l'union ne combine pas dans la fusion car l'algorithme SoA dépend de pouvoir basculer l'état premier basé sur le nombre de solutions quadratiques trouvées pour une valeur particulière. Le code est plus lent que tree fusion des entreprises D'État en raison d'environ trois fois le nombre d'opérations aériennes traitant environ trois fois le nombre de séquences un peu plus complexes, mais il ya probablement une gamme de tamisage de très grands nombres où il passera L'entreprise d'État en raison de son avantage théorique de performance.

le code suivant est fidèle à la formulation de la SoA, utilise des types de flux Coinductifs plutôt que Lazylist's ou des séquences que la memoization n'est pas nécessaire et la performance est meilleure, aussi ne pas utiliser de Victimes de Syndicats et évite de filtrage pour des raisons de performances:

#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 la roue arbre pliant optimisé SoE comme affiché dans une autre réponse à environ 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 a trouvé la performance que le meilleur des autres solutions fonctionnelles pures. D'autres optimisations qui pourraient être essayées sont que les séquences de nombres premiers carrés n'utilisent pas factorisation de roue pour éliminer les 357 multiples des carrés ou même utiliser seulement les multiples premiers des carrés premiers pour réduire le nombre de valeurs dans les flux de séquence de carrés et peut-être d'autres optimisations liées aux flux de séquence d'expression d'équation quadratique.

EDIT_ADD: j'ai pris un peu de temps pour examiner les optimisations modulo SoA et voir qu'en plus des optimisations "squarefree" ci-dessus, qui probablement cela ne fera pas beaucoup de différence, que les séquences quadratiques ont un motif modulo sur chaque 15 éléments qui permettrait à beaucoup de valeurs de test composite passées toggled d'être pré-criblé et éliminerait le besoin pour les opérations spécifiques modulo 12 pour chaque nombre composite. Toutes ces optimisations sont susceptibles d'entraîner une réduction de calcul de travail soumis à l'arbre de pliage jusqu'à environ 50% pour faire un peu plus version optimisée de la SoA exécuter proche ou légèrement c'est mieux que le meilleur arbre fusionnant SoE. Je ne sais pas où 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 vont en effet augmenter la performance d'environ un facteur de deux, je vois pourquoi la performance empirique actuelle avec l'augmentation de n n'est pas aussi bonne que SoE: alors que le SoE est particulièrement adapté à l'arbre opérations de repliement en ce que les premières séquences sont plus denses et se répètent plus souvent avec des séquences plus tardives beaucoup moins denses, les séquences "4x" SoA sont plus denses pour les séquences plus tardives lorsqu'elles sont ajoutées et alors que les séquences "3X" commencent moins denses, elles deviennent plus denses à mesure que y approche zéro, puis deviennent moins denses à nouveau; cela signifie que les séquences appel/retour ne sont pas maintenues à une profondeur minimale comme pour SoE mais que cette profondeur augmente au-delà d'être proportionnelle à la plage de nombres. En utilisant des Solutions le repliement n'est pas joli car on peut mettre en œuvre le repliement gauche pour les séquences qui augmentent en densité avec le temps, mais cela laisse encore les parties négatives des séquences "3X" mal optimisées, comme le fait de briser les séquences "3X" EN parties positives et négatives. La solution la plus simple est de sauvegarder toutes les séquences sur une carte, ce qui signifie que le temps d'accès augmentera de quelque chose comme le log de la racine carrée de la gamme, mais ce sera mieux pour une plus grande gamme de nombres que l'arbre actuel. pliant. 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 de manière impérative au code purement fonctionnel en F#. Il fournit des exemples d'utilisation de continuations comme dans les flux Coinductifs pour mettre en œuvre la paresse sans utiliser le type paresseux, mettre en œuvre (queue) boucles récursives pour éviter toute exigence de mutabilité, filetage d'un accumulateur (tgl) par des appels récursifs à obtenir un résultat (le nombre de fois que les équations quadratiques "a frappé" le nombre testé), présentant les solutions aux équations comme (paresseux) séquences (ou flux dans ce cas), etcetera.

pour ceux qui voudraient jouer plus loin avec ce code, même sans un système de développement basé sur Windows, je l'ai aussi posté sur tryfsharp.org et Ideone.com , bien qu'il fonctionne plus lentement sur ces deux plates-formes, avec tryfsharp à la fois proportionnel à la vitesse de la machine client locale et plus lent en raison de l'exécution sous Silverlight, et Ideone tournant sur le serveur Linux CPU sous Mono-projet 2.0, qui est notoirement lent dans la mise en œuvre et en particulier en ce qui concerne la collecte des ordures.

Je ne pense pas que cette question soit complète en ne considérant que des algorithmes purement fonctionnels qui cachent un État soit dans une file D'attente Map ou Priority dans le cas de quelques réponses ou un arbre de fusion plié dans le cas de l'une de mes autres réponses en ce que l'une de ces réponses est assez limitée quant à l'utilisabilité pour de grandes gammes de nombres premiers en raison de leur performance approximative O(N^1.2) ('^'signifie soulevé à la puissance de où n est le nombre supérieur dans la séquence) ainsi que leurs frais généraux de abattage. Cela signifie que même pour la gamme 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 n'est pas quelque chose qui est très utilisable.

il y a eu plusieurs réponses présentant des solutions utilisant divers degrés de mutabilité, mais soit ils n'ont pas tiré pleinement profit de leur mutabilité et ont été inefficaces ou ont été des traductions très simplistes de code impératif et laid fonctionnellement. Il me semble que le tableau F# mutable n'est qu'une autre forme d'état mutable caché à l'intérieur d'une structure de données, et qu'un algorithme efficace peut être développé qui a aucune autre mutabilité utilisée autre que le tableau de tampon mutable utilisé pour l'élimination efficace de nombres composés par segments de tampon en pagination, le reste du code étant écrit dans le style fonctionnel pur.

le le code suivant a été développé après avoir vu le code par Jon Harrop , et améliore ces idées comme suit:

1. le code de Jon échoue en termes de grande utilisation de la mémoire (sauve tous les nombres premiers générés au lieu de seulement les nombres premiers à la racine carrée du premier candidat le plus élevé, et régénère continuellement des matrices de tampon de taille toujours croissante (égale à la taille du dernier nombre premier trouvé) indépendamment des tailles de cache CPU.

2. de même, son code tel qu'il est présenté ne comprend pas de séquence génératrice.

3. de plus, le code tel qu'il est présenté n'a pas les optimisations nécessaires pour traiter au moins les nombres impairs et encore moins pour ne pas envisager l'utilisation de la factorisation des roues.

si le code de Jon a été utilisé pour générer la gamme de nombres premiers à la gamme de 32 bits de quatre de plus, il aurait un besoin de mémoire de gigaoctets pour les nombres premiers sauvegardés dans la structure de liste et un autre multi-gigaoctets pour le tampon tamis, bien qu'il n'y ait aucune raison réelle que ce dernier ne puisse pas être d'une taille fixe plus petite. Une fois que le tampon tamis dépasse la taille des tailles de cache CPU, les performances vont rapidement se détériorer dans le "cache thrashing", avec une perte de performance croissante, d'abord la L1, puis la L2, et finalement les tailles L3 (si présentes) sont dépassées.

C'est pourquoi le code de Jon ne calculera des nombres premiers jusqu'à environ 25 millions, même sur ma machine 64 bits avec huit gigaoctets de mémoire, avant de générer une exception hors mémoire, et il explique aussi pourquoi il y a une baisse de plus en plus grande de la performance relative alors que les gammes augmentent avec une performance O(N^1.4) avec une portée croissante et n'est que quelque peu sauvegardée parce qu'elle a une si faible complexité de calcul pour commencer.

le code suivant adresse toutes ces limites, en ce qu'il ne mémoize les nombres premiers de base jusqu'à la racine carrée du nombre maximum dans la gamme qui sont calculés comme nécessaire (seulement quelques kilo-octets dans le cas de la gamme de nombre de 32 bits) et n'utilise que de très petits tampons d'environ seize kilo-octets pour chacun des générateurs de nombres premiers de base et le filtre à tamis segmenté de la page principale (plus petit que la taille de cache L1 de la plupart des CPU modernes), ainsi qu'en incluant le code de séquence génératrice et (actuellement) être optimisé pour seulement tamisage pour les nombres impairs, ce qui signifie que la mémoire est utilisée plus efficacement. En outre, un tableau de bits empaquetés est utilisé pour améliorer encore l'efficacité de la mémoire; son coût de calcul est la plupart du temps compensé en moins de calculs devant être effectués lors de la numérisation 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"

ce nouveau code calcule les 203,280,221 nombres premiers dans la gamme de 32 bits dans environ ajouté/ corrigé: 25,4 secondes avec temps de fonctionnement pour les premiers 100000, un million, 10 millions, et 100 millions testés comme 0,01, 0,088, 0,94, et 11,25 secondes, respectivement sur un ordinateur de bureau rapide (i7-2700K @ 3,5 GHz), et peut fonctionner sur tryfsharp.org et ideone.com , mais dans une moindre mesure pour ce dernier en raison de 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 son augmentation la complexité du calcul, mais le dépasse très rapidement pour de plus grandes gammes en raison de son meilleur algorithme de performance qui compense cette complexité de sorte qu'il est environ cinq fois plus rapide pour le 10 millionième prime et environ sept fois plus rapide juste avant que le code de Jon explose à environ le 25 millionième prime.

du temps total d'exécution, plus de la moitié est dépensée dans l'énumération de la séquence de base et ne serait donc pas aidé dans une grande mesure en exécutant le nombre composite les opérations d'abattage en tant qu'opérations de fond, bien que les optimisations de la factorisation de la roue en combinaison aideraient (bien que plus intensives en calcul, cette complexité se déroulerait en arrière-plan) en ce sens qu'elles réduisent le nombre d'opérations d'essai de tampon nécessaires au dénombrement. D'autres optimisations pourraient être faites si l'ordre des séquences n'avait pas besoin d'être préservé (comme dans le simple comptage du nombre de nombres premiers ou dans la sommation des nombres premiers) car les séquences pourraient être écrites pour supporter les les interfaces d'énumération parallèles ou le code pourrait être écrit comme une classe de sorte que les méthodes de membre pourraient faire le calcul sans énumération. Ce code pourrait facilement être ajusté pour offrir le même type de performance que le code C# mais plus succinctement exprimé, bien qu'il n'atteindra jamais la performance du code natif c++ optimisé tel que PrimeSieve qui a été optimisé principalement pour la tâche de compter le nombre de nombres premiers sur une gamme et peut calculer le nombre de nombres premiers dans la gamme de nombre de 32 bits est une petite fraction d'une seconde (0.25 secondes sur le i7-2700K).

ainsi, d'autres optimisations seraient concentrées autour d'un empilage de bits supplémentaires le réseau de tamisage en utilisant la factorisation de roue pour minimiser le travail fait dans l'élimination des nombres composés, en essayant d'augmenter l'efficacité de l'énumération des nombres premiers résultants, et reléguer toute élimination composite aux fils de fond où un processeur de quatre à huit noyaux pourrait masquer la complexité de calcul supplémentaire nécessaire.

et c'est surtout du pur code fonctionnel, juste qu'il utilise un réseau mutable pour fusionner l'abattage composite....

Je ne suis pas très familier avec Haskell multimaps, mais le F# Power Pack a une classe HashMultiMap, dont le résumé xmldoc est: "tables de hachage, par défaut basé sur F# structural" hachage " et (=) fonctions. La table peut mapper une seule clé pour plusieurs reliures."Peut-être cela peut vous aider?

2 * 10^6 en 1 sec 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)