Черпаем простые числа решетом Аткина

Программирование - Практика программирования

Предложен способ увеличения быстродействия запроса поиска простых чисел с помощью "решета Аткина".

При обсуждении статьи "Как получить список простых чисел в запросе" среди прочих методов был предложен способ,известный как Решето Аткина. Данный алгоритм хорош тем, что позволяет получить список простых чисел в диапазоне от 2 до N в квадрате, используя начальное множество в диапазоне 1...N. То есть, чтобы получить все простые числа меньшие 10 000, нам надо взять исходный набор в диапазоне от 1 до 100 и на его основе получить  числа кандидаты, а затем  внутри построенного набора выбрать искомые простые числа. Получение чисел кандидатов состоит в вычислении квадратичных форм и определении остатков от деления на 12. Приведем эти формы:

Квадратичная форма Остаток от деления на 12
4*X*X+Y*Y 1,5
3*X*X+Y*Y 7
3*X*X-Y*Y 11

Основная идея построения запроса состоит в следующем. Формируем исходное множество чисел в диапазоне от 1 до N. Для  каждой пары чисел из данного набора получаем значение квадратичной формы и вычисляем для него остаток от деления на 12. Если остаток соответствует указанным условиям, то включаем найденное число в список кандидатов. Среди чисел кандидатов оставляем те, которые встречались нечетное число раз. Затем полученный набор проверяем на делимость на числа, не превышающие корень квадратный из тестируемого числа. Ниже приведен текст запроса Сергея (ildarovich)

 

ВЫБРАТЬ
    0 КАК Х
ПОМЕСТИТЬ Регистр1

ОБЪЕДИНИТЬ

ВЫБРАТЬ
    1
;
ВЫБРАТЬ
    Младшие.Х + 2 * Старшие.Х КАК Х
ПОМЕСТИТЬ Регистр2
ИЗ
    Регистр1 КАК Младшие,
    Регистр1 КАК Старшие
;
ВЫБРАТЬ
    Младшие.Х + 4 * Старшие.Х КАК Х
ПОМЕСТИТЬ Регистр4
ИЗ
    Регистр2 КАК Младшие,
    Регистр2 КАК Старшие
;
ВЫБРАТЬ
    Младшие.Х + 16 * Старшие.Х КАК Х
ПОМЕСТИТЬ Регистр8
ИЗ
    Регистр4 КАК Младшие,
    Регистр4 КАК Старшие
;
ВЫБРАТЬ
    1 + Младшие.Х + 256 * Старшие.Х КАК Х
ПОМЕСТИТЬ Ряд
ИЗ
    Регистр8 КАК Младшие,
    Регистр2 КАК Старшие
ГДЕ
    1 + Младшие.Х + 256 * Старшие.Х <= &КореньИзГраницы
;

////////////////////////////////////////////////////////////­////////////////////
ВЫБРАТЬ
    КвадратичныеФормы.Х,
    КвадратичныеФормы.Х * КвадратичныеФормы.Х КАК ХХ
ПОМЕСТИТЬ Кандидаты
ИЗ
    (ВЫБРАТЬ
        Случай1.Х КАК Х
    ИЗ
        (ВЫБРАТЬ
            4 * А.Х * А.Х + Б.Х * Б.Х КАК Х
        ИЗ
            Ряд КАК А,
            Ряд КАК Б) КАК Случай1
    ГДЕ
        СЕКУНДА(ДОБАВИТЬКДАТЕ(ДАТАВРЕМЯ(1, 1, 1), СЕКУНДА, Случай1.Х * 5)) В (5, 25)
    
    ОБЪЕДИНИТЬ ВСЕ
    
    ВЫБРАТЬ
        Случай2.Х
    ИЗ
        (ВЫБРАТЬ
            3 * А.Х * А.Х + Б.Х * Б.Х КАК Х
        ИЗ
            Ряд КАК А,
            Ряд КАК Б) КАК Случай2
    ГДЕ
        СЕКУНДА(ДОБАВИТЬКДАТЕ(ДАТАВРЕМЯ(1, 1, 1), СЕКУНДА, Случай2.Х * 5)) = 35
    
    ОБЪЕДИНИТЬ ВСЕ
    
    ВЫБРАТЬ
        Случай3.Х
    ИЗ
        (ВЫБРАТЬ
            3 * А.Х * А.Х - Б.Х * Б.Х КАК Х
        ИЗ
            Ряд КАК А
                ВНУТРЕННЕЕ СОЕДИНЕНИЕ Ряд КАК Б
                ПО А.Х > Б.Х) КАК Случай3
    ГДЕ
        СЕКУНДА(ДОБАВИТЬКДАТЕ(ДАТАВРЕМЯ(1, 1, 1), СЕКУНДА, Случай3.Х * 5)) = 55) КАК КвадратичныеФормы
