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

[Dmitriy40]

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

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

Ну в общем вот код с модульной арифметикой вместо forprime.
...

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

 \\ вычисляем положение числа которое делится на это простое относительно начала цепочки
 d=len-(pmods_n0[i]+2*k*pmods_m[i]+len-1)%pr[i];

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

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

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

 n=n0+2*k*m+d-1; \\ вычислим текущее число 
 \\Если уже перебор множителей, то отбрасываем всю цепочку
 if(nf[d]==nu[d],q1++;if(!ispseudoprime(n/(v[d]*nn[d])), next(2)));

Здесь вычисление n лучше внести под первый if, к q1++ или сразу в ispseudoprime, оно ведь понадобится лишь для ispseudoprime, иначе в этом цикле не нужно.
Кроме того, здесь потеряли второй параметр 1 у ispseudoprime - а это нехорошо, это её ускоряет процентов на 10-15 судя по тестам.

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

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

if(nf[i]<nu[i] && ispseudoprime((n0+2*k*m+i-1)/(v[i]*nn[i])),q1++;next(2)));

Непонятно что здесь считается в q1, выше это был счётчик вызовов ispseudoprime, здесь же счётчик сколько отброшено условием nf<nu. Не стоит смешивать разные вещи в одном счётчике, заведите лишний, это не больно.

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

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

