algorithm - metodo - numeros primos del 1 al 100
El tamiz de Eratóstenes en F# (14)
Al leer ese artículo, se me ocurrió una idea que no requiere un multimap. Maneja las teclas del mapa que colisionan moviendo la tecla colisionante hacia adelante por su valor principal una y otra vez hasta que alcanza una tecla que no está en el mapa. Debajo de primes hay un mapa con las claves del siguiente valor del iterador y valores que son primos.
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
Aquí está el algoritmo basado en cola de prioridad de ese documento sin la optimización cuadrada. Coloqué las funciones de cola de prioridad genérica en la parte superior. Usé una tupla para representar los iteradores de la lista perezosa.
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)))
Aquí está el algoritmo basado en la cola de prioridad con la optimización cuadrada. Para facilitar el uso de primos de adición diferidos en la tabla de búsqueda, los desplazamientos de rueda debían devolverse junto con los valores primos. Esta versión del algoritmo tiene un uso de memoria O (sqrt (n)) donde el ninguno optimizado es 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)))
Aquí está mi programa de prueba.
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
Estoy interesado en una implementación del tamiz de eratostenes en F # puramente funcional. Me interesa una implementación del tamiz real, no la implementación funcional ingenua que no es realmente el tamiz , por lo que no es algo como esto:
let rec PseudoSieve list =
match list with
| hd::tl -> hd :: (PseudoSieve <| List.filter (fun x -> x % hd <> 0) tl)
| [] -> []
El segundo enlace arriba describe brevemente un algoritmo que requeriría el uso de un multimap, que no está disponible en F # hasta donde yo sé. La implementación de Haskell proporcionada utiliza un mapa que admite un método insertWith , que no he visto disponible en el mapa funcional F # .
¿Alguien sabe una forma de traducir el código del mapa Haskell dado a F #, o tal vez conoce métodos alternativos de implementación o algoritmos de cribado que son tan eficientes y más adecuados para una implementación funcional o F #?
Aunque ha habido una respuesta que proporciona un algoritmo utilizando una cola de prioridad (PQ) como en un SkewBinomialHeap , quizás no sea la PQ correcta para el trabajo. Lo que requiere el tamiz incremental de Eratóstenes (iEoS) es un PQ que tiene un rendimiento excelente para obtener el valor mínimo y reinsertar los valores en su mayor parte en la cola, pero no necesita el máximo rendimiento para agregar nuevos valores ya que iSoE solo agrega como nuevo valores un total de los números primos hasta la raíz cuadrada del rango (que es una pequeña fracción del número de reinserciones que ocurren una vez por reducción). El SkewBinomialHeap PQ en realidad no da mucho más que usar el Mapa incorporado que usa un árbol de búsqueda binaria equilibrado - todas las operaciones O (log n) - aparte de que cambia la ponderación de las operaciones ligeramente a favor de los requisitos del SoE. Sin embargo, SkewBinaryHeap aún requiere muchas operaciones O (log n) por reducción.
Un PQ implementado como un Heap en particular como un Binary Heap y aún más particularmente como un MinHeap que satisface los requisitos de iSoE con O (1) desempeño para obtener el rendimiento mínimo y O (log n) para reinserciones y agregar nuevas entradas , aunque el rendimiento es en realidad una fracción de O (log n) ya que la mayoría de las reinserciones ocurren cerca de la parte superior de la cola y la mayoría de las adiciones de nuevos valores (que no importan ya que son poco frecuentes) ocurren cerca del final de la cola donde estas operaciones son más eficientes. Además, el MinHeap PQ puede implementar eficientemente la función de eliminar e insertar en una (por lo general, una fracción) de un pase O (log n). Entonces, en lugar del Mapa (que se implementa como un árbol AVL ) donde hay una operación O (log n) con un rango de ''log n'' completo debido al valor mínimo, necesitamos estar en la última hoja de la izquierda el árbol, generalmente estamos agregando y eliminando el mínimo en la raíz e insertándolo en el promedio de algunos niveles en una pasada. Por lo tanto, el MinHeap PQ puede usarse con solo una fracción de la operación O (log n) por reducción de eliminación selectiva en lugar de múltiples operaciones de fracción más grande O (log n).
El MinHeap PQ se puede implementar con código funcional puro (sin implementar "removeMin" ya que el iSoE no lo requiere, pero hay una función de "ajuste" para usar en la segmentación), como se muestra a continuación:
[<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
Utilizando el módulo anterior, el iSoE se puede escribir con las optimizaciones de factorización de rueda y el uso de corrientes co-inductivas eficientes (CIS) de la siguiente manera:
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 }
El código anterior calcula los primeros 100,000 primos en aproximadamente 0.077 segundos, los primeros 1,000,000 primos en 0.977 segundos, los primeros 10,000,000 primos en aproximadamente 14,33 segundos, y los primeros 100,000,000 cebadores en aproximadamente 221.87 segundos, todo en un i7-2700K (3.5GHz) como código de 64 bits. Este código puramente funcional es ligeramente más rápido que el código mutable del diccionario de Dustin Cambell con las optimizaciones comunes agregadas de factorización de rueda, adición diferida de primos de base y uso de los CID más eficientes agregados ( tryfsharp e ideone ), pero sigue siendo funcional puro código donde el uso de la clase Dictionary no es . Sin embargo, para rangos primos mayores de alrededor de dos mil millones (aproximadamente 100 millones de primos), el código que utiliza el diccionario basado en tabla hash será más rápido ya que las operaciones del diccionario no tienen un factor O (log n) y esta ganancia supera al computacional complejidad del uso de tablas hash de diccionario.
El programa anterior tiene la característica adicional de que la rueda de factorización está parametrizada de modo que, por ejemplo, uno puede usar una rueda extremadamente grande configurando WHLPRMS en [| 2u; 3u; 5u; 7u; 11u; 13u; 17u; 19u |] y FSTPRM a 23u para obtener un tiempo de ejecución de aproximadamente dos tercios para intervalos grandes a aproximadamente 9.34 segundos para diez millones de primos, aunque tenga en cuenta que lleva varios segundos calcule el WHLPTRN antes de que el programa comience a ejecutarse, que es una sobrecarga constante sin importar el rango principal.
Análisis comparativo : en comparación con la implementación de plegamiento incremental de árbol puramente funcional, este algoritmo es ligeramente más rápido porque la altura promedio utilizada del árbol MinHeap es menor en un factor de dos que la profundidad del árbol plegado, pero eso se compensa con un equivalente pérdida de factor constante en eficiencia en la capacidad de atravesar los niveles de árbol PQ debido a que se basa en un montón binario que requiere el procesamiento de las hojas derecha e izquierda para cada nivel de pila y una rama en ambos sentidos en lugar de una comparación individual por nivel para el árbol plegable con generalmente la rama menos profunda la tomada. En comparación con otros algoritmos funcionales basados en PQ y Map, las mejoras son generalmente por un factor constante en la reducción del número de operaciones O (log n) al atravesar cada nivel de las estructuras de árbol respectivas.
El MinHeap generalmente se implementa como un montón binario de matriz mutable después de un modelo de árbol genealógico inventado por Michael Eytzinger hace más de 400 años. Sé que la pregunta decía que no había interés en el código mutable no funcional, pero si uno debe evitar todo el código secundario que usa la mutabilidad, entonces no podríamos usar list o LazyList, que usan la mutabilidad "por debajo" por razones de rendimiento. Por lo tanto, imagine que la siguiente versión mutable alternativa de MinHeap PQ es la que proporciona una biblioteca y disfruta de otro factor de más de dos para rangos principales de mayor rendimiento:
[<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
Nota de Geek: en realidad había esperado que la versión mutable ofreciera una mejor relación de rendimiento mejorada, pero se atasca en las reinserciones debido a la estructura de código anidada if-then-else y al comportamiento aleatorio de los valores de descarte principales que significan la predicción de la bifurcación de la CPU falla para una gran proporción de las bifurcaciones, lo que da como resultado muchos 10 ciclos adicionales de reloj de la CPU por reducción de eliminación para reconstruir la memoria caché previa a la búsqueda.
Las únicas otras ganancias de rendimiento de factor constante en este algoritmo serían la segmentación y el uso de tareas múltiples para obtener una ganancia de rendimiento proporcional al número de núcleos de CPU; sin embargo, tal como está, este es el algoritmo SoE funcional más rápido hasta la fecha, e incluso la forma funcional pura que utiliza el funcional MinHeap supera implementaciones imperativas simplistas como el código de Jon Harrop o Sieve of Atkin de Johan Kullbom (que está equivocado en su tiempo ya que solo calculó los números primos a 10 millones en lugar de los 10 millones del primo ), pero esos algoritmos serían aproximadamente cinco veces más rápidos si se utilizaran mejores optimizaciones. Esa relación de aproximadamente cinco entre código funcional e imperativo se reducirá de alguna manera cuando agreguemos multi-threading de mayor factorización de rueda ya que la complejidad computacional del código imperativo aumenta más rápido que el código funcional y multi-threading ayuda al código funcional más lento que el código imperativo más rápido a medida que este último se acerca al límite base del tiempo requerido para enumerar a través de los primos encontrados.
EDIT_ADD: aunque se podría optar por seguir utilizando la versión funcional pura de MinHeap, agregar una segmentación eficiente en preparación para multihebras podría "romper" ligeramente la "pureza" del código funcional de la siguiente manera: 1) La forma más eficiente de la transferencia de una representación de primos compuestos-seleccionados es una matriz de bits empaquetados del tamaño del segmento, 2) Si bien se conoce el tamaño de la matriz, usar una comprensión de matriz para inicializarla de una manera funcional no es eficiente ya que utiliza " ResizeArray "debajo de las cubiertas que necesita copiarse para cada x adiciones (creo que ''x'' es ocho para la implementación actual) y el uso de Array.init no funciona ya que muchos valores en índices particulares se saltan, 3) Por lo tanto, el La forma más fácil de llenar el conjunto compuesto-entresacado es hacer cero el tamaño correcto y luego ejecutar una función de inicialización que podría escribir en cada índice de matriz mutable no más de una vez. Aunque esto no es estrictamente "funcional", está cerca de que la matriz se inicialice y luego nunca se modifique nuevamente.
El código con segmentación agregada, multi-threading, circunferencia factorial de rueda programable y muchos ajustes de rendimiento es el siguiente (aparte de algunas nuevas constantes adicionales, el código extra ajustado para implementar la segmentación y multi-threading es la parte inferior aproximadamente la mitad del código comenzando en la función "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)) }
Tenga en cuenta que los módulos MinHeap, tanto funcionales como basados en arreglos, han agregado una función de "ajuste" para permitir el ajuste del estado de exclusión de la versión de cada hilo de la PQ al comienzo de cada página de segmento nuevo. También tenga en cuenta que fue posible ajustar el código para que la mayor parte del cálculo se realice usando rangos de 32 bits con la salida de secuencia final como uint64 a bajo costo en tiempo de cálculo, de modo que actualmente el rango teórico es algo más de 100 billones (diez elevado a el poder de catorce) si uno estaba dispuesto a esperar los tres o cuatro meses necesarios para calcular ese rango. Las comprobaciones de rango numérico se eliminaron, ya que es poco probable que alguien use este algoritmo para calcular hasta ese rango y mucho menos pasarlo.
Utilizando MinHeap funcional puro y factorización de rueda 2,3,5,7, el programa anterior calcula los primeros cien mil, un millón, diez millones y cien millones de primos en 0.062, 0.629, 10.53 y 195.62 segundos, respectivamente. El uso de MinHeap basado en arreglos acelera esto hasta 0.097, 0.276, 3.48 y 51.60 segundos, respectivamente. Usando la rueda 2,3,5,7,11,13,17 cambiando WHLPRMS a "[| 2u; 3u; 5u; 7u; 11u; 13u; 17u |]" y FSTPRM a 19u acelera eso un poco más a 0.181, 0.308, 2.49 y 36.58 segundos, respectivamente (para una mejora constante del factor con una sobrecarga constante). Este ajuste más rápido calcula los 203,280,221 primos en el rango numérico de 32 bits en aproximadamente 88,37 segundos. La constante "BFSZ" se puede ajustar con compensaciones entre tiempos más lentos para rangos más pequeños versión tiempos más rápidos para rangos más grandes,con un valor de "1 <<< 14" recomendado para probar para los rangos más grandes. Esta constante solo establece el tamaño de búfer mínimo, con el programa ajustando el tamaño del búfer por encima de ese tamaño automáticamente para rangos más grandes, de modo que el búfer sea suficiente para que el mayor básico base requerido para el rango de página siempre "ataque" cada página al menos una vez ; esto significa que no se requiere la complejidad y la sobrecarga de un "tamiz de cubeta" adicional. Esta última versión completamente optimizada puede calcular los números primos de hasta 10 y 100 mil millones en aproximadamente 256.8 y 3617.4 segundos (poco más de una hora para los 100 mil millones) como se prueba utilizando "primesPQOWSE () |> Seq.takeWhile ((> =) 100000000000UL) |> Seq.fold (diversión sp -> s + 1UL) 0UL "para salida.Aquí es donde provienen las estimaciones de aproximadamente medio día para el conteo de números primos a un billón, una semana por hasta diez billones y de tres a cuatro meses por hasta cien billones.
No creo que sea posible hacer código funcional o casi funcional usando el algoritmo SoE incremental para correr mucho más rápido que esto. Como se puede ver al mirar el código, la optimización del algoritmo incremental básico ha aumentado en gran medida la complejidad del código de tal manera que es ligeramente más complejo que el código optimizado de manera equivalente basado en la eliminación de matriz directa con ese código capaz de ejecutarse aproximadamente diez veces más rápido esto y sin el exponente adicional en el rendimiento, lo que significa que este código incremental funcional tiene una sobrecarga de porcentaje adicional cada vez mayor.
Entonces, ¿esto es útil más que desde un punto de vista teórico e intelectual interesante? Probablemente no lo sea Para rangos más pequeños de números primos de hasta diez millones, los mejores SoE funcionales incrementales no optimizados básicos son probablemente adecuados y bastante simples de escribir o tienen menos memoria RAM que los SoE imperativos más simples. Sin embargo, son mucho más lentos que el código más imperativo al usar una matriz, por lo que se "agotan" para rangos superiores a ese. Si bien se ha demostrado aquí que el código se puede acelerar mediante la optimización, todavía es 10 veces más lento que una versión más imperativa basada en arreglos puros, pero se ha agregado complejidad para ser al menos tan complejo como ese código con optimizaciones equivalentes. ,e incluso ese código bajo F # en DotNet es aproximadamente cuatro veces más lento que usar un lenguaje como C ++ compilado directamente a código nativo; si uno realmente quisiera investigar grandes rangos de primos, uno probablemente usaría uno de esos otros lenguajes y técnicas dondeprimesieve puede calcular el número de números primos en el rango de cien billones en menos de cuatro horas en lugar de los tres meses requeridos para este código. END_EDIT_ADD
2 * 10 ^ 6 en 1 segundo en 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)
Aquí está mi intento de una traducción razonablemente fiel del código Haskell a 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)
Como esta pregunta específicamente solicita otros algoritmos, proporciono la siguiente implementación:
o quizás conozca métodos de implementación alternativos o algoritmos de tamizado
Ningún envío de varios algoritmos de Sieve of Eratóstenes (SoE) es realmente completo sin mencionar el Sieve of Atkin (SoA), que de hecho es una variación de SoE que usa las soluciones de un conjunto de ecuaciones cuadráticas para implementar el sacrificio compuesto así como eliminar todos los múltiplos de los cuadrados de los primos base (primos menos o igual a la raíz cuadrada del número más alto probado para primalidad). Teóricamente, el SoA es más eficiente que el SoE en cuanto a que hay operaciones ligeramente menores en el rango, por lo que debería tener un 20% menos de complejidad en el rango de 10 a 100 millones, pero en general es más lento debido a la carga computacional de la complejidad de resolver varias ecuaciones cuadráticas. Aunque, el altamente optimizadoLa implementación C de Daniel J. Bernstein es más rápida que la implementación SoE contra la cual lo probó para ese rango particular de números de prueba , la implementación SoE contra la cual probó no era la más óptima y las versiones más optimizadas de SoE directo son aún más rápidas. Este parece ser el caso aquí, aunque admito que puede haber más optimizaciones que me he perdido.
Dado que O''Neill en su artículo sobre el SoE utilizando incrementales Tamices ilimitados se propuso principalmente para mostrar que el Turner Sieve no es SoE tanto en cuanto a algoritmo como en cuanto a rendimiento, no consideró muchas otras variaciones del SoE como SoA. Al hacer una búsqueda rápida de literatura, no puedo encontrar ninguna aplicación de SoA a las secuencias incrementales ilimitadas que estamos discutiendo aquí, así que la adapté yo mismo como en el siguiente código.
Del mismo modo que se puede considerar que el caso puro ilimitado de SoE tiene como números compuestos una secuencia ilimitada de secuencias de primos múltiples, el SoA considera tener como potenciales primos la secuencia ilimitada de las secuencias ilimitadas de todas las expresiones de las ecuaciones cuadráticas con una de las dos variables libres, ''x'' o ''y'' fijadas a un valor inicial y con una secuencia de "eliminación" separada de las secuencias de todos los múltiplos de los primos base, que es muy similar a las secuencias de eliminación compuestas de secuencias para SoE, excepto que las secuencias avanzan más rápidamente por el cuadrado de los números primos en lugar de por un (menor) múltiplo de los números primos. He tratado de reducir el número de secuencias de ecuaciones cuadráticas expresadas al reconocer que, a los fines de un tamiz incremental, el "3 * x ^ 2 + y ^ 2 "y las secuencias" 3 * x ^ 2 - y ^ 2 "son en realidad lo mismo, excepto por el signo del segundo término y la eliminación de todas las soluciones que no son extrañas, además de aplicar 2357 factorización de rueda (aunque el SoA ya tiene factorización de rueda 235 inherente). Utiliza el eficiente algoritmo de fusión / combinación de plegado de árbol como en la fusión de árbol SoE para tratar cada secuencia de secuencias pero con una simplificación que el operador de unión no combina al fusionarse. el algoritmo SoA depende de poder alternar el estado principal en función del número de soluciones cuadráticas para un valor determinado. El código es más lento que el árbol que fusiona SoE debido a aproximadamente tres veces el número de operaciones generales que trata de aproximadamente tres veces el número de secuencias algo más complejas,pero es probable que haya un rango de tamizado de números muy grandes en los que pasará SoE debido a su ventaja de rendimiento teórico.
El siguiente código es fiel a la formulación del SoA, utiliza los tipos CoInductive Stream en lugar de las secuencias o secuencias de LazyList ya que no se requiere memoria y el rendimiento es mejor, tampoco usa Sindicatos discriminados y evita la coincidencia de patrones por motivos de rendimiento:
#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 }
Como se indicó, el código es más lento que SoE optimizado para Tree Folding Wheel como se publicó en otra respuesta a aproximadamente medio segundo para los primeros 100,000 primos, y tiene aproximadamente el mismo O empírico (n ^ 1.2) para primos encontró desempeño como el mejor de otras soluciones funcionales puras Algunas optimizaciones adicionales que podrían probarse son que las secuencias cuadradas de primos no usan la factorización de rueda para eliminar los 357 múltiplos de los cuadrados o incluso usan solo los múltiplos principales de los cuadrados principales para reducir el número de valores en las secuencias de secuencia de cuadrados y quizás otras optimizaciones relacionadas con las secuencias de secuencia de expresión de ecuación cuadrática.
EDIT_ADD: me he tomado un poco de tiempo para analizar las optimizaciones del módulo SoA y ver que, además de las optimizaciones "sin cuadrados" anteriores, que probablemente no harán mucha diferencia, las secuencias cuadráticas tienen un patrón de módulo sobre cada 15 elementos que permitiría que muchos de los valores de prueba compuestos pasados conmutados sean preseleccionados y eliminaría la necesidad de las operaciones específicas del módulo 12 para cada número compuesto. Todas estas optimizaciones probablemente den como resultado una reducción en el trabajo computacional enviado al árbol plegado de hasta aproximadamente el 50% para hacer una versión ligeramente más optimizada de la ejecución de SoA cercana o ligeramente mejor que la mejor fusión de dobleces de árboles SoE. No sé cuándo podría encontrar el tiempo para hacer estos pocos días de investigación para determinar el resultado.END_EDIT_ADD
EDIT_ADD2:Al trabajar en las optimizaciones anteriores que de hecho aumentarán el rendimiento en aproximadamente un factor de dos, veo por qué el rendimiento empírico actual al aumentar n no es tan bueno como SoE: mientras que el SoE es particularmente adecuado para las operaciones de plegado de árboles en que el primero las secuencias son más densas y repiten más a menudo con secuencias posteriores mucho menos densas, las secuencias SoA "4X" son más densas para secuencias posteriores cuando se agregan y mientras que las secuencias "3X" comienzan menos densas, se vuelven más densas a medida que y se aproxima a cero , luego vuelve a ser menos denso; esto significa que las secuencias de llamada / retorno no se mantienen a una profundidad mínima como para SoE, pero que esta profundidad aumenta más allá de ser proporcional al rango numérico. Soluciones que usan plegado aren ''Es bonito como se puede implementar el plegado a la izquierda para las secuencias que aumentan en densidad con el tiempo, pero eso aún deja las porciones negativas de las secuencias "3X" pobremente optimizadas, al igual que dividir las secuencias "3X" en porciones positivas y negativas. La solución más sencilla es probable que guarde todas las secuencias en un Mapa, lo que significa que el tiempo de acceso aumentará en algo así como el registro de la raíz cuadrada del rango, pero será mejor para un mayor rango numérico que el plegado actual del árbol.pero eso será mejor para un mayor rango numérico que el actual plegado de árboles.pero eso será mejor para un mayor rango numérico que el actual plegado de árboles.END_EDIT_ADD2
Aunque más lento, presento esta solución aquí para mostrar cómo se puede desarrollar el código para expresar ideas originalmente concebidas imperativamente a código funcional puro en F #. Proporciona ejemplos del uso de continuaciones como en CoInductive Streams para implementar la pereza sin utilizar el tipo Lazy, implementando bucles recursivos (tail) para evitar cualquier requisito de mutabilidad, enhebrando un acumulador (tgl) a través de llamadas recursivas para obtener un resultado (el número de veces las ecuaciones cuadráticas "golpearon" el número probado), presentando las soluciones a las ecuaciones como secuencias (perezosas) (o corrientes en este caso), etcétera.
Para aquellos que quieran jugar aún más con este código incluso sin un sistema de desarrollo basado en Windows, también lo he publicado en tryfsharp.org y Ideone.com , aunque se ejecuta más lento en ambas plataformas, con tryfsharp ambos proporcionales a la velocidad de la máquina cliente local y más lenta debido a que se ejecuta bajo Silverlight, y Ideone se ejecuta en la CPU del servidor Linux bajo Mono-project 2.0, que es notoriamente lenta tanto en la implementación como en particular en cuanto a las recolecciones de basura.
I''m not very familiar with Haskell multimaps, but the F# Power Pack has a HashMultiMap class, whose xmldoc summary is: "Hash tables, by default based on F# structural "hash" and (=) functions. The table may map a single key to multiple bindings." Perhaps this might help you?
Tamiz Prime implementado con procesadores de buzones:
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
Aquí están mis dos centavos, aunque no estoy seguro de que cumpla con el criterio del OP de ser verdaderamente el tamiz de Eratosthenes. No utiliza división modular e implementa una optimización del documento citado por el OP. Solo funciona para listas finitas, pero me parece que está en el espíritu de cómo se describió originalmente el tamiz. Como un aparte, el documento habla sobre la complejidad en términos del número de marcas y el número de divisiones. Parece que, como tenemos que atravesar una lista enlazada, esto quizás ignora algunos aspectos clave de los diversos algoritmos en términos de rendimiento. En general, aunque la división modular con computadoras es una operación costosa.
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]
EDITAR: corrigió un error en el código de mi publicación original. Intenté seguir la lógica original del tamiz con algunas modificaciones. A saber, comenzar con el primer elemento y tachar los múltiplos de ese elemento del conjunto. Este algoritmo literalmente busca el siguiente ítem que es un múltiplo del primo en lugar de hacer una división modular en cada número del conjunto. Una optimización del documento es que comienza a buscar múltiplos de la prima mayor que p ^ 2.
La parte en la función auxiliar con el multinivel trata con la posibilidad de que el siguiente múltiplo del primo ya se haya eliminado de la lista. Entonces, por ejemplo, con el prime 5, tratará de eliminar el número 30, pero nunca lo encontrará porque ya fue eliminado por el primer 3. Espero que aclare la lógica del algoritmo.
Aquí hay otro método más para lograr el tamiz incremental de Eratóstenes (SoE) usando solo código F # funcional puro. Es una adaptación del código Haskell desarrollado como "Esta idea se debe a Dave Bayer, aunque usó una formulación compleja y una estructura equilibrada de árboles ternarios, profundizando progresivamente de manera uniforme (formulación simplificada y una estructura de árbol binario sesgada y más profunda introducida por Heinrich Apfelmus, simplificado aún más por Will Ness). Idea de producción por etapas debida a M. O''Neill "según el siguiente enlace: código optimizado de plegado del árbol utilizando una rueda factorial en Haskell .
El siguiente código tiene varias optimizaciones que lo hacen más adecuado para la ejecución en F #, de la siguiente manera:
El código usa flujos coinductivos en lugar de LazyList, ya que este algoritmo no tiene (o poca) necesidad de memoria de LazyList y mis flujos coinductivos son más eficientes que los LazyLists (del FSharp.PowerPack) o las secuencias integradas. Otra ventaja es que mi código se puede ejecutar en tryFSharp.org y ideone.com sin tener que copiar y pegar en el código fuente de Microsoft.FSharp.PowerPack Core para el tipo y módulo LazyList (junto con el aviso de copyright)
Se descubrió que hay bastante sobrecarga para la coincidencia de patrones de F # en los parámetros de función, por lo que el tipo de unión discriminada anterior más legible utilizando tuplas se sacrificó a favor de la estructura de valor by (o clase como se ejecuta más rápido en algunas plataformas) tipos de un factor de dos o más aceleración.
Las optimizaciones de Ness van desde el plegado lineal de árboles hasta el plegado bilateral y el plegado multidireccional y las mejoras utilizando la factorización de ruedas son tan efectivas ratiométricamente para F # como lo son para Haskell, con la diferencia principal entre los dos idiomas que Haskell puede compilarse código nativo y tiene un compilador más altamente optimizado, mientras que F # tiene más gastos generales bajo el sistema 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: Esto se ha corregido ya que el código original no manejaba adecuadamente la cola de la secuencia y pasaba la cola de la secuencia de parámetros a la función de pares para la función joinT3 en lugar de la cola que sigue a la secuencia zs. El tiempo que se muestra a continuación también se corrigió, con una aceleración extra del 30%. Los códigos de enlace tryFSharp e ideone también se han corregido. END_EDIT_ADD
El programa anterior funciona a aproximadamente O (n ^ 1.1) rendimiento con n el primo máximo calculado o aproximadamente O (n ^ 1.18) cuando n es el número de primos calculados, y tarda aproximadamente 2,16 segundos para calcular el primer millón de primos (aproximadamente 0,14 segundo para los primeros 100.000 números primos) en una computadora rápida que ejecuta código de 64 bits utilizando tipos de estructura en lugar de clases (parece que algunas implementaciones encapsulan y unbox las estructuras de by-value al formar cierres). Considero que se trata del rango práctico máximo para cualquiera de estos algoritmos primarios funcionales puros. Probablemente este sea el más rápido con el que se puede ejecutar un algoritmo SoE funcional puro aparte de algunos ajustes menores para reducir los factores constantes.
Además de combinar la segmentación y el multihilo para compartir el cálculo entre múltiples núcleos de CPU, la mayoría de los "ajustes" que se pueden realizar a este algoritmo se basan en aumentar la circunferencia de la factorización de rueda para obtener una ganancia de hasta 40% en rendimiento y ganancias menores debido a ajustes en el uso de estructuras, clases, tuplas o más parámetros individuales directos en el paso del estado entre funciones.
EDIT_ADD2: He realizado las optimizaciones anteriores con el resultado de que el código ahora es casi el doble de rápido debido a las optimizaciones de la estructura con la ventaja añadida de utilizar opcionalmente circunferencias de factorización de rueda más grandes para la reducción más pequeña agregada. Tenga en cuenta que el siguiente código evita el uso de continuaciones en el bucle de generación de secuencia principal y solo las utiliza cuando sea necesario para las secuencias de primos de base y las secuencias de eliminación de números compuestos posteriores derivadas de esos primos de base. El nuevo código es el siguiente:
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 }
El código anterior toma alrededor de 0.07, 1.02 y 14.58 segundos para enumerar los primeros cien mil, un millón y diez millones de primos, respectivamente, todos en la máquina de referencia Intel i7-2700K (3.5 GHz) en el modo de 64 bits. Esto no es mucho más lento que la implementación de Haskell de referencia de la que se derivó este código, aunque es un poco más lento en tryfsharp e ideone debido a estar en modo de 32 bits para tryfsharp bajo Silverlight (aproximadamente la mitad más lenta) y funcionando en una máquina más lenta bajo Mono 2.0 (que es inherentemente mucho más lento para F #) en ideone, por lo que es hasta aproximadamente cinco veces más lento que la máquina de referencia. Tenga en cuenta que el tiempo de ejecución informado por ideone incluye el tiempo de inicialización para las matrices de la tabla de búsqueda incorporada, cuyo tiempo debe tenerse en cuenta.
El programa anterior tiene la característica adicional de que la rueda de factorización está parametrizada de modo que, por ejemplo, uno puede usar una rueda extremadamente grande configurando WHLPRMS en [| 2u; 3u; 5u; 7u; 11u; 13u; 17u; 19u |] y FSTPRM a 23u para obtener un tiempo de ejecución de aproximadamente dos tercios para intervalos grandes a aproximadamente 10.02 segundos para diez millones de primos, aunque tenga en cuenta que lleva varios segundos compute WHLPTRN antes de que el programa comience a ejecutarse.
Nota de Geek: no he implementado un "combinador de punto fijo no compartido para la producción de primo multietapa telescópico" según el código Haskell de referencia, aunque intenté hacer esto, porque uno necesita tener algo así como la lista perezosa de Haskell para que esto funcione sin correr de distancia en un bucle infinito y desbordamiento de la pila. Aunque mis Co-Inductive Streams (CIS) tienen algunas propiedades de holgazanería, no son formalmente listas perezosas o secuencias en caché (se convierten en secuencias no almacenadas en caché y pueden ser almacenadas en caché cuando pasan de modo que una función como el "genseq" proporciono para la secuencia de salida final). No quería usar la implementación de PowerPack de LazyList porque no es muy eficiente y requeriría que copie ese código fuente en tryfsharp e ideone, que no proporcionan módulos importados.El uso de las secuencias incorporadas (incluso almacenadas en caché) es muy ineficiente cuando uno quiere usar las operaciones de cabeza / cola como se requiere para este algoritmo ya que la única forma de obtener la cola de una secuencia es usar "Seq.skip 1" que en los usos múltiples producen una nueva secuencia basada en la secuencia original repetidamente omitida muchas veces. Podría implementar mi propia clase eficaz de LazyList basada en CIS, pero no parece que valga la pena demostrar un punto cuando los objetos recursivos "initcmpsts" y "baseprimes" toman muy poco código. Además, al pasar un LazyList a una función para producir extensiones a ese LazyList cuya función solo utiliza valores cercanos al comienzo del Lazylist, se requiere que casi todo el LazyList sea memorado para una reducción en la eficiencia de la memoria:un pase para los primeros 10 millones de primos requerirá una LazyList en memoria con casi 180 millones de elementos. Así que aprobé esto.
Tenga en cuenta que para rangos más grandes (10 millones de primos o más), este código puramente funcional es aproximadamente la misma velocidad que muchas implementaciones imperativas simplistas del Tamiz de Eratóstenes o Atkins, aunque eso se debe a la falta de optimización de esos algoritmos imperativos; una implementación más imperativa que esta usando optimizaciones equivalentes y matrices de cribado segmentadas seguirá siendo unas diez veces más rápida que esta según mi respuesta "casi funcional".
También tenga en cuenta que si bien es posible implementar el cribado segmentado utilizando plegamiento de árboles, es más difícil ya que los algoritmos de avance de descarte están enterrados dentro de las continuaciones usadas para el operador "union - ^" y trabajar en esto significaría un rango de descarte que avanza continuamente necesita ser usado; esto es diferente a otros algoritmos donde el estado de la variable de exclusión se puede restablecer para cada nueva página, incluyendo la reducción de su rango, de modo que si se utilizan rangos mayores que 32 bits, el rango interno seleccionado aún puede reiniciarse para operar dentro de los 32 intervalo de bits incluso cuando se determina un rango de primos de 64 bits a bajo costo en el tiempo de ejecución por primo. END_EDIT_ADD2
Aquí hay un algoritmo incremental (y recursivo) basado en algoritmo de Eratóstenes con secuencias optimizado, ya que no hay necesidad de memorizar valores de secuencia previos (aparte de que hay una ligera ventaja en el almacenamiento en memoria caché de los valores primos base usando Seq. caché), con las principales optimizaciones es que utiliza la factorización de rueda para la secuencia de entrada y que utiliza múltiples (recursivas) secuencias para mantener los primos base que son menores que la raíz cuadrada del último número que se tamiza, de la siguiente manera:
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) }
Encuentra 100.000 primos hasta 1.299.721 en aproximadamente 0.445 segundos, pero no siendo un algoritmo EoS imperativo adecuado, no escala casi linealmente con un mayor número de primos, tarda 7.775 segundos para encontrar los 1.000.000 de prima hasta 15.485.867 para un rendimiento sobre este rango de aproximadamente O (n ^ 1.2) donde n es el máximo principal encontrado.
Hay un poco más de ajuste que se puede hacer, pero probablemente no va a hacer mucha diferencia en cuanto a un gran porcentaje en un mejor rendimiento de la siguiente manera:
Como la biblioteca de secuencias F # es notablemente lenta, se podría usar un tipo autodefinido que implemente IEnumerable para reducir el tiempo que se pasa en la secuencia interna, pero como las operaciones de secuencia solo toman alrededor del 20% del tiempo total, incluso si éstas se redujeron a tiempo cero el resultado sería solo una reducción al 80% del tiempo.
Se podrían probar otras formas de almacenamiento de mapas, como una cola de prioridad mencionada por O''Neil o SkewBinomialHeap tal como la utiliza @gradbot, pero al menos para SkewBinomialHeap, la mejora en el rendimiento es de solo un pequeño porcentaje. Parece que al elegir diferentes implementaciones de mapas, uno solo está negociando una mejor respuesta para encontrar y eliminar elementos que están cerca del comienzo de la lista en comparación con el tiempo invertido en insertar nuevas entradas para disfrutar de esos beneficios, por lo que la ganancia neta es casi un lavado y todavía tiene un rendimiento O (log n) con entradas crecientes en el mapa. La optimización anterior utilizando múltiples flujos de entradas solo para la raíz cuadrada reduce el número de entradas en el mapa y, por lo tanto, las mejoras no tienen mucha importancia.
EDIT_ADD: Hice un poco más de optimización y el rendimiento mejoró algo más de lo esperado, probablemente debido a la forma mejorada de eliminar el Seq.skip como una forma de avanzar a través de la secuencia de primos base. Esta optimización utiliza un reemplazo para la generación de secuencia interna como una tupla de valor entero y una función de continuación utilizada para avanzar al siguiente valor en la secuencia, con la secuencia F # final generada por una operación de desplegado general. El código es el siguiente:
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 }
Los tiempos requeridos para encontrar los primos 100,000 y 1,000,000 son aproximadamente 0,31 y 5,1 segundos, respectivamente, por lo que hay un porcentaje considerable de ganancia para este pequeño cambio. Probé mi propia implementación de las interfaces IEnumerable / IEnumerator que son la base de las secuencias, y aunque son más rápidas que las versiones usadas por el módulo Seq, apenas hacen una diferencia más en este algoritmo donde la mayor parte del tiempo se usa en el Funciones del mapa END_EDIT_ADD
Además de las implementaciones de EoS incrementales basadas en mapas, existe otra implementación "puramente funcional" que usa Tree Folding que se dice que es ligeramente más rápida, pero como todavía tiene un término O (log n) en el plegado del árbol, sospecho que principalmente será más rápido (si es así) debido a cómo se implementa el algoritmo en cuanto a los números de operaciones de la computadora en comparación con el uso de un mapa. Si las personas están interesadas, desarrollaré esa versión también.
Al final, uno debe aceptar que ninguna implementación funcional pura de la EoS incremental se acercará alguna vez a la velocidad de procesamiento sin procesar de una implementación imperativa buena para rangos numéricos más grandes. Sin embargo, se podría llegar a un enfoque en el que todo el código es puramente funcional excepto por el cribado segmentado de números compuestos en un rango que utiliza una matriz (mutable) que se acercaría al rendimiento de O (n) y en el uso práctico sería de cincuenta a cien veces más rápido que los algoritmos funcionales para rangos grandes, como los primeros 200,000,000 primos. Esto ha sido hecho por @Jon Harrop en fsharpnews.blogspot.com/2010/02/sieve-of-eratosthenes.html , pero esto podría ajustarse aún más con muy poco código adicional.
En realidad, traté de hacer lo mismo, probé primero la misma implementación ingenua que en cuestión, pero fue demasiado lento. Luego encontré esta página YAPES: Problem Seven, Part 2 , donde usó Sieve of Eratosthenes real basado en Melissa E. O''Neill. Tomé el código de allí, solo un poco lo modifiqué, porque F # cambió un poco desde la publicación.
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
No creo que esta pregunta sea completa solo considerando algoritmos puramente funcionales que ocultan estado en un mapa o en una cola de prioridad en el caso de unas pocas respuestas o en un árbol de fusión doblado en el caso de una de mis otras respuestas en que cualquiera de estos son bastante limitados en cuanto a la usabilidad para grandes rangos de primos debido a su rendimiento aproximado de O (n ^ 1.2) (''^'' significa elevado a la potencia de donde n es el número superior en la secuencia) así como su sobrecarga computacional por operación de descarte. Esto significa que incluso para el rango numérico de 32 bits, estos algoritmos tomarán algo en el rango de al menos muchos minutos para generar los números primos de hasta cuatro mil millones más, lo cual no es algo muy útil.
Ha habido varias respuestas que presentan soluciones con diversos grados de mutabilidad, pero o bien no han aprovechado al máximo su mutabilidad y han sido ineficientes o han sido traducciones muy simplistas del código imperativo y feo funcionalmente. Me parece que la matriz mutable F # es simplemente otra forma de ocultar el estado mutable dentro de una estructura de datos, y que se puede desarrollar un algoritmo eficiente que no tenga otra mutabilidad utilizada que no sea la matriz de buffer mutable utilizada para el sacrificio eficiente de números compuestos por segmentos de búfer paginado, con el resto del código escrito en estilo funcional puro.
El siguiente código fue desarrollado después de ver el código de Jon Harrop , y mejora esas ideas de la siguiente manera:
El código de Jon falla en términos de uso de memoria alta (guarda todos los números primos generados en lugar de solo los números primos en la raíz cuadrada del primo candidato más alto, y regenera continuamente matrices de almacenamiento intermedio de tamaño cada vez mayor (igual al tamaño del último valor primo encontrado) independientemente de los tamaños de caché de la CPU.
Además, su código como se presenta no incluye una secuencia generadora.
Además, el código tal como se presenta no tiene las optimizaciones para, al menos, tratar con los números impares y mucho menos sin considerar el uso del uso de la factorización de las ruedas.
Si el código de Jon se usara para generar el rango de números primos en el rango numérico de 32 bits de más de cuatro mil millones, tendría un requisito de memoria de Gigabytes para los primos guardados en la estructura de lista y otro multi-Gigabytes para el tampón de criba, aunque no hay una razón real para que este último no sea de un tamaño fijo más pequeño. Una vez que el buffer tamiz excede el tamaño de los tamaños de caché de CPU, el rendimiento se deteriorará rápidamente en "cache thrashing", con una mayor pérdida de rendimiento ya que primero se superan los tamaños L1, L2 y finalmente L3 (si está presente).
Esta es la razón por la cual el código de Jon solo calculará números primos de hasta aproximadamente 25 millones, incluso en mi máquina de 64 bits con ocho gigabytes de memoria antes de generar una excepción de memoria insuficiente y también explica por qué hay una disminución mayor en el rendimiento a medida que los rangos aumentan con un rendimiento aproximado de O (n ^ 1,4) con un rango creciente y solo se salva de alguna manera porque tiene una complejidad computacional tan baja para empezar.
El siguiente código aborda todas estas limitaciones, ya que solo memoriza los primos base hasta la raíz cuadrada del número máximo en el rango que se calculan según sea necesario (solo unos pocos kilobytes en el caso del rango numérico de 32 bits) y solo utiliza búferes muy pequeños de aproximadamente dieciséis Kilobytes para cada uno de los generadores de primos base y el filtro de tamiz segmentado de la página principal (más pequeño que el tamaño de caché L1 de la mayoría de las CPU modernas), además de incluir el código de secuencia generadora y (actualmente) algo optimizado para solo cribar números impares, lo que significa que la memoria se usa de manera más eficiente. Además, se usa una matriz de bits empaquetados para mejorar aún más la eficiencia de la memoria; su costo de cálculo está compuesto principalmente por menos cálculos que deben hacerse para escanear el búfer.
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"
Este nuevo código calcula los 203,280,221 primos en el rango numérico de 32 bits en aproximadamente AGREGADO / CORREGIDO: 25,4 segundos con tiempos de funcionamiento para los primeros 100000, un millón, 10 millones y 100 millones probados como 0.01, 0.088, 0.94 y 11.25 segundos , respectivamente, en una computadora de escritorio rápida (i7-2700K @ 3.5 GHz), y puede ejecutarse en tryfsharp.org y ideone.com, aunque en un rango menor para este último debido a restricciones de tiempo de ejecución. Tiene un peor rendimiento que el código de Jon Harrop para rangos pequeños de unos pocos miles de primos debido a su mayor complejidad computacional, pero lo pasa muy rápidamente para rangos más grandes debido a su mejor algoritmo de rendimiento que compensa esa complejidad de manera que es aproximadamente cinco veces más rápido para la primo 10 millonésima y aproximadamente siete veces más rápido justo antes de que el código de Jon explote alrededor de la primo 25 millonésima.
Del tiempo total de ejecución, más de la mitad se gasta en la enumeración de secuencia básica y, por lo tanto, no se ayudaría en gran medida ejecutando las operaciones de selección de números compuestos como operaciones de fondo, aunque las optimizaciones de factorización de rueda en combinación ayudarían (aunque más computacionalmente intensivo, esa complejidad se ejecutaría en segundo plano) ya que reducen el número de operaciones de prueba de memoria intermedia requeridas en la enumeración. Se podrían realizar optimizaciones adicionales si el orden de las secuencias no necesitara ser preservado (como al contar el número de primos o al sumar los números primos) ya que las secuencias podrían escribirse para soportar las interfaces de enumeración paralelas o el código podría ser escrito como una clase para que los métodos miembros puedan hacer el cálculo sin enumeración.Este código se puede sintonizar fácilmente para ofrecer el mismo tipo de rendimiento que el código C # pero se expresa de forma más concisa, aunque nunca alcanzará el rendimiento del código nativo optimizado de C ++ comoprimesieve que se ha optimizado principalmente para la tarea de contar el número de primos en un rango y puede calcular el número de primos en el rango numérico de 32 bits, es una pequeña fracción de segundo (0,25 segundos en el i7-2700K).
Por lo tanto, optimizaciones adicionales se concentrarían alrededor de la agrupación de tamizado usando factorización de rueda para minimizar el trabajo realizado en el sacrificio de los números compuestos, tratando de aumentar la eficiencia de la enumeración de los primos resultantes, y relegando todos los descartes compuestos a los hilos de fondo cuando procesador central de cuatro a ocho podría ocultar la complejidad computacional adicional requerida.
Y es sobre todo código funcional puro, solo que usa una matriz mutable para combinar el sacrificio compuesto ...
Por lo que vale, esto no es un tamiz de Erathothenes, pero es muy rápido:
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
Encuentra la prima número 100.000 en 1.25 segundos en mi máquina (AMD Phenom II, 3.2GHZ quadcore).
Sé que usted declaró explícitamente que estaba interesado en una implementación de criba puramente funcional, por lo que esperé a presentar mi criba hasta ahora. Pero al volver a leer el documento al que hizo referencia, veo que el algoritmo de criba incremental presentado allí es esencialmente el mismo que el mío, la única diferencia son los detalles de implementación del uso de técnicas puramente funcionales versus técnicas decididamente imperativas. Así que creo que al menos califico a medias para satisfacer su curiosidad. Además, argumentaría que el uso de técnicas imperativas cuando se pueden lograr mejoras significativas en el rendimiento pero que están ocultas por interfaces funcionales es una de las técnicas más poderosas que se fomenta en la programación F #, en oposición a la cultura pura de Haskell. Primero publiqué esta implementación en mi proyecto Euler para el blog F # unpero vuelva a publicar aquí con el código de prerrequisito sustituido de nuevo y la eliminación estructural. primespuede calcular los primeros 100,000 primos en 0.248 segundos y los primeros 1,000,000 primos en 4.8 segundos en mi computadora (tenga en cuenta que primesalmacena en caché sus resultados, por lo que tendrá que volver a evaluarlo cada vez que realice un punto de referencia).
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)
}