ГДЕ
    СЕКУНДА(ДОБАВИТЬКДАТЕ(ДАТАВРЕМЯ(1, 1, 1), СЕКУНДА, КвадратичныеФормы.Х * 12)) > 0
    И КвадратичныеФормы.Х <= &Граница

СГРУППИРОВАТЬ ПО
    КвадратичныеФормы.Х

ИМЕЮЩИЕ
    СЕКУНДА(ДОБАВИТЬКДАТЕ(ДАТАВРЕМЯ(1, 1, 1), СЕКУНДА, КОЛИЧЕСТВО(*) * 30)) = 30
;

////////////////////////////////////////////////////////////­////////////////////
ВЫБРАТЬ
    Делимое.Х КАК Х
ИЗ
    Кандидаты КАК Делимое
        ЛЕВОЕ СОЕДИНЕНИЕ Кандидаты КАК Делитель
        ПО (Делитель.ХХ <= Делимое.Х)
            И ((ВЫРАЗИТЬ(Делимое.Х / Делитель.Х КАК ЧИСЛО(15, 0))) = Делимое.Х / Делитель.Х)
ГДЕ
    Делитель.Х ЕСТЬ NULL

ОБЪЕДИНИТЬ

ВЫБРАТЬ
    2

ОБЪЕДИНИТЬ

ВЫБРАТЬ
    3

ОБЪЕДИНИТЬ

ВЫБРАТЬ
    5

УПОРЯДОЧИТЬ ПО
    Делимое.Х
R03;

 

Первое, что бросается в глаза это использование функции СЕКУНДА для проверки делимости на 12. Действительно,  если число  записать в виде N=12*А+B, где B < 12, тогда 5*N=60*A+5*B. Следовательно остаток от деления числа 5*N на 60, равен остатку от деления числа N на 12 умноженному на 5. Отсюда ,например, условие, что N mod 12 = 7  равносильно условию СЕКУНДА( ДОБАВИТЬКДАТЕ( ДАТАВРЕМЯ(1,1,1),5*N))=35. Данная конструкция имеет свои ограничения и  об этом Сергей пишет в обсуждении данного способа.

Да, но и с секундами будет когда-нибудь ругаться. Там предел 2^32 - 1 или 2^31 - 1. Это примерно миллиард. Надеюсь, никто не додумается из миллиарда простые числа в запросе определять. Так как будет уже не с временем (оно почти линейно растет, поэтому миллиард за 16000 секунд должен посчитаться), а с памятью проблемы.

Зададимся вопросом, как можно ускорить данный алгоритм и уйти от ограничений ? Предлагается следующий подход.

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

ВЫБРАТЬ
	0 КАК Х,
	0 КАК Остаток
ПОМЕСТИТЬ Регистр1

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	1,
	1

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	2,
	2

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	3,
	3

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	4,
	4

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	5,
	5

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	6,
	6

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	7,
	7

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	8,
	8

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	9,
	9

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	10,
	10

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	11,
	11
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	Младшие.Х + 12 * Старшие.Х КАК Х,
	Младшие.Остаток КАК Остаток
ПОМЕСТИТЬ Регистр2
ИЗ
	Регистр1 КАК Младшие,
	Регистр1 КАК Старшие
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	Младшие.Х + 144 * Старшие.Х КАК Х,
	Младшие.Остаток КАК Остаток
ПОМЕСТИТЬ Регистр4
ИЗ
	Регистр2 КАК Младшие,
	Регистр2 КАК Старшие
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	Младшие.Х + 1728 * Старшие.Х КАК Н,
	Младшие.Остаток КАК Остаток
ПОМЕСТИТЬ Числа
ИЗ
	Регистр4 КАК Младшие,
	Регистр1 КАК Старшие
