Как писать быстрые программы

[Yadryara]

VAL в сообщении #1703736 писал(а):

решил отписаться

Правильно. А то я уж подумал, что решили с grisа пример взять...

Меня спрашивают — я отвечаю. Разве уж на совсем личные вопросы могу не ответить.

VAL в сообщении #1703736 писал(а):

Василий исследование (а также школу, МФТИ и магистратуру) закончил,

Да, это в корне меняет дело... Я же ведь именно об этом и спрашивал, не правда ли.

VAL в сообщении #1703736 писал(а):

Эти пятерки я искал на медленном компе, используя медленный maple. А времени потратил, разумеется, не помню сколько.

Это и несущественно сколько потратили на медленном. Желательно запускать проги на одном и том же компе и сравнивать.

Если тот код у Вас сохранился и Вы его уже на PARI переписали, давайте запустим и попробуем сравнить.

Если на повестке дня 24 и 48 делителей, то с какого перепуга я занимался 4 и 6 делителями?

[gris]

Yadryara, я не прячусь под кроватью от стыда, но совсем мало уделяю внимания интереснейшей теме. Программку написал, какую хотели. Простенькую, без ляпов. Но и без ухищрений. Брутфорс. Хоть я не понимаю этого термина. Ну медленная, а кто сроки поставил? Я всегда вас читаю.

Dmitriy40, спасибо за красивые образцы. Воспользовался. Но мне показалось, что ускорения не получается. Конечно, рядом с numdiv ожидать трудно. Но я подумал, что это PARI оптимизирует, что без ошибок. Хотя ошибки же не исправляет... Во всяком случае даже чисто эстетически образцы хороши и полезны!

[Dmitriy40]

Yadryara в сообщении #1703725 писал(а):

интереснее с другими тремя паттернами разобраться

Смотрите, уравнение $2p^2+1=3q^2$ образует серию в натуральных $p,q$
$p_0=1,p_1=11,p_{i+1}=10p_i-p_{i-1}$
и имеет следующее решение в (псевдо)простых $p,q$:
p=109,q=89
Всё, больше решений нет вплоть до 4100 цифр в $p$ (это примерно 4120-е $p$, из которых лишь 18шт оказались псевдопростыми).

С уравнением $2p^2+3=5q^2$ чуть сложнее, оно распадается на две серии
$p_0=+1,p_1=49,p_{i+1}=38p_i-p_{i-1}$
$p_0=-1,p_1=11,p_{i+1}=38p_i-p_{i-1}$
В первой есть решение
p=3869693101,q=2447408809
и больше их нет вплоть до 4100 цифр в $p$ (это примерно 2600-е $p$, из которых лишь 8шт псевдопростое).
Во второй есть решение
p=11,q=7
и больше их нет вплоть до 4100 цифр в $p$ (это примерно 2600-е $p$, из которых лишь 6шт псевдопростое).

Итого, есть такие цепочки:
241: [2, 6, 6, 6, 6], len=4, p=11
23761: [2, 6, 6, 12, 12], len=2, p=109
29949049391853992401: [8, 6, 32, 36, 6], len=1, p=3869693101
Как видите цепочки из 5-ти нету. Вплоть до 4-х тысяч знаков в $p$. Потому ожидания найти таковую до 1e24 ... мягко говоря не оправданы.

(Программа)

Пример программы для третьей цепочки:
p0=-1; p=11; n=1; while(logint(p,10)<4100,[p,p0]=[p*38-p0,p]; n++; if(!ispseudoprime(p,1), next); print(n,". #p=",logint(p,10)+1); x=2*p^2+3; q=0; if(x%5==0 && issquare(x/5,&q) && ispseudoprime(q,1), print("p=",p,",q=",q)))
Каждая серия считалась около 7-12 минут.

PS. Все серии выделены ручным анализом из примеров решений каждого уравнения вольфрамальфой. Его формулы решений ужасны (хотя сразу дают $p_i,q_i$ по $i$), проще оказалось вручную подобрать рекурсивные выражения (так как помню что такие пелле-подобные уравнения часто имеют простую рекурсивную формулу).

-- 30.09.2025, 00:32 --

