غربال أتكين

في الرياضيات، غربال أتكين هو خوارزمية سريعة وعصرية لإيجاد جميع الأعداد الأولية الأصغر من عدد ما. هو صيغة متطورة للغربال القديم والمعروف بغربال إراتوستينس.

نظرت هاته الخوارزمية من طرف العالمين أرثور أوليفر لوندال أتكين ودانييل يوليوس بيرنشتاين.

نسخة الخوارزمية على لغة الباسكال

program atkin;
var is_prime:array[1..10001] of boolean; jj:int64;
procedure dosieve(limit:int64);
var i,k,x,y:int64; n:int64;
begin
  for i:=5 to limit do
    is_prime[i]:=false;
  for x:=1 to trunc(sqrt(limit)) do
    for y:=1 to trunc(sqrt(limit)) do
    begin
      n:=4*sqr(x)+sqr(y);
      if (n<=limit) and ((n mod 12 = 1) or (n mod 12 = 5)) then
        is_prime[n]:=not is_prime[n];
      n:=n-sqr(x);
      if (n<=limit) and (n mod 12 = 7) then
        is_prime[n]:=not is_prime[n];
      n:=n-2*sqr(y);
      if (x>y) and (n<=limit) and (n mod 12 = 11) then
        is_prime[n]:=not is_prime[n];
    end;
  for i:=5 to trunc(sqrt(limit)) do
  begin
    if is_prime[i] then
    begin
      k:=sqr(i);
      n:=k;
      while n<=limit do
      begin
        is_prime[n]:=false;
        n:=n+k;
      end;
    end;
  end;
  is_prime[2]:=true;
  is_prime[3]:=true;
end;
begin
  dosieve(10000);
  for jj:=1 to 10000 do
    if is_prime[jj] then
      writeln(jj);
end.

النسخة المبسطة للخوارزمية

int limit = 1000;
int sqr_lim;
bool is_prime[1001];
int x2, y2;
int i, j;
int n;

// Инициализация решета
sqr_lim = (int)sqrt((long double)limit);
for (i = 0; i <= limit; i++) is_prime[i] = false;
is_prime[2] = true;
is_prime[3] = true;

// Предположительно простые - это целые с нечетным числом
// представлений в данных квадратных формах.
// x2 и y2 - это квадраты i и j (оптимизация).
x2 = 0;
for (i = 1; i <= sqr_lim; i++) {
    x2 += 2 * i - 1;
    y2 = 0;
    for (j = 1; j <= sqr_lim; j++) {
        y2 += 2 * j - 1;
        
        n = 4 * x2 + y2;
        if ((n <= limit) && (n % 12 == 1 || n % 12 == 5))
            is_prime[n] = !is_prime[n];
        
        // n = 3 * x2 + y2; 
        n -= x2; // Оптимизация
        if ((n <= limit) && (n % 12 == 7))
            is_prime[n] = !is_prime[n];
        
        // n = 3 * x2 - y2;
        n -= 2 * y2; // Оптимизация
        if ((i> j) && (n <= limit) && (n % 12 == 11))
            is_prime[n] = !is_prime[n];
    }
}

// Отсеиваем кратные квадратам простых чисел в интервале [5, sqrt(limit)].
// (основной этап не может их отсеять)
for (i = 5; i <= sqr_lim; i++) {
    if (is_prime[i]) {
        n = i * i;
        for (j = n; j <= limit; j += n) {
            is_prime[j] = false;
        }
    }
}

// Вывод списка простых чисел в консоль.
printf("2, 3, 5"); 
for (i = 6; i <= limit; i++) {  // добавлена проверка делимости на 3 и 5. В оригинальной версии алгоритма потребности в ней нет.
    if ((is_prime[i]) && (i % 3 != 0) && (i % 5 !=  0)){
       printf(", %d", i);
    }
}