ГДЕ
	Младшие.Х + 1728 * Старшие.Х <= &КореньИзГраницы
	И Младшие.Х + Старшие.Х <> 0

ИНДЕКСИРОВАТЬ ПО
	Н
;

////////////////////////////////////////////////////////////////////////////////
УНИЧТОЖИТЬ Регистр1
;

////////////////////////////////////////////////////////////////////////////////
УНИЧТОЖИТЬ Регистр2
;

////////////////////////////////////////////////////////////////////////////////
УНИЧТОЖИТЬ Регистр4

Хорошо. Мы сформировали список чисел и для каждого числа знаем его значение по модулю 12. А как получить остаток для квадратичной формы ? И здесь тоже все очень просто.  Вычислим значения квадратичной формы только для чисел в диапазоне  0...11. Для каждого из значений найдем остаток от деления на 12 и оставим только те, которые удовлетворяют нашим условиям. Остаток от деления на 12  будем вычислять по формуле

КвадратичнаяФорма- (ВЫРАЗИТЬ(КвадратичнаяФорма/12 как ЧИСЛО(5,0)))*12

При проверке условия на значение остатка учитывается возможность округления как в меньшую, так и большую сторону. Поэтому, например, условие, что остаток по модулю 12 равен 7, записывается как ОСТАТОК В (-5,7) . Ниже приводится текст запроса.

ВЫБРАТЬ
	0 КАК Х,
	0 КАК КвадратХ
ПОМЕСТИТЬ мОстатки

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	1,
	1

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	2,
	4

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	3,
	9

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	4,
	16

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	5,
	25

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	6,
	36

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	7,
	49

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	8,
	64

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	9,
	81

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	10,
	100

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	11,
	121
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	4 * А.КвадратХ + Б.КвадратХ - (ВЫРАЗИТЬ((4 * А.КвадратХ + Б.КвадратХ) / 12 КАК ЧИСЛО(5, 0))) * 12 КАК Остатки,
	А.Х КАК А,
	Б.Х КАК Б,
	12 * А.Х + Б.Х КАК АБ
ПОМЕСТИТЬ мОстаткиТ1
ИЗ
	мОстатки КАК А,
	мОстатки КАК Б
ГДЕ
	4 * А.КвадратХ + Б.КвадратХ - (ВЫРАЗИТЬ((4 * А.КвадратХ + Б.КвадратХ) / 12 КАК ЧИСЛО(5, 0))) * 12 В (1, -11, 5, -7)

ИНДЕКСИРОВАТЬ ПО
	АБ
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	3 * А.КвадратХ + Б.КвадратХ - (ВЫРАЗИТЬ((3 * А.КвадратХ + Б.КвадратХ) / 12 КАК ЧИСЛО(5, 0))) * 12 КАК Остатки,
	А.Х КАК А,
	Б.Х КАК Б,
	12 * А.Х + Б.Х КАК АБ
ПОМЕСТИТЬ мОстаткиТ2
ИЗ
	мОстатки КАК А,
	мОстатки КАК Б
ГДЕ
	3 * А.КвадратХ + Б.КвадратХ - (ВЫРАЗИТЬ((3 * А.КвадратХ + Б.КвадратХ) / 12 КАК ЧИСЛО(5, 0))) * 12 В (-5, 7)

ИНДЕКСИРОВАТЬ ПО
	АБ
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	3 * А.КвадратХ - Б.КвадратХ - (ВЫРАЗИТЬ((3 * А.КвадратХ - Б.КвадратХ) / 12 КАК ЧИСЛО(5, 0))) * 12 КАК Остатки,
	А.Х КАК А,
	Б.Х КАК Б,
	12 * А.Х + Б.Х КАК АБ
ПОМЕСТИТЬ мОстаткиТ3
ИЗ
	мОстатки КАК А,
	мОстатки КАК Б
ГДЕ
	3 * А.КвадратХ - Б.КвадратХ - (ВЫРАЗИТЬ((3 * А.КвадратХ - Б.КвадратХ) / 12 КАК ЧИСЛО(5, 0))) * 12 В (-1, 11)

ИНДЕКСИРОВАТЬ ПО
	АБ
;

////////////////////////////////////////////////////////////////////////////////
УНИЧТОЖИТЬ мОстатки

