HuzHugo, похоже я так и не ответил Вам об устройстве ускорителей и возможности их применения для других паттернов (patterns). Отвечу здесь, возможно это будет интересно и другим.
Предварительные замечания.
Задача в нахождении цепочки чисел с постоянным количеством делителей. Поделим числа в цепочке на два класса: проверяемые (checked) и непроверяемые (nonchecked), в первых осталось найти только одно большое простое в первой степени, во вторых простых может быть два и более. Ускорять внешней программой будем лишь поиск чисел первого класса, остальную часть работы оставим за
PARI/GP, который заодно и перепроверит и найденные числа первого класса.
В соответствии с
КТО первое число в цепочке не может быть произвольным, а должно соответствовать форме
, где
это некое начальное значение в зависимости от конкретного расположения малых простых в паттерне цепочки, а
это модуль/шаг. Таким образом достаточно перебирать числа
для всех
и проверять количество делителей чисел
.
Для проверяемых (checked) чисел количество делителей однозначно определяется простотой большого числа, равного
, где
это коэффициент с ним, а
это номер проверяемого числа в цепочке (с нуля). Таким образом поиск кандидата в цепочку свёлся к поиску некоторого количества больших простых чисел
в строгом соотношении друг с другом. Найденный кандидат в цепочку потом будет допроверен в PARI/GP и только там станет (или не станет) найденной цепочкой.
Что делают мои ускорители.
В PARI/GP, из-за того что это интерпретатор, все операции довольно медленны, потому основной идеей ускорителей было замена внутреннего цикла перебора
на внешнюю программу, которая возвращает вероятно подходящие
, относительно быстро допроверяемые потом в PARI/GP.
Для этого достаточно передать внешней программе величины
и
, а она найдёт вероятно подходящие
. Для деления работы на кусочки можно передавать и диапазон
в котором искать вероятно подходящие
. В текущей реализации ускорителей величины
и
жёстко зашиты в исходный код программы и передаются только
.
Так как все
являются гарантированно целыми (это гарантировала КТО), то при переходе от
к
величины
изменятся всегда на одно и то же целое значение (каждая на своё). Потому вместо вычисления
на каждом шаге и для каждого проверяемого числа можно инициализировать вектор проверяемых чисел начальными значениями
и потом лишь прибавлять к нему величину
. Это эффективно реализуется в SSE/AVX одной командой векторного сложения.
Далее, для проверки чисел
на простоту используется крайне простой алгоритм проверки делимости на малые простые. Но вместо делений
(где
это малые простые до некоего порога) используется арифметика по модулую
, т.е. вместо
хранятся
и при переходе от
к
они все пересчитываются по формуле
. Как только ни одно из чисел
не стало равным
(равенство означает что число делится на соответствующее
) — значит нашли кандидата в цепочку. Операция суммирования по модулю реализуется тремя векторными командами AVX (в SSE шестью).
Итак, поиск кандидата в цепочку свёлся к правильной инициализации нескольких таблиц и потом циклу суммирования чисел по модулю малых простых с проверкой обнуления каждого из чисел. Так как простые малы, то достаточно отвести под каждое одно слово (word, 16 битов), что позволит проверять делимость кандидатов в цепочки на простые до
(почему не
поясню ниже). До этого порога 3512 простых, или меньше 220 слов AVX по 256 битов каждое. Но это только на одно проверяемое число в цепочке, для 11-ти проверяемых чисел слов AVX будет уже
или всего 75КБ.
Основной цикл перехода от
к
(порты выполнения и времена указаны для архитектуры Haswell):
; Ports Time
.add0: mov EDI,0 ;0156 1/0.25
.add: vmovdqu ymm15,[save_ymm+EDI] ;23 3/0.5
vpaddw ymm15,ymm15,[steps+EDI] ;15 1/0.5
vpsubw ymm14,ymm15,[primes+EDI] ;15 1/0.5
vpminuw ymm15,ymm14,ymm15 ;15 1/0.5 текущее значение остатков по модулю простых
vmovdqu [save_ymm+EDI],ymm15 ;4 1/1
add EDI,32 ;0156 1/0.25
cmp EDI,stop1 ;0156 1/0.25
jc .add ;06 1/0.5
В save_ymm сидят текущие значения
, без разделения к какому из проверяемых чисел и какому
оно относится, это деление выполняется внешней программой подготовки таблиц, для ускорителя это деление по индексам совершенно неважно, он проверяет факт обнуления любого из этих чисел.
В steps сидят величины
, т.е. приращения чисел
при переходе от
к
.
В primes сидят значения малых простых, по модулю которых и производятся вычисления.
После сложения с приращением вычитаем из результата величину модуля (малого простого) и берём минимум из них, при этом если сумма превышает простое, то будет взято уменьшенное значение (т.е. как раз по модулю простого), если же сумма не превышает простое, то результат вычитания станет больше суммы и беззнаковый (unsigned) минимум возьмёт сумму. Именно в этом месте появляется ограничение на величину малых простых, чтобы сумма не переполняла отведённое поле и правильно работала команда беззнакового (unsigned) минимума.
Можно было бы сразу здесь же проверить результат на ноль в любой позиции и если он там есть, то переходить к следующему
, но из-за принятого решения заменить умножения
на сложения
необходимо всегда подсчитывать все остатки
, даже если уже понятно что текущее
не подходит.
Это обходится частичным разворачиванием цикла (loop unrolling) по
по длине
и дублированием таблиц в памяти
раз. При этом можно сразу выкинуть из таблиц те
, которые точно не подойдут, например по модулю некоторых простых (если взять
таким простым), что позволит ещё немного ускорить вычисления. При этом самый внутренний цикл проверяет лишь неравенство нулю чисел
и при обнаружении нуля в любой позиции сразу же переходит к следующему
, не выполняя увеличения
, главное правильно инициализировать все таблицы, что выполняет внешняя программа, ускоритель же пользуется уже готовыми таблицами. В результате полное вычисление всех
выполняется в
раз реже чем проверка их на нули. В таблицах при этом сидят значения для
(если вдруг они все допустимы).
Итак, самый внутренний цикл проверки
на нули:
.next: inc R15 ;0156 1/0.25 k++
add RSI,const1 ;0156 1/0.25
cmp RSI,stop2 ;0156 1/0.25
jnc .add0 ;06 1/0.5
mov EDI,0 ;0156 1/0.25
.check: vmovdqu ymm15,[save_ymm+EDI] ;23 3/0.5
vpcmpeqw ymm15,ymm15,[RSI+EDI] ;15 1/0.5
vptest ymm15,ymm15 ;0+5 2/1
jnz .next ;06 1/0.5 если нули есть, то что-то где-то разделилось и переходим к следующему
add EDI,32 ;0156 1/0.25
cmp EDI,stop1 ;0156 1/0.25
jc .check ;06 1/0.5
На самом деле значения
проверяются не на нули, а на значения дающие ноль после суммирования с
, так как после разворачивания цикла (по
на длину
, в RSI начало таблицы для соответствующего
) при каждом значении RSI числа в
должны быть разными чтобы при суммировании давать ноль.
Полный кусок цикла поиска подходящего
:
; Ports Time
.add0: mov EDI,0 ;0156 1/0.25
.add: vmovdqu ymm15,[save_ymm+EDI] ;23 3/0.5
vpaddw ymm15,ymm15,[steps+EDI] ;15 1/0.5
vpsubw ymm14,ymm15,[primes+EDI] ;15 1/0.5
vpminuw ymm15,ymm14,ymm15 ;15 1/0.5 текущее значение остатков по модулю простых
vmovdqu [save_ymm+EDI],ymm15 ;4 1/1
add EDI,32 ;0156 1/0.25
cmp EDI,stop1 ;0156 1/0.25
jc .add ;06 1/0.5
mov RSI,table1 ;0156 1/0.25
jmp .check ;06 1/0.5
.next: inc R15 ;0156 1/0.25 k++
add RSI,size_table1 ;0156 1/0.25
cmp RSI,stop2 ;0156 1/0.25 а не достигли ли конца разворота внутреннего цикла
jnc .add0 ;06 1/0.5 если достигли, уходим на пересчёт всех x_i,j
mov EDI,0 ;0156 1/0.25
.check: vmovdqu ymm15,[save_ymm+EDI] ;23 3/0.5
vpcmpeqw ymm15,ymm15,[RSI+EDI] ;15 1/0.5 проверка на (-km/v_i,j)%p_j - т.е. даст ли ноль после суммирования для текущего k
vptest ymm15,ymm15 ;0+5 2/1
jnz .next ;06 1/0.5 если нули есть, то что-то где-то разделилось и переходим к следующему k
add EDI,32 ;0156 1/0.25
cmp EDI,stop1 ;0156 1/0.25
jc .check ;06 1/0.5 цикл по j проверки делимости на p_j
.ok: ;здесь должен быть вывод найденного k
jmp .next
В table1 сидят значения
, т.е. те, которые дадут ноль при суммировании для текущего
.
Порядок расположения кусков кода специально подобран для лучшего предсказания условных переходов при наиболее вероятном пути выполнения.
Собственно это всё, остальное вынесено в генератор таблиц для асм кода.
Ну и пару мелких оптимизаций (кэширование в регистрах ymm первых значений save_ymm, разворачивание первого десятка итераций цикла .check, объединение нескольких vpcmpeqw в группу с одним vptest+jnz после неё, использование 8-ми битных простых вместо 16-ти битных) не показал так как они существенно усложняют код не давая принципиального выигрыша.
PS. Асм код правильный, а вот мои пояснения к нему могут быть не везде точными, особенно про индексы, кажется я в них запутался.
PPS. По скорости: наиболее часто выполняются несколько (1-2-3-5) итераций цикла .check, потом уход на .next, и лишь один из
раз происходит уход на вычисление всех
в .add0. В итоге, для 11-ти проверяемых чисел и порога малых простых до 3584 вместо 32768 средняя скорость составила 1.5 варианта
за такт или примерно 5 млрд
в секунду в одном потоке (thread).