Плюс сейчас как-то немного за кадром у меня остались предфильтры которые есть у вас в эталоне,
if((n+#v-1)%3^6==21-19, next);\\Такие n не подходят
я там не понял что к чему, а их можно дополнить и применять не только к началу цепочки а и ко всем членам?
Я их у себя объединил в два:
if((n+#v-1)%3^6<#v, next);
forprime(p=5,53, if((n+#v-1)%p^3<#v, next(2)); );А делают они следующее: вот есть в паттерне v[] число 3^5, если при каком-то n на него попадёт ещё одна тройка, то число станет 3^6, что даст множитель в numdiv равный 7, но 48 не делится на 7 и потому такие n не подходят. Вероятность попасть сразу двум тройкам (т.е. 3^2) в место с 3^5 чтобы дать 3^7 и допустимый множитель 8 в numdiv пренебрежимо мала.
Но тройка (точнее 3^6) не может попасть ни на какие позиции v[] кроме как с 3^5, значит проверив просто что 3^6 в паттерн не попала мы полностью исключаем её и все высшие степени 3 в любой позиции. И их проверять уже не надо.
Аналогично и для простых 5...53, только они есть лишь в квадратах, потому проверяем не превратились ли они в кубы (и более старшие степени).
Где-то куб может оказаться и допустимым, но тогда там не хватит квадрата, а квадрат простого намного реже первой степени и на такие n тоже плюём ради скорости (задача поиска строго наименьшего решения не ставилась).
На самом деле эти длинные деления можно сделать короткими пересчитав
((n=n0+m*2*i)+#v-1)%p<#v в
i%p<>C[p] (ведь остальное известные константы) где C[p] будет зависеть от места квадрата данного простого p в v[] (как и было в эталоне, видите эти 21-19, это оно, позиция в v[]), но мне лень, и на дальнейшую факторизацию это не влияет. В рабочей программе это сделано. Я даже проверил, на скорость замена в пределах погрешности не влияет.
-- 07.12.2025, 03:20 --Чтобы ещё меньше делить, например. Ну может предрассчитать остатки сразу скажем для начала 10 следующих цепочек.
Чтобы меньше делить надо полностью перейти к арифметике по модулю, тогда останутся только сложения/вычитания вместо умножений и делений. Но придётся складывать по модулю два вектора. Я же писал про это вчера днём (5 декабря):
Другое дело что да, можно заменить деление на сложение по модулю, но тогда вместо одного 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, но это вероятно медленнее).
...
И эта операция прекрасно ложится в SSE/AVX2. Например в AVX2 она реализуется тремя командами над регистром длиной 256 битов или 8 штук int32 и все три команды выполняются в среднем за полтора такта, плюс две команды накопления флага равенства нулю любого из остатков за такт, т.е. скорость такого суммирования длинных векторов составляет примерно 3.3 элементов int32 за такт. Или 6.6 элементов int16 за такт. Деление 192 битового числа (до 57 знаков) занимает 100-300 тактов в самом лучшем случае. Разница на три порядка.
Сможет ли gcc выдать столь же быстрый код - очень сомневаюсь, но это не мои проблемы, у меня есть ассемблер.
И ранее:
Умножение тоже не нужно: переход от предыдущего 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 это массивы. И я об этом уже писал где-то выше. Именно так работают мои программе на асме.
Да, проблема в том что их нельзя взять и оборвать досрочно если где-то остаток обнулился (т.е. разделилось на простое), надо просуммировать до конца - для следующих n. Зато весь (ну или начальная часть) цикл forprime с длинным делением внутри заменяется на сложение векторов, которое относительно быстрое. Как его превратить в быстрое сложение по модулю это вопрос для исследования. Можно как показал с циклом и if, можно типа
n=[t=n[j]+s[j];if(t<p[j], t, t-p[j]) | j<-[1..#n]], можно объявить их (n[], s[]) с Mod(x,p) и тогда простое сложение прокатит, но вот проверка на наличие нуля вместо простого vecmin(n)==0 превратится снова в цикл, может ещё что придумать ...
-- 07.12.2025, 03:27 --Если вдруг gcc в убунте поддерживает встроенный ассемблер или интриксы для вставки AVX команд в С код, то можно прямо на C поправить цикл для быстрого суммирования по модулю векторов. А их инициализацию и дальнейшие действия оставить на PARI. Но с этим надо разбираться. А у меня ни убунты ни gcc нету.
07.12.2025, 02:15 