Отдельный интересный вопрос может ли первая серия пересекаться с любой из второй или третьей кроме как при $p=11$.
Вольфрамальфа говорит что уравнение $2p^2=3q^2-1=5r^2-3$ имеет лишь одно решение в целых $p=\pm1$, т.е. якобы таких решений вовсе нет, но веры ему нет так как $p=11,q=9,r=7$ он не нашёл (и вообще отвратно ищет такие двойные уравнения с тремя неизвестными).
Впрочем, добавив ему условие $p>1$ он решение находит (если самому последить за знаками, он глючит). Но ровно одно указанное. Верить ему в этом или нет - вопрос открытый.

[Yadryara]

Dmitriy40 в сообщении #1703806 писал(а):

Смотрите, уравнение $2p^2+1=3q^2$ образует серию в натуральных $p,q$
$p_0=1,p_1=11,p_{i+1}=10p_i-p_{i-1}$

Да, я как раз возился с этим уравнением, причём без Альфы и вышел на эту серию.

Dmitriy40 в сообщении #1703806 писал(а):

Вольфрамальфа говорит
что уравнение $2p^2=3q^2-1=5r^2-3$ имеет лишь одно решение в целых $p=\pm1$, т.е. якобы таких решений вовсе нет, но веры ему нет так как $p=11,q=9,r=7$ он не нашёл

Не нашла и правильно сделала. Вы, видимо, не учли что $9=3^2$, а Альфа учла. Она, видимо, подразумевает что в уравнении $2p^2=3q^2-1=5r^2-3$ кэфы попарно не равны переменным: $p \neq 2^k$, $q \neq 3^k$, $r \neq 5^k$, где $k$ — целое положительное.

То есть при таком $q$ имеется $3q^2 = 3^5$, а это другой случай и он в ММ77 рассмотрен.

Так что вопрос, похоже, закрыт: до 1e25 все 3493 пятёрки из 4-го паттерна.

gris
Я же не про кровать спрашивал. Понятны ли вам идеи ускорения, вот в чём вопрос. Откуда такой большой проигрыш берётся?

Ключевой момент: надо сначала перебирать и проверять кандидатов, а потом уже возводить в квадрат или сначала возвести, а потом уже перебирать и проверять?

[Dmitriy40]

Yadryara в сообщении #1703822 писал(а):

Так что вопрос, похоже, закрыт: до 1e25 все 3493 пятёрки из 4-го паттерна.

Я бы сказал до 1e4100.

Yadryara в сообщении #1703822 писал(а):

Ключевой момент: надо сначала перебирать и проверять кандидатов, а потом уже возводить в квадрат или сначала возвести, а потом уже перебирать и проверять?

Очень непонятный вопрос (почти как "перестали пить водку по утрам?"): в программе gris возведения в квадрат нету, однако программа рабочая, значит возводить в квадрат не обязательно и обе альтернативы вопроса могут быть неверны и должна быть как минимум третья.
Хотя смысл вопроса конечно понятен: что лучше, линейный или квадратичный перебор. И ответ на него очевиден.
Просто, и я снова это повторю, чтобы придумать квадратичный перебор надо разобраться как устроен паттерн (что там есть $p^2$), а для линейного ни с чем разбираться не надо, достаточно знать numdiv и в общем всё (шаг проверки равный k достаточно очевиден). И если не ставилась задача максимально ускорить поиск, то и разбираться никто не будет и получит работающую программу линейного перебора (так любимый заказчиком брутфорс).

-- 30.09.2025, 11:46 --

Yadryara в сообщении #1703822 писал(а):

Не нашла и правильно сделала. Вы, видимо, не учли что $9=3^2$, а Альфа учла. Она, видимо, подразумевает что в уравнении $2p^2=3q^2-1=5r^2-3$ кэфы попарно не равны переменным: $p \neq 2^k$, $q \neq 3^k$, $r \neq 5^k$, где $k$ — целое положительное.
То есть при таком $q$ имеется $3q^2 = 3^5$, а это другой случай и он в ММ77 рассмотрен.

Этого вообще не понял: вольфрамальфа вполне себе находит решение $p=\pm11,q=\pm9$ для уравнения $2p^2+1=3q^2$. Однако для двойного уравнения такие $q$ уже не находит.
А я этого учитывать и не собирался, меня устроят любые целые $q$, я их потом отдельно допроверю на простоту и всё.

