C ++ Сито Аткина пропускает несколько простых чисел - PullRequest
3 голосов
/ 27 января 2010

Недавно я работал над генератором простых чисел C ++, который использует Sieve of Atkin (http://en.wikipedia.org/wiki/Sieve_of_atkin) для генерации простых чисел. Моя цель - создать 32-битное число. Я буду использовать его в основном для задач проекта Эйлера. в основном это просто летний проект.

Программа использует битовую доску для хранения первичности: то есть серии единиц и нулей, где, например, 11-й бит будет 1, 12-й 0 и 13-й 1 и т. Д. Для эффективного использования памяти, на самом деле это массив символов, каждый из которых содержит 8 бит. Я использую флаги и побитовые операторы для установки и получения битов. Суть алгоритма проста: сделать первый проход, используя некоторые уравнения, которые я не претендую понять, чтобы определить, считается ли число «простым» или нет. В большинстве случаев это даст правильные ответы, но пара не простых чисел будет помечена как простые. Следовательно, при переборе по списку вы устанавливаете все кратные простому числу, которое вы только что нашли, как «не простое». Это имеет удобное преимущество: требуется меньше процессорного времени, чем больше простое число.

У меня это завершено на 90%, с одним уловом: некоторые простые числа отсутствуют.

Благодаря проверке битборда я установил, что эти простые числа опущены во время первого прохода, что в основном переключает число для каждого решения, которое оно имеет для ряда уравнений (см. Статью в Википедии). Я перебирал этот кусок кода снова и снова. Я даже пытался увеличить границы показанных в статьях Википедии, что менее эффективно, но я подумал, что может попасть в несколько цифр, которые я как-то пропустил. Ничего не сработало. Эти числа просто оценивают как не простые. Большая часть моего теста была на всех простых числах до 128. Из этого диапазона, эти простые числа опущены:

23 и 59.

Я не сомневаюсь, что в более высоком диапазоне больше будет отсутствовать (просто не хочу считать их всех). Я не знаю, почему они отсутствуют, но они есть. Есть ли что-то особенное в этих двух простых числах? Я дважды и трижды проверил, нахожу и исправляю ошибки, но все же, наверное, что-то глупое, что мне не хватает.

В любом случае, вот мой код:

#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;
    }
}

Я сделал все возможное, чтобы очистить его и сделать его читабельным. Я не профессиональный программист, поэтому, пожалуйста, будь милостив.

Вот вывод, который я получаю, когда я инициализирую объект PrimeGen с именем pgen, распечатываю его начальную битовую доску с помощью pgen.printBoard () (обратите внимание, что 23 и 59 отсутствуют до итерации next ()), а затем перебираете следующую () и выведите все возвращенные простые числа:

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.

1 Ответ

5 голосов
/ 27 января 2010

Эврика !!!

Как и ожидалось, это была глупая ошибка с моей стороны.

У уравнения 3x ^ 2 - y ^ 2 есть небольшая оговорка, которую я пропустил: x> y. Учитывая это, я переключался на 23 и 59 слишком много раз, что приводило к их сбою.

Спасибо за помощь, ребята. Спас мой бекон.

...