Начну пожалуй про модулярную арифметику в приложении к поиску последовательностей простых по заданному паттерну (расстоянию между соседними простыми). И для примера возьму цепочку из 4-х простых с интервалами 2,4,2, или в другой записи 0,2,6,8 - две максимально близкие пары простых близнецов. Как
Yadryara объяснил в другой теме, из каждых 30 последовательных чисел подходит лишь одно, с остатком 11 от деления на 30 (только с такого числа может начинаться искомая цепочка). Руками не считал, нагрузил расчётом PARI:
? chinese([Mod(1,2),Mod(2,3),Mod(1,5)])
%1 = Mod(11, 30)Первые решения:
C:\>primesieve.x64con.exe 0 -d1000 -p4
(5, 7, 11, 13)
(11, 13, 17, 19)
(101, 103, 107, 109)
(191, 193, 197, 199)
(821, 823, 827, 829)
Первая строка отпадает - это исключение, остаток по модулю 30 другой, но такие исключения могут быть только до числа 30, не дальше, на это плюнем.
Так как между искомыми простыми никаких других простых быть не может в принципе, то можно просто идти по всем натуральным числам с 11 с шагом 30 и проверять их на простоту. Если простое - проверять остальные три числа на простоту. Для чисел около 4млрд (доступных нам пока что для проверок на простоту) простым является примерно каждое 11 нечётное число (
), а мы будем проверять лишь каждое 15-е из них, то в среднем первое число цепочки окажется простым один раз из 165. Разумеется это намного медленнее решета Эратосфена, но это же чисто модельный пример.
Все пояснения про цепочки, паттерны, допустимые остатки, модули - всё в другой теме (например указанной
Yadryara выше), здесь займёмся приложением математики к ассемблеру.
Вспомним что значит простота числа: это значит что число не делится на простые меньшие него (корня, не корня - не суть). Попробуем отказаться от проверки делимости первого числа цепочки на 7 (на меньшие простые не делится гарантированно и их проверять не надо). Число делится на простое если остаток от деления равен 0. Но остаток 11%7=4 мы знаем (или можем легко посчитать один раз). Раз он не 0, то и 11 на 7 не делится (кто бы сомневался).
А теперь фокус: а будет ли делиться на 7 следующее проверяемое число? А какое оно? Ну так 11+30=41. Конечно всем с 3-го класса очевидно что не делится. А компьютер в школе не учили и ему не очевидно! И делить мы не хотим. Сделаем так: 41%7=(11+30)%7=(11%7+30%7)%7. Казалось бы что в этом хорошего, одно деление превратили в три? Но ведь не в три, 11%7 мы и так уже знали на предыдущей итерации перебора, это 4. 30%7=2 мы тоже знаем (или можем легко посчитать), и оно не будет меняться ни для каких проверяемых нами чисел, значит его можно вычислить один раз в начале. Остаётся последнее взятие остатка. Но смотрим внимательно, в скобках оба слагаемых меньше модуля, их сумма уж всяко меньше удвоенного модуля, значит делить не обязательно, достаточно сравнить с 7 и если превышает, то вычесть 7 чтобы результат стал всегда меньше 7 (модуля), и это гораздо проще чем делить.
При повторении прибавления шага 30 снова сработает ровно тот же фокус.
Приходим к следующей проверке делимости первого числа цепочки на 7:
x=11; m=30;
x7=x%7; m7=m%7;
for(i=1; i<33; ++i) {//Все числа <30 пропустим и начнём сразу дальше
x7=x7+m7;
if(x7>=7) x7=x7-7;
if(x7!=0) printf("First number is not div by 7 at i= %u",i);
}
Ну, начальная инициализация и цикл нам не интересны, как и printf, речь про две строки вычисления x7 без делений.
Вариант на асме вычисления x7 (в ESI), с условным переходом должен быть понятен всем:
;Вариант с условным переходом
add ESI,m7 ;x7=x7+m7
cmp ESI,7
jc .nsub
sub ESI,7
.nsub:
;Вариант без перехода
add ESI,m7 ;x7=x7+m7
mov EAX,ESI
sub EAX,7
cmovnc ESI,EAX
А вот второй вариант интереснее, есть лишняя mov (в таком коде по идее совершенно бесплатная) и непонятная
cmovnc - а она скопирует EAX в ESI только если перенос будет =0 (после cmp, а значит ESI>=7), иначе не сделает ничего. Тут не только не нужна метка (программисты все лентяи), но и работать должно быстрее - любой условный переход заметно тормозит программу если его нельзя предсказать более чем в 90% случаев (будет он выполнен или нет), а тут его нет. В данном конкретном случае у нас m7=2 и переход будет выполняться с вероятностью 5/7 (лишь два остатка из 7 возможных после прибавления двойки дадут число больше или равное 7 и переход выполнен не будет). С такой вероятностью предсказание перехода будет ошибаться раза 4 из 7, т.е. даже больше половины, и каждый раз будет тормозить программу. А команда cmov программу не тормозит. На моём процессоре (от 2013г) этот код должен выполняться за 1 такт, все 4 команды. Вместо тактов 10 для команды div.
Фактически связка команд mov (бесплатная) + cmp + cmov превращает предыдущее сложение в сложение по модулю.
Ещё причина избавиться от условного перехода: его нет в MMX/SSE/AVX, а мы целимся их применить для ускорения, ну если получится.
Итак, при последовательном переборе с неким постоянным шагом мы можем вычислять остатки по модулю простых достаточно быстро (и без делений).Дальше мы можем с остатком делать что угодно: хоть проверить на 0 и если да, то остальные числа в цепочке не проверять (первое же кратно 7 и следовательно не простое), хоть командами bt + jnc взять признак разрешённого конкретно такого остатка по модулю простого для всей цепочки (их 3шт: 2,3,4) и остальные числа на делимость на 7 вообще не проверять, хоть ещё что-нибудь забавное.
А в SSE4.1 и AVX есть команда pminu, для беззнаковых чисел, выдающая минимум из двух чисел. У нас же при вычитании sub EAX,7 может получиться отрицательное число (если ESI<7 после сложения), а отрицательные числа имеют 1 в старшем бите и уж точно больше (беззнаково!) любого положительного числа, даже удвоенного модуля (в данном случае 14). Т.е. команду cmovnc ESI,EAX можно заменить на команду pminu ESI,EAX (не прямо такую, но очень похожую). Если в ESI число больше или равно 7, то EAX положителен и меньше ESI (он же на 7 меньше) и тогда надо взять именно EAX (минимум), если же EAX<0 (и старший бит 1), то ESI<7 и в беззнаковой форме уже ESI меньше беззнакового EAX и команда минимума выдаст именно ESI.
Зачем менять одну команду на другую? Так в AVX нет команды cmov (pblendvb только байтовая и вчетверо медленнее pminu, фтопку), зато каждая команда может работать сразу с 16-ю числами (для модулей <32768)! А для модулей до 128 и сразу с 32-я числами! Т.е. одновременно можно вычислять сразу 32 новых остатка по малому простому модулю. Все три (mov в AVX не нужна) команды за 1.5 такта. Сразу 10 (20 для модулей до 128) чисел за такт, в 10 (20) раз быстрее показанного выше кода! И кстати в 100 (200) раз быстрее варианта с div.
А целочисленных делений в SSE/AVX нет. Связываться же с плавающей точкой: 7 тактов на 4 числа в SSE и 14 тактов на 8 чисел в AVX, для чисел одинарной точности, т.е. для целых лишь до примерно 16млн - вообще неактуально. С двойной точностью вдвое меньше чисел и деление вдвое медленнее: 14 тактов на 2 числа в SSE и 28 тактов на 4 числа в AVX, для целых чисел примерно до 9e15 - ещё более неактуально.
-- 24.01.2024, 02:08 --Ещё одно огромное преимущество данного подхода: делить длинное/большое исходное число на модуль надо лишь один раз при инициализации. И сколько это займёт времени не так важно, конечно если потом нам надо не тысячи итераций, а сильно поболее (ради тысячи проще на чём угодно посчитать, хоть на калькуляторе), миллиарды и триллионы.
А раз делить нам внутри циклов нигде не надо, и вообще нигде при проверке на простоту не появляются исходные длинные/большие числа и все вычисления идут лишь по модулю небольших простых, то и в каком диапазоне работает проверка
в принципе вообще не важно! Хоть около 4млрд, хоть около 97e9999, главное правильно всё инициализировать и потом спокойно считать миллиарды итераций, с одинаковой скоростью. Фактически счёт в описанном варианте идёт просто по модулю 210 (и кроме шага 30 все числа никогда не превышают 7*2=14), а уж где именно в натуральном ряду расположен этот кусок совершенно не важно (для вычислений остатков на следующей итерации из известных на текущей).
Да, ещё. Вычисление x=x+30 из цикла специально убрал: пока не нашли решение он нам и не нужен, работаем только с остатками, а когда нашли - вот тогда и посчитаем, например как z=11+30*(i-1) или z=x+m*(i-1).
, но ... реально это мало кому нужно. А когда понадобится - тогда и изучат, по мере необходимости и в разумном объёме.
И кое-что нашёл, может поэтому и теоретическая скорость поднялась, надо дописать и померить реальную.
. Т.е. деление самое медленное, умножение лишь немногим медленнее почти всех прочих команд. По мере накопления опыта можно будет уточнять скорость каждой команды.