[Yadryara]

Dmitriy40 в сообщении #1703872 писал(а):

в программе gris возведения в квадрат нету,

И напрасно. В первую очередь из-за этого она в итоге медленнее в триллион раз.

Dmitriy40 в сообщении #1703872 писал(а):

Хотя смысл вопроса конечно понятен: что лучше, линейный или квадратичный перебор.

Это он понятен тем, кто знаком с этими терминами. Я бы сказал по-другому: при равном шаге обычно намного быстрей перебирать не квадраты, а подквадратные числа.

Dmitriy40 в сообщении #1703872 писал(а):

И если не ставилась задача максимально ускорить поиск,

Автором темы именно такая задача ставилась.

И когда уже есть прога, то зачем же не до, а именно после неё постить ту, которая работает в миллиарды раз медленнее??

Видимо, gris опять не разобрался в моей проге. Уж не знаю, смотрел ли её вообще.

Dmitriy40 в сообщении #1703872 писал(а):

Этого вообще не понял: вольфрамальфа вполне себе находит решение $p=\pm11,q=\pm9$ для уравнения $2p^2+1=3q^2$. Однако для двойного уравнения такие $q$ уже не находит.

Значит я тоже не понимаю как здесь работает Альфа. Может это даже вопрос для отдельной темы. И вроде бы Вы уже эту тему поднимали.

[Dmitriy40]

Yadryara в сообщении #1703884 писал(а):

В первую очередь из-за этого она в итоге медленнее в триллион раз.

Не в триллион, а в зависимости от диапазона проверки. Может и в сто, и в секстиллион.
Нормально это называется сложностью и указывается как $O(n)$ или $O(\sqrt n)$. Отношение которых и зависит от величины $n$.

(Оффтоп)

Yadryara в сообщении #1703884 писал(а):

Я бы сказал по-другому: при равном шаге обычно намного быстрей перебирать не квадраты, а подквадратные числа.

Вот если бы сразу так и сказали - я бы промолчал, это выражение вполне корректно и понятно.
Хотя само наличие подквадратных чисел (т.е. неизвестных чисел в квадрате) совсем не очевидно (если не смотреть ни Вашу прогу, ни паттерн), тем более смотрите далее:

Yadryara в сообщении #1703884 писал(а):

Автором темы именно такая задача ставилась.

В данном месте я говорил не про Вас как автора данной темы, а про заказчика проги gris-а (которого мы все знаем). А там задача ставилась ускорить работу вполне конкретной проги (и эта задача очевидно успешно выполнена, сами посмотрите какое ускорение достигнуто, это было опубликовано), а не придумать другой быстрый алгоритм.
Вот зачем эту прогу показывать здесь я тоже не вполне понимаю. Но может просто не сравнивали скорость с Вашей прогой ...

(Оффтоп)

Yadryara в сообщении #1703884 писал(а):

Значит я тоже не понимаю как здесь работает Альфа. Может это даже вопрос для отдельной темы. И вроде бы Вы уже эту тему поднимали.

Ага, была такая тема, но толку особо не было.

[Yadryara]

Dmitriy40 в сообщении #1703891 писал(а):

Не в триллион, а в зависимости от диапазона проверки. Может и в сто, и в секстиллион.

Вы это уже говорили. И я не спорил. Чего же второй раз. Интервал я обозначил. Так что на нём в триллион.

Dmitriy40 в сообщении #1703891 писал(а):

В данном месте я говорил не про Вас как автора данной темы, а про заказчика проги

Думаете я не понял Только gris, увы, опять не сформулировал зачем было выкладывать столь медленную, когда уже показана быстрая.

Dmitriy40 в сообщении #1703891 писал(а):

Вот зачем эту прогу показывать здесь я тоже не вполне понимаю.

Конечно. Вот он, ключевой момент. Выделил крупно. Конечно я рад, что тема не заброшена, но зачем же такую прогу, gris? Вы правда не понимали, что и для D24 квадратов в паттернах полным-полно?

[Yadryara]

Yadryara в сообщении #1703953 писал(а):

