C++ Tamis d'Atkin surplombe quelques nombres premiers

Question

Récemment, j'ai travaillé sur un générateur C++ prime qui utilise le tamis d'Atkin (http://en.wikipedia.org/wiki/Sieve_of_atkin) pour générer ses nombres premiers. Mon objectif est de pouvoir générer n'importe quel nombre de 32 bits. Je vais l'utiliser surtout pour des problèmes de projet euler. surtout c'est juste un projet d'été.

Le programme utilise un bitboard pour stocker la primalité: c'est-à-dire une série de uns et de zéros où par exemple le 11ème bit serait un 1, le 12 a 0 et le 13 a 1, etc. , c'est en fait un tableau de caractères, chaque caractère contenant 8 bits. J'utilise des drapeaux et des opérateurs au niveau du bit pour définir et récupérer des bits. Le gist de l'algorithme est simple: faire un premier passage en utilisant des équations que je ne prétends pas comprendre pour définir si un nombre est considéré comme «premier» ou non. Ce sera pour la plupart obtenir les bonnes réponses, mais quelques numéros nonprime seront marqués comme premier. Par conséquent, lorsque vous parcourez la liste, vous définissez tous les multiples du nombre premier que vous venez de trouver comme "non premiers". Ceci a l'avantage pratique de nécessiter moins de temps de traitement, plus le premier devient important.

Je l'ai terminé à 90%, avec un crochet: certains nombres premiers sont manquants. Par l'inspection du tableau, j'ai vérifié que ces nombres premiers sont omis lors du premier passage, ce qui change fondamentalement un nombre pour chaque solution qu'il a pour un certain nombre d'équations (voir entrée wikipedia). J'ai parcouru ce morceau de code maintes et maintes fois. J'ai même essayé d'augmenter les limites de ce qui est montré dans les articles de Wikipédia, ce qui est moins efficace mais j'ai pensé que je pourrais frapper quelques chiffres que j'ai omis d'une certaine manière. Rien n'a fonctionné. Ces chiffres évaluent simplement à ne pas amorcer. La plupart de mon test a été sur tous les nombres premiers en 128. De cette gamme, ce sont les nombres premiers qui sont omis:

23 et 59.

Je ne doute pas que sur une gamme plus élevée, plus il manquerait (Je ne veux pas les compter tous). Je ne sais pas pourquoi ils manquent, mais ils le sont. Y a-t-il quelque chose de spécial à propos de ces deux nombres premiers? J'ai vérifié deux fois et trois fois, trouvé et réparé des erreurs, mais c'est probablement encore quelque chose de stupide qui me manque.

de toute façon, voici mon code:

#include <iostream> 
#include <limits.h> 
#include <math.h> 

using namespace std; 

const unsigned short DWORD_BITS = 8; 

unsigned char flag(const unsigned char); 
void printBinary(unsigned char); 


class PrimeGen 
{ 
    public: 
     unsigned char* sieve; 
     unsigned sievelen; 
     unsigned limit; 
     unsigned bookmark; 


     PrimeGen(const unsigned); 

     void firstPass(); 
     unsigned next(); 

     bool getBit(const unsigned); 
     void onBit(const unsigned); 
     void offBit(const unsigned); 
     void switchBit(const unsigned); 

     void printBoard(); 
}; 


PrimeGen::PrimeGen(const unsigned max_num) 
{ 
    limit = max_num; 
    sievelen = limit/DWORD_BITS + 1; 
    bookmark = 0; 

    sieve = (unsigned char*) malloc(sievelen); 
    for (unsigned i = 0; i < sievelen; i++) {sieve[i] = 0;} 

    firstPass(); 
} 


inline bool PrimeGen::getBit(const unsigned index) 
{ 
    return sieve[index/DWORD_BITS] & flag(index%DWORD_BITS); 
} 


inline void PrimeGen::onBit(const unsigned index) 
{ 
    sieve[index/DWORD_BITS] |= flag(index%DWORD_BITS); 
} 


inline void PrimeGen::offBit(const unsigned index) 
{ 
    sieve[index/DWORD_BITS] &= ~flag(index%DWORD_BITS); 
} 


inline void PrimeGen::switchBit(const unsigned index) 
{ 
    sieve[index/DWORD_BITS] ^= flag(index%DWORD_BITS); 
} 


void PrimeGen::firstPass() 
{ 
    unsigned nmod,n,x,y,xroof, yroof; 

    //n = 4x^2 + y^2 
    xroof = (unsigned) sqrt(((double)(limit - 1))/4); 
    for(x = 1; x <= xroof; x++){ 
     yroof = (unsigned) sqrt((double)(limit - 4 * x * x)); 
     for(y = 1; y <= yroof; y++){ 
      n = (4 * x * x) + (y * y); 
      nmod = n % 12; 
      if (nmod == 1 || nmod == 5){ 
       switchBit(n); 
      } 
     } 
    } 

    xroof = (unsigned) sqrt(((double)(limit - 1))/3); 
    for(x = 1; x <= xroof; x++){ 
     yroof = (unsigned) sqrt((double)(limit - 3 * x * x)); 
     for(y = 1; y <= yroof; y++){ 
      n = (3 * x * x) + (y * y); 
      nmod = n % 12; 
      if (nmod == 7){ 
       switchBit(n); 
      } 
     } 
    } 

    xroof = (unsigned) sqrt(((double)(limit + 1))/3); 
    for(x = 1; x <= xroof; x++){ 
     yroof = (unsigned) sqrt((double)(3 * x * x - 1)); 
     for(y = 1; y <= yroof; y++){ 
      n = (3 * x * x) - (y * y); 
      nmod = n % 12; 
      if (nmod == 11){ 
       switchBit(n); 
      } 
     } 
    } 
} 


unsigned PrimeGen::next() 
{ 
    while (bookmark <= limit) 
    { 
     bookmark++; 

     if (getBit(bookmark)) 
     { 
      unsigned out = bookmark; 

      for(unsigned num = bookmark * 2; num <= limit; num += bookmark) 
      { 
       offBit(num); 
      } 

      return out; 
     } 
    } 

    return 0; 
} 


inline void PrimeGen::printBoard() 
{ 
    for(unsigned i = 0; i < sievelen; i++) 
    { 
     if (i % 4 == 0) 
      cout << endl; 

     printBinary(sieve[i]); 
     cout << " "; 
    } 
} 


inline unsigned char flag(const unsigned char bit_index) 
{ 
    return ((unsigned char) 128) >> bit_index; 
} 


inline void printBinary(unsigned char byte) 
{ 
    unsigned int i = 1 << (sizeof(byte) * 8 - 1); 

    while (i > 0) { 
     if (byte & i) 
      cout << "1"; 
     else 
      cout << "0"; 
     i >>= 1; 
    } 
} 

Je l'ai fait de mon mieux pour le nettoyer et le rendre lisible. Je ne suis pas un programmeur professionnel, alors s'il vous plaît soyez miséricordieux.

Voici la sortie que j'obtiens, lorsque j'initialise un objet PrimeGen nommé pgen, j'imprime son bitboard initial avec pgen.printBoard() (notez que 23 et 59 sont manquants avant l'itération next()), puis passez en revue next() et imprimer tous les nombres premiers retournés:

00000101 00010100 01010000 01000101 
00000100 01010001 00000100 00000100 
00010001 01000001 00010000 01000000 
01000101 00010100 01000000 00000001 

5 
7 
11 
13 
17 
19 
29 
31 
37 
41 
43 
47 
53 
61 
67 
71 
73 
79 
83 
89 
97 
101 
103 
107 
109 
113 
127 

DONE 

Process returned 0 (0x0) execution time : 0.064 s 
Press any key to continue. 

Answer 1

Eureka !!!

Comme prévu, c'était une erreur stupide de ma part.

L'équation 3x^2 - y^2 a une petite mise en garde que j'ai ignorée: x> y. Avec cela pris en compte, je changeais trop souvent de 23 et 59, ce qui les conduisait à l'échec.

Merci pour toute l'aide que vous avez. Sauvé mon bacon.