for(k=k_start,k_start+k_len, \\ идём по цепочкам

Здесь должно быть +k_len-1.

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

Показатели такие, на миллион цепочек.

А каковы показатели эталона? Пусть только в один поток, многопоточность допилим отдельно.

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

К сожалению, как запустить многопоток "эталона" пока не очень представляю

Да так же как и этот - переделать в функцию и сравнить что фильтрация на 10млн и время (с точностью до погрешности) не изменились.
Но я не понял как Вы запускаете многопоток, что, несколько потоков вызывают функцию с разными параметрами?
А почему бы не сделать многопоточным цикл for(k=k_start,k_start+k_len, ...)? Локальных переменных у него немного, k,nf,nn,d,n,q,q1,q2,c,i и они все или скалярные или короткие векторные, остальные объявить в export. Придётся извратиться с передачей счётчиков qX наверх, но это не смертельно.

[wrest]

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

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

Увы, но и от длинной отказаться не удалось. Я пробовал типизировать переменные, которые там участвуют, но что-то идёт не так и там падает с segmentstion fault. Так что там по-прежнему длинная арифметика.

-- 06.12.2025, 12:45 --

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

Здесь вычисление n лучше внести под первый if, к q1++ или сразу в ispseudoprime, оно ведь понадобится лишь для ispseudoprime, иначе в этом цикле не нужно.

Я в какой-то из моментов замерял это вычисление, но оно оказалось очень быстрым и не влияющим на общее время ( проценты).

-- 06.12.2025, 12:48 --

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

Кроме того, здесь потеряли второй параметр 1 у ispseudoprime - а это нехорошо, это её ускоряет процентов на 10-15 судя по тестам.

А, я обратил внимание что у вас с параметром вызывается, но не понял зачем, всё хотел спросить, и убрал.

-- 06.12.2025, 13:00 --

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

Да так же как и этот - переделать в функцию и сравнить что фильтрация на 10млн и время (с точностью до погрешности) не изменились.

Да, я даже думаю что код с forprime немного быстрее. Идея с отказом от длинного деления пока проваливается -- оно все равно там осталось. Так что этот мой последний код с модулями только как proof of concept.

-- 06.12.2025, 13:10 --

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

Но я не понял как Вы запускаете многопоток, что, несколько потоков вызывают функцию с разными параметрами?

Parfor и каждой функции по своему миллиону цепочек (третий аргумент k_len равен 10^6)

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

А почему бы не сделать многопоточным цикл for(k=k_start,k_start+k_len, ...)? Локальных переменных у него немного, k,nf,nn,d,n,q,q1,q2,c,i

Так ведь параллелить стоит то что долго работает, а тут у нас в среднем от 50 до 100 микросекунд на цепочку.

[Dmitriy40]

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

А, я обратил внимание что у вас с параметром вызывается, но не понял зачем, всё хотел спросить,

Я уже объяснял:

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

Можно наверное немного ускорить заменив nf==4 && !ispseudoprime(x) на nf==4 && !ispseudoprime(x,1) (что почти эквивалентно nf==4 && lift(Mod(2,x)^(x-1))<>1), так не будет точной проверки на простоту, но она именно тут и не нужна, тут нужна точная проверка на составное, а обе конструкции именно это гарантированно и делают. Какая из них быстрее после компиляции мне неведомо.

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

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

Смысл в том что ispseudoprime без 1 выполняет сначала тест Ферма (Миллера-Рабина) с основанием 2, а потом ещё и тест Люка. Если же указать второй параметр, то ограничится лишь указанным числом тестов Ферма (со случайным основанием), чего в 99.9999% случаев достаточно для определения что число составное (что простое - недостаточно). Но именно в этом месте нам надо знать как раз что число составное и потому можно ограничиться одной итерацией теста Ферма (с любым основанием).
В проверке nf<nu так делать нельзя, там надо быть уверенным в простоте числа, а не что оно составное (для ispseudoprime это не взаимоисключающие вещи).

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

Так ведь параллелить стоит то что долго работает, а тут у нас в среднем от 50 до 100 микросекунд на цепочку.

Это да.
Тогда разбивать этот for на два и параллелить внешний. Ну как Вы собственно и делаете.

[wrest]

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

Непонятно что здесь считается в q1, выше это был счётчик вызовов ispseudoprime, здесь же счётчик сколько отброшено условием nf<nu.

Поскольку там условие nf<nu AND составное, то счётчик не срабатывает если простота не проверялась вовсе, но не учитывает случаи когда проверялась но результат оказался составным, тут вы правы, надо было последовательно двумя if-ами.

[Yadryara]

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

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

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

Кстати, почему за эталон взяли такой пример. Ну да, из-за рекордной 21-ки.

Но сейчас поиск уже весь идёт только среди 52-значных чисел. Может лучше было бы взять недавно найденную рекордную (для А-файла https://oeis.org/A292580/a292580_12.txt)) 20-ку:

1343688799116140536185862959273033940943830770546841

Это было бы гораздо ближе к телу.

Ждал что wrest прокомментит мой сбой, но увы. Возможно, что это из-за установки обновления (без моего участия) ?

Перезапустил. В новом варианте время заметно нелинейно. Так что считая повыше, удалось достичь ускорения в 2.4 раза в сравнении со старым способом, близким к эталону. Или к этанолу...

[wrest]

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

Ждал что wrest прокомментит мой сбой, но увы. Возможно, что это из-за обновления?

У меня такого никогда не случалось. Но у меня wsl2 а не wsl1, и я не гоняю wsl с утра до вечера :)
Заморочиться имеет смысл если сбой будет повторяться.

[Dmitriy40]

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

Кстати, почему за эталон взяли такой пример.

Просто первым попался на глаза. И по нему больше вариантов программы было сделано для разных тестов и они остались лежать в легко доступном месте, вот и попались первыми.

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

Но сейчас поиск уже весь идёт только среди 52-значных чисел.

В чём проблема поменять диапазон i в цикле ...

[wrest]

Dmitriy40
Я таки победил в борьбе за короткую арифметику.

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

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

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

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

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

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

print(myname, " приступил. Задание: Старт: ",k_start," длина:",k_len," Таблица простых до:", lim," Подготовился за:",strtime(getwalltime()-t0));

\\ основной цикл
t0=getwalltime();\\Замер времени отсеивания
for(k:small=k_start,k_start+k_len-1, \\ идём по цепочкам
\\Здесь начинается проверка конкретного n
   nf=vectorsmall(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((n0+2*k*m+d-1)/(v[d]*nn[d]),1), 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))); 
        if(nf[i]<nu[i],q1++;
            if(ispseudoprime((n0+2*k*m+i-1)/(v[i]*nn[i])),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,". Cкорость счета: ", 1000*k_len\(getwalltime()-t0)," цеп/сек. Время счёта:",strtime(getwalltime()-t0));
return(q2);
}

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

(Оффтоп)

