Решето́ Сундара́ма — детерминированный алгоритм нахождения всех простых чисел до некоторого целого числа $n$ .

Haskell

 sundaram n = 2 : [2*m+1 | m <- [1..k], m `notElem` zs]
   where
     zs = [i + j + 2*i*j | i <- [1..k`div`3], j <- [i..(k-i)`div`(2*i + 1)]]
     k  = (n-1) `div` 2

C++

void sundaram(int n){
    int *a = new int [n], i, j, k;
    memset(a, 0, sizeof(int) * n);
    for(i = 1; 3*i+1 < n; i++){
 	  for(j = 1; (k = i+j+2*i*j) < n && j <= i; j++){
 		a[k] = 1;
 	  }
    }
    for(i = 1; i < n; i++){
 	  if(a[i] == 0) cout << (2 * i + 1) << " ";
    }
    delete [] a;
 }

Ассемблер

                К.Королев
                О соотношении Сундарама

Так называемое соотношение Сундарама

i _{m n} = 2mn + m + n , m, n - натуральные (1)

используется как решето простых чисел по принципу «как есть» и это соотношение действительно есть готовый оператор для любого языка программирования высокого уровня. Однако, если рассматривать это соотношение как симметричную матрицу, то решето можно реализовать более эффективно. Но сначала обратимся ко 2-му тому культовой монографии Д.Кнута «Искусство программирования», точнее, к упр.8 раздела 4.5.4. «Разложение на простые множители». В упражнении предлагается найти алгоритм высеивания, использующий только операции сложения и вычитания, а в ответе на него приведены соотношения, непосредственно вытекающие из (1): в оригинальных обозначениях монографии это p = 2j + 1 и q = 2j + 2j² соответственно. Сам же Д.Кнут происхождение этих выражений не объясняет, возможно, по причине весьма скромного рейтинга упражнения (задача на 15 – 20 минут), которое имеет, хотя тематически и не совсем логичное, продолжение в упражнении 15 раздела 5.2.3 3-го тома. Приведенное в ответе, цитата, «решение несколько замысловато: в нем пропускаются все числа, кратные 3». Решение предложено Б.Чартресом.

Однако, используя (1), можно построить более простой и эффективный алгоритм, пропускающий все числа, кратные 3. Вернее, с точки зрения матрицы (1) этот алгоритм пропускает все строки и столбцы матрицы, кратные 3-м. Приведем этот алгоритм, реализованный на ассемблере. Смысл используемых переменных следующий:

NOdd+1 – количество нечетных чисел, из которых надо отсеять составные;
Prime – однобайтовый массив такой, что если Prime[J – 1] = 1, то соответствующее простое

число вычисляется как 2J+1. Понятно, что отсутствует элемент массива, соответствующий тривиальной двойке.

SS_33-процедура с исключением строк и столбцов, «кратных» 3-ем

;Решето Сундарама: отсеивание из последовательности нечетных чисел всех составных чисел.
;Инициализация массива Prime, нулевой элемент которого, соответствующий 
;простому числу 3, равен 1. При этом отсеиваются все числа, кратные 3-ем. 
    mov    Byte Ptr  Prime[0],1       ; простое
    mov    eax,1
@@INIT:
    mov    Word Ptr  Prime[eax],0101h ; возможно, простые числа
    add    eax,2
    mov    Byte Ptr  Prime[eax],0 ; составное число, кратное 3-ём (исключение столбца)
    inc    eax
    cmp    eax,DWord Ptr NOdd
    jb     @@INIT
;====================== 
;Отсеивание остальных чисел
    mov    eax,12-1   ; второй диагональный элемент матрицы минус единица
    mov    ecx,5      ; период второй строки   
    mov    edx,10     ; для исключения столбца второй строки, «кратного» 3-ем
    mov    bl,2       ; счетчик для исключения  строк, «кратных» 3-ем
   ;=======  Движение по строке  =======
;Из-за «переопределенности» (1) следующие пересылки выполняются 
;многократно по одному и тому же адресу
@@S0:
    mov    edi,eax 	; сохранение диагонального элемента
@@S1:
    mov    Byte Ptr Prime[eax],0      ; составное число
    add    eax,ecx 
    cmp    eax,DWord Ptr NOdd
    jnb    @@S3 	; eax ≥ NOdd
@@S2:
    mov    Byte Ptr Prime[eax],0    ; составное число
    add    eax,edx                  ; исключение столбца, «кратного» 3-ем
    cmp    eax, DWord Ptr NOdd
    jb     @@S1	; eax < NOdd
@@S3:
    ja     @@S4	; eax > NOdd
    mov    Byte Ptr Prime[eax],0    ; последнее число в строке составное 
@@S4:
; ======= Переход к следующей строке =======
    mov    eax,edi	    ; восстановление диагонального элемента
    add    eax,ecx
    add    ecx,2     	    ; период следующей строки
    add    eax,ecx          ; следующий диагональный элемент минус единица
    mov    edx,ecx		     
    shl    edx,1            ; для исключения столбца очередной строки, «кратного» 3-ем
    dec    bl 
    jnz    @@S5            ; строка не «кратна» 3-ем
; Исключение строки, «кратной» 3-ем
    mov     bl,2                                                     
    add     eax,ecx
    add     ecx,2     	    ; период следующей строки
    add     eax,ecx         ; следующий диагональный элемент минус единица               
    mov     edx,ecx 
    shl     edx,1           ; для исключения столбца очередной строки, «кратного» 3-ем
    cmp     eax, DWord Ptr NOdd
    jb      @@S0                 ; на столбец, не «кратный» 3-ем
    ja      @@End 
    jz      @@S6
@@S5:
    mov     edi,eax 
    cmp     eax, DWord Ptr NOdd
    jb      @@S2                 ; исключить столбец, «кратный» 3-ем
    ja      @@End 
@@S6:
    mov     Byte Ptr Prime[eax],0    ; последнее нечетное число составное
@@End:

Приведенный алгоритм работает в 3-4 раза быстрее (по времени), например, алгоритма Аткина-Бернстейна. Кроме того, последний, так же, кстати, как и решение Чартреса, требует в два раза больше памяти, т.к. работает на всем множестве натуральных чисел. Несколько сложнее выглядят процедуры, исключающие дополнительно строки, «кратные» 5-ти и 7-ми соответственно, но они имеют еще большее преимущество.

Кстати, это не единственный результат, который можно получить из (1). Так, за 15-20 минут из (1) невозможно не получить бинарный алгоритм вычисления НОД (наибольшего общего делителя).

См. также

Решето Аткина

Решето Эратосфена