Список простых чисел с использованием метода Sieve с использованием битовой маски

Я написал следующий код, чтобы перечислить все простые числа до 2 миллиардов, используя метод Sieve. Я использовал битмаскирование для целей маркировки. Хотя я могу правильно получить простые числа, несколько простых чисел в начале каждый раз отсутствуют. Помогите пожалуйста найти ошибку в программе.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>

#define MAX 2000000000

char* listPrimes(){
int block = sqrt(MAX);
char* mark = calloc((MAX/8),sizeof(char));
int i = 2;
int j;
char mask[8];
for(j=0;j<8;j++)
    mask[j] = 0;
mask[7] = 1;
mask[6] |= mask[7] << 1;
mask[5] |= mask[7] << 2;
mask[4] |= mask[7] << 3;
mask[3] |= mask[7] << 4;
mask[2] |= mask[7] << 5;
mask[1] |= mask[7] << 6;
mask[0] |= mask[7] << 7;

for(j=0;j<8;j++)
    printf("%d ",mask[j]);
mark[0] |= mask[0];
mark[0] |= mask[1];

while (i < block){

        for (j = 2; i*j <= block; j++)
                mark[(i*j) / 8] |= mask[((i*j) % 8 )];
        i++;
    }
printf("\n");
printf("The block size is\t:\t%d\n",block);


j = 2;
while(j<=block){
    if((mark[j / 8] & mask[j]) == 0 ){
        for(i = 2;i <= MAX; i++){
            if((i%j) == 0){
                mark[i / 8] |= mask[(i % 8)];
            }
        }
    }
while((mark[++j / 8] & mask[j % 8]) != 0);
}


for(j=0;j<=MAX;j++)
        if((mark[j / 8] & mask[(j % 8)]) == 0)
            printf("%d\n", ((8*(j / 8)) + (j % 8)));

return mark;
}   

int main(int argc,char* argv[]){

listPrimes();

return 0;
}

c arrays bitmask primes sieve-of-eratosthenes
person Krishna    schedule 08.01.2013    source источник

Ответы (3)


arrow_upward
1
arrow_downward

Как сказал ArunMK, во втором цикле while вы отмечаете простое число j как кратное j. И, как сказал Ли Меадор, вам нужно взять модуль j по модулю 8 для индекса mask, иначе вы получите доступ выходит за пределы и вызывает неопределенное поведение.

Еще один момент, когда вы вызываете неопределенное поведение,

while((mark[++j / 8] & mask[j % 8]) != 0);

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

do {
    ++j;
}while((mark[j/8] & mask[j%8]) != 0);

или, если вы настаиваете на цикле while с пустым телом

while(++j, (mark[j/8] & mask[j%8]) != 0);

вы можете использовать оператор запятой.

Более неопределенное поведение при доступе к mark[MAX/8], который не выделен в

for(i = 2;i <= MAX; i++){

и

for(j=0;j<=MAX;j++)

Кроме того, если char имеет знак и имеет ширину восемь бит,

mask[0] |= mask[7] << 7;

определяется реализацией (и может вызвать сигнал, определяемый реализацией), поскольку результат

mask[0] | (mask[7] << 7)

(int 128) не может быть представлено как char.

Но почему вы делите каждое число на все простые числа, не превышающие квадратный корень из оценки во втором цикле while?

    for(i = 2;i <= MAX; i++){
        if((i%j) == 0){

Это делает ваш алгоритм не решетом Эратосфена, а пробным делением.

Почему бы вам не использовать и здесь технику из первого while цикла? (А зачем вообще две петли?)

while (i <= block){
    if ((mark[i/8] & mask[i%8]) == 0) {
        for (j = 2; i*j < MAX; j++) {
            mark[(i*j) / 8] |= mask[((i*j) % 8 )];
        }
    }
    i++;
}

не будет переполняться (для заданного значения MAX, если это можно представить как int) и будет производить правильные выходные порядки быстрее.

person Daniel Fischer    schedule 08.01.2013
comment
Определено, что ++j произойдет до того, как значение j будет использовано в любом месте, поскольку оценка гарантированно выполняется слева направо. Нет необходимости менять «пока». Однако предложенный вами код легче читать и понимать, особенно в 2012 году, когда использование оператора ++ таким образом не одобряется, и новые разработчики могут найти его незнакомым. - person Lee Meador; 09.01.2013
comment
поскольку оценка гарантированно выполняется слева направо. Нет, это не так. Порядок оценки операндов & не указан, и не гарантируется, что сохранение увеличенного значения произошло до следующей точки последовательности. Таким образом, даже если mark[++j/8] оценивается первым, mask[j%8] все еще может считывать неинкрементированное значение j. Скомпилируйте с -Wsequence-point (и оптимизацией), и ваш gcc предупредит вас об этом. - person Daniel Fischer; 09.01.2013
comment
Ты прав. Я не прав. Только && и || гарантируют слева направо, и они не гарантируют, что выражение справа вообще вычисляется. Есть и другие очевидные, такие как , и ?: которые имеют порядок. «запятая» предназначена специально для выполнения действий по порядку, а ?: должен оценить условие, чтобы указать, какое из двух других выражений следует оценивать. - person Lee Meador; 09.01.2013
comment
Что этот второй цикл должен выполнить? Все числа, кратные простым числам, уже должны быть отмечены в 1-м цикле. Даже ваше предложение проверить текущее i, уже помеченное как кратное чему-то меньшему, не заставляет его работать, а просто делает его более эффективным, пропуская, когда мы знаем, что все кратные уже установлены. - person Lee Meador; 09.01.2013
comment
В оригинале первая петля помечает только квадратный корень из предела, поэтому вторая необходима для пометки композитов > block = sqrt(MAX). Я предлагаю использовать только один цикл и отмечать только кратные простые числа. Оно работает. Я опустил очевидные улучшения (начните с i*i, не вычисляйте i*j, а увеличивайте j на i [2*i для нечетных i], ...), чтобы оставаться ближе к коду OP. - person Daniel Fischer; 09.01.2013
comment
Я пропустил ваше изменение в конце внутреннего цикла с «блока» на «МАКС». Вы также можете вдвое уменьшить размер массива или сократить время (путем предварительного ИЛИ с 0xCC), отметив, что биты 0, 2, 4 и 6 установлены во всех байтах, кроме 1-го, но это не то, что спросил OP . - person Lee Meador; 09.01.2013
comment
Да, удаление четных чисел из сита — одно из очевидных улучшений, о которых я не упомянул. Для данного MAX это нужно сделать, вероятно, также следует удалить числа, кратные 3. - person Daniel Fischer; 09.01.2013

arrow_upward
1
arrow_downward

Измените средний цикл, чтобы добавить модуль:

j = 2;
while(j<=block){
    if((mark[j / 8] & mask[j % 8]) == 0 ){
        for(i = 2;i <= MAX; i++){
            if((i%j) == 0){
                mark[i / 8] |= mask[(i % 8)];
            }
        }
    }
}
person Lee Meador    schedule 08.01.2013

arrow_upward
1
arrow_downward

Во втором цикле while вы перебираете i, начиная с 2, и выполняете if (i%j == 0). Это будет верно и для i, когда это простое число. Вам нужно проверить (i != j). Также по модулю, как указано выше. Отсюда получается: if ((i%j == 0) { if (i!=j) mark[i/j] |= mask[i%j]; }

person user1952500    schedule 08.01.2013
comment
Или измените этот цикл for(i = 2;i <= MAX; i++), чтобы он начинался с 2*j - person Lee Meador; 09.01.2013