C #: Attuazione del crivello di Atkin
https://stackoverflow.com/questions/1569127
-
21-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianDomanda
Mi chiedevo se qualcuno qui hanno una buona implementazione del crivello di Atkin che vorrebbero condividere.
Sto cercando di attuarlo, ma non riesco a avvolgere la mia testa intorno ad esso. Ecco quello che ho finora.
public class Atkin : IEnumerable<ulong>
{
private readonly List<ulong> primes;
private readonly ulong limit;
public Atkin(ulong limit)
{
this.limit = limit;
primes = new List<ulong>();
}
private void FindPrimes()
{
var isPrime = new bool[limit + 1];
var sqrt = Math.Sqrt(limit);
for (ulong x = 1; x <= sqrt; x++)
for (ulong y = 1; y <= sqrt; y++)
{
var n = 4*x*x + y*y;
if (n <= limit && (n % 12 == 1 || n % 12 == 5))
isPrime[n] ^= true;
n = 3*x*x + y*y;
if (n <= limit && n % 12 == 7)
isPrime[n] ^= true;
n = 3*x*x - y*y;
if (x > y && n <= limit && n % 12 == 11)
isPrime[n] ^= true;
}
for (ulong n = 5; n <= sqrt; n++)
if (isPrime[n])
for (ulong k = n*n; k <= limit; k *= k)
isPrime[k] = false;
primes.Add(2);
primes.Add(3);
for (ulong n = 5; n <= limit; n++)
if (isPrime[n])
primes.Add(n);
}
public IEnumerator<ulong> GetEnumerator()
{
if (!primes.Any())
FindPrimes();
foreach (var p in primes)
yield return p;
}
IEnumerator IEnumerable.GetEnumerator()
{
return GetEnumerator();
}
}
Ho praticamente solo cercato di "tradurre" il pseudocodice elencato a Wikipedia , ma non funziona correttamente. Quindi, o ho frainteso qualcosa o semplicemente fatto qualcosa di sbagliato. O più probabilmente sia ...
Avere un elenco dei primi 500 numeri primi che uso come test e la mia realizzazione non riesce al numero 40 (o 41?).
I valori differiscono di indice [40]
Previsto: 179
Ma era: 175
Sei in grado di trovare il mio errore, avete un'implementazione che in giro che si potrebbe condividere, o entrambi?
Il test esatto che sto usando assomiglia a questo:
public abstract class AtkinTests
{
[Test]
public void GetEnumerator_FirstFiveHundredNumbers_AreCorrect()
{
var sequence = new Atkin(2000000);
var actual = sequence.Take(500).ToArray();
var expected = First500;
CollectionAssert.AreEqual(expected, actual);
}
private static readonly ulong[] First500 = new ulong[]
{
2, 3, 5, 7, 11, 13, 17, ...
};
}
Soluzione
Il codice:
for (ulong k = n*n; k <= limit; k *= k)
isPrime[k] = false;
non sembra essere una traduzione fedele di questo pseudocodice:
is_prime(k) ← false, k ∈ {n², 2n², 3n², ..., limit}
Il tuo codice sembra che avrà una durata di n * n, n ^ 4, n ^ 8, ecc vale a dire quadratura ogni volta che invece di aggiungere n-squared ogni volta. Prova questo:
ulong nSquared = n * n;
for (ulong k = nSquared; k <= limit; k += nSquared)
isPrime[k] = false;
Altri suggerimenti
Il ultima risposta da Aaron Mugatroyd come dalla sorgente Python tradotto per un Crivello di Atkin (SOA) non è troppo male, ma può essere migliorato sotto diversi aspetti, come manca alcune ottimizzazioni importanti, come segue:
-
La sua risposta non usa la piena modulo 60 versione originale Atkin e Bernstein della Sieve, ma piuttosto un po 'migliorato la variazione del codice pseudo alle voci di Wikipedia in modo utilizza circa un fattore di 0,36 della gamma vaglio numerica operazioni snodata / riforma combinati; il mio codice di seguito utilizza la non-tabella codici segmento pseudo ragionevolmente efficiente come da i miei commenti in una risposta commentando il Crivello di Atkin che usa un fattore di circa 0,26 volte l'intervallo numerico per ridurre la quantità di lavoro svolto per circa un fattore di circa due settimi.
-
Il suo codice riduce la dimensione del buffer avendo solo rappresentazioni dispari, mentre il mio codice ulteriore bit dei pacchetti per eliminare qualsiasi rappresentazione dei numeri divisibili per tre e cinque nonché quelli divisibile per due implicito "odds sola "; ciò riduce il requisito di memoria di un ulteriore fattore di quasi la metà (a 8/15) e contribuisce a rendere migliore utilizzo delle cache CPU per un ulteriore aumento della velocità a causa del ridotto tempo medio di accesso alla memoria.
-
Il mio codice conta il numero di numeri primi utilizzando una tabella veloce Look Up (LUT) tecnica del conteggio delle pop a prendere pochissimo tempo per contare come rispetto al secondo circa un utilizzando la tecnica bit-per-bit che usa; tuttavia, in questo codice di esempio è preso anche il tempo che piccolo fuori della porzione temporizzata del codice.
-
Infine, il codice di ottimizzare le operazioni di manipolazione bit per un minimo di codice per ogni ciclo interno. Ad esempio, non utilizza spostamento continuo destra di uno per produrre l'indice dispari e rappresentazione di fatto poco spostamento affatto scrivendo tutti i cicli interni come costante modulo (uguale posizione di bit) operazioni. Come pure, il codice tradotto di Aaron è molto inefficiente in operazioni come ad esempio nel primo abbattimento di Piazza si aggiunge il quadrato del primo all'indice poi controlli per un risultato strano, piuttosto che solo l'aggiunta di due volte il quadrato in modo da non richiedere l'assegno ; quindi rende anche il superfluo controllo spostando il giusto numero di uno (dividendo per due) prima di iniziare l'operazione di abbattimento nel ciclo interno, così come avviene per tutti i loop. Questo codice inefficiente non farà molta differenza nel tempo di esecuzione per grandi campi di utilizzo di questa tecnica "grande setaccio buffer di matrice", come la maggior parte del tempo per operazione viene utilizzato in accesso alla memoria RAM (circa 37 cicli di clock CPU o più di una gamma di un miliardo), ma renderà il tempo di esecuzione molto più lento di quello che deve essere per intervalli più piccoli che si inseriscono le cache CPU; in altre parole si imposta un limite troppo elevata disponibilità nella velocità di esecuzione per operazione.
Il codice è il seguente:
//Sieve of Atkin based on full non page segmented modulo 60 implementation...
using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
namespace NonPagedSoA {
//implements the non-paged Sieve of Atkin (full modulo 60 version)...
class SoA : IEnumerable<ulong> {
private ushort[] buf = null;
private long cnt = 0;
private long opcnt = 0;
private static byte[] modPRMS = { 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59, 61 };
private static ushort[] modLUT;
private static byte[] cntLUT;
//initialize the private LUT's...
static SoA() {
modLUT = new ushort[60];
for (int i = 0, m = 0; i < modLUT.Length; ++i) {
if ((i & 1) != 0 || (i + 7) % 3 == 0 || (i + 7) % 5 == 0) modLUT[i] = 0;
else modLUT[i] = (ushort)(1 << (m++));
}
cntLUT = new byte[65536];
for (int i = 0; i < cntLUT.Length; ++i) {
var c = 0;
for (int j = i; j > 0; j >>= 1) c += j & 1;
cntLUT[i] = (byte)c;
}
}
//initialization and all the work producing the prime bit array done in the constructor...
public SoA(ulong range) {
this.opcnt = 0;
if (range < 7) {
if (range > 1) {
cnt = 1;
if (range > 2) this.cnt += (long)(range - 1) / 2;
}
this.buf = new ushort[0];
}
else {
this.cnt = 3;
var nrng = range - 7; var lmtw = nrng / 60;
//initialize sufficient wheels to non-prime
this.buf = new ushort[lmtw + 1];
//Put in candidate primes:
//for the 4 * x ^ 2 + y ^ 2 quadratic solution toggles - all x odd y...
ulong n = 6; // equivalent to 13 - 7 = 6...
for (uint x = 1, y = 3; n <= nrng; n += (x << 3) + 4, ++x, y = 1) {
var cb = n; if (x <= 1) n -= 8; //cancel the effect of skipping the first one...
for (uint i = 0; i < 15 && cb <= range; cb += (y << 2) + 4, y += 2, ++i) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if (cm != 0)
for (uint c = (uint)cbd, my = y + 15; c < buf.Length; c += my, my += 30) {
buf[c] ^= cm; // ++this.opcnt;
}
}
}
//for the 3 * x ^ 2 + y ^ 2 quadratic solution toggles - x odd y even...
n = 0; // equivalent to 7 - 7 = 0...
for (uint x = 1, y = 2; n <= nrng; n += ((x + x + x) << 2) + 12, x += 2, y = 2) {
var cb = n;
for (var i = 0; i < 15 && cb <= range; cb += (y << 2) + 4, y += 2, ++i) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if (cm != 0)
for (uint c = (uint)cbd, my = y + 15; c < buf.Length; c += my, my += 30) {
buf[c] ^= cm; // ++this.opcnt;
}
}
}
//for the 3 * x ^ 2 - y ^ 2 quadratic solution toggles all x and opposite y = x - 1...
n = 4; // equivalent to 11 - 7 = 4...
for (uint x = 2, y = x - 1; n <= nrng; n += (ulong)(x << 2) + 4, y = x, ++x) {
var cb = n; int i = 0;
for ( ; y > 1 && i < 15 && cb <= nrng; cb += (ulong)(y << 2) - 4, y -= 2, ++i) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if (cm != 0) {
uint c = (uint)cbd, my = y;
for ( ; my >= 30 && c < buf.Length; c += my - 15, my -= 30) {
buf[c] ^= cm; // ++this.opcnt;
}
if (my > 0 && c < buf.Length) { buf[c] ^= cm; /* ++this.opcnt; */ }
}
}
if (y == 1 && i < 15) {
var cbd = cb / 60; var cm = modLUT[cb % 60];
if ((cm & 0x4822) != 0 && cbd < (ulong)buf.Length) { buf[cbd] ^= cm; /* ++this.opcnt; */ }
}
}
//Eliminate squares of base primes, only for those on the wheel:
for (uint i = 0, w = 0, pd = 0, pn = 0, msk = 1; w < this.buf.Length ; ++i) {
uint p = pd + modPRMS[pn];
ulong sqr = (ulong)p * (ulong)p; //to handle ranges above UInt32.MaxValue
if (sqr > range) break;
if ((this.buf[w] & msk) != 0) { //found base prime, square free it...
ulong s = sqr - 7;
for (int j = 0; s <= nrng && j < modPRMS.Length; s = sqr * modPRMS[j] - 7, ++j) {
var cd = s / 60; var cm = (ushort)(modLUT[s % 60] ^ 0xFFFF);
//may need ulong loop index for ranges larger than two billion
//but buf length only good to about 2^31 * 60 = 120 million anyway,
//even with large array setting and half that with 32-bit...
for (ulong c = cd; c < (ulong)this.buf.Length; c += sqr) {
this.buf[c] &= cm; // ++this.opcnt;
}
}
}
if (msk >= 0x8000) { msk = 1; pn = 0; ++w; pd += 60; }
else { msk <<= 1; ++pn; }
}
//clear any overflow primes in the excess space in the last wheel/word:
var ndx = nrng % 60; //clear any primes beyond the range
for (; modLUT[ndx] == 0; --ndx) ;
this.buf[lmtw] &= (ushort)((modLUT[ndx] << 1) - 1);
}
}
//uses a fast pop count Look Up Table to return the total number of primes...
public long Count {
get {
long cnt = this.cnt;
for (int i = 0; i < this.buf.Length; ++i) cnt += cntLUT[this.buf[i]];
return cnt;
}
}
//returns the number of toggle/cull operations used to sieve the prime bit array...
public long Ops {
get {
return this.opcnt;
}
}
//generate the enumeration of primes...
public IEnumerator<ulong> GetEnumerator() {
yield return 2; yield return 3; yield return 5;
ulong pd = 0;
for (uint i = 0, w = 0, pn = 0, msk = 1; w < this.buf.Length; ++i) {
if ((this.buf[w] & msk) != 0) //found a prime bit...
yield return pd + modPRMS[pn]; //add it to the list
if (msk >= 0x8000) { msk = 1; pn = 0; ++w; pd += 60; }
else { msk <<= 1; ++pn; }
}
}
//required for the above enumeration...
IEnumerator IEnumerable.GetEnumerator() {
return this.GetEnumerator();
}
}
class Program {
static void Main(string[] args) {
Console.WriteLine("This program calculates primes by a simple full version of the Sieve of Atkin.\r\n");
const ulong n = 1000000000;
var elpsd = -DateTime.Now.Ticks;
var gen = new SoA(n);
elpsd += DateTime.Now.Ticks;
Console.WriteLine("{0} primes found to {1} using {2} operations in {3} milliseconds.", gen.Count, n, gen.Ops, elpsd / 10000);
//Output prime list for testing...
//Console.WriteLine();
//foreach (var p in gen) {
// Console.Write(p + " ");
//}
//Console.WriteLine();
//Test options showing what one can do with the enumeration, although more slowly...
// Console.WriteLine("\r\nThere are {0} primes with the last one {1} and the sum {2}.",gen.Count(),gen.Last(),gen.Sum(x => (long)x));
Console.Write("\r\nPress any key to exit:");
Console.ReadKey(true);
Console.WriteLine();
}
}
}
Questo codice viene eseguito circa due volte più veloce come codice di Aaron (circa 2,7 secondi con modalità a 64 bit o 32-bit su un i7-2700K (3,5 GHz) con il buffer di circa 16,5 MB e circa 0.258 milioni a levetta combinato / Prime Square operazioni libera da riforma (visualizzabili decommentando il "++ this.opcnt" dichiarazioni) per possibilità vaglio di un miliardo, rispetto a 5,4 / 6,2 secondi (32-bit / 64-bit) per il suo codice senza il conteggio tempo e quasi il doppio l'utilizzo di memoria utilizzando circa 0,359 miliardi operazioni di commutazione / da riforma combinate per setacciatura fino a un miliardo.
Anche se è più veloce rispetto le sue probabilità di sola ingenui più ottimizzata implementazione della Sieve non di paging di Eratostene (SOE), che non rendono il crivello di Atkin più veloce rispetto al crivello di Eratostene , come se si applica tecniche simili a quelle utilizzate nella realizzazione di affidabilità sopra alla RSA Plus utilizza fattorizzazione ruota massima, il SOE circa the stessa velocità di questo.
Analisi: Sebbene il numero di operazioni per il RSA completamente ottimizzata sono circa lo stesso del numero di operazioni per il SOA per una gamma vaglio di un miliardo, il collo di bottiglia principale per questi non di paging implementazioni è accesso alla memoria una volta che la dimensione del buffer vaglio supera le dimensioni della cache CPU (32 kilobyte di cache L1 ad un accesso di clock, 256 di cache Kilobyte L2 a circa quattro cicli di clock tempo di accesso e 8 megabyte di cache L3 a circa 20 cicli di clock tempo di accesso per il mio i7), dopo di che l'accesso alla memoria può superare un centinaio di cicli di clock.
Ora entrambi hanno un fattore di circa otto miglioramento nella velocità di accesso di memoria quando si adegua gli algoritmi alla pagina segmentazione così si può setacciare intervalli che altrimenti non inserirsi nella memoria disponibile. Tuttavia, SOE continua a guadagnare oltre la dichiarazione di affidabilità come la gamma vaglio inizia a diventare molto grande a causa di difficoltà di attuazione del "numeri primi quadrato libero" parte dell'algoritmo a causa delle enormi passi avanti nella abbattimento scansioni che crescono rapidamente a molte centinaia di volte la dimensione dei buffer di pagina. Come pure, e forse più grave, diventa molto memoria e / o intensità di calcolo per calcolare il nuovo punto di partenza per ogni valore di 'x' per il valore di 'y' alla rappresentazione più basso di ciascun buffer di pagina per un ulteriore piuttosto grande perdita di efficienza del paging di affidabilità comparaed alla RSA come la gamma cresce.
EDIT_ADD: Le probabilità per soli SoE come usati da Aaron Murgatroyd utilizza circa 1,026 miliardi le operazioni da riforma per una serie vaglio di un miliardo in modo da circa quattro volte il numero di operazioni come la SOA e quindi dovrebbe funzionare su quattro volte più lenta, ma la dichiarazione di affidabilità anche se qui implementata ha un più complesso ciclo interno e soprattutto a causa di una proporzione molto maggiore delle quote soli abbattimenti SOE hanno un passo molto più breve nelle scansioni abbattimento rispetto ai progressi della SOA odds naive -solo SoE ha molto migliori tempi di accesso alla memoria media nonostante buffer vaglio supera notevolmente le dimensioni della cache CPU (migliore utilizzo delle cache di associatività). Questo spiega il motivo per cui quanto sopra SOA è solo circa due volte più veloce le probabilità per soli SoE anche se sarebbe teoricamente sembrano fare solo un quarto del lavoro.
Se si dovesse utilizzare un algoritmo simile utilizzando costanti cicli interni rollover come per SOA alto e attuato la stessa fattorizzazione 2/3/5 ruota, il RSA potrebbe ridurre il numero di operazioni di riforma a circa 0,405 miliardi di operazioni così solo 50% più operazioni rispetto alla SOA e sarebbe teoricamente funzionare solo leggermente più lento rispetto alla dichiarazione di affidabilità, ma possono funzionare a circa la stessa velocità a causa dei passi da riforma ancora essere un po 'più piccola per la dichiarazione di affidabilità sulla media di questo grande memoria "naive" uso dei buffer. Aumentando la fattorizzazione ruota per ruota 2/3/5/7: le operazioni SoE da riforma sono ridotti a circa 0,314 per una serie abbattimento di un miliardo e possono rendere quella versione della corsa RSA circa la stessa velocità per questo algoritmo.
ulteriore utilizzo di fattorizzazione ruota può essere effettuato tramite pre-culling matrice setaccio (copia in un modello) per i 2/3/5/7/11/13/17/19 fattori primi praticamente senza alcuna perdita nel tempo di esecuzione per ridurre il numero totale di operazioni di riforma a circa 0,251 miliardi una gamma vaglio di un miliardo e RSA sarà più veloce o circa la stessa velocità che anche questa versione ottimizzata di SOA, anche per queste grandi versioni di memoria buffer, con il SoE continua ad avere molto meno complessità del codice di cui sopra.
Così, si può vedere che il numero di operazioni per il RSA può essere notevolmente ridotto da un anche quote sola o 2/3/5 versione fattorizzazione ruota naive o tale che il numero di operazioni è circa lo stesso come per la dichiarazione di affidabilità mentre allo stesso tempo il tempo per operazione può effettivamente essere meno dovuto sia cicli interni meno complessi e accesso alla memoria più efficiente. END_EDIT_ADD
EDIT_ADD2: qui aggiungoil codice per un RSA usando una tecnica simile posizione costante modulo / bit per i cicli più interni come per la dichiarazione di affidabilità sopra in base al codice pseudo più in basso la risposta, come legato soprattutto . Il codice è un po 'meno complessa SOA sopra, pur avendo alta fattorizzazione ruota e pre-abbattimento applicati tale che il numero totale di operazioni di riforma sono in realtà meno operazioni snodata / riforma combinati per la dichiarazione di affidabilità fino ad un rang setacciatura di circa due miliardi. Il codice nel seguente modo:
EDIT_FINAL codice corretto qui sotto e commenti relativi ad esso END_EDIT_FINAL
//Sieve of Eratosthenes based on maximum wheel factorization and pre-culling implementation...
using System;
using System.Collections;
using System.Collections.Generic;
namespace NonPagedSoE {
//implements the non-paged Sieve of Eratosthenes (full modulo 210 version with preculling)...
class SoE : IEnumerable<ulong> {
private ushort[] buf = null;
private long cnt = 0;
private long opcnt = 0;
private static byte[] basePRMS = { 2, 3, 5, 7, 11, 13, 17, 19 };
private static byte[] modPRMS = { 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, //positions + 23
97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143, 149, 151, 157, 163,
167, 169, 173, 179, 181 ,187 ,191 ,193, 197, 199, 209, 211, 221, 223, 227, 229 };
private static byte[] gapsPRMS = { 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, 2, 4, 2, 4 };
private static ulong[] modLUT;
private static byte[] cntLUT;
//initialize the private LUT's...
static SoE() {
modLUT = new ulong[210];
for (int i = 0, m = 0; i < modLUT.Length; ++i) {
if ((i & 1) != 0 || (i + 23) % 3 == 0 || (i + 23) % 5 == 0 || (i + 23) % 7 == 0) modLUT[i] = 0;
else modLUT[i] = 1UL << (m++);
}
cntLUT = new byte[65536];
for (int i = 0; i < cntLUT.Length; ++i) {
var c = 0;
for (int j = i ^ 0xFFFF; j > 0; j >>= 1) c += j & 1; //reverse logic; 0 is prime; 1 is composite
cntLUT[i] = (byte)c;
}
}
//initialization and all the work producing the prime bit array done in the constructor...
public SoE(ulong range) {
this.opcnt = 0;
if (range < 23) {
if (range > 1) {
for (int i = 0; i < modPRMS.Length; ++i) if (modPRMS[i] <= range) this.cnt++; else break;
}
this.buf = new ushort[0];
}
else {
this.cnt = 8;
var nrng = range - 23; var lmtw = nrng / 210; var lmtwt3 = lmtw * 3;
//initialize sufficient wheels to prime
this.buf = new ushort[lmtwt3 + 3]; //initial state of all zero's is all potential prime.
//initialize array to account for preculling the primes of 11, 13, 17, and 19;
//(2, 3, 5, and 7 already eliminated by the bit packing to residues).
for (int pn = modPRMS.Length - 4; pn < modPRMS.Length; ++pn) {
uint p = modPRMS[pn] - 210u; ulong pt3 = p * 3;
ulong s = p * p - 23;
ulong xrng = Math.Min(9699709, nrng); // only do for the repeating master pattern size
ulong nwrds = (ulong)Math.Min(138567, this.buf.Length);
for (int j = 0; s <= xrng && j < modPRMS.Length; s += p * gapsPRMS[(pn + j++) % 48]) {
var sm = modLUT[s % 210];
var si = (sm < (1UL << 16)) ? 0UL : ((sm < (1UL << 32)) ? 1UL : 2UL);
var cd = s / 210 * 3 + si; var cm = (ushort)(sm >> (int)(si << 4));
for (ulong c = cd; c < nwrds; c += pt3) { //tight culling loop for size of master pattern
this.buf[c] |= cm; // ++this.opcnt; //reverse logic; mark composites with ones.
}
}
}
//Now copy the master pattern so it repeats across the main buffer, allow for overflow...
for (long i = 138567; i < this.buf.Length; i += 138567)
if (i + 138567 <= this.buf.Length)
Array.Copy(this.buf, 0, this.buf, i, 138567);
else Array.Copy(this.buf, 0, this.buf, i, this.buf.Length - i);
//Eliminate all composites which are factors of base primes, only for those on the wheel:
for (uint i = 0, w = 0, wi = 0, pd = 0, pn = 0, msk = 1; w < this.buf.Length; ++i) {
uint p = pd + modPRMS[pn];
ulong sqr = (ulong)p * (ulong)p;
if (sqr > range) break;
if ((this.buf[w] & msk) == 0) { //found base prime, mark its composites...
ulong s = sqr - 23; ulong pt3 = p * 3;
for (int j = 0; s <= nrng && j < modPRMS.Length; s += p * gapsPRMS[(pn + j++) % 48]) {
var sm = modLUT[s % 210];
var si = (sm < (1UL << 16)) ? 0UL : ((sm < (1UL << 32)) ? 1UL : 2UL);
var cd = s / 210 * 3 + si; var cm = (ushort)(sm >> (int)(si << 4));
for (ulong c = cd; c < (ulong)this.buf.Length; c += pt3) { //tight culling loop
this.buf[c] |= cm; // ++this.opcnt; //reverse logic; mark composites with ones.
}
}
}
++pn;
if (msk >= 0x8000) { msk = 1; ++w; ++wi; if (wi == 3) { wi = 0; pn = 0; pd += 210; } }
else msk <<= 1;
}
//clear any overflow primes in the excess space in the last wheel/word:
var ndx = nrng % 210; //clear any primes beyond the range
for (; modLUT[ndx] == 0; --ndx) ;
var cmsk = (~(modLUT[ndx] - 1)) << 1; //force all bits above to be composite ones.
this.buf[lmtwt3++] |= (ushort)cmsk;
this.buf[lmtwt3++] |= (ushort)(cmsk >> 16);
this.buf[lmtwt3] |= (ushort)(cmsk >> 32);
}
}
//uses a fast pop count Look Up Table to return the total number of primes...
public long Count {
get {
long cnt = this.cnt;
for (int i = 0; i < this.buf.Length; ++i) cnt += cntLUT[this.buf[i]];
return cnt;
}
}
//returns the number of cull operations used to sieve the prime bit array...
public long Ops {
get {
return this.opcnt;
}
}
//generate the enumeration of primes...
public IEnumerator<ulong> GetEnumerator() {
yield return 2; yield return 3; yield return 5; yield return 7;
yield return 11; yield return 13; yield return 17; yield return 19;
ulong pd = 0;
for (uint i = 0, w = 0, wi = 0, pn = 0, msk = 1; w < this.buf.Length; ++i) {
if ((this.buf[w] & msk) == 0) //found a prime bit...
yield return pd + modPRMS[pn];
++pn;
if (msk >= 0x8000) { msk = 1; ++w; ++wi; if (wi == 3) { wi = 0; pn = 0; pd += 210; } }
else msk <<= 1;
}
}
//required for the above enumeration...
IEnumerator IEnumerable.GetEnumerator() {
return this.GetEnumerator();
}
}
class Program {
static void Main(string[] args) {
Console.WriteLine("This program calculates primes by a simple maximually wheel factorized version of the Sieve of Eratosthenes.\r\n");
const ulong n = 1000000000;
var elpsd = -DateTime.Now.Ticks;
var gen = new SoE(n);
elpsd += DateTime.Now.Ticks;
Console.WriteLine("{0} primes found to {1} using {2} operations in {3} milliseconds.", gen.Count, n, gen.Ops, elpsd / 10000);
// Console.WriteLine();
// foreach (var p in gen) {
// Console.Write(p + " ");
// }
// Console.WriteLine();
// Console.WriteLine("\r\nThere are {0} primes with the last one {1} and the sum {2}.",gen.Count(),gen.Last(),gen.Sum(x => (long)x));
Console.Write("\r\nPress any key to exit:");
Console.ReadKey(true);
Console.WriteLine();
}
}
}
Questo codice viene eseguito in realtà una piccola percentuale più veloce di SOA di cui sopra come dovrebbe in quanto vi sono un po 'meno le operazioni e il collo di bottiglia principale per questa grande dimensione della matrice per una gamma di un miliardo è il tempo di accesso alla memoria di qualcosa come 40 a oltre 100 cicli di clock della CPU a seconda delle specifiche della CPU e della memoria; questo significa che ottimizzazioni del codice (diversi dalla riduzione del numero totale di operazioni) sono inefficaci come la maggior parte del tempo è spendere attesa di accesso alla memoria. In ogni caso, utilizzando un buffer di memoria enorme non è il modo più efficiente per setacciare grandi intervalli, con un fattore fino a circa otto volte miglioramento per il RSA utilizzando la pagina segmentazione con la stessa fattorizzazione massimo della ruota (che costituisce anche la base per multi-processing).
È nell'attuazione pagina segmentazione e multi-processing che la SOA è davvero carente per intervalli molto al di sopra quattro miliardi rispetto al RSA come gli utili a causa della ridotta complessità asintotica del SoA rapidamente mangiati dalla pagina elaborazione overhead fattori legati al primo trattamento di piazza e calcolando il numero molto maggiore di indirizzi iniziali pagina; in alternativa, si supera questo memorizzando marcatori nella memoria RAM ad un costo enorme consumo di memoria ed ulteriori inefficienze nell'accesso queste strutture deposito marcatore. END_EDIT_ADD2
In breve, la dichiarazione di affidabilità non è realmente un setaccio pratico rispetto al la ruota completamente fattorizzata SoE poiché come l'aumento di complessità asintotica comincia ad avvicinarla prestazioni al RSA completamente ottimizzata, comincia a perdere efficienza causa i particolari di realizzazione pratica come in volta accesso alla memoria e relativa complessità pagina segmentazione nonché essendo generalmente più complesso e difficile da scrivere. A mio parere è più di un concetto intellettuale interessante e l'esercizio mentale di un setaccio pratico rispetto alla RSA.
Un giorno mi adatterò queste tecniche per una pagina multi-threaded segmentato Crivello di Eratostene di essere circa veloce in C # come Atkin e Bernstein di implementazione "PrimeGen" della dichiarazione di affidabilità in 'C' e vi soffiare fuori dall'acqua per le grandi catene di sopra di circa quattro miliardi anche a thread singolo, con una spinta in più in termini di velocità di fino a circa quattro anni quando multi-threading sul mio i7 (otto core compreso Hyper Threading).
Ecco un'altra implementazione. Esso utilizza BitArray per risparmiare memoria. Il Parallel.For ha bisogno di .NET Framework 4.
static List<int> FindPrimesBySieveOfAtkins(int max)
{
// var isPrime = new BitArray((int)max+1, false);
// Can't use BitArray because of threading issues.
var isPrime = new bool[max + 1];
var sqrt = (int)Math.Sqrt(max);
Parallel.For(1, sqrt, x =>
{
var xx = x * x;
for (int y = 1; y <= sqrt; y++)
{
var yy = y * y;
var n = 4 * xx + yy;
if (n <= max && (n % 12 == 1 || n % 12 == 5))
isPrime[n] ^= true;
n = 3 * xx + yy;
if (n <= max && n % 12 == 7)
isPrime[n] ^= true;
n = 3 * xx - yy;
if (x > y && n <= max && n % 12 == 11)
isPrime[n] ^= true;
}
});
var primes = new List<int>() { 2, 3 };
for (int n = 5; n <= sqrt; n++)
{
if (isPrime[n])
{
primes.Add(n);
int nn = n * n;
for (int k = nn; k <= max; k += nn)
isPrime[k] = false;
}
}
for (int n = sqrt + 1; n <= max; n++)
if (isPrime[n])
primes.Add(n);
return primes;
}
Ecco una più rapida attuazione del crivello di Atkin, ho rubato l'algoritmo da questo script Python qui (mi prendo nessun credito per l'algoritmo):
http://programmingpraxis.com/2010/02/ 19 / setaccio-of-Atkin-migliorato /
using System;
using System.Collections;
using System.Collections.Generic;
namespace PrimeGenerator
{
// The block element type for the bit array,
// use any unsigned value. WARNING: UInt64 is
// slower even on x64 architectures.
using BitArrayType = System.UInt32;
// This should never be any bigger than 256 bits - leave as is.
using BitsPerBlockType = System.Byte;
// The prime data type, this can be any unsigned value, the limit
// of this type determines the limit of Prime value that can be
// found. WARNING: UInt64 is slower even on x64 architectures.
using PrimeType = System.Int32;
/// <summary>
/// Calculates prime number using the Sieve of Eratosthenes method.
/// </summary>
/// <example>
/// <code>
/// var lpPrimes = new Eratosthenes(1e7);
/// foreach (UInt32 luiPrime in lpPrimes)
/// Console.WriteLine(luiPrime);
/// </example>
public class Atkin : IEnumerable<PrimeType>
{
#region Constants
/// <summary>
/// Constant for number of bits per block, calculated based on size of BitArrayType.
/// </summary>
const BitsPerBlockType cbBitsPerBlock = sizeof(BitArrayType) * 8;
#endregion
#region Protected Locals
/// <summary>
/// The limit for the maximum prime value to find.
/// </summary>
protected readonly PrimeType mpLimit;
/// <summary>
/// The number of primes calculated or null if not calculated yet.
/// </summary>
protected PrimeType? mpCount = null;
/// <summary>
/// The current bit array where a set bit means
/// the odd value at that location has been determined
/// to not be prime.
/// </summary>
protected BitArrayType[] mbaOddPrime;
#endregion
#region Initialisation
/// <summary>
/// Create Sieve of Atkin generator.
/// </summary>
/// <param name="limit">The limit for the maximum prime value to find.</param>
public Atkin(PrimeType limit)
{
// Check limit range
if (limit > PrimeType.MaxValue - (PrimeType)Math.Sqrt(PrimeType.MaxValue))
throw new ArgumentOutOfRangeException();
mpLimit = limit;
FindPrimes();
}
#endregion
#region Private Methods
/// <summary>
/// Finds the prime number within range.
/// </summary>
private unsafe void FindPrimes()
{
// Allocate bit array.
mbaOddPrime = new BitArrayType[(((mpLimit >> 1) + 1) / cbBitsPerBlock) + 1];
PrimeType lpYLimit, lpN, lpXX3, lpXX4, lpDXX, lpDN, lpDXX4, lpXX, lpX, lpYY, lpMinY, lpS, lpK;
fixed (BitArrayType* lpbOddPrime = &mbaOddPrime[0])
{
// n = 3x^2 + y^2 section
lpXX3 = 3;
for (lpDXX = 0; lpDXX < 12 * SQRT((mpLimit - 1) / 3); lpDXX += 24)
{
lpXX3 += lpDXX;
lpYLimit = (12 * SQRT(mpLimit - lpXX3)) - 36;
lpN = lpXX3 + 16;
for (lpDN = -12; lpDN < lpYLimit + 1; lpDN += 72)
{
lpN += lpDN;
lpbOddPrime[(lpN >> 1) / cbBitsPerBlock] ^=
(BitArrayType)((BitArrayType)1 << (int)((lpN >> 1) % cbBitsPerBlock));
}
lpN = lpXX3 + 4;
for (lpDN = 12; lpDN < lpYLimit + 1; lpDN += 72)
{
lpN += lpDN;
lpbOddPrime[(lpN >> 1) / cbBitsPerBlock] ^=
(BitArrayType)((BitArrayType)1 << (int)((lpN >> 1) % cbBitsPerBlock));
}
}
// # n = 4x^2 + y^2 section
lpXX4 = 0;
for (lpDXX4 = 4; lpDXX4 < 8 * SQRT((mpLimit - 1) / 4) + 4; lpDXX4 += 8)
{
lpXX4 += lpDXX4;
lpN = lpXX4 + 1;
if ((lpXX4 % 3) != 0)
{
for (lpDN = 0; lpDN < (4 * SQRT(mpLimit - lpXX4)) - 3; lpDN += 8)
{
lpN += lpDN;
lpbOddPrime[(lpN >> 1) / cbBitsPerBlock] ^=
(BitArrayType)((BitArrayType)1 << (int)((lpN >> 1) % cbBitsPerBlock));
}
}
else
{
lpYLimit = (12 * SQRT(mpLimit - lpXX4)) - 36;
lpN = lpXX4 + 25;
for (lpDN = -24; lpDN < lpYLimit + 1; lpDN += 72)
{
lpN += lpDN;
lpbOddPrime[(lpN >> 1) / cbBitsPerBlock] ^=
(BitArrayType)((BitArrayType)1 << (int)((lpN >> 1) % cbBitsPerBlock));
}
lpN = lpXX4 + 1;
for (lpDN = 24; lpDN < lpYLimit + 1; lpDN += 72)
{
lpN += lpDN;
lpbOddPrime[(lpN >> 1) / cbBitsPerBlock] ^=
(BitArrayType)((BitArrayType)1 << (int)((lpN >> 1) % cbBitsPerBlock));
}
}
}
// # n = 3x^2 - y^2 section
lpXX = 1;
for (lpX = 3; lpX < SQRT(mpLimit / 2) + 1; lpX += 2)
{
lpXX += 4 * lpX - 4;
lpN = 3 * lpXX;
if (lpN > mpLimit)
{
lpMinY = ((SQRT(lpN - mpLimit) >> 2) << 2);
lpYY = lpMinY * lpMinY;
lpN -= lpYY;
lpS = 4 * lpMinY + 4;
}
else
lpS = 4;
for (lpDN = lpS; lpDN < 4 * lpX; lpDN += 8)
{
lpN -= lpDN;
if (lpN <= mpLimit && lpN % 12 == 11)
lpbOddPrime[(lpN >> 1) / cbBitsPerBlock] ^=
(BitArrayType)((BitArrayType)1 << (int)((lpN >> 1) % cbBitsPerBlock));
}
}
// xx = 0
lpXX = 0;
for (lpX = 2; lpX < SQRT(mpLimit / 2) + 1; lpX += 2)
{
lpXX += 4*lpX - 4;
lpN = 3*lpXX;
if (lpN > mpLimit)
{
lpMinY = ((SQRT(lpN - mpLimit) >> 2) << 2) - 1;
lpYY = lpMinY * lpMinY;
lpN -= lpYY;
lpS = 4*lpMinY + 4;
}
else
{
lpN -= 1;
lpS = 0;
}
for (lpDN = lpS; lpDN < 4 * lpX; lpDN += 8)
{
lpN -= lpDN;
if (lpN <= mpLimit && lpN % 12 == 11)
lpbOddPrime[(lpN>>1) / cbBitsPerBlock] ^=
(BitArrayType)((BitArrayType)1 << (int)((lpN>>1) % cbBitsPerBlock));
}
}
// # eliminate squares
for (lpN = 5; lpN < SQRT(mpLimit) + 1; lpN += 2)
if ((lpbOddPrime[(lpN >> 1) / cbBitsPerBlock] & ((BitArrayType)1 << (int)((lpN >> 1) % cbBitsPerBlock))) != 0)
for (lpK = lpN * lpN; lpK < mpLimit; lpK += lpN * lpN)
if ((lpK & 1) == 1)
lpbOddPrime[(lpK >> 1) / cbBitsPerBlock] &=
(BitArrayType)~((BitArrayType)1 << (int)((lpK >> 1) % cbBitsPerBlock));
}
}
/// <summary>
/// Calculates the truncated square root for a number.
/// </summary>
/// <param name="value">The value to get the square root for.</param>
/// <returns>The truncated sqrt of the value.</returns>
private unsafe PrimeType SQRT(PrimeType value)
{
return (PrimeType)Math.Sqrt(value);
}
/// <summary>
/// Gets a bit value by index.
/// </summary>
/// <param name="bits">The blocks containing the bits.</param>
/// <param name="index">The index of the bit.</param>
/// <returns>True if bit is set, false if cleared.</returns>
private bool GetBitSafe(BitArrayType[] bits, PrimeType index)
{
if ((index & 1) == 1)
return (bits[(index >> 1) / cbBitsPerBlock] & ((BitArrayType)1 << (int)((index >> 1) % cbBitsPerBlock))) != 0;
else
return false;
}
#endregion
#region Public Properties
/// <summary>
/// Get the limit for the maximum prime value to find.
/// </summary>
public PrimeType Limit
{
get
{
return mpLimit;
}
}
/// <summary>
/// Returns the number of primes found in the range.
/// </summary>
public PrimeType Count
{
get
{
if (!mpCount.HasValue)
{
PrimeType lpCount = 0;
foreach (PrimeType liPrime in this) lpCount++;
mpCount = lpCount;
}
return mpCount.Value;
}
}
/// <summary>
/// Determines if a value in range is prime or not.
/// </summary>
/// <param name="test">The value to test for primality.</param>
/// <returns>True if the value is prime, false otherwise.</returns>
public bool this[PrimeType test]
{
get
{
if (test > mpLimit) throw new ArgumentOutOfRangeException();
if (test <= 1) return false;
if (test == 2) return true;
if ((test & 1) == 0) return false;
return !GetBitSafe(mbaOddPrime, test >> 1);
}
}
#endregion
#region Public Methods
/// <summary>
/// Gets the enumerator for the primes.
/// </summary>
/// <returns>The enumerator of the primes.</returns>
public IEnumerator<PrimeType> GetEnumerator()
{
// return [2,3] + filter(primes.__getitem__, xrange(5,limit,2))
// Two & Three always prime.
yield return 2;
yield return 3;
// Start at first block, third MSB (5).
int liBlock = 0;
byte lbBit = 2;
BitArrayType lbaCurrent = mbaOddPrime[0] >> lbBit;
// For each value in range stepping in incrments of two for odd values.
for (PrimeType lpN = 5; lpN <= mpLimit; lpN += 2)
{
// If current bit not set then value is prime.
if ((lbaCurrent & 1) == 1)
yield return lpN;
// Move to NSB.
lbaCurrent >>= 1;
// Increment bit value.
lbBit++;
// If block is finished.
if (lbBit == cbBitsPerBlock)
{
lbBit = 0;
lbaCurrent = mbaOddPrime[++liBlock];
//// Move to first bit of next block skipping full blocks.
while (lbaCurrent == 0)
{
lpN += ((PrimeType)cbBitsPerBlock) << 1;
if (lpN <= mpLimit)
lbaCurrent = mbaOddPrime[++liBlock];
else
break;
}
}
}
}
#endregion
#region IEnumerable<PrimeType> Implementation
/// <summary>
/// Gets the enumerator for the primes.
/// </summary>
/// <returns></returns>
IEnumerator IEnumerable.GetEnumerator()
{
return GetEnumerator();
}
#endregion
}
}
Il suo vicino di velocità alla mia versione più ottimizzata del crivello di Eratostene, ma la sua ancora più lento di circa il 20%, si può trovare qui:
Heres mio, utilizza una classe denominata CompartmentalisedParallel che consente di eseguire in parallelo cicli for ma di controllare il numero di fili in modo che gli indici sono raggruppati in su. Tuttavia, a causa dei problemi di threading che è necessario sia bloccare il BitArray ogni volta che viene alterato o creare un BitArray separato per ciascun thread e poi XOR insieme alla fine, la prima opzione è stato piuttosto lento perché la quantità di serrature, il secondo opzione sembrava più veloce per me!
using System;
using System.Collections;
using System.Collections.Generic;
using System.Threading.Tasks;
namespace PrimeGenerator
{
public class Atkin : Primes
{
protected BitArray mbaPrimes;
protected bool mbThreaded = true;
public Atkin(int limit)
: this(limit, true)
{
}
public Atkin(int limit, bool threaded)
: base(limit)
{
mbThreaded = threaded;
if (mbaPrimes == null) FindPrimes();
}
public bool Threaded
{
get
{
return mbThreaded;
}
}
public override IEnumerator<int> GetEnumerator()
{
yield return 2;
yield return 3;
for (int lsN = 5; lsN <= msLimit; lsN += 2)
if (mbaPrimes[lsN]) yield return lsN;
}
private void FindPrimes()
{
mbaPrimes = new BitArray(msLimit + 1, false);
int lsSQRT = (int)Math.Sqrt(msLimit);
int[] lsSquares = new int[lsSQRT + 1];
for (int lsN = 0; lsN <= lsSQRT; lsN++)
lsSquares[lsN] = lsN * lsN;
if (Threaded)
{
CompartmentalisedParallel.For<BitArray>(
1, lsSQRT + 1, new ParallelOptions(),
(start, finish) => { return new BitArray(msLimit + 1, false); },
(lsX, lsState, lbaLocal) =>
{
int lsX2 = lsSquares[lsX];
for (int lsY = 1; lsY <= lsSQRT; lsY++)
{
int lsY2 = lsSquares[lsY];
int lsN = 4 * lsX2 + lsY2;
if (lsN <= msLimit && (lsN % 12 == 1 || lsN % 12 == 5))
lbaLocal[lsN] ^= true;
lsN -= lsX2;
if (lsN <= msLimit && lsN % 12 == 7)
lbaLocal[lsN] ^= true;
if (lsX > lsY)
{
lsN -= lsY2 * 2;
if (lsN <= msLimit && lsN % 12 == 11)
lbaLocal[lsN] ^= true;
}
}
return lbaLocal;
},
(lbaResult, start, finish) =>
{
lock (mbaPrimes)
mbaPrimes.Xor(lbaResult);
},
-1
);
}
else
{
for (int lsX = 1; lsX <= lsSQRT; lsX++)
{
int lsX2 = lsSquares[lsX];
for (int lsY = 1; lsY <= lsSQRT; lsY++)
{
int lsY2 = lsSquares[lsY];
int lsN = 4 * lsX2 + lsY2;
if (lsN <= msLimit && (lsN % 12 == 1 || lsN % 12 == 5))
mbaPrimes[lsN] ^= true;
lsN -= lsX2;
if (lsN <= msLimit && lsN % 12 == 7)
mbaPrimes[lsN] ^= true;
if (lsX > lsY)
{
lsN -= lsY2 * 2;
if (lsN <= msLimit && lsN % 12 == 11)
mbaPrimes[lsN] ^= true;
}
}
}
}
for (int lsN = 5; lsN < lsSQRT; lsN += 2)
if (mbaPrimes[lsN])
{
var lsS = lsSquares[lsN];
for (int lsK = lsS; lsK <= msLimit; lsK += lsS)
mbaPrimes[lsK] = false;
}
}
}
}
E la classe CompartmentalisedParallel:
using System;
using System.Threading.Tasks;
namespace PrimeGenerator
{
public static class CompartmentalisedParallel
{
#region Int
private static int[] CalculateCompartments(int startInclusive, int endExclusive, ref int threads)
{
if (threads == 0) threads = 1;
if (threads == -1) threads = Environment.ProcessorCount;
if (threads > endExclusive - startInclusive) threads = endExclusive - startInclusive;
int[] liThreadIndexes = new int[threads + 1];
liThreadIndexes[threads] = endExclusive - 1;
int liIndexesPerThread = (endExclusive - startInclusive) / threads;
for (int liCount = 0; liCount < threads; liCount++)
liThreadIndexes[liCount] = liCount * liIndexesPerThread;
return liThreadIndexes;
}
public static void For<TLocal>(
int startInclusive, int endExclusive,
ParallelOptions parallelOptions,
Func<int, int, TLocal> localInit,
Func<int, ParallelLoopState, TLocal, TLocal> body,
Action<TLocal, int, int> localFinally,
int threads)
{
int[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
if (threads > 1)
Parallel.For(
0, threads, parallelOptions,
(liThread, lsState) =>
{
TLocal llLocal = localInit(liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
for (int liCounter = liThreadIndexes[liThread]; liCounter < liThreadIndexes[liThread + 1]; liCounter++)
body(liCounter, lsState, llLocal);
localFinally(llLocal, liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
}
);
else
{
TLocal llLocal = localInit(startInclusive, endExclusive);
for (int liCounter = startInclusive; liCounter < endExclusive; liCounter++)
body(liCounter, null, llLocal);
localFinally(llLocal, startInclusive, endExclusive);
}
}
public static void For(
int startInclusive, int endExclusive,
ParallelOptions parallelOptions,
Action<int, ParallelLoopState> body,
int threads)
{
int[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
if (threads > 1)
Parallel.For(
0, threads, parallelOptions,
(liThread, lsState) =>
{
for (int liCounter = liThreadIndexes[liThread]; liCounter < liThreadIndexes[liThread + 1]; liCounter++)
body(liCounter, lsState);
}
);
else
for (int liCounter = startInclusive; liCounter < endExclusive; liCounter++)
body(liCounter, null);
}
public static void For(
int startInclusive, int endExclusive,
ParallelOptions parallelOptions,
Action<int> body,
int threads)
{
int[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
if (threads > 1)
Parallel.For(
0, threads, parallelOptions,
(liThread) =>
{
for (int liCounter = liThreadIndexes[liThread]; liCounter < liThreadIndexes[liThread + 1]; liCounter++)
body(liCounter);
}
);
else
for (int liCounter = startInclusive; liCounter < endExclusive; liCounter++)
body(liCounter);
}
public static void For(
int startInclusive, int endExclusive,
Action<int, ParallelLoopState> body,
int threads)
{
For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
}
public static void For(
int startInclusive, int endExclusive,
Action<int> body,
int threads)
{
For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
}
public static void For<TLocal>(
int startInclusive, int endExclusive,
Func<int, int, TLocal> localInit,
Func<int, ParallelLoopState, TLocal, TLocal> body,
Action<TLocal, int, int> localFinally,
int threads)
{
For<TLocal>(startInclusive, endExclusive, new ParallelOptions(), localInit, body, localFinally, threads);
}
#endregion
#region Long
private static long[] CalculateCompartments(long startInclusive, long endExclusive, ref long threads)
{
if (threads == 0) threads = 1;
if (threads == -1) threads = Environment.ProcessorCount;
if (threads > endExclusive - startInclusive) threads = endExclusive - startInclusive;
long[] liThreadIndexes = new long[threads + 1];
liThreadIndexes[threads] = endExclusive - 1;
long liIndexesPerThread = (endExclusive - startInclusive) / threads;
for (long liCount = 0; liCount < threads; liCount++)
liThreadIndexes[liCount] = liCount * liIndexesPerThread;
return liThreadIndexes;
}
public static void For<TLocal>(
long startInclusive, long endExclusive,
ParallelOptions parallelOptions,
Func<long, long, TLocal> localInit,
Func<long, ParallelLoopState, TLocal, TLocal> body,
Action<TLocal, long, long> localFinally,
long threads)
{
long[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
if (threads > 1)
Parallel.For(
0, threads, parallelOptions,
(liThread, lsState) =>
{
TLocal llLocal = localInit(liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
for (long liCounter = liThreadIndexes[liThread]; liCounter < liThreadIndexes[liThread + 1]; liCounter++)
body(liCounter, lsState, llLocal);
localFinally(llLocal, liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
}
);
else
{
TLocal llLocal = localInit(startInclusive, endExclusive);
for (long liCounter = startInclusive; liCounter < endExclusive; liCounter++)
body(liCounter, null, llLocal);
localFinally(llLocal, startInclusive, endExclusive);
}
}
public static void For(
long startInclusive, long endExclusive,
ParallelOptions parallelOptions,
Action<long, ParallelLoopState> body,
long threads)
{
long[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
if (threads > 1)
Parallel.For(
0, threads, parallelOptions,
(liThread, lsState) =>
{
for (long liCounter = liThreadIndexes[liThread]; liCounter < liThreadIndexes[liThread + 1]; liCounter++)
body(liCounter, lsState);
}
);
else
for (long liCounter = startInclusive; liCounter < endExclusive; liCounter++)
body(liCounter, null);
}
public static void For(
long startInclusive, long endExclusive,
ParallelOptions parallelOptions,
Action<long> body,
long threads)
{
long[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
if (threads > 1)
Parallel.For(
0, threads, parallelOptions,
(liThread) =>
{
for (long liCounter = liThreadIndexes[liThread]; liCounter < liThreadIndexes[liThread + 1]; liCounter++)
body(liCounter);
}
);
else
for (long liCounter = startInclusive; liCounter < endExclusive; liCounter++)
body(liCounter);
}
public static void For(
long startInclusive, long endExclusive,
Action<long, ParallelLoopState> body,
long threads)
{
For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
}
public static void For(
long startInclusive, long endExclusive,
Action<long> body,
long threads)
{
For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
}
public static void For<TLocal>(
long startInclusive, long endExclusive,
Func<long, long, TLocal> localInit,
Func<long, ParallelLoopState, TLocal, TLocal> body,
Action<TLocal, long, long> localFinally,
long threads)
{
For<TLocal>(startInclusive, endExclusive, new ParallelOptions(), localInit, body, localFinally, threads);
}
#endregion
}
}
classe di numeri primi di base:
using System.Collections;
using System.Collections.Generic;
namespace PrimeGenerator
{
public abstract class Primes : IEnumerable<int>
{
protected readonly int msLimit;
public Primes(int limit)
{
msLimit = limit;
}
public int Limit
{
get
{
return msLimit;
}
}
public int Count
{
get
{
int liCount = 0;
foreach (int liPrime in this)
liCount++;
return liCount;
}
}
public abstract IEnumerator<int> GetEnumerator();
IEnumerator IEnumerable.GetEnumerator()
{
return GetEnumerator();
}
}
}
Si usa nel modo seguente:
var lpPrimes = new Atkin(count, true);
Console.WriteLine(lpPrimes.Count);
Console.WriteLine(s.ElapsedMilliseconds);
Tuttavia, ho trovato le Eratostene di essere più veloce in ogni caso, anche con una CPU quattro core funziona in modalità multithread per l'Atkin:
using System;
using System.Collections;
using System.Collections.Generic;
namespace PrimeGenerator
{
public class Eratosthenes : Primes
{
protected BitArray mbaOddEliminated;
public Eratosthenes(int limit)
: base(limit)
{
if (mbaOddEliminated == null) FindPrimes();
}
public override IEnumerator<int> GetEnumerator()
{
yield return 2;
for (int lsN = 3; lsN <= msLimit; lsN+=2)
if (!mbaOddEliminated[lsN>>1]) yield return lsN;
}
private void FindPrimes()
{
mbaOddEliminated = new BitArray((msLimit>>1) + 1);
int lsSQRT = (int)Math.Sqrt(msLimit);
for (int lsN = 3; lsN < lsSQRT + 1; lsN += 2)
if (!mbaOddEliminated[lsN>>1])
for (int lsM = lsN*lsN; lsM <= msLimit; lsM += lsN<<1)
mbaOddEliminated[lsM>>1] = true;
}
}
}
Se si ottiene l'Atkin a correre più veloce, per favore fatemelo sapere!
Ecco un miglioramento del crivello di Eratostene utilizzando FixBitArrays personalizzati e codice non sicuro per i risultati di velocità, questo è circa 225% più veloce rispetto al mio algoritmo precedente Eratostene, e la classe è standalone (questo non è multithreaded - Eratostene non è compatibile con il multi threading), su un AMD Phenom II X4 965 Processor posso calcolare Primes a 1.000.000.000 limite in 9.261 ms:
using System;
using System.Collections;
using System.Collections.Generic;
namespace PrimeGenerator
{
// The block element type for the bit array,
// use any unsigned value. WARNING: UInt64 is
// slower even on x64 architectures.
using BitArrayType = System.UInt32;
// This should never be any bigger than 256 bits - leave as is.
using BitsPerBlockType = System.Byte;
// The prime data type, this can be any unsigned value, the limit
// of this type determines the limit of Prime value that can be
// found. WARNING: UInt64 is slower even on x64 architectures.
using PrimeType = System.UInt32;
/// <summary>
/// Calculates prime number using the Sieve of Eratosthenes method.
/// </summary>
/// <example>
/// <code>
/// var lpPrimes = new Eratosthenes(1e7);
/// foreach (UInt32 luiPrime in lpPrimes)
/// Console.WriteLine(luiPrime);
/// </example>
public class Eratosthenes : IEnumerable<PrimeType>
{
#region Constants
/// <summary>
/// Constant for number of bits per block, calculated based on size of BitArrayType.
/// </summary>
const BitsPerBlockType cbBitsPerBlock = sizeof(BitArrayType) * 8;
#endregion
#region Protected Locals
/// <summary>
/// The limit for the maximum prime value to find.
/// </summary>
protected readonly PrimeType mpLimit;
/// <summary>
/// The current bit array where a set bit means
/// the odd value at that location has been determined
/// to not be prime.
/// </summary>
protected BitArrayType[] mbaOddNotPrime;
#endregion
#region Initialisation
/// <summary>
/// Create Sieve of Eratosthenes generator.
/// </summary>
/// <param name="limit">The limit for the maximum prime value to find.</param>
public Eratosthenes(PrimeType limit)
{
// Check limit range
if (limit > PrimeType.MaxValue - (PrimeType)Math.Sqrt(PrimeType.MaxValue))
throw new ArgumentOutOfRangeException();
mpLimit = limit;
FindPrimes();
}
#endregion
#region Private Methods
/// <summary>
/// Finds the prime number within range.
/// </summary>
private unsafe void FindPrimes()
{
// Allocate bit array.
mbaOddNotPrime = new BitArrayType[(((mpLimit >> 1) + 1) / cbBitsPerBlock) + 1];
// Cache Sqrt of limit.
PrimeType lpSQRT = (PrimeType)Math.Sqrt(mpLimit);
// Fix the bit array for pointer access
fixed (BitArrayType* lpbOddNotPrime = &mbaOddNotPrime[0])
// Scan primes up to lpSQRT
for (PrimeType lpN = 3; lpN <= lpSQRT; lpN += 2)
// If the current bit value for index lpN is cleared (prime)
if (
(
lpbOddNotPrime[(lpN >> 1) / cbBitsPerBlock] &
((BitArrayType)1 << (BitsPerBlockType)((lpN >> 1) % cbBitsPerBlock))
) == 0
)
// Leave it cleared (prime) and mark all multiples of lpN*2 from lpN*lpN as not prime
for (PrimeType lpM = lpN * lpN; lpM <= mpLimit; lpM += lpN << 1)
// Set as not prime
lpbOddNotPrime[(lpM >> 1) / cbBitsPerBlock] |=
(BitArrayType)((BitArrayType)1 << (BitsPerBlockType)((lpM >> 1) % cbBitsPerBlock));
}
/// <summary>
/// Gets a bit value by index.
/// </summary>
/// <param name="bits">The blocks containing the bits.</param>
/// <param name="index">The index of the bit.</param>
/// <returns>True if bit is set, false if cleared.</returns>
private bool GetBitSafe(BitArrayType[] bits, PrimeType index)
{
return (bits[index / cbBitsPerBlock] & ((BitArrayType)1 << (BitsPerBlockType)(index % cbBitsPerBlock))) != 0;
}
#endregion
#region Public Properties
/// <summary>
/// Get the limit for the maximum prime value to find.
/// </summary>
public PrimeType Limit
{
get
{
return mpLimit;
}
}
/// <summary>
/// Returns the number of primes found in the range.
/// </summary>
public PrimeType Count
{
get
{
PrimeType lptCount = 0;
foreach (PrimeType liPrime in this)
lptCount++;
return lptCount;
}
}
/// <summary>
/// Determines if a value in range is prime or not.
/// </summary>
/// <param name="test">The value to test for primality.</param>
/// <returns>True if the value is prime, false otherwise.</returns>
public bool this[PrimeType test]
{
get
{
if (test > mpLimit) throw new ArgumentOutOfRangeException();
if (test <= 1) return false;
if (test == 2) return true;
if ((test & 1) == 0) return false;
return !GetBitSafe(mbaOddNotPrime, test >> 1);
}
}
#endregion
#region Public Methods
/// <summary>
/// Gets the enumerator for the primes.
/// </summary>
/// <returns>The enumerator of the primes.</returns>
public IEnumerator<PrimeType> GetEnumerator()
{
// Two always prime.
yield return 2;
// Start at first block, second MSB.
int liBlock = 0;
byte lbBit = 1;
BitArrayType lbaCurrent = mbaOddNotPrime[0] >> 1;
// For each value in range stepping in incrments of two for odd values.
for (PrimeType lpN = 3; lpN <= mpLimit; lpN += 2)
{
// If current bit not set then value is prime.
if ((lbaCurrent & 1) == 0)
yield return lpN;
// Move to NSB.
lbaCurrent >>= 1;
// Increment bit value.
lbBit++;
// If block is finished.
if (lbBit == cbBitsPerBlock)
{
// Move to first bit of next block.
lbBit = 0;
liBlock++;
lbaCurrent = mbaOddNotPrime[liBlock];
}
}
}
#endregion
#region IEnumerable<PrimeType> Implementation
/// <summary>
/// Gets the enumerator for the primes.
/// </summary>
/// <returns>The enumerator for the prime numbers.</returns>
IEnumerator IEnumerable.GetEnumerator()
{
return GetEnumerator();
}
#endregion
}
}
Primes trovati nelle miliardo: 50.847.534 a 9.261 ms
