Реализации алгоритмов/Решето Сундарама
Решето́ Сундара́ма — детерминированный алгоритм нахождения всех простых чисел до некоторого целого числа .
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
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 + 2j2 соответственно. Сам же Д.Кнут происхождение этих выражений не объясняет, возможно, по причине весьма скромного рейтинга упражнения (задача на 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) невозможно не получить бинарный алгоритм вычисления НОД (наибольшего общего делителя).