C #: La mise en œuvre de la Crible d'Atkin
https://stackoverflow.com/questions/1569127
-
21-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Je me demandais si quelqu'un ici a une bonne mise en œuvre du Crible d'Atkin qu'ils aimeraient partager.
Je suis en train de le mettre en œuvre, mais ne peut pas tout à fait envelopper la tête autour d'elle. Voici ce que j'ai à ce jour.
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();
}
}
Je suis à peu près juste essayé de « traduire » le pseudo-code figurant à Wikipedia , mais il ne fonctionne pas correctement. Donc, soit j'ai mal compris quelque chose ou juste fait quelque chose de mal. Ou plus probablement les deux ...
Avoir une liste des 500 premiers nombres premiers que j'utilise comme un test et ma mise en œuvre échoue au numéro 40 (ou 41?).
Les valeurs diffèrent à l'index [40]
Prévu: 179
Mais était: 175
Êtes-vous capable de trouver mon erreur, avoir une mise en œuvre autour de la pose que vous pourriez partager vous, ou les deux?
Le test exact que je utilise ressemble à ceci:
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, ...
};
}
La solution
Ce code:
for (ulong k = n*n; k <= limit; k *= k)
isPrime[k] = false;
ne semble pas être une traduction fidèle de ce pseudo-code:
is_prime(k) ← false, k ∈ {n², 2n², 3n², ..., limit}
Votre code ressemble il fonctionnera pour n * n, n ^ 4, n ^ 8, etc. dire à chaque fois élevant au carré au lieu d'ajouter n-carré à chaque fois. Essayez ceci:
ulong nSquared = n * n;
for (ulong k = nSquared; k <= limit; k += nSquared)
isPrime[k] = false;
Autres conseils
Le dernière réponse par Aaron Mugatroyd à partir de la source Python traduit pour une Crible d'Atkin (SoA) n'est pas trop mal, mais il peut être amélioré à plusieurs égards car il manque quelques optimisations importantes, comme suit:
-
Sa réponse n'utilise pas le plein modulo version 60 Atkin et Bernstein original de la Sieve, mais plutôt une variation légèrement améliorée du code pseudo de l'article Wikipedia utilise donc d'un facteur de 0,36 de la gamme de tamis numérique opérations bascule / réforme combinée; mon code ci-dessous utilise le code pseudo segment non paginé raisonnablement efficace par mes commentaires dans une réponse commentant le Crible d'Atkin qui utilise un facteur d'environ 0,26 fois la plage numérique pour réduire la quantité de travail effectué à environ un facteur d'environ deux septièmes.
-
son code réduit la taille de la mémoire tampon en ayant seulement des représentations de nombre impair, alors que mon code supplémentaire emballe peu pour éliminer toute représentation des nombres divisibles par trois et cinq, ainsi que ceux divisible par deux sous-entendus par « odds seulement « ; ce qui réduit l'exigence de mémoire par un facteur supplémentaire de près de la moitié (à 8/15) et aider à mieux utiliser des caches CPU pour une nouvelle augmentation de la vitesse en raison du temps d'accès à la mémoire moyenne réduite.
-
Mon code compte le nombre de nombres premiers en utilisant une technique de comptage pop Regardez rapide Up Table (LUT) de prendre presque pas le temps de compter par rapport à l'environ une seconde en utilisant la technique bit par bit, il utilise; Cependant, dans cet exemple de code même que peu de temps est retiré de la partie chronométrée du code.
-
Enfin, mon code optimise les opérations de manipulation de bits pour un minimum de code par boucle interne. Par exemple, il ne se sert pas d'un changement continuel droit de produire l'indice de représentation impairs et en fait peu de décalage du tout en écrivant toutes les boucles intérieures que les opérations modulo constante (égale position de bit). De plus, le code traduit d'Aaron est tout à fait inefficace dans les opérations comme par exemple dans l'abattage gratuit premier carré, il ajoute la place du premier à l'indice vérifie ensuite pour un résultat étrange plutôt que d'ajouter deux fois la place afin de ne pas exiger la vérification ; il fait même la redondance de contrôle en déplaçant le bon nombre par un (en divisant par deux) avant de faire l'opération d'abattage dans la boucle intérieure, comme il le fait pour toutes les boucles. Ce code inefficace ne fera pas beaucoup de différence en temps d'exécution pour les grandes plages en utilisant cette technique « large éventail tampon tamis », comme la plupart du temps par opération est utilisé dans l'accès à la mémoire RAM (environ 37 cycles d'horloge du CPU ou plus pour un gamme d'un milliard), mais fera le temps d'exécution beaucoup plus lent que ce doit être pour les gammes plus petites qui entrent dans les caches CPU; en d'autres termes, il fixe une limite trop élevée la plus faible vitesse d'exécution par opération.
Le code est le suivant:
//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();
}
}
}
Ce code fonctionne deux fois plus vite que le code d'Aaron (environ 2,7 secondes en utilisant le mode 64 bits ou 32 bits sur un i7-2700K (3,5 GHz) avec le tampon d'environ 16,5 méga-octets et environ 0258000000 bascule combiné / carré premier les opérations de réforme libre (qui peuvent être représentés par la décommentant « ++ » this.opcnt états) pour une gamme de tamis d'un milliard, comparativement à 5,4 / 6,2 secondes (32 bits / 64 bits) pour son code sans le décompte temps et presque deux fois l'utilisation de la mémoire en utilisant environ 0359000000 opérations bascule / combinés pour tamiser réforme jusqu'à un milliard.
Bien qu'il soit plus rapide que ses chances seulement naïves les plus optimisées de mise en œuvre de la Sieve d'Eratosthène (SOE) non paginée, qui ne font pas le Crible d'Atkin plus vite que le Crible d'Eratosthène , comme si l'on applique des techniques similaires à celles utilisées dans la mise en œuvre de SoA ci-dessus pour Soe, plus utilise factorisation de roue maximale, le SOE sur les ee même vitesse que cela.
Analyse: Bien que le nombre d'opérations pour l'état de l'environnement entièrement optimisé sont environ le même que le nombre d'opérations pour le SoA pour une plage de tamis d'un milliard, le principal goulot d'étranglement pour ces non paginée implémentations est accès à la mémoire une fois que la taille du tampon de tamisage est supérieure à la taille de la mémoire cache du processeur (32 kilo-octets de mémoire cache L1 à un accès de cycle d'horloge, 256 cache kilo-octets L2 à environ quatre cycles d'horloge du temps d'accès et de 8 Mo de cache L3 au moment de l'accès à environ 20 cycles d'horloge pour mon i7), après quoi l'accès mémoire peut dépasser une centaine de cycles d'horloge.
Maintenant, les deux ont un facteur d'environ huit amélioration des vitesses d'accès mémoire lorsque l'on adapte les algorithmes de segmentation de la page, donc on peut tamiser des plages qui ne répondraient pas autrement dans la mémoire disponible. Cependant, l'entreprise publique continue à gagner les SoA que la gamme de tamis commence à être très grand en raison de difficultés dans la mise en œuvre de la partie de l'algorithme « de nombres premiers sans carré » en raison des énormes progrès dans l'abattage des analyses qui poussent rapidement à plusieurs centaines de fois la taille des tampons de page. De plus, et peut-être plus grave, il est très mémoire et / ou de calculs pour calculer le nouveau point de départ pour chaque valeur de « x » à la valeur de « y » à la plus faible représentation de chaque tampon de page pour une autre tout à fait grande perte d'efficacité de l'paginée SoA comparaed à Soe que la gamme se développe.
Edit_add: Les odds que SOE utilisé par Aaron Murgatroyd utilise environ 1,026 milliards d'opérations de réforme pour une gamme de tamis d'un milliard donc environ quatre fois plus d'opérations que l'SoA et devrait donc courir quatre fois plus lente, mais la SoA même mis en œuvre ici a une boucle intérieure plus complexe et surtout en raison d'une proportion beaucoup plus élevée des cotes-seulement abattages SoE ont un pas beaucoup plus courte dans les analyses d'abattage que les progrès de la SoA les odds naïves -Seulement SoE a bien meilleurs temps d'accès mémoire moyenne malgré le tampon de tamis dépassant largement la taille du cache du processeur (meilleure utilisation du cache associativité). Ceci explique pourquoi le SoA ci-dessus est seulement deux fois plus vite que les chances ne-SoE même si il semblerait théoriquement être faire seulement un quart du travail.
Si l'on devait utiliser un algorithme similaire en utilisant constante modulo boucles internes que pour les SoA ci-dessus et mis en œuvre le même 2/3/5 roue factorisation, Soe réduirait le nombre d'opérations de réforme à environ 0405000000 opérations si seulement 50% plus d'opérations que la SoA et courraient théoriquement juste un peu plus lent que le SoA, mais peuvent fonctionner à peu près la même vitesse en raison des progrès de réforme étant encore un peu plus faible que pour la SoA en moyenne pour cette « naïve » grande mémoire utilisation de la mémoire tampon. L'augmentation de la roue factorisation au moyen de roue 2/3/5/7 les opérations de réforme des entreprises publiques sont réduites à environ 0,314 pour une gamme d'abattage d'un milliard et peut rendre cette version de la course soe la même vitesse pour cet algorithme.
En outre l'utilisation de la roue factorisation peut être faite par pré-abattage du tableau de tamis (copie dans un motif) pour les 2/3/5/7/11/13/17/19 facteurs premiers à presque aucun coût en temps d'exécution de réduire le nombre total d'opérations de réforme à environ 0,251 milliards pour une gamme de tamis d'un milliard et Soe tournera plus vite vers la même vitesse que même cette version optimisée du SoA, même pour ces grandes versions de mémoire tampon, avec Soe ayant encore beaucoup moins que la complexité du code ci-dessus.
Ainsi, on peut voir que le nombre d'opérations pour l'entreprise publique peut être considérablement réduit d'un naïf ou même odds seule ou 2/3/5 roue version factorisation de telle sorte que le nombre d'opérations sont du même ordre que pour SoA tout en même temps le temps par opération peut effectivement être moins en raison de deux boucles internes moins complexes et un accès mémoire plus efficace. END_EDIT_ADD
EDIT_ADD2: ajouter icile code d'un état de l'environnement en utilisant une technique de position modulo / binaire constant semblable pour les boucles les plus internes que pour le SoA ci-dessus selon le pseudo-code plus bas la réponse ci-dessus comme liée . Le code est un peu moins complexe que les SoA ci-dessus en dépit d'avoir factorisation haut de la roue et de pré-abattage appliquée de telle sorte que le nombre total d'opérations de réforme sont en fait moins que les opérations de basculement / réforme combinées pour l'SoA jusqu'à un rang de tamisage d'environ deux milliards. Le code comme suit:
EDIT_FINAL code corrigé ci-dessous et les commentaires qui s'y rapporte 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();
}
}
}
Ce code fonctionne en fait quelques pour cent plus rapidement que les SoA ci-dessus comme il devrait comme il y a des opérations un peu moins et le principal goulot d'étranglement pour cette grande taille du tableau pour une gamme d'un milliard est temps d'accès à la mémoire de quelque chose comme 40 à plus de 100 cycles d'horloge du processeur en fonction des spécifications du processeur et de la mémoire; cela signifie que les optimisations de code (autres que la réduction du nombre total d'opérations) sont inefficaces, car la plupart du temps est passé à attendre sur l'accès mémoire. En tout cas, en utilisant un énorme tampon de mémoire n'est pas le plus efficace moyen de tamis de grandes gammes, avec un facteur allant jusqu'à environ huit fois l'amélioration de la segmentation de la page à l'aide SoE avec le même factorisation de roue maximale (qui ouvre la voie à multi-traitement).
Il est dans l'application de la page segmentation et multi-traitement que le SoA est vraiment déficiente pour les gammes bien au-dessus de quatre milliards par rapport à l'état de l'environnement comme des gains en raison de la complexité asymptotique réduite du SoA se rapidement mangées par la page de traitement frais généraux facteurs liés au traitement gratuit premier carré et le calcul du plus grand nombre d'adresses de début de la page; En variante, on surmonte cette marqueurs en stockant dans la mémoire RAM à un coût énorme consommation de mémoire et d'autres inefficacités dans l'accès à ces structures de magasin de marqueur. END_EDIT_ADD2
En bref, le SoA est pas vraiment un tamis pratique par rapport au Soe factorisée entièrement roue depuis tout comme le gain de complexité asymptotique commence à la rapprocher des performances de l'entreprise publique entièrement optimisé, il commence à perdre de leur efficacité en raison des détails de mise en œuvre pratique dans le temps d'accès à la mémoire relative et la complexité de la segmentation de la page ainsi que généralement être plus complexe et difficile à écrire. À mon avis, il est plus d'un concept intellectuel intéressant et exercice mental qu'un tamis pratique par rapport à l'état de l'environnement.
Un jour, j'adapter ces techniques à une page multi-thread segmentée Sieve d'Ératosthène pour être à peu près aussi rapide en C # en tant que mise en œuvre « primegen » de Atkin et Bernstein du SoA en « C » et soufflera hors de l'eau pour les grandes gammes supérieures à environ quatre milliards, même un seul thread, avec un coup de pouce supplémentaire de la vitesse allant jusqu'à environ quatre ans quand multithread sur mon i7 (huit cœurs, y compris Hyper Threading).
Voici une autre mise en œuvre. Il utilise BitArray pour économiser la mémoire. Parallel.For a besoin .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;
}
Voici une rapide mise en œuvre du Crible d'Atkin, j'ai volé l'algorithme de ce script Python ici (je prends pas de crédit pour l'algorithme):
http://programmingpraxis.com/2010/02/ 19 / tamis-of-Atkin amélioré /
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
}
}
Sa vitesse proche de ma version la plus optimisée de la Sieve d'Eratosthène, mais il est encore plus lente d'environ 20%, il se trouve ici:
le mien Heres, il utilise une classe appelée CompartmentalisedParallel qui vous permet d'effectuer en parallèle pour les boucles, mais contrôler le nombre de threads de telle sorte que les indices sont regroupés vers le haut. Cependant, en raison des problèmes de filetage, vous devez soit verrouiller le BitArray chaque fois qu'il est modifié ou créer un BitArray pour chaque fil, puis les XOR ensemble à la fin, la première option était assez lent car la quantité de serrures, la deuxième l'option semblait plus rapide pour moi!
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;
}
}
}
}
Et 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 de base nombres premiers:
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();
}
}
}
Utilisez-en procédant comme suit:
var lpPrimes = new Atkin(count, true);
Console.WriteLine(lpPrimes.Count);
Console.WriteLine(s.ElapsedMilliseconds);
Cependant, je trouve les Ératosthène pour être plus rapide dans tous les cas, même avec un processeur quatre coeur fonctionnant en mode multithread pour le 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;
}
}
}
Si vous obtenez le Atkin de courir plus vite, s'il vous plaît laissez-moi savoir!
Heres une amélioration de la Sieve d'Eratosthène en utilisant FixBitArrays personnalisés et le code dangereux pour les résultats de la vitesse, cela est environ 225% plus rapide que mon précédent algorithme Eratosthène et la classe est autonome (ce n'est pas multithread - Eratosthène n'est pas compatible avec le multi filetage), sur un processeur AMD Phenom II X4 965 Processor Je peux calculer PRIMES pour l'1.000.000.000 limite 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
}
}
trouvés dans 1000000000 nombres premiers: 50847534 en 9,261 ms
