Симметричные кортежи из последовательных простых чисел

[Yadryara]

Dmitriy40, единственное что я пока понял, это то, что новости хорошие. А можно на простейших числовых примерах, для дятлов ?

[Dmitriy40]

Yadryara
Можно, глядишь и сам лучше разберусь.
Для примера возьмём паттерн 0,2,6,8 (из темы про асм). Он имеет следующие допустимые остатки:

? v=[0,2,6,8]; forprime(p=2,15,m=setminus(vector(p,i,i-1),Set(-v%p)); print(p,": ",m,", len=",#m);)
2: [1], len=1
3: [2], len=1
5: [1], len=1
7: [2, 3, 4], len=3
11: [1, 2, 4, 6, 7, 8, 10], len=7
13: [1, 2, 3, 4, 6, 8, 9, 10, 12], len=9

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

Узнаем какой единственный остаток по модулю 2*3*5=30:

? chinese([Mod(1,2),Mod(2,3),Mod(1,5)])
%1 = Mod(11, 30)

Т.е. в таблице будет пока единственное число 11.
Посмотрим что за 3 остатка по модулю 30*7=210 нам надо получить:

? foreach([2,3,4],y,print1(chinese(Mod(11,30),Mod(y,7)),", "))
Mod(191, 210), Mod(101, 210), Mod(11, 210),

Т.е. в таблице должны получиться числа 191,101,11 в любом порядке.

Это пока были образцовые данные, те что и надо получить.
Вспоминаем формулу r[]=x[]+Mx*(((y[]-x[])*d)%My). Здесь x[]=[11], Mx=30, y[]=[2,3,4], My=7, вычисляем d=Mod(1/Mx,My)=Mod(4,7)=4. Пробуем получить:

? 11+30*(((2-11)*4)%7)
%1 = 191
? 11+30*(((3-11)*4)%7)
%2 = 101
? 11+30*(((4-11)*4)%7)
%3 = 11

Получили, КТО по такой формуле работает (ну никто и не сомневался, и в вики и вообще везде она и приведена под названием алгоритм Гарнера).

А теперь попробую (самому интересно) получить то же самое через A[]=Mx*((y[]*d)%My) и B[]=Mx*((x[]*d)%My), так как все y[] меньше единственного x[], то условие y[]<x[] выполняется всегда и значит четвертое слагаемое должно быть всегда:

? 11+30*((2*4)%7)-30*((11*4)%7)+30*7
%1 = 191
? 11+30*((3*4)%7)-30*((11*4)%7)+30*7
%2 = 311
? 11+30*((4*4)%7)-30*((11*4)%7)+30*7
%3 = 221

А вот и глюк. Условие y[]<x[] вовсе не равнозначно условию A[]<B[]. Хотя по модулю 210 получили ровно те же числа, почти правильные, надо лишь с условным слагаемым разобраться.
Посмотрим чему равны A[] и B[]:

? B=30*(([11]*4)%7)
%1 = [60]
? A=30*(([2,3,4]*4)%7)
%2 = [30, 150, 60]

И действительно, условие A[]<B[] выполняется лишь для первого числа и не выполняется для двух остальных, которые и получились на Mx*My больше правильных.
Т.е. для правильного вычисления КТО по такой формуле, с A[] и B[], надо использовать именно условие A[]<B[], как у меня выше в формуле и было написано. Проверим:

? foreach([2,3,4],y, A=30*((y*4)%7); B=30*((11*4)%7); print(11+A-B+if(A<B,30*7)); )
191
101
11

Вуаля, с КТО всё получилось.


Дальше речь уже не про КТО, а про мою программу, которой эти числа в общем и не нужны, а нужны остатки от них по модулю например 71:

? [191,101,11]%71
%1 = [49, 30, 11]

Именно этот вектор надо получить имея на руках x[]%71=[11]%71=[11] и y[]%71=[2,3,4]%71=[2,3,4].
Предполагаем что x[] и B[] и Mx*My и все прочие константы можем посчитать вне цикла, в цикле же будем перебирать лишь y[]. Сразу видим что сам y[] нам нафиг не сдался, в формуле x[]+A[]-B[]+if(A[]<B[],Mx*My) вместо него участвует только A[]=Mx*((y[]*d)%My) и логично именно его и хранить как A[]=[30, 150, 60].
Далее, x[] и B[] внутри цикла по y[] (A[]) не меняются, значит их сумму вычисляем заранее.
А раз нас интересует только остаток всего этого дела по модулю 71, то и сумму и A[]%71=A1[] вычислим заранее.
Приходим к следующей программе:

? x=11; B=30*((x*4)%7); S=(x+B)%71
%1 = 0
? (30*7)%71
%2 = 68
? A=[30,150,60]; A1=A%71
%3 = [30, 8, 60]
? for(i=1,#A, r=S+A1[i]+if(A[i]<B,68,0); print(r); )
98
8
60

Вообще ни разу не совпадают с требуемыми [49, 30, 11], даже по модулю 71.
Чешу репу, вот же засада, а я уже кинулся код на асме писать ...
Но формула в принципе правильная:

? for(i=1,#A, r=11%71+A[i]%71-B%71+if(A[i]<B,(30*7)%71,0); print(r); )
49
-41
11
? for(i=1,#A, r=11%71+A[i]%71-B%71+if(A[i]<B,(30*7)%71,0); if(r<0, r+=71); print(r); )
49
30
11

Осталось понять какое из сложений вызывает глюк.
Оказалось объединение x[]-B[] по общему модулю 71:

? for(i=1,#A, r=(11-B)%71+A1[i]+if(i<2,(30*7)%71,0); if(r<0, r+=71); print(r); )
120
30
82

Но при этом цифры почти правильные - по модулю 71 это именно то что нужно. Т.е. достаточно все два сложения (с A1[] и с if) выполнить по модулю 71 и получим искомое. Развернув цикл и учтя какие индексы дают именно A[]<B[] получим:

? ((11-60+68)%71+A1[1])%71
%1 = 49
? ((11-60)%71+A1[2])%71
%2 = 30
? ((11-60)%71+A1[3])%71
%3 = 11

Ура! Именно то что нужно.
При этом разумеется левое слагаемое вычисляем заранее, не в цикле перебора A1[] (y[]). А так как мы всегда можем отсортировать A1[], то будет точная граница индекса, до которой слева будет одно слагаемое, а после неё другое, это легко реализуется.
Остаётся ровно одно сложение по модулю, что как уже известно реализуется тремя командами AVX за 1 такт (4 обычных команды c cmov). И именно его и делает текущий выполняемый код. Надо только добавить ему не одно левое слагаемое, а два разных, до границы и после. Границу можно вычислять не в самом цикле, а снаружи, поиском по A[], это быстро.
Заметим что вычисление левого слагаемого по имеющемуся x[] достаточно просто, два деления, два умножения, несколько сложений, поиск в массиве. Всё это займёт ну пусть полтыщи тактов, но выбрав длину y[] (A[], A1[]) достаточно большой чтобы его обработка занимала хотя бы полсотни тысяч тактов, накладные расходы на вычисление левого слагаемого уйдут ниже 1% и ими можно будет пренебречь.
А далее, вовсе не обязательно x[] хранить, его тоже можно получать по КТО из остатков по любым простым кроме используемых в y[],My. КТО же даже в исходном виде x[]+Mx*(((y[]-x[])*d)%My) требует всего одно деление и два умножения, выбирая всегда чтобы и Mx и My на каждом шаге (для каждой пары простых или накопленного вектора и нового произведения простых) помещались в регистры (x64 разумеется) это вычисление довольно быстрое. А если даже и не слишком быстрое, то всегда можно увеличить y[] для внутреннего цикла настолько чтобы его обработка занимала более 99% общего времени, для этого достаточно длин порядка сотен тысяч, в много миллионов уходить не надо.
На этом для меня всё стало ясно, дальше уже оптимизация (ассемблерного) кода.

-- 26.01.2024, 16:13 --

Уф, получив по модулю чушь реально подумал всё, кирдык, ещё одна прекрасная идея улетела в туман, нифига не работает ...
Стал судорожно усложнять исходную очевидно работающую формулу КТО буквально по мелким шажкам (большинство из которых показывать уже не стал) и таки обнаружил затык. И при его лечении оказалось я был таки прав изначально, всё получается как надо, код внутреннего цикла укладывается в одно сложение по модулю, как и сейчас. С дополнительным условием, но которое можно вынести и за цикл, скорость оно не испортит. Теперь как это всё реализовать, в коде полно своих заморочек ...

[Dmitriy40]

И кстати, никто не мешает применить этот же подход и к PARI, не парясь про какие-то то там остатки (по 71 в примере): перебирать блоки с таблицей скажем 37# внешними циклами не линейно, а формируя правильную "базу" применением КТО (даже стандартной chinese без заоморочек) для любых не задействованных простых вместо индекса множимого на 37#, в самом же цикле прохода по таблице 37# для вычисления КТО достаточно условного сложения, при том что обычное там и так есть. Вот наверняка вычисление r=A[j]+if(A[j]<B,xBQ,xB) не столь уж тормознее существующего r=m[j]+i*37# (A[] посчитан из m[] 1-в-1 показанным выше про y[]->A[] способом, xB=x-B и xBQ=x-B+37#). При этом x, B, xB, xBQ надо перевычислять лишь один раз на всю таблицу 37#, в обрамляющем цикле, это вообще ничтожное время займёт (в процентах от времени обхода таблицы 37#).
Т.е. скажем до 59#=2e21 перебирать всю таблицу 37# не 59#/37#=41*43*47*53*59 раз, а лишь 24*24*28*34*40 раз (для паттерна 19-252), почти в 12 раз быстрее. Предел перебора можно гибко регулировать выбирая какие простые учитывать (никто же не заставляет прям исключительно только наименьшие подряд, можно любые, просто для меньших выигрыш скорости больше). Вот только ограничить его снизу не получится, хотя лишний if(r<W, next) во внутреннем цикле тоже не сильно вредит (по сравнению с ispseudoprime там дальше).
При большом желании один-два самых верхних цикла (по самым большим простым) можно оставить и линейными, всё равно там выигрыш самый маленький (хоть и есть всегда), тогда можно перебирать к примеру кусками по 59#=2e21, не получая выигрыша 61/42=1.45 раза (для 19-252), зато выигрыш почти 12 раз (для таблицы 37#) никуда не денется. В 1.45 раза медленнее, зато 61 блок по 2e21 подряд.

И чем длиннее паттерн, тем больше будет выигрыш, ведь для каждого (достаточно большого, для 19-252 от 43) простого p он ровен p/(p-len) раз.

[Yadryara]

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

Или Вам получение m[] непонятно?

Понятно, я сам недавно приводил способ:

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

Второе, 3-е, 4-е и 5-е числа больше начального соответственно на 24, 30, 36 и 60. Посмотрим на ближайшие к ним сверху числа кратные 7 и вычитанием определим запрещённый остаток:

$28-24=4$

$35-30=5$

$42-36=6$

$63-60=3$

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

Дальнейшее ещё не всё внимательно посмотрел.

[Yadryara]

Вроде разобрался, круто!

Ещё и приятную для Вас ошибку нашёл:

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

Но зато приятные следствия - в повышении скорости перебора: для таблицы 37# с 6.636e6 строками до полных 67#=7.858e24 надо перебрать не 41*43*47*53*59*61*67 блоков, а лишь 24*24*2834*40*42*48, в 8 раз меньше (т.е. быстрее).

Не в 8, а почти в 24 раза меньше (быстрее) !

[Dmitriy40]

И правда, знак умножения пропал ...
Дочитайте что и в PARI это (то что касается вычисления КТО) в общем применимо, для поиска по любым паттернам. И тоже не требует сверхбольших таблиц (думаю 37# или даже 31# для 19-252 будет достаточно). Те же Ваши кэфы наверняка так можно считать заметно быстрее.

-- 26.01.2024, 21:11 --

Как всем заинтересованным известно моя программа ищет вовсе не КПППЧ, а более широкий класс цепочек, в который включена и искомая КПППЧ19d252 (19-252), так что она пропущена не будет. Но вот что при этом найдётся (и что не найдётся) другого - вопрос (с известным ответом разумеется). В частности решил проверить не нашлись ли какие-то незамеченные 17-ки с центром в центре искомой 19-ки, но с любым паттерном. Перепроверил логи и да, нашёл две новые (не показанные ранее как именно 17-ки):
309436831159079027513491: [0, 6, 12, 30, 42, 72, 90, 96, 126, 156, 162, 180, 210, 222, 240, 246, 252], n=17
709403968751236831236377: [0, 6, 12, 30, 42, 72, 90, 96, 126, 156, 162, 180, 210, 222, 240, 246, 252], n=17
Первая была показана раньше, вторая на форуме не появлялась (поиск не находит).
Обе не входят в искомую 19-ку (имеют две дырки сразу вокруг центра +126).

[Yadryara]

Так это уникальный результат: 17-252. В TBEG и SPT вместе 35 17-к. Самый маленький диаметр из этих 35-ти — 312.

[Dmitriy40]

Так то с решетом, Врублевский на какой-то из прошлых конкурсов нашёл три 17-ки с диаметром 252 (первую чуть выше 4e20, все с разными паттернами) и ещё штук 25 других, все до 3.3e22. Мне же ни одна из этих 60 с лишним 17-ек не попалась.

[Dmitriy40]

Придумал как заставить КТО выдавать числа по возрастанию и заодно ограничить проверяемый диапазон и снизу и сверху.
В формуле x+A-B+if(A<B,M) мы всегда можем отсортировать A по возрастанию. Тогда цикл по A разбивается на два цикла подряд: пока A<B и остальное (один из них может быть и пустым). При этом A только прибавляется и потому в каждом из двух этих новых циклов числа будут по возрастанию. Дублируем вектора x и B в новые удвоенного размера x1 и B1, в половине x1=x-B, в другой половине x1=x-B+M, в B1=B в обоих половинках, создаём вектор f1 флагов добавленного M тоже двойного размера (первая половина f1=0, вторая f1=1), сортируем x1 по возрастанию, переставляя в том же порядке и f1 и B1. Перебираем x1 в порядке возрастания, для каждого берём значение B из B1 и из f1 признак какую часть A из двух перебирать (меньше B или нет). Всё, получили все значения в порядке возрастания. Ограничить диапазон теперь уже тривиально.
Внутри циклов никаких сортировок нет, они все снаружи и однократно и на каждый вектор отдельно (отсортировать два вектора размерами 2*X и Y сильно быстрее сортировки вектора размером X*Y).

Но на три и более векторов не обобщается.

[Yadryara]

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

Дочитайте что и в PARI это (то что касается вычисления КТО) в общем применимо, для поиска по любым паттернам. И тоже не требует сверхбольших таблиц (думаю 37# или даже 31# для 19-252 будет достаточно). Те же Ваши кэфы наверняка так можно считать заметно быстрее.

Дочитал, но может лучше на примере 7-ки и 23# покажете? Самая быстрая прога сейчас такая

(PARI)

allocatemem(2^27);
{print();
tz0=getwalltime;
v=[0,6,30,36,42,66,72];
kp=0;k3=0;k5=0;k7=0;k9=0;
start=358705*10^11;fin=start+10*10^9;
w=2*3*5*7*11*13*17*19*23; \\*29;

ko=1;m=vector(23,i,[]);
forprime(p=2,#m, m[p]=setminus(vector(p,i,i-1),Set(-v%p));
printf("%d: x %d: %d\n",p,#m[p],m[p]);
ko=ko*#m[p]);o=vector(ko);print();

print(w);print();
print(ko);

x0=Mod(1,2);
foreach(m[23],m23,
foreach(m[19],m19,
foreach(m[17],m17,
foreach(m[13],m13,
foreach(m[11],m11,
foreach(m[7],m7,
foreach(m[5],m5,
foreach(m[3],m3,
x=lift(chinese([x0,Mod(m23,23),Mod(m19,19),Mod(m17,17),Mod(m13,13),Mod(m11,11),Mod(m7,7),Mod(m5,5),Mod(m3,3)]));
no++;o[no]=x))))))));

print(no);print();

tz1=getwalltime;
t1=tz1-tz0;
print(floor(t1/1000),".",floor((t1-floor(t1/1000)*1000)/10), " seconds");
print();

mm=vector(200,i,[]); forprime(p=28,#mm, mm[p]=Set(-v%p);
\\print(p, "   ", #mm[p], "   ", mm[p]);print();
);


for(y=floor(start/w), floor(fin/w)+1,
for(i=1,ko,
kan=y*w+o[i];
if(kan<start || kan>=fin, next);
kkan++;

forprime(p=28,#mm, if(setsearch(mm[p], kan%p), next(2)));

if(ispseudoprime(kan)
&& nextprime(kan+1)-kan==v[2]
&& nextprime(kan+v[2]+1)-kan==v[3]
&& nextprime(kan+v[3]+1)-kan==v[4]
&& nextprime(kan+v[4]+1)-kan==v[5]
&& nextprime(kan+v[5]+1)-kan==v[6]
&& nextprime(kan+v[6]+1)-kan==v[7],

di=nextprime(kan+v[6]+1)-kan;

\\if(kan>=start && kan<=fin,
k7++;
print1(k7,"   ");
forprime(p=kan,kan+v[1],print1(p,"   "));
print(di);

)));
print();print(kkan);print();
print(k7);

tz2=getwalltime;
t2=tz2-tz1;
print();
print(floor(t2/1000),".",floor((t2-floor(t2/1000)*1000)/10), " seconds");
print();


}quit;

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

Врублевский на какой-то из прошлых конкурсов нашёл три 17-ки с диаметром 252 (первую чуть выше 4e20, все с разными паттернами) и ещё штук 25 других, все до 3.3e22.

Да, я забыл, что он не только 17-240 нашёл.

[Dmitriy40]

Отчего же не показать, глядишь ещё какую ошибку у себя найду.
Перебирать в тесте буду полный диапазон от 0 до 31#=2e11 - ну вот так проще совместить два метода перебора.
Выигрыш от таблицы 23#, равный $\frac{31}{31-7}\frac{29}{29-7}=1.7$, как-то маловат, потому ограничу таблицу 19#, тогда выигрыш составит $\frac{31}{31-7}\frac{29}{29-7}\frac{23}{23-7}=2.45$ раза, уже что-то.
Вот сокращённая (не люблю простыни в логах) и несколько оптимизированная выдача вашей программы (под x64):

19#: 34560/9699690, time: 57 ms
714597120/3546, time: 19min, 53,144 ms

Теперь меняем строки

for(y=floor(start/w), floor(fin/w),
   for(i=1,ko,
      kan=y*w+o[i];

На новые

Mx=23*29*31; Q=Mx*w; d=lift(Mod(1/Mx,w)); o=vecsort(Mx*((o*d)%w));
m23=setminus(vector(23,i,i-1),Set(-v%23));
m29=setminus(vector(29,i,i-1),Set(-v%29));
m31=setminus(vector(31,i,i-1),Set(-v%31));
xx=vector(Mx); n=0; foreach(m23,a, foreach(m29,b, foreach(m31,c, xx[n++]=lift(chinese([Mod(a,23),Mod(b,29),Mod(c,31)])); ))); xx=xx[1..n];
for(ix=1,#xx, x=xx[ix];
   B=Mx*((x*d)%w); xB=x-B; xBQ=xB+Q;\\Специально под терминологию формулы КТО в моих сообщениях выше
   for(i=1,ko,
      kan=o[i]+if(o[i]<B,xBQ,xB);

Обратите внимание на пересчёт o[], теперь в нём нужны A[] (в терминологии формул КТО у меня выше). Сортировка по возрастанию - лишняя, но на 35тыс элементов почти мгновенна.
В строке чуть ниже

forprime(p=20,#mm, if(setsearch(mm[p], kan%p), next(2)));

число 20 меняем на 32 (меньшие простые проверены циклом по xx[]).

Запускаем:

19#: 34560/9699690, time: 56 ms
291962880/3546, time: 10min, 12,178 ms

Количество найденных цепочек ровно то же, я их не сравнивал, верю что те же.

Расчётное ускорение за счёт уменьшение количества итераций внутреннего цикла $\frac{714597120}{291962880}=2.45$, реальное ускорение $\frac{19m53s}{10m12s}=1.95$. Принцип рабочий.
Почему столь велика разница (26%) разбираться с PARI что-то влом, слишком он непредсказуем в плане скорости.
Плюс я не совсем уверен что запускаемая прога на 100% соответствует вашей - например в цикле for(y=floor(start/w), floor(fin/w), не вижу смысла добавлять 1, более того, и fin тут (или в самом начале) надо бы сделать на 1 меньше, иначе будет перебираться на 1 блок 19# больше необходимого (пусть и лишь до первого if, однако и это время). Но никаких глобальных или принципиальных отличий быть не должно.

[Yadryara]

Вот ещё момент, где можно выиграть в скорости.

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

число 20 меняем на 32 (меньшие простые проверены циклом по xx[]).

Можно ведь и 36 и даже 37 взять.

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

Естественно проверку по setsearch переносил: 24-200 и 20-200 соответственно.

А можно и 29-200 и 23-200 брать. Или 199? Замена 24-200 на 28-200 (которая видна в моём коде) дала ускорение не меньше 15%.

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

Почему столь велика разница (26%) разбираться с PARI что-то влом,

Одну из причин вроде нашёл. См. выше.

[Dmitriy40]

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

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

число 20 меняем на 32 (меньшие простые проверены циклом по xx[]).

Можно ведь и 36 и даже 37 взять.

Там forprime, потому 32=36=37. Думаю скорость определения что 32-36 не простые настолько высока, что её хрен измеришь средствами PARI:

? n=0; for(i=1,10^6, forprime(p=37,1000, n++););
time = 20,748 ms.
? n=0; for(i=1,10^6, forprime(p=32,1000, n++););
time = 20,530 ms.

Случайно оказалось даже быстрее, чего ну никак быть не должно - явный признак недостоверности измерений.
И даже если бы вдруг 0.2с оказались в правильную сторону, то даже на 30млн итераций это всего лишь 6с или около 1% от 612с общего времени - просто ни о чём, сравнимо с погрешностью.

[Yadryara]

Да, закрыл Хром, стал менять с 24 на 29 и перепроверять. 29 выигрывает, но только на доли процента.

[Yadryara]

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

Расчётное ускорение за счёт уменьшение количества итераций внутреннего цикла $\frac{714597120}{291962880}=2.45$, реальное ускорение $\frac{19m53s}{10m12s}=1.95$. Принцип рабочий.

Проверил. У меня при печати всех 3546 начальных чисел 7-к отработала за 28.1 минут, а при печати только итоговых количеств — за 27.1 минут. Спасибо.

Так что те кто считает на PARI, а таких нынче несколько человек, пользуйтесь добротой Dmitriy40. Берите алгоритм и считайте.

За счёт перехода на более высокие периоды (вплоть до 41#) вы уже ускорили счёт в десятки раз, надо ускорить ещё в сотни. А иначе даже 17-ку за год вам не найти.