D(48,20).

Интервал               Консоль               PARI
10  млн               44.8 sec           49.4 sec
100 млн         7min, 19.0 sec     7min, 26.1 sec

Попробовал ещё такой приём:

Вместо i%5 !=3 && i%7 !=0 сделал проверку 24 разрешённых остатков сразу по модулю 35:

bittest(0x5EF5B9E76, i%35)

Это позволило сбросить считанные секунды (до 7min, 11.7 sec), но зато это устойчивый выигрыш, то есть вроде можно с этим объединением остатков ещё поиграться.

Например, дальше есть setsearch(T[j],i%P[j]). И здесь тоже модуль именно простой, а не произведение простых...

[Dmitriy40]

Можно немного ускорить заменив в одной строке три конструкции
if(ispseudoprime(...), ... )
на конструкции
if(!ispseudoprime(...,1), next); ...
Это ускоряет с 60.5с до 52.5с для интервала 10млн, на 13%.
Здесь главное добавка в !ispseudoprime() вторым параметром 1 - с ! условие выполняется для составных чисел, а для них ispseudoprime никогда не ошибается, даже если проводит тест и не до конца (только одну итерацию проверки как указано). Ведь выгоднее быстрее отбросить неподходящие составные числа, чем проверять точно ли оно простое.
Без ! в ispseudoprime добавлять 1 нежелательно, это ухудшит фильтрацию (за счёт очень редких ошибок) и придётся эти числа обязательно допроверять.
Собственно допроверить придётся и эти три числа, которые теперь уже не обязательно простые. Но это можно сделать в самом конце, такие случае очень редки.

-- 01.10.2025, 02:07 --

Ещё 4с до 48.4с можно выиграть продублировав проверку (весь цикл) с factor(...,3000000) с меньшим порогом, например factor(...,30000). Всё же такие малые делители встречаются довольно часто и позволяют быстрее отбросить неподходящие кандидаты.
Итого выигрыш 60.5с -> 48.4с, уже 20%.

-- 01.10.2025, 02:20 --