? parfor(i=0,7,em3(i,i*10^6,10^6,2^20))
2 приступил. Задание: Старт: 2000000 длина:1000000 Таблица простых до:1048576 Подготовился за:15 ms
1 приступил. Задание: Старт: 1000000 длина:1000000 Таблица простых до:1048576 Подготовился за:16 ms
4 приступил. Задание: Старт: 4000000 длина:1000000 Таблица простых до:1048576 Подготовился за:17 ms
0 приступил. Задание: Старт: 0 длина:1000000 Таблица простых до:1048576 Подготовился за:19 ms
5 приступил. Задание: Старт: 5000000 длина:1000000 Таблица простых до:1048576 Подготовился за:20 ms
6 приступил. Задание: Старт: 6000000 длина:1000000 Таблица простых до:1048576 Подготовился за:20 ms
7 приступил. Задание: Старт: 7000000 длина:1000000 Таблица простых до:1048576 Подготовился за:21 ms
3 приступил. Задание: Старт: 3000000 длина:1000000 Таблица простых до:1048576 Подготовился за:32 ms
1 готов. Найдено:0. Прошло:0. Проверено псевдопростых:1098511. Cкорость счета: 33040 цеп/сек. Время счёта:30,266 ms
2 готов. Найдено:0. Прошло:1. Проверено псевдопростых:1098296. Cкорость счета: 31093 цеп/сек. Время счёта:32,161 ms
6 готов. Найдено:0. Прошло:1. Проверено псевдопростых:1097228. Cкорость счета: 30142 цеп/сек. Время счёта:33,176 ms
3 готов. Найдено:0. Прошло:1. Проверено псевдопростых:1097853. Cкорость счета: 30109 цеп/сек. Время счёта:33,212 ms
0 готов. Найдено:1. Прошло:2. Проверено псевдопростых:1100107. Cкорость счета: 29754 цеп/сек. Время счёта:33,608 ms
5 готов. Найдено:0. Прошло:2. Проверено псевдопростых:1097272. Cкорость счета: 29183 цеп/сек. Время счёта:34,266 ms
7 готов. Найдено:0. Прошло:2. Проверено псевдопростых:1096437. Cкорость счета: 28946 цеп/сек. Время счёта:34,547 ms
4 готов. Найдено:0. Прошло:3. Проверено псевдопростых:1097228. Cкорость счета: 28452 цеп/сек. Время счёта:35,146 ms
cpu time = 4min, 26,196 ms, real time = 35,176 ms.
?

В общем, результат теперь, после компиляции, в многопоточном режиме 230 тыс цепочек в секунду в сумме. 8 миллионов за 35 секунд.
Оптимальным оказался предел по таблице простых 2^20, поскольку я факторизую numdiv-ом все прошедшие через фильтр числа. До типизации и избавления от длинных делений рекорд (по модульному алгоритму) был 95 тыс цепочек в секунду, т.е. рост в скорости более чем двукратный! До типизации, длинная арифметика:

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

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

Думаю, что подход уйти от длинных делений в циклах forprime был все-таки верным.
Результат достойный, как мне кажется.

Теперь хотелось бы, наверное, по максимуму избавиться от проверки на [псевдо]простоту, через которую проходит, в среднем, одно число в цепочке. И при этом мы знаем, что в проверяемых числах нет простых множителей меньше 2^20.

[Dmitriy40]

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

Результат достойный, как мне кажется.

Более двухкратного роста - это прекрасно!

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

Теперь хотелось бы, наверное, по максимуму избавиться от проверки на [псевдо]простоту, через которую проходит, в среднем, одно число в цепочке. И при этом мы знаем, что в проверяемых числах нет простых множителей меньше 2^20.

Можно вместо !ispseudoprime(t,1) использовать lift(Mod(2,t)^(t-1))<>1 как говорил выше, это отличается собственно лишь отсутствием предварительных делений на малые простые (мы же сами их уже исключили), но в интерпретаторе я заметного ускорения не увидел, попробуйте скомпилить.
Никаких более быстрых способов проверки простоты я не знаю, они в общем все основаны на малой теореме Ферма $a^{p-1}=1\pmod p$ или используют её как первый шаг, а это и есть возведение длинного числа (оно сначала короткое, но по мере возведения в степень становится длинным) в степень длинного числа по модулю длинного числа, кошмар-кошмар.

-- 06.12.2025, 18:30 --

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

через которую проходит, в среднем, одно число в цепочке.

Это кстати отдельно интересно, выходит около 5% чисел имеют достаточно много малых делителей чтобы накопить nf==nu.

[wrest]

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

Никаких более быстрых способов проверки простоты я не знаю, они в общем все основаны на малой теореме Ферма $a^{p-1}=1\pmod p$ или используют её как первый шаг, а это и есть возведение длинного числа (оно сначала короткое, но по мере возведения в степень становится длинным) в степень длинного числа по модулю длинного числа, кошмар-кошмар.

