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

[Dmitriy40]

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

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

Ну так уберите цепочки с valids<18, например. С экрана, в лог пусть пишутся.

[Yadryara]

Вот как выглядят сейчас мои действия и отклик на них:

yadryara@DESKTOP-QPP2F5P:~/D48-22$ gp2c-run -g Rab_60_67_1.gp
Rab_60_67_1.gp.c: In function ‘init_Rab_60_67_1’:
Rab_60_67_1.gp.c:1132:3: warning: implicit declaration of function ‘gp_quit’ [-Wimplicit-function-declaration]
 1132 |   gp_quit(0);
      |   ^~~~~~~
Reading GPRC: /etc/gprc
GPRC Done.

                                          GP/PARI CALCULATOR Version 2.15.4 (released)
                                  amd64 running linux (x86-64/GMP-6.3.0 kernel) 64-bit version
                              compiled: Apr  1 2024, gcc version 13.2.0 (Ubuntu 13.2.0-23ubuntu3)
                                                   threading engine: pthread
                                         (readline v8.2 enabled, extended help enabled)

                                             Copyright (C) 2000-2022 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?18 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 500000, nbthreads = 12
? allocatemem(2^26)
  ***   Warning: new stack size = 67108864 (64.000 Mbytes).
? init_Rab_60_67_1()



bolv    = [13, 2, 3, 28, 5, 18, 1, 32, 3, 50, 7, 12, 11, 338, 45, 8, 1, 102, 1, 20, 3, 2]     22

onlyp   = [18]                                                                                 1

pq      = [16]                                                                                 1

pqr     = [1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 20, 21, 22]                        17

pqrs    = [7, 17, 19]                                                                          3

mkv     = [1, 2, 3, 5, 7, 9, 11, 13, 16, 17, 18, 19, 21, 22]                                  14

rkp     = [7, 11, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61]                             14


prter = [5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61]



Potok 4


60.   [18, 2, 1, 5]

1938056648387506790981431134219620047372044882450841             11  11111 1  1111111        15
1952310394628263019166481192813446835466617941263641              111111 111 1 111111        16
1028311535443738759614215960146102814009831834447641                11 1  1 11 11111         11
1955937841294378999893294305625093262693190775878041              111  1 1 11 1 1111  1      13

То есть сначала команда на компиляцию, затем задаю память, затем инициализация.

Программа работает вроде нормально.

bolv — болванка. Это частично заполненный паттерн. Далее он очищается и вновь заполняется.

onlyp — место, где надо найти только простое число. По нему и идёт перебор.

pq — место, где надо найти число с 4 делителями. Как правило, это произведение двух простых чисел p и q.

pqr — 17 мест, где надо найти число с 8 делителями. Как правило, это произведение трёх простых чисел p, q и r.

pqrs — 3 места, где надо найти число с 16 делителями. Как правило, это произведение четырёх простых чисел p, q, r и s.

mkv — 14 мест для расстановки квадратов.

rkp — 14 расставляемых (квадратов) простых

prter — простые для терапевтики.

Номер и множитель потока.

Номер комплекта и расположение 4-х характерных квадратов простых.

Ну и дальше найденные цепочки с valids по всему полю шириной 31.

[Yadryara]

Ну вот какая-то ошибка возникла. Причём во всех 9-ти открытых окнах :

1715572919914207030518870899895376379783886899716441              11  1111 1111 1111111      17
Error: 0xd00002fe
Код ошибки: Wsl/Service/0xd00002fe
Press any key to continue...

И при нажатии на пробел окно закрывается.

И при нажатии Win на самом верху появилась надпись

Недавно добавленные

WSL Settings
Система

[wrest]

wrest в сообщении #1711691 писал(а):

Ещё появились кое-какие мысли, но надо потестить.