На этом шаге мы формируем временные таблицы мОстаткиТ1,мОстаткиТ2,мОстаткиТ3 каждая  для своей квадратичной формы. Теперь когда проведена подготовительна работа получим список чисел кандидатов. Здесь надо остановиться на следующем моменте. При вычислении значения квадратичной формы мы одновременно вычисляем значение индекса 12*(X mod 12)+ (Y mod 12) , в запросе это поле  12 * X.Остаток + Y.Остаток КАК АБ.  В дальнейшем мы сравниваем этот индекс с индексом из подготовленных таблиц и ,если они совпадают, то значение квадратичной формы  включается в  список чисел кандидатов.
 

ВЫБРАТЬ
	Числа.Н КАК Н,
	Числа.Остаток КАК Остаток,
	Числа.Н * Числа.Н КАК КвадратН
ПОМЕСТИТЬ ВТ
ИЗ
	Числа КАК Числа
ГДЕ
	Числа.Н <= &КореньИзГраницы

ИНДЕКСИРОВАТЬ ПО
	Н
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	КвадратичнаяФорма.Н КАК Н,
	КОЛИЧЕСТВО(*) КАК КолВо
ПОМЕСТИТЬ Т
ИЗ
	(ВЫБРАТЬ
		4 * X.КвадратН + Y.КвадратН КАК Н,
		12 * X.Остаток + Y.Остаток КАК АБ
	ИЗ
		ВТ КАК X
			ВНУТРЕННЕЕ СОЕДИНЕНИЕ ВТ КАК Y
			ПО (4 * X.КвадратН + Y.КвадратН <= &Граница)) КАК КвадратичнаяФорма
		ВНУТРЕННЕЕ СОЕДИНЕНИЕ мОстаткиТ1 КАК мОстаткиТ1
		ПО КвадратичнаяФорма.АБ = мОстаткиТ1.АБ

СГРУППИРОВАТЬ ПО
	КвадратичнаяФорма.Н

ИМЕЮЩИЕ
	КОЛИЧЕСТВО(*) В (1, 3, 5, 7, 9, 11)

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	КвадратичнаяФорма.Н,
	КОЛИЧЕСТВО(*)
ИЗ
	(ВЫБРАТЬ
		3 * X.КвадратН + Y.КвадратН КАК Н,
		12 * X.Остаток + Y.Остаток КАК АБ
	ИЗ
		ВТ КАК X
			ВНУТРЕННЕЕ СОЕДИНЕНИЕ ВТ КАК Y
			ПО (3 * X.КвадратН + Y.КвадратН < &Граница)) КАК КвадратичнаяФорма
		ВНУТРЕННЕЕ СОЕДИНЕНИЕ мОстаткиТ2 КАК мОстаткиТ2
		ПО КвадратичнаяФорма.АБ = мОстаткиТ2.АБ

СГРУППИРОВАТЬ ПО
	КвадратичнаяФорма.Н

ИМЕЮЩИЕ
	КОЛИЧЕСТВО(*) В (1, 3, 5, 7, 9, 11)

ОБЪЕДИНИТЬ

ВЫБРАТЬ
	КвадратичнаяФорма.Н,
	КОЛИЧЕСТВО(*)
ИЗ
	(ВЫБРАТЬ
		3 * X.КвадратН - Y.КвадратН КАК Н,
		12 * X.Остаток + Y.Остаток КАК АБ
	ИЗ
		ВТ КАК X
			ВНУТРЕННЕЕ СОЕДИНЕНИЕ ВТ КАК Y
			ПО X.Н >= Y.Н
				И (3 * X.КвадратН - Y.КвадратН < &Граница)) КАК КвадратичнаяФорма
		ВНУТРЕННЕЕ СОЕДИНЕНИЕ мОстаткиТ3 КАК мОстаткиТ3
		ПО КвадратичнаяФорма.АБ = мОстаткиТ3.АБ

СГРУППИРОВАТЬ ПО
	КвадратичнаяФорма.Н

ИМЕЮЩИЕ
	КОЛИЧЕСТВО(*) В (1, 3, 5, 7, 9, 11)

 

И наконец, последний участок кода, в котором находится количество простых чисел не превышающих заданную границу.

ВЫБРАТЬ
	Т.Н КАК Н
ПОМЕСТИТЬ мДляПроверки
ИЗ
	Т КАК Т