Ну так речь не о том чтобы ускорить саму проверку, а том чтобы поменьше проверять, или проверять числа поменьше если есть выбор и т.п.
Т.е. нужен анализ а почему мы должны проверять, нет ли ещё каких дефектов у цепочки.

[wrest]

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

Это кстати отдельно интересно, выходит около 5% чисел имеют достаточно много малых делителей чтобы накопить nf==nu.

Ну это как минимум, мы же делаем ранний выход. Может их и больше, несложно посчитать через numdiv :). Вообще вся эта проверка же и основана на том, что в среднем, в каждой цепочке как минимум одно число избыточно по количеству делителей, причём малых. Это легко проверить, что мы и делаем. Цепочки с недостаточными числами мы выявляем позже, на последнем этапе проверки, и их оказывается немного...

[wrest]

Посчитал такую статистику. Простое число, на котором в среднем прекращается проверка цепочки, равно 593. Вычисляется как сумма простых, на которых проверка прекратилась и принято решение цепочку отбросить, делить на количество цепочек, отброшенных по перебору делителей (матожидание). А простых чисел между 59 и 593 меньше сотни.
Можно посчитать статистику по каждому простому до 2^20 -- сколько раз на него что-нибудь разделилось, сколько раз оно прервало проверку и т.п. Исходя из этого, можно заранее поменять порядок следования простых в таблице простых и делить сперва на самые ходовые, обеспечивая ещё более ранний выход.
Только нужно ли...
В принципе, написав прогу-бенчмарк паттерна и прогнав через неё скажем первый миллион или миллион случайных цепочек этого паттерна, можно вычислить подстроечные параметры (тот же порядок следования в таблице простых), которые ускорят проверку цепочек этого паттерна.
А если согласовать набор подстроечных параметров не зависящий от задачи (т.е. набор тот же и для D48 и для D24: искуемая длина цепочки, параметры прогрессии цепочки, вектор встроенных в цепочку множителей, порядок следования пробных простых, лимит по простым, имя и параметры логирования/логфайла и т.п.), можно достигнуть универсальности и считать задачи с разные количествами делителей с одинаковой эффективностью.
Можно наверное придумать и какой-то анализатор паттернов, который посчитает разложения случайного миллиона цепочек и подскажет "хороший" ли получился паттерн, посчитает для него оптимальные параметры перебора и т.п.
...
...
...
profit!

[Dmitriy40]

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

Можно посчитать статистику по каждому простому до 2^20 -- сколько раз на него что-нибудь разделилось, сколько раз оно прервало проверку и т.п. Исходя из этого, можно заранее поменять порядок следования простых в таблице простых и делить сперва на самые ходовые, обеспечивая ещё более ранний выход.

Теория говорит что получится список возрастающих простых - как уже и есть.

Потому что m взаимно простое с любым простым больше 53 и значит остатки по таким простым от n распределены равномерно и значит остаток 0 встречается не чаще других и значит вероятность его обнаружить равна 1/p для простого p. Практика прекрасно согласуется с теорией (штуки на миллион):

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