Dmitriy40
А мысли такие. Поскольку тестируемые кандидаты (начальные члены цепочек) у нас из арфметической прогрессии $n_k=n_0+2km$ с шагом 2m то длинное деление на малые простые нам по сути нужно только для определение остатков от деления n0 и 2m на простые. Так же как мы вычисляем остатки членов цепочки по остаткам первого члена, мы можем вычислять и остатки первых членов один раз вычислив остатки от деления n0 и 2m на малые простые, и перейти от длинной арифметики в 64-битную.
Получается же, что любое проверяемое нами число вычисляется по формуле $n_{kj}=n_0+k\cdot 2m + j, k=1...10^6, j=0..20$ где k - номер кортежа, j - номер числа в кортеже.
Тогда, длинное деление с остатком для решения сравнения $n_{k0} \equiv x \mod p$ мы можем заменить сложением/умножением.

-- 05.12.2025, 10:35 --

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

А v[] и nu[] не используются что ли?

Используются (как глобальные переменные), то есть в логике всё, кроме предфильтров типа
if((n+#v-1)%3^6==21-19, next);\\Такие n не подходят
и двух проверок (ранний выход по явному переполнению и проверки на бесквадратность) полностью соответствует. Но есть и особенности реализации - я вынес проверки в отдельную функцию, там инициализируются аккумуляторы
nf=vector(#v,d,1); nn=vector(#v,d,1);\\Вектор количества делителей (omega) и вектор накопления их для каждой позиции
Так что у меня ещё есть накладные расходы на вызов функции проверки, и передачу в неё длинного числа n
Похоже, от функции тоже надо отказаться...

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

Судя по данным ставить порог выше 2^18 нет смысла, несколько штук проще допроверить прямым numdiv/factor.

Ну не знаю, время скакнуло при переходе порога с 2^19 на 2^20, так что наверное порог нужен 2^19
Это надо бы перетестировать, пожалуй.
Ну и ещё - надо бы посмотреть сколько отъедает последняя проверка случаев недостатка делителей, когда число разложилось полностью и 48 делителей не набралось.

[Yadryara]

wrest в сообщении #1711712 писал(а):

Поскольку тестируемые кандидаты (начальные члены цепочек) у нас из арфметической прогрессии $n_k=n_0+2km$ с шагом 2m

О, как! У меня шаг уже давно 6-кратный. То есть 6m или 6*lcm(v). Я сделал автозатачивание под 6-кратный шаг.

А у Дмитрия шаг двукратный, потому что m = lcm([63, 3610, 3179, 12, 17797, 3362, 75, 392, 841, 18, 2209, 20, 1587, 242, 6727, 96, 14045, 338, 243, 68, 35131]) = 229403502054304075118105309394590566946400

То есть перебор здесь втрое менее эффективен чем надо бы. Хотя... 3^6 ведь проверяется, так что шаг видимо всё-таки 6-кратный, просто достигается это по-другому.

[Dmitriy40]

wrest в сообщении #1711712 писал(а):

и перейти от длинной арифметики в 64-битную.

Даже в 32-х битную, все простые у нас меньше 2^32.

wrest в сообщении #1711712 писал(а):

Тогда, длинное деление с остатком для решения сравнения $n_{k0} \equiv x \mod p$ мы можем заменить сложением/умножением.

Умножение тоже не нужно: переход от предыдущего n к следующему производится просто прибавлением 2m. Соответственно переход от n%p к следующему n%p производится командой (n%p+m%p)%p (тут и далее двойка уже учтена в m) - оба слагаемых уже известны, достаточно сложить и проверить не превысило ли p, и если да, то просто вычесть p, ведь больше 2p стать не может: np=np+mp;if(np>=p,np-=p);. Разумеется np и mp это массивы. И я об этом уже писал где-то выше. Именно так работают мои программе на асме.
В PARI есть специальный тип Mod(x,p), их можно складывать, даже в виде матриц n[1900,22]+s[1900,22] (матрицы указанного размера, в каждом элементе обеих матриц по Mod(x,p), остаток по одному из 1900 простых на каждом из 22 мест цепочки и увеличение этого остатка при добавлении m), но тяжело проверить матрицу на появление Mod(0,p), т.е. что где-то разделилось. Насколько это быстро по сравнению с ручным вычислением сложения по модулю, тем более после компиляции, я не знаю.

wrest в сообщении #1711712 писал(а):

Так что у меня ещё есть накладные расходы на вызов функции проверки, и передачу в неё длинного числа n

Длинное число передаётся реально по ссылке, т.е. это просто регистр. Вы думаете так уж долго передать ссылку? По сравнению с тысячами делений и проверкой простоты?
Сколько занимает вызов функции после компиляции я опять же не знаю, может и сравнимо с парой проверок (деление всё же весьма медленно).

wrest в сообщении #1711712 писал(а):

Ну и ещё - надо бы посмотреть сколько отъедает последняя проверка случаев недостатка делителей, когда число разложилось полностью и 48 делителей не набралось.

Дык она же выполняется всего 22 раза, не для каждого простого, только если кандидат не отбросился, и занимает около 2% времени. А на самом деле ещё реже, ведь много чисел будет отброшено выкинутыми Вами предпроверками по простым 3...19 или 3...53.
Ну и скорость таки надо сравнивать на диапазон i, а не kpod, последнее полезнее для оценки скорости именно факторизации и коэффициента фильтрации, но в реальной программе важна скорость перебора i. Чтобы набрать 1млн kpod с предпроверками по простым 3...53 надо перебрать 7348061/2 i.
И я перебираю i так: forstep(i=1,oo,2, n=n0+m*i; ...), это лишь для теста Вам сделал for вместо forstep и 2m вместо m, и зря.

-- 05.12.2025, 13:27 --

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

О, как! У меня шаг уже давно 6-кратный. То есть 6m или 6*lcm(v). Я сделал автозатачивание под 6-кратный шаг.

А у Дмитрия шаг двукратный, потому что

... паттерны разные, у меня без места p, у Вас с ним, вот поэтому Вы можете идти с шагом 6, а мне приходится идти с шагом 2 и выкидывать лишь 1 из 3 остатков по модулю 3 (проверкой на 3^6), а не 2, как надо для 6 кратного шага.
Так что кандидатов перебирается вдвое больше, не втрое. Третий быстро откидывается проверкой по 3^6.
И не забываем, это не рабочая программа, это тестовая для камланий с факторизацией.

-- 05.12.2025, 13:39 --

wrest в сообщении #1711712 писал(а):

Ну не знаю, время скакнуло при переходе порога с 2^19 на 2^20, так что наверное порог нужен 2^19

Вы неправильно ставите условие. Надо проверить что дольше, допроверить 19-11 цепочек или увеличить порог с 2^16 до 2^17. И аналогично дальше. Может оказаться что допроверка большего количества цепочек быстрее замедления времени от увеличения шага и тогда шаг лучше не увеличивать.
Допроверка в самом примитивном виде выглядит так: nd=[numdiv(n+d-1)|d<-[1..#v]] (просто numdiv на все числа цепочки). Или так: for(d=1,#v, if(numdiv(n+d-1)<>48, next(2))) (но так не будет никаких приближений, только само решение). Numdiv будет тормозить только на недоразложенных местах, разложенные пролетят мгновенно из-за factor_add_primes=1.

[Yadryara]

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

паттерны разные, у меня без места p, у Вас с ним, вот поэтому Вы можете идти с шагом 6, а мне приходится идти с шагом 2 и выкидывать лишь 1 из 3 остатков по модулю 3

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

Вот как у меня сделано:

for(j = 1, #v, if(v[j]%32==0, m32 = j));
for(j = 1, 9, if(v[j]%9 ==0, m9 = j));

aa = v[ip+1]; m = 6*lcm(v)/aa;

n0 = (lift(chinese([chinese([Mod(-j+1, v[j]) | j<-[1..#v]]),
Mod(32-m32+1, 64), Mod(9-m9+1, 27)])) + ip) / aa;

Таким образом, если нет перебираемого места, то aa = 1; ip = 0;

То есть формулы становятся проще, но шаг-то всё равно 6-кратный.

[Dmitriy40]

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

Вот как у меня сделано:

Только в показанном мною паттерне не 3^2, а 3^5.

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

То есть формула становится проще, но шаг-то всё равно 6-кратный.

Про формулу не скажу, она похожа на правильную для 6-ти кратного шага, но точнее можно определить шаг по количеству допустимых остатков:

{v=[63, 3610, 3179, 12, 17797, 3362, 75, 392, 841, 18, 2209, 20, 1587, 242, 6727, 96, 14045, 338, 243, 68, 35131];
print("v=",v);
print("m= ",m=lcm(v));
n0=lift(chinese([Mod(-j+1, v[j]) | j<-[1..#v]]));
for(d=0,5,
   n=n0+m*d; print1("d=",d,": ",y=[((n+j-1)/v[j])%6 | j<-[1..#v]]);
   if(#select(t->(t%2==0||t==3),y)==0, print(" OK, n0=",n) , print);
);}

Запуск:
v=[63, 3610, 3179, 12, 17797, 3362, 75, 392, 841, 18, 2209, 20, 1587, 242, 6727, 96, 14045, 338, 243, 68, 35131]
m= 229403502054304075118105309394590566946400
d=0: [1, 1, 1, 1, 1, 1, 5, 5, 5, 1, 1, 1, 1, 5, 5, 4, 5, 1, 1, 5, 5]
d=1: [1, 1, 1, 1, 1, 1, 5, 5, 5, 1, 1, 1, 1, 5, 5, 1, 5, 1, 3, 5, 5]
d=2: [1, 1, 1, 1, 1, 1, 5, 5, 5, 1, 1, 1, 1, 5, 5, 4, 5, 1, 5, 5, 5]
d=3: [1, 1, 1, 1, 1, 1, 5, 5, 5, 1, 1, 1, 1, 5, 5, 1, 5, 1, 1, 5, 5] OK, n0=742745703976917864157851884948204505669969
d=4: [1, 1, 1, 1, 1, 1, 5, 5, 5, 1, 1, 1, 1, 5, 5, 4, 5, 1, 3, 5, 5]
d=5: [1, 1, 1, 1, 1, 1, 5, 5, 5, 1, 1, 1, 1, 5, 5, 1, 5, 1, 5, 5, 5] OK, n0=1201552708085526014394062503737385639562769
Надеюсь два допустимых остатка по модулю 6 хорошо видны?

Исправление 5 на 2 в for и %6 на %2 в следующей строке даёт такие допустимые остатки по модулю 2:
d=0: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1]
d=1: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] OK, n0=283938699868309713921641266159023371777169
Именно такое n0 я и указал в тестовой программе.

[wrest]

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

В PARI есть специальный тип Mod(x,p), их можно складывать, даже в виде матриц n[1900,22]+s[1900,22] (матрицы указанного размера, в каждом элементе обеих матриц по Mod(x,p), остаток по одному из 1900 простых на каждом из 22 мест цепочки и увеличение этого остатка при добавлении m),

Они все наверняка длинные, как принято в pari.

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

но тяжело проверить матрицу на появление Mod(0,p), т.е. что где-то разделилось. Насколько это быстро по сравнению с ручным вычислением сложения по модулю, тем более после компиляции, я не знаю.

Не, ну так мы будем вычислять слишком много. В текущем варианте же мы останавливаемся сразу как только поняли что цепочку надо отсеять. Этот подход надо сохранить. Наверное об этом вы и пишете говоря "тяжело проверить матрицу" -- нам не нужна матрица. Нам (как мне кажется) главное не повторять длинных делений на одно и то же малое простое, если мы их уже делали, вот то что жирным выделено.
forprime(p=59,2^20,\\Проверяем делимость только на такие простые
d=#v-x%p;if(d<1, next); nf[d]++; nn[d]*=p;\\Получение остатка, проверка что он (и простое p) попадает в кортеж, накопление числа и самих делителей

-- 05.12.2025, 15:45 --

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

Вы неправильно ставите условие. Надо проверить что дольше, допроверить 19-11 цепочек или увеличить порог с 2^16 до 2^17.

Я согласен. Ну то что показали мои измерения (надо повторить) -- скачок по времени, 10 секундный, произошёл при переходе порога от 2^19 к 2^20
lim=512 passed=11007 time: 27,136 ms
lim=1024 passed=4068 time: 26,546 ms
lim=2048 passed=1510 time: 24,978 ms
lim=4096 passed=542 time: 25,053 ms
lim=8192 passed=228 time: 25,473 ms
lim=16384 passed=102 time: 27,477 ms
lim=32768 passed=40 time: 27,170 ms
lim=65536 passed=19 time: 28,314 ms
lim=131072 passed=11 time: 29,373 ms
lim=262144 passed=5 time: 30,858 ms
lim=524288 passed=3 time: 32,138 ms
lim=1048576 passed=2 time: 43,510 ms
ну за 10 секунд то уж всяко можно проверить одну цепочку.
Остальные изменения по времени между разными порогами на уровне погрешности измерений, я им не верю, возможно процессор нагрелся или венда пошла проверять не пора ли обновиться. Типа ускорения на 10% (3 секунды) при переходе с порога 2^9 к порогу 2^11 Да я и этом скачку с 32 до 42 секунд не верю. Надо больше статистики. Там же проверялся первый миллион цепочек конкретной структуры (паттерна), надо проверить например пятый и восьмой миллионы...

[Dmitriy40]

wrest в сообщении #1711742 писал(а):

Нам (как мне кажется) главное не повторять длинных делений на одно и то же малое простое, если мы их уже делали, вот то что жирным выделено.

Где Вы видите повторы? В цикле forprime меняется p, потому деление будет каждый раз на разные числа.
Между циклами меняется x (которое вычисляется по n), так что тоже все деления будут другими.

Другое дело что да, можно заменить деление на сложение по модулю, но тогда вместо одного n или x придётся хранить и вычислять много разных n[] или x[], для каждого p своё. Про матрицы был не прав, достаточно вектора длиной с количество простых.
И да, сложение по модулю двух не слишком больших векторов вполне может быть быстрее деления длинного числа на короткое и тогда выгоднее посчитать больше сложений чем одно деление. Даже на PARI. Пробовать надо.
Например вот сложение двух векторов n[] и s[] по модулю третьего p[]: n+=s; t=n-p; for(j=1,#n, if(t[j]>0, n[j]=t[j]);, n,s,t,p - вектора длиной с количество проверяемых простых. И все они могут быть int32 (или даже int16 для порога до 2^14, но это вероятно медленнее).

wrest в сообщении #1711742 писал(а):

надо проверить например пятый и восьмой миллионы...

Тут Вы абсолютно правы.
Но вся разница будет связана исключительно с ispseudoprime, простых будет меньше, будет чаще срабатывать условие nf=nu и реже условие nf<nu. Что пересилит кстати так за глаза непонятно, фильтрация может и улучшиться, и ухудшиться. А вот частота появления малых делителей (те что в forprime) не изменится (пока берём интервал по kpod больше чем порог для p, который 2^20), как соответственно и частота вызова ispseudoprime для nf=nu.

-- 05.12.2025, 17:04 --

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

Например вот сложение двух векторов n[] и s[] по модулю третьего p[]: n+=s; t=n-p; for(j=1,#n, if(t[j]>0, n[j]=t[j]);, n,s,t,p - вектора длиной с количество проверяемых простых.

И эта операция прекрасно ложится в SSE/AVX2. Например в AVX2 она реализуется тремя командами над регистром длиной 256 битов или 8 штук int32 и все три команды выполняются в среднем за полтора такта, плюс две команды накопления флага равенства нулю любого из остатков за такт, т.е. скорость такого суммирования длинных векторов составляет примерно 3.3 элементов int32 за такт. Или 6.6 элементов int16 за такт. Деление 192 битового числа (до 57 знаков) занимает 100-300 тактов в самом лучшем случае. Разница на три порядка.
Сможет ли gcc выдать столь же быстрый код - очень сомневаюсь, но это не мои проблемы, у меня есть ассемблер.

[wrest]

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

Длинное число передаётся реально по ссылке, т.е. это просто регистр.

Не, в pari передаётся значение. Во встроенных функциях иногда можно по ссылке передавать. User Guide к 2.17.2, стр. 42

Цитата:

Important technical note. Built-in functions are a special case since they are read-only (you
cannot overwrite their default meaning), and they use features not available to user functions,
in particular pointer arguments. In the present version 2.17.2, it is possible to assign a built-in
function to a variable, or to use a built-in function name to create an anonymous function, but
some special argument combinations may not be available:
? issquare(9, &e)
%1 = 1
? e
%2 = 3
? g = issquare;
? g(9)
%4 = 1
? g(9, &e) \\ pointers are not implemented for user functions
*** unexpected &: g(9,&e)
*** ^---

-- 05.12.2025, 17:30 --

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

Другое дело что да, можно заменить деление на сложение по модулю, но тогда вместо одного n или x придётся хранить и вычислять много разных n[] или x[], для каждого p своё.

Да, вот за оптимизацию тут и надо подумать. Вопрос-то конечно не в вычислять (это и так делаем), а в хранить и иметь быстрый доступ (мемоизация и вот это всё).

[Dmitriy40]

wrest в сообщении #1711753 писал(а):

Не, в pari передаётся значение.

Ну-ну. Как Вы собрались передать по значению многогигабайтный вектор? А если таких векторов несколько передаются одной функции?
Кроме того, насколько мне известно, С в принципе не позволяет передавать аргументы переменного размера, единственное исключение - последний аргумент, точнее один список последних аргументов в функциях типа printf. Всё что не помещается в регистр - передаётся по ссылке в память. А PARI написан на C, так что наследует его ограничения.
То что в синтаксисе PARI передача аргумента выглядит как по значению - ни о чём, физически передаётся ссылка в регистре (ну или стеке).
UPD. А, понял, Вы копирование вектора и передачу ссылки на копию назвали передачей по значению, т.е. типа не ссылка на уже существующий объект, а ссылка на новый, старый же не портится. Да, это не быстро, но уж точно быстрее даже одного деления, скопировать 3-5 слов вообще не проблема, дольше память из кучи выделить. Но от передачи n или n0 с m никуда не деться же.

wrest в сообщении #1711753 писал(а):

в хранить и иметь быстрый доступ (мемоизация и вот это всё).

Я повторю вопрос: где в x%p Вы нашли повторные вычисления?!

[wrest]

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

Ну-ну. Как Вы собрались передать по значению многогигабайтный вектор?

Так я и не собираюсь, поэтому. Но вот кстати как раз векторы вроде теперь можно ссылкой передавать (то ли с 2.15 то ли позже), но позабыл я как именно. Кажется, там тильду надо писать перед именем переменной. Т.е. во встроенных функциях передача по ссылке обозначается как &varname а в пользовательских -- как ~varname но там тоже какие-то заморочки, типа так передавать можно только векторы или какие-то такие (контейнерные) переменные.

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

Кроме того, насколько мне известно, С в принципе не позволяет передавать аргументы переменного размера

Так мы же про pari а не Си

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

Я повторю вопрос: где в x%p Вы нашли повторные вычисления?!

Сейчас мы для каждого начала цепочки их вычисляем (конечно не по всем простым а по небольшой их части). Миллион раз для миллиона цепочек. А могли бы и не.
Но с модульной арифметикой (т.е. с вычислением остатков от деления на простые только для n0 и m) я обнаружил ещё одну засаду. Нам надо проверить нефакторизованный остаток (частное) на простоту. А для этого надо вычислить текущий член цепочки и поделить его на накопленный делитель, вот это:
ispseudoprime((n+d-1)/v[d]/nn[d])
то есть сперва длинное умножение, затем длинное деление, и затем проверка на простоту.
Я попробовал с Mod но что-то у меня не получается быстро. Сами Mod-ы предвычисляются по всем простым от 59 до 2^20 для n0 и m довольно быстро (30мс или около того). Но потом быстрее чем сейчас (с forprime по каждой цепочке) пока не выходит. Соотношение такое: 36 сек на forprime vs 57 сек с подготовкой простых и их Mod-ов. Видимо от Mod-ов надо перейти к обычной арифметике. Ну или эти длинные деления не особо то и мешают...

[Dmitriy40]

wrest в сообщении #1711786 писал(а):

Нам надо проверить нефакторизованный остаток (частное) на простоту.

Не всегда: может набраться nf>nu делителей и тогда на простоту проверять не придётся (потому я и добавил проверку nf>nu, по инерции). Можете закомментировать проверку nf==nu и посмотреть насколько плохо фильтруется без неё, только по nf>nu (ну и квадратам). У меня фильтрация упала с 0.1млн:16 (без условия nf<nu, но с проверкой квадратов) до 0.1млн:2932, при этом время выросло с 4.7с до 134с.
Так что обычно да, проверять простоту надо. А это весьма не быстро.

wrest в сообщении #1711786 писал(а):

вот это:
ispseudoprime((n+d-1)/v[d]/nn[d])
то есть сперва длинное умножение, затем длинное деление, и затем проверка на простоту.

Длинное умножение и деление - сущие мелочи по сравнению с возведением числа в степень длинного числа по модулю длинного числа (т.е. Mod(a,x)^(x-1)) в ispseudoprime ... Если что, такое возведение в степень требует умножений двух длинных чисел по модулю третьего длинного числа, и таких умножений примерно в полтора раза больше длины чисел в битах (полтора умножения на бит длины). Это если сделано правильно и без деления, а то ещё и делений двух длинных чисел столько же. Так что при трёх сотнях умножений и возможно ещё и делений пара лишних умножений с делениями погоды вообще не сделают.

[wrest]

Ну в общем вот код с модульной арифметикой вместо forprime.
Сделан в виде функции, для удобства компиляции и многопоточного запуска.
К сожалению, расставить типы не удалось, постоянно segmentation fault.

Используется синтаксис C

em3(myname,k_start,k_len,lim)={
/* myname - будет печататься в отчёте
* k_start - начало диапазона
* k_len - количество цепочек для просчёта
* lim - порог для пробных простых, оптимальный будет 2^18
*/

my(t0=getwalltime());\\Замер времени подготовки

\\ переменные специфичные для паттерна
my(v=[63, 3610, 3179, 12, 17797, 3362, 75, 392, 841, 18, 2209, 20, 1587, 242, 6727, 96, 14045, 338, 243, 68, 35131]);
my(n0=283938699868309713921641266159023371777169);
my(m=229403502054304075118105309394590566946400);
my(nu=[3, 2, 3, 3, 3, 3, 3, 2, 4, 3, 4, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3]);\\Сколько разных простых делителей должно быть в каждой позиции
my(len=#v);

\\ вычисляем остатки по всем малым простым диапазона
my(pr=primes([59,lim])); \\ вектор с простыми
my(pr_q=#pr); \\ и сколько их там оказалось
my(pmods_n0=vector(#pr,i,n0%pr[i])); \\вектор с остатками по простым по n0
my(pmods_m=vector(#pr,i,m%pr[i])); \\ вектор с остатками по простым по m

\\ вспомогательные переменные
my(q=0,q1=0,q2=0,c=0);\\Счётчики
my(d); \\ индекс числа от начала цепочки
my(nf=vector(len,d,1)); \\Вектор количества делителей (omega) 
my(nn=vector(len,d,1)); \\ и вектор накопления их для каждой позиции
my(n);

\\ print(myname, " приступил. Подготовился за ",strtime(getwalltime()-t0));

\\ основной цикл
t0=getwalltime();\\Замер времени отсеивания
for(k=k_start,k_start+k_len, \\ идём по цепочкам
\\Здесь начинается проверка конкретного n
   nf=vector(len,d,1); nn=vector(len,d,1); \\ инициализируем аккумуляторы для цепочки
   for(i=1,pr_q,\\ идём по простым для каждой цепочки
      \\ вычисляем положение числа которое делится на это простое относительно начала цепочки
      d=len-(pmods_n0[i]+2*k*pmods_m[i]+len-1)%pr[i]; 
      \\ на i-е простое ни одно число цепочки не разделилось, выходим на проверку следующего
      if(d<1, next);
      nf[d]++; \\ на i-е простое разделился d-й член цепочки накопим по нему омегу 
      nn[d]*=pr[i]; \\ и накопим по нему частное
      n=n0+2*k*m+d-1; \\ вычислим текущее число 
      \\Если уже перебор множителей, то отбрасываем всю цепочку
      if(nf[d]==nu[d],q1++;if(!ispseudoprime(n/(v[d]*nn[d])), next(2)));
   );
   \\Проверка что разложение закончено, а делителей меньше требуемого
   for(i=1,len, 
        if(nf[i]<nu[i] && ispseudoprime((n0+2*k*m+i-1)/(v[i]*nn[i])),q1++;next(2))); 
   \\Проверка закончена, что осталось то осталось, подсчитаем сколько их для отчёта о проделанной работе

  q++; \\ посчитаем сколько прошло через фильтр
  \\ прошедшие полностью факторизуем, не будем откладывать.
  n=n0+2*k*m; c=0; for(i=1,len,if(numdiv(n+i-1)==48,c++));if(c==len,q2++); \\ сюда добавить логирование и фанфары если найшлось
 \\ print1(n0+2*k*m," ");for(i=0,len-1,if(numdiv(n+i)==48,print1(1);c++,print1(0)));if(c==len,q2++);print1(" ",c);print();
 
);
\\ отчитаемся о работе
print(myname," готов,",q2," найдено,",q," прошло,",q1+q," проверено псевдопростых. скорость счета ", 1000*k_len\(getwalltime()-t0)," цеп/сек, время счёта:",strtime(getwalltime()-t0));
return(q2);
}

Показатели такие, на миллион цепочек. Проверяется всё полностью (все прошедшие дофакторизуются, от 2 до 8 на миллион).
Интерпретация, один поток ~23 тыс/сек
Интерпретация, 4 ядра/8 потоков ~47 тыс/сек
Компиляция, один поток ~47 тыс/сек
Компиляция, 4 ядра/8 потоков ~95 тыс/сек

Так что хоть и не на порядок, но в 4 раза поднять скорость компиляцией + многопоточностью на "кортежных" вычислениях вполне можно. Это с учётом того, что там много проверки на простоту -- чуть больше единицы на цепочку (1.1) то есть на миллион цепочек - 1.1 миллиона проверок на простоту нефакторизованного остатка (частного). Ну и с учётом того что с хорошей типизацией я не справился.

К сожалению, как запустить многопоток "эталона" пока не очень представляю -- надо конвертировать в функцию (но тогда он не будет эталоном) или запускать в разных сессиях линукса.