Замена setsearch(T[j],i%P[j]) на bittest(bT[j],i%P[j]) ускоряет с 48.4с до 46.1с.
Массив вычисляется например так: bT=[sum(i=1,#j, 2^j[i]) | j<-T];
Объединять здесь модули смысла не вижу, слишком они уже большие.

-- 01.10.2025, 03:22 --

А переставив проверку с factor(...,30000) на условие выше, сразу за двойным ispseudoprime, можно ускорить до 44.3с.

[Yadryara]

Это всё равно круто, что даже на PARI всё ещё удаётся ускорять. Суммарно на 35%.

Теперь с целью ускорения, я хочу небольшой кусочек программы написать на Асме, скомпилировать и интегрировать в прогу. Какой-то самый простой, пусть даже 1-2 операции, чтоб я понял что делается.

[Dmitriy40]

Во-первых на С писать гораздо проще. Пусть результат и не столь же быстр.

Во-вторых, одним маленьким кусочком не обойтись если хотите получить ускорение вместо замедления. Всегда же есть накладные расходы на запуск внешней проги, в PARI это сотые доли секунды, значит внешняя прога должна полезно работать хотя бы десятые секунды. Тест накладных расходов в PARI:

{t0=getwalltime(); n=0; shag=10^2;
while((tt=getwalltime()-t0)<10000,\\Считаем сколько раз за примерно 10с сможем запустить внешнюю прогу
   q=0; for(i=1,shag, q+=#externstr("cmd.exe /c cls")); n+=shag;
);
print("Speed: ",round(n/tt*1e3),"/c (",n," iterations)");

У меня выводит около 180 раз в секунду, значит накладные расходы примерно 5мс.

Перенести длинный if и цикл с setsearch/bittest в прогу будет недостаточно, они явно выполняются быстрее 10мс (у меня 17с на 10млн вместе с циклом по i, т.е. 1.7мкс). Разве что для тренировки. Но ускорения так не получить, скорее замедление в тысячи раз.
Переносить три isprime в прогу - трудно, числа большие, писать самому нормальный тест простоты откровенно влом, да и не факт что получится сильно быстрее isprime в PARI (помнится когда с Hugo занимались минимизацией цепочек, у меня простой тест простоты работал лишь немногим быстрее PARI и только для чисел меньше $2^{64}$).
В итоге в прогу стоит переносить сразу цикл по i с проверкой длинным if и bittest (это одна команда, а setsearch целая процедура) и возвратом списка кандидатов в том же 10млн отрезке i. Их будет 1.5млн из 10млн, это слишком много, фильтрация 7:1 вообще ни о чём. А выполняться внешняя прога будет точно дольше 1.5млн*54модуля*150тактов/3.5ГГц=3.5с (только три деления числа 192 битов на мелкое простое плюс 0.1с всё остальное), на самом деле время будет даже больше, раза в два. Экономия всего 13с (17с-4с) из 44с.
Выходит что isprime придётся таки делать в проге в каком-то виде, иначе про существенное ускорение говорить не приходится. Проверка делимости на простые $59...2^{17}$ оставляет 35% исходных чисел, для трёх проверяемых чисел это даёт $0.35^3=4.3\%$. Вот 23:1 уже неплохой коэффициент фильтрации, но 3*12250*3 делений (для трёх проверяемых чисел длиной по 192 бита) займут 1.5мс на каждый кандидат из 1/7 исходных. Пусть даже в 15 раз меньше, типа много составных не будут доходить до конца делений, всё равно это 0.1мс на кандидата и 150с на 1.5млн кандидатов после первой фильтрации 7:1. Замедление в несколько раз.

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

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

Yadryara
Так что Вы хотите перенести в асм (или С)?
Компилятор С кстати есть уникально маленький, для простых программ достаточно. Но как к нему нормально подключить асм+AVX я пока не разобрался.

[Dmitriy40]

Dmitriy40 в сообщении #1704123 писал(а):

Переносить три isprime в прогу - трудно, числа большие, писать самому нормальный тест простоты откровенно влом,

Ещё раз подумав, прикинув затраты, выходит не так уж и трудно. Муторно, но принципиально понятно. Тест Ферма гарантированно обнаруживает все составные, что нам и нужно, а его вычисление идеологически несложно. Для чисел длиной 150-190 битов конечно ещё та морока, но принципиальных трудностей нет. Даже прикинул на глазок скорость, получается порядка 100мкс на число (10000 чисел в секунду). PARI кстати раз в 10 быстрее, 12мкс на число. Ну может у меня "глазок" кривой ... я ж код не писал и не мерил. Или в PARI тест Ферма хорошо оптимизирован, лучше чем я сходу втупую прикинул.

[Yadryara]

Да, я сразу же дал понятную ссылку на тему про Асм. И там решето Эратосфена есть, но только вроде для 32-х разрядов.

Кстати, сумел ещё ускорить прогу — c 317 до 309-310 секунд на 100млн.

Похоже, что цикл forstep в PARI намного медленнее чем for. Может нужно как раз по-другому организовать цикл? Но я через while не люблю — каждый раз забываю как он работает.

[Dmitriy40]

Yadryara в сообщении #1704177 писал(а):

Похоже, что цикл forstep в PARI намного медленнее чем for. Может нужно как раз по-другому организовать цикл? Но я через while не люблю — каждый раз забываю как он работает.

Да не сказал бы что намного медленнее, всего 10%:

? for(i=10^9-10^8,10^9, i%5!=3&&i%7!=0&&i%11!=9)
time = 18,799 ms.
? forstep(i=10^9-3*10^8,10^9,3, i%5!=3&&i%7!=0&&i%11!=9)
time = 20,592 ms.
? i=10^9-10^8-1; while((i++)<=10^9, i%5!=3&&i%7!=0&&i%11!=9)
time = 42,652 ms.
? i=10^9-10^8-1; while(i<10^9, i++; i%5!=3&&i%7!=0&&i%11!=9)
time = 43,056 ms.
? i=10^9-3*10^8-3; while(i<10^9, i+=3; i%5!=3&&i%7!=0&&i%11!=9)
time = 45,037 ms.

Если у Вас разница существенно больше, то вероятно влияют какие-то ещё команды, которые есть с forstep и нет с for.
А вот с while аж вдвое тормознее.