...
frq=matrix(100,100);
for(i=0,oo,
...
        forprime(p=59,100,
                frq[p,(n+#v-1)%p+1]++;
        );
);
forprime(p=59,100,printf("%2u: %6u\n",p,frq[p,1..p]););

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

59: [ 16943, 16960, 16948, 16946, 16944, 16959, 16944, 16947, 16946, 16939, 16967, 16953, 16957, 16945, 16943, 16947, 16943, 16951, 16958, 16954, 16951, 16952, 16947, 16955, 16950, 16961, 16933, 16949, 16945, 16933, 16960, 16952, 16951, 16968, 16949, 16935, 16951, 16940, 16943, 16956, 16933, 16950, 16952, 16933, 16946, 16951, 16948, 16956, 16940, 16934, 16942, 16946, 16953, 16959, 16958, 16954, 16963, 16961, 16946]
61: [ 16363, 16401, 16399, 16375, 16394, 16389, 16390, 16393, 16393, 16402, 16383, 16395, 16394, 16383, 16397, 16395, 16395, 16394, 16385, 16384, 16398, 16388, 16378, 16395, 16386, 16397, 16395, 16396, 16395, 16388, 16395, 16396, 16389, 16401, 16387, 16397, 16388, 16401, 16402, 16400, 16393, 16385, 16397, 16399, 16407, 16393, 16392, 16385, 16390, 16390, 16411, 16405, 16392, 16384, 16392, 16406, 16391, 16401, 16395, 16410, 16406]
67: [ 14929, 14919, 14926, 14919, 14939, 14927, 14918, 14926, 14926, 14938, 14907, 14911, 14931, 14922, 14919, 14929, 14916, 14927, 14936, 14924, 14924, 14918, 14937, 14917, 14924, 14924, 14926, 14923, 14915, 14926, 14932, 14921, 14932, 14925, 14921, 14919, 14942, 14944, 14931, 14919, 14921, 14927, 14925, 14933, 14917, 14923, 14922, 14928, 14912, 14926, 14934, 14929, 14921, 14926, 14916, 14920, 14932, 14929, 14931, 14930, 14920, 14932, 14938, 14916, 14935, 14926, 14922]
71: [ 14077, 14088, 14093, 14083, 14081, 14080, 14080, 14085, 14075, 14091, 14092, 14081, 14087, 14078, 14083, 14083, 14090, 14096, 14073, 14083, 14080, 14075, 14092, 14076, 14071, 14092, 14096, 14077, 14079, 14088, 14080, 14086, 14084, 14096, 14085, 14104, 14092, 14069, 14074, 14070, 14087, 14065, 14095, 14079, 14095, 14085, 14085, 14073, 14071, 14096, 14089, 14086, 14088, 14098, 14084, 14078, 14086, 14095, 14097, 14088, 14072, 14082, 14101, 14099, 14081, 14085, 14079, 14084, 14086, 14079, 14087]
73: [ 13697, 13695, 13706, 13703, 13690, 13689, 13690, 13693, 13691, 13686, 13708, 13703, 13692, 13694, 13690, 13694, 13694, 13712, 13699, 13696, 13711, 13695, 13696, 13704, 13693, 13695, 13700, 13705, 13700, 13711, 13711, 13703, 13703, 13690, 13692, 13696, 13699, 13695, 13705, 13702, 13697, 13695, 13700, 13707, 13705, 13698, 13691, 13713, 13690, 13698, 13693, 13701, 13720, 13690, 13693, 13695, 13710, 13704, 13706, 13701, 13694, 13708, 13688, 13697, 13694, 13706, 13689, 13711, 13687, 13717, 13694, 13684, 13696]
79: [ 12673, 12659, 12658, 12663, 12650, 12653, 12659, 12658, 12671, 12670, 12664, 12652, 12661, 12657, 12653, 12646, 12648, 12659, 12662, 12652, 12661, 12665, 12662, 12656, 12642, 12663, 12671, 12666, 12640, 12661, 12649, 12648, 12672, 12649, 12666, 12657, 12658, 12664, 12665, 12659, 12641, 12650, 12662, 12653, 12662, 12676, 12668, 12651, 12665, 12658, 12652, 12649, 12650, 12650, 12655, 12644, 12672, 12648, 12651, 12651, 12652, 12648, 12660, 12671, 12649, 12657, 12643, 12661, 12658, 12651, 12670, 12658, 12663, 12678, 12669, 12669, 12657, 12674, 12662]
83: [ 12051, 12058, 12052, 12047, 12055, 12050, 12044, 12056, 12039, 12045, 12047, 12051, 12044, 12035, 12049, 12044, 12035, 12063, 12032, 12032, 12045, 12059, 12058, 12055, 12053, 12042, 12045, 12042, 12050, 12043, 12065, 12042, 12064, 12042, 12051, 12042, 12060, 12031, 12041, 12036, 12049, 12044, 12051, 12046, 12053, 12037, 12053, 12050, 12046, 12054, 12051, 12040, 12043, 12043, 12052, 12062, 12065, 12053, 12038, 12063, 12055, 12050, 12054, 12044, 12067, 12060, 12045, 12056, 12052, 12046, 12034, 12062, 12047, 12043, 12056, 12051, 12044, 12035, 12028, 12058, 12038, 12044, 12043]
89: [ 11242, 11229, 11232, 11237, 11242, 11233, 11232, 11233, 11242, 11249, 11240, 11249, 11234, 11221, 11231, 11232, 11242, 11233, 11243, 11224, 11230, 11243, 11245, 11227, 11237, 11245, 11246, 11234, 11244, 11232, 11232, 11228, 11243, 11233, 11230, 11234, 11232, 11237, 11239, 11243, 11241, 11243, 11236, 11241, 11235, 11249, 11241, 11244, 11236, 11234, 11235, 11228, 11238, 11248, 11236, 11227, 11232, 11236, 11237, 11253, 11240, 11230, 11234, 11243, 11235, 11224, 11222, 11220, 11236, 11227, 11221, 11228, 11253, 11232, 11251, 11250, 11225, 11223, 11226, 11227, 11238, 11232, 11226, 11249, 11228, 11242, 11245, 11239, 11230]
97: [ 10301, 10303, 10321, 10322, 10298, 10307, 10301, 10301, 10318, 10305, 10315, 10319, 10316, 10311, 10302, 10307, 10306, 10311, 10333, 10319, 10293, 10300, 10293, 10315, 10304, 10307, 10307, 10297, 10315, 10317, 10302, 10317, 10316, 10314, 10321, 10309, 10323, 10310, 10305, 10309, 10304, 10309, 10313, 10310, 10309, 10310, 10326, 10310, 10310, 10314, 10301, 10297, 10310, 10310, 10304, 10307, 10313, 10315, 10307, 10316, 10321, 10312, 10321, 10313, 10320, 10312, 10307, 10315, 10310, 10317, 10304, 10301, 10306, 10319, 10301, 10291, 10293, 10319, 10303, 10301, 10307, 10309, 10303, 10302, 10314, 10291, 10314, 10324, 10314, 10312, 10318, 10310, 10299, 10305, 10298, 10304, 10305]

Соответственно на интервалах превышающих произведение проверяемых простых все комбинации остатков по всем проверяемым простым равновероятны. И если бы все вероятности разделиться на простое были одинаковы, то равновероятно на каком из них наберётся nf==nu, но их вероятность 1/p вместо одинаковой, потому меньшие простые выгоднее, они чаще делят n.
На меньших интервалах, а $\prod\limits_{p\in P}^{2^{20}}p\approx10^{454814}$ офигенно велико, распределение комбинаций остатков будет неким образом зависеть от проверяемых чисел, но будет своим для каждого интервала n (и длины и положения), плюс отличия будут малы (вот это стоит конечно проверить), плюс определить их будет не проще проверки всего интервала, КТО неплохо всё перемешивает. Т.е. зависимость будет сильнее не от паттерна (что само собой будет, но ничем не помогает в силу следующего), а от позиции и длины интервала по n внутри того офигенно огромного.

Проверка какое простое сколько раз отбросило кандидата:

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

...
frq=vectorsmall(2^20);
for(i=0,oo,
...
        x=n+#v-1;
        forprime(p=59,2^20,
                d=#v-x%p; if(d<1, next); nf[d]++; nn[d]*=p;
                if(nf[d]>nu[d], frq[p]++;next(2));
                if(x%p^2<#v, frq[p]++;next(2));
                if(nf[d]==nu[d] && !ispseudoprime((n+d-1)/v[d]/nn[d],1), frq[p]++;next(2));
        );
);
forprime(p=59,2^20, if(frq[p]>=100, print1(p,":",frq[p],",  "));); print;

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

59:52292,  61:51568,  67:47241,  71:44496,  73:43001,  79:39037,  83:36610,  89:33441,  97:30072,  101:28250,  103:27163,  107:25230,  109:24368,  113:22893,  127:19761,  131:18598,  137:17390,  139:16769,  149:15038,  151:14564,  157:13608,  163:12861,  167:12187,  173:11450,  179:10913,  181:10510,  191:9547,  193:9410,  197:9002,  199:8791,  211:8026,  223:7352,   227:7031,  229:7077,  233:6797,  239:6388,  241:6288,  251:5742,  257:5447,  263:5319,  269:5059,  271:4966,  277:4610,  281:4613,  283:4583,  293:4190,  307:4152,  311:3869,  313:3723, 317:3636,  331:3385,  337:3309,  347:3138,  349:3078,  353:2944,  359:2934,  367:2740,  373:2750,  379:2609,  383:2556,  389:2454,  397:2354,  401:2399,  409:2325,  419:2158,  421:2109,  431:2129,  433:1934,  439:2035,  443:1909,  449:1870,  457:1761,  461:1789,  463:1728,  467:1743,  479:1627,  487:1582,  491:1573,  499:1490,  503:1471,  509:1457,  521:1395,  523:1456,  541:1307,  547:1298,  557:1223,  563:1184,  569:1285,  571:1150,  577:1167,  587:1132,  593:1086,  599:1148,  601:1132,  607:1028,  613:997,  617:1038,  619:1005,  631:953,  641:948, 643:933, 647:919,  653:896,  659:838,  661:890,  673:832,  677:841,  683:790,  691:839,  701:780,  709:763,  719:712,  727:741,  733:694,  739:742,  743:719,  751:621,  757:622,  761:681,  769:679,  773:604,  787:598,  797:635,  809:653,  811:562,  821:576,  823:533,  827:545,  829:546,  839:583,  853:541,  857:486,  859:532,  863:516,  877:497,  881:516,  883:543,  887:472,  907:492,  911:500,  919:449,  929:494,  937:410,  941:431,  947:441,  953:415,  967:382,  971:420,  977:403,  983:369,  991:395,  997:395,  1009:354,  1013:381,  1019:359,  1021:395,  1031:417,  1033:364,  1039:364,  1049:338,  1051:405,  1061:325,  1063:346,  1069:353,  1087:333,  1091:362,  1093:304,  1097:329,  1103:298,  1109:305,  1117:314,  1123:301,  1129:304,  1151:298,  1153:284,  1163:297,  1171:277,  1181:283,  1187:273,  1193:296,  1201:265,  1213:265,  1217:262,  1223:264,  1229:268,  1231:271,  1237:236,  1249:241,  1259:239,  1277:235,  1279:251,  1283:223,  1289:246,  1291:246,  1297:231,  1301:234,  1303:236,  1307:247,  1319:246,  1321:218,  1327:214,  1361:243,  1367:232,  1373:195,  1381:209,  1399:208,  1409:206,  1423:192,  1427:195,  1429:203,  1433:195,  1439:210,  1447:189,  1451:191,  1453:187,  1459:183,  1471:187,  1481:168,  1483:168,  1487:180,  1489:176,  1493:187,  1499:187,  1511:171,  1523:162,  1531:148,  1543:177,  1549:156,  1553:156,  1559:136,  1567:139,  1571:160,  1579:177,  1583:166,  1597:144,  1601:123,  1607:150,  1609:155,  1613:169,  1619:146,  1621:164,  1627:141,  1637:155, 1657:146,  1663:154,  1667:133,  1669:150,  1693:130,  1697:124,  1699:150,  1709:121,  1721:125,  1723:119,  1733:121,  1741:151,  1747:141,  1753:125,  1759:127,  1777:119,  1783:131,  1787:133,  1789:113,  1801:128,  1811:131,  1823:130,  1831:127,  1847:109,  1861:129,  1867:111,  1871:110,  1873:112,  1879:120,  1901:109,  1907:104,  1913:107,  1931:100,  1949:119,  1973:114,  1987:112,  2003:116,  2011:113,  2027:103,  2039:111,  2087:108,

Меньше 100 раз там ещё тысячи простых, но они уже точно не интересны.

Ну и всё, список простых и так уже оптимален.

[wrest]

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

Ну и всё, список простых и так уже оптимален.

Да, видно хорошо, спасибо. То есть порядок следования простых для попыток факторизации от паттерна никак не зависит, они должны просто идти по возрастанию, и в этом смысле уже всё хорошо. Ясно. Профит отменяется.

Да кстати по поводу isdseudoprime, я как-то неясно там выразился и не то считал.
Ситуация тут такая. Хотя бы один из членов каждой цепочки проходит через проверку простоты. Вернее, не он сам, а известный к тому моменту его остаток (частное). Решение по отсеву 100% цепочек принимается в результате именно проверки простоты остатка (частного). Так выясняется и перебор и недобор множителей. Так вот, помимо миллиона проверок на простоту, который приводит к отсеву, около 100 тыс проверок к отсеву цепочки не приводят.
Но попытка попробовать переполнять числа малыми делителями (ждать пока станет nf>nu) приводит к увеличению времени на пару порядков. То есть, таких чисел в цепочках очень мало.

[Yadryara]

Ну вот стоило поспать 7 часов, а тут такие события. Я конечно буду постепенно вникать в то, что вы написали.

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