ИНДЕКСИРОВАТЬ ПО
	Н
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	мДляПроверки.Н
ПОМЕСТИТЬ мПростыеЧисла
ИЗ
	мДляПроверки КАК мДляПроверки
		ВНУТРЕННЕЕ СОЕДИНЕНИЕ ВТ КАК мДелители
		ПО (мДелители.КвадратН <= мДляПроверки.Н)
			И ((ВЫРАЗИТЬ(мДляПроверки.Н / мДелители.Н КАК ЧИСЛО(19, 0))) = мДляПроверки.Н / мДелители.Н)

СГРУППИРОВАТЬ ПО
	мДляПроверки.Н

ИМЕЮЩИЕ
	КОЛИЧЕСТВО(мДляПроверки.Н) = 1
;

////////////////////////////////////////////////////////////////////////////////
ВЫБРАТЬ
	КОЛИЧЕСТВО(мПростыеЧисла.Н) + 2 КАК НатуральноеЧисло
ИЗ
	мПростыеЧисла КАК мПростыеЧисла

Теперь вполне законный вопрос, а стоит ли овчинка выделки ? Дает ли предложенный подход выигрыш в скорости ? В следующей таблице приведены результаты замера времени выполнения запросов  (в миллисекундах)   для различных значений верхней границы области поиска. 

верхняя граница 1000 10 000 50 000 100 000 300 000 1 000 000
Валерий (scientes) 423 599 1 569 5 040 14 447 49 434
Сергей (ildarovich) 85 473 2 978 7 069 40 517 348 492

При верхней границе диапазона поиска примерно равной 30 000 авторский вариант начинает работать быстрее и по мере роста границы разрыв увеличивается. Но это, как видно, достигается за счет увеличения объема кода и подготовительных вычислений.  К статье прилагается обработка (управляемая форма), которую я использовал для измерения скорости выполнения различных вариантов запросов. 

Поскольку в комментариях, как обычно, много вопросов "Зачем это нужно и где может пригодиться ?", то повторю старую мысль. Главное не инструмент, а идеи и алгоритмы. Например, мы видим ,как постоянно появляются новые языки программирования и оболочки. Если стараться изучить все это многообразие, то жизнь, скорее всего, пройдет мимо. Быть лидером в этой гонке невозможно. Поэтому надо знать и овладевать тем, что не подвержено  столь быстрым изменениям и что лежит в основе предметной области, которой мы занимаемся. А это идеи и алгоритмы. Кстати, в статье на Википедии приведен код простой программы на C и Pascal для получения списка простых чисел методом решета Аткина. Данный код  может быть оптимизирован путем применения критерия отнесения значения квадратичной формы к числам кандидатам, как это предложено в статье. Впрочем, я полагаю, что  в других реализациях это уже сделано.

Скачать файлы

Наименование Файл Версия Размер
SieveOfAtkin
.epf 20,04Kb
10.01.18
0
.epf 20,04Kb Скачать

См. также

Комментарии
1. Сергей Ожерельев (Поручик) 4066 09.01.18 21:19 Сейчас в теме
Только я один с ужасом читаю статьи по мотивам статей ildarovich'а?
корум; fancy; Evil Beaver; maxdmt; +4 Ответить
2. Петр Малыгин (pm74) 61 09.01.18 21:56 Сейчас в теме
(1) Всегда читал его статьи с интересом. Кстати куда он пропал ?
Somebody1; bulpi; Evil Beaver; +3 Ответить
3. Андрей Овсянкин (Evil Beaver) 4911 10.01.18 01:38 Сейчас в теме
4. Сергей Блинов (BlinOFF) 10.01.18 07:38 Сейчас в теме
я может быть фигню спрошу, но зачем в 1ске простые числа? Просто из спортивного интереса?
5. Александр Лаптев (SaschaL) 10.01.18 07:58 Сейчас в теме
Статься написана интересно. Но вот пока не могу сообразить куда её применить на практике
корум; +1 Ответить
6. Павел Князев (viking7) 10.01.18 13:15 Сейчас в теме
Интересная статья! Классика жанра!)
7. Максим Сагайдак (maximsagaydak) 7 11.01.18 18:25 Сейчас в теме
Интересные изыскания, но не вполне ясно для чего! :)
Оставьте свое сообщение