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

[wrest]

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

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

Ну вот к моему удивлению, если и быстрее то незаметно.

Но с точки зрения чтения текста, менее понятно. Так что я у себя вернул код типа
if(n%3^7>len && n%3^6<len , next);- тут n длинное и оно вычисляется до того
а с индексом цикла k у меня получался такой
if(k%3==0 && k%9!=3, next()); - тут всё короткое и n не вычисляется до этого оператора

-- 26.03.2026, 18:46 --

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

А вот для неразмещённых простых придётся остаток проверять по списку, а не просто <len, и как из него быстро получить правильный d для каждого простого не очень понятно.

Там это и не нужно по крайней мере пока. Для неразмещённых главное ускорение должно быть в сокращении использования ispseudoprime, всё остальное там быстрое.
Но как я писал выше, высших "плохих" степеней неразмещённых простых относительно немного, и по ним я от проверки простоты остатка (частного) уже избавился.

-- 26.03.2026, 19:05 --

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

Если бы Вы посмотрели на любой код VAL, то там так и сделано.

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

my(pat_prime_lim:small=nextprime(vecmax(factor(vecprod(v))[,1])+1));

Например, "терапевтика". С ней надо бы подразобраться. Например: а почему она такая и какая будет для других паттернов? Ну по тройкам она такая потому, что в одно из чисел цепочки встроено простое 3^5, и через два на третий там будет появляться 3^6, а поскольку 48 не делится на 6+1=7, можно откниуть практически треть цепочек "не глядя". Но если поглядеть, то из этой трети надо вернуть в оборот тоже треть которая будет давать допустимое 3^7 т.к. 48 делится на 7+1=8, и таким образом откинуть надо не 3/9, а 2/9. И это должно стать понятно программе, чтобы программно определить параметры "терапевтики".

Почему кубы встроенных простых выкидываем на втором шаге терапевтики?
Ну потому что куб даёт 3+1=4 делителя и из этого делается пока для меня загадочный вывод, что для успеха такого числа в нём должен засесть квадрат ещё одного простого, который даст 2+1=3 делителя, а вместе они дадут 3*4=12 делителей. Вот из структуры паттерна программа, зная что дальше начальные числа будут считаться как n=n0+k*m, должна бы вычислить, что именно надо подвергнуть "терапевтике".
Все предварительные вычисления, включая таблицу простых например до 2^24 это меньше секунды.

[Dmitriy40]

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

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

Задача не стояла найти все цепочки, задача стояла найти любую цепочку - и это позволяет пропускать сильно менее вероятные комбинации в пользу простоты кода и скорости перебора.
А вероятность обнаружить в разложении куб простого плюс ещё и квадрат простого вместо нескольких простых в первой степени (а квадрат первого простого мы уже и так сунули) существенно ниже, на порядки. И потому игнорируется.
Как и возможность получить 3^7, её проверка добавит лишь 1/9=11% цепочек в проверку, зато потребует лишнего деления (вычисления остатка). Плюс к ней ведь ещё понадобится простое в квадрате, что слишком маловероятно.
Помнится некоторое время назад собирали статистику вероятностей степеней, правда среди всех нечётных чисел.

[wrest]

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

Как и возможность получить 3^7, её проверка добавит лишь 1/9=11% цепочек в проверку, зато потребует лишнего деления (вычисления остатка).

Одно лишнее деление на цепочку это мелкие расходы. Лишний ispseudoprime -- средние, лишний numdiv - крупные.

-- 26.03.2026, 21:16 --

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

её проверка добавит лишь 1/9=11% цепочек в проверку,

Да, увеличивает нагрузку на табличную проверку с 27% до 38% от общего количества цепочек.

[Dmitriy40]

wrest
Не просто одно лишнее деление на цепочку, а лишнее деление на 11% цепочек ради выигрыша тысячных процента цепочек (и 3^7 и ещё p^2). ИМХО скорость перебора +10% важнее обнаружения на 0.01% больше подходящих цепочек.

[wrest]

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

ИМХО скорость перебора +10% важнее обнаружения на 0.01% больше подходящих цепочек.

Хорошо. Ну раз вы вероятности разные считали, то наверное прикинули и на каких местах наименее вероятно получение нужного количестства делителей? Хотя эти ваши ласточки и дротики с дырками в разных местах...
Надо наверное посчитать какие места выпали сколько раз из таблицы Это паттерн-специфичная ситуация, очевидно.

Посчитал. Ну всё банально и ровно. Первое число - количество делителей в паттерне, второе - относительная частота отсева этого места. Больше делителей в паттерне - чаще отсев
3;12
6;100
12;760

[DemISdx]

К сведению.

(Оффтоп)

Yadryara сообщил, по другим каналам связи, что у него форум перестал открываться.
От слова совсем.
Остальные сайты вроде работают, по крайней мере частично.

[wrest]

DemISdx в сообщении #1721072 писал(а):

Yadryara сообщил, по другим каналам связи, что у него форум перестал открываться.
От слова совсем

Подождём...

-- 29.03.2026, 18:56 --

В общем мне приходили и были отброшены такие идеи по улучшению табличного фильтра.
Проверять не осталось ли найти три делителя (numdiv остатка-частного) и проверять не куб ли остаток (частное). Отброшено по причине редкости кубов.
Если осталось найти 4 делителя (numdiv остатка частного должен быть 4, т.еи pq), то подождать, вдруг найдётся ещё делитель в этом числе или переполнится в соседнем. Это очень сильно увеличивает проверки по таблице (сейчас в среднем цепочка отбрасывается на 85-м простом при таббшице от 59 до 2^20).

Заполнять таблицу 20х10^6 Тут ещё можно побороться, но реализации через компиляцию с помощью gp2c плохо чистят память, в итоге память разрастается до гигабайта - неприемлемо, можно было бы вручную поправить Си-код по сборке мусора но это не наш метод. Порог наступаео именно на таблице из 10^6 цепочек, на 10^5 всё норм, а при 10^6 время увеличиваетс я в 15 раз (а должноьменьше чем в 10). Оценка по количеству операций собсно простая: сумма обратных чисел к простым в таблице простых от 59 до 2^20 примерно 1,21, анализируемых чисел ну скажем 20 млн. Итого 24 млн инкрементов элементов массива. Немного, но надо предвычислять инверсии и что-то ещё, в общем похоже что этот способ может быть немного быстрее, но требует проработки.

Похоже, кстати, причина того, что у Yadryara растёт стек примерно понятна: неэффективные процедуры сбора мусора при длинных и запутанных циклах большой вложенности. Надо конечно смотреть на Си-код, но я у вот натолкнулся на эту проблему, как говорится где не ждали.

В общем, пока радикальных идей нет.

[wrest]

Да, ещё мелькала такая мысль, что раз мы по таблице проверяем до 2^20 (или скажем до 2^24), то потом, если мы получим остаток (частное) меньше 2^40 (или 2^48), то он заведомо простой. Но тоже как-то не пригодилось. Даже не стальпроверять попадаются ли остатки (частные) меньшие 2^40...2^48

-- 30.03.2026, 00:14 --

О подсчёте всей таблицы разом (результат -- омега для каждого числа во всех цепочках).
Текст:

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

divt(k_len:small, lim_min:small, lim_max:small)={
    my(t0=getwalltime());\\Замер времени подготовки
    my(n0=283938699868309713921641266159023371777169):int;
    my(m=229403502054304075118105309394590566946400*2):int;
    my(len=21):small;
    my(pr:vecsmall=Vecsmall(primes([lim_min,lim_max])));
    my(pr_q=#pr):small;
    my(pmods_n0:vecsmall=vectorsmall(pr_q,i,n0%pr[i]));
    my(pmods_m:vecsmall=vectorsmall(pr_q,i,m%pr[i]));
    
    my(chunk=vector(k_len,i,vector(len,j,0)));
    my(inv_mods:vecsmall=vectorsmall(pr_q,i,if(pmods_m[i]==0, 0, (1/pmods_m[i])%pr[i])));
    print("Таблицы: модули простых для n0, m и обратные m : ",pr_q," каждая." );
    print("Таблица со счётчиками: ",len,"x",k_len," ~ 10^",logint(len*k_len,10)+1);
    print("Время подготовки: ",strtime(getwalltime()-t0));
    print("+++++ Старт вычислений +++++");
    my(t1=getwalltime());
    for(col:small=1,len,
        for(i:small=1,pr_q,
            \\ my(p:small=pr[i]);
            \\ my(r_n0:small=pmods_n0[i]);
            \\ my(r_m:small=pmods_m[i]);
            
            if(pmods_m[i] == 0,
                \\ Если m кратно p, то Бинго! проверяем только первый элемент
                if(((pmods_n0[i] + (col-1)) % pr[i]) == 0,
                    for(row=1,k_len,
                        chunk[row][col]++
                    )
                )
            ,
                \\ Иначе решаем: (r_n0 + (row-1)*r_m + (col-1)) % p == 0
                \\ (row-1)*r_m ≡ -r_n0 - (col-1) (mod p)
                my(first_row:small  = 1+ ((-pmods_n0[i] - (col-1)) % pr[i] * inv_mods[i]) % pr[i]);
                forstep(row:small=first_row, k_len, pr[i],
                    chunk[row][col]++
                )
            )
        )
    );
    print("Счёт закончен, скорость: ",1000*k_len\(getwalltime()-t1)," цеп/сек (только таблица)");
    print("Счёт закончен, скорость: ",1000*k_len\(getwalltime()-t0)," цеп/сек (включая подготовку)");
    
    return(chunk);
}

И вот такое странное поведение после компиляции через gp2c:

? divt(10^4,59,2^20);
Таблицы: модули простых для n0, m и обратные m : 82009 каждая. 
Таблица со счётчиками: 21x10000 ~ 10^6 
Время подготовки: 39 ms
+++++ Старт вычислений +++++ 
Счёт закончен, скорость: 333333 цеп/сек (только таблица)
Счёт закончен, скорость: 144927 цеп/сек (включая подготовку) time = 69 ms.
? divt(10^5,59,2^20);
Таблицы: модули простых для n0, m и обратные m : 82009 каждая. 
Таблица со счётчиками: 21x100000 ~ 10^7
Время подготовки: 34 ms
+++++ Старт вычислений +++++ 
Счёт закончен, скорость: 529100 цеп/сек (только таблица)
Счёт закончен, скорость: 448430 цеп/сек (включая подготовку)
time = 441 ms. 
? divt(10^6,59,2^20);
Таблицы: модули простых для n0, m и обратные m : 82009 каждая. 
Таблица со счётчиками: 21x1000000 ~ 10^8 
Время подготовки: 86 ms
+++++ Старт вычислений +++++
*** divt: Warning: increasing stack size to 800000000. 
Счёт закончен, скорость: 70397 цеп/сек (только таблица) 
Счёт закончен, скорость: 69974 цеп/сек (включая подготовку)
time = 12,607 ms.
?

Первый параметр - количество цепочек, следующие два параметра -- диапазон табличных простых.

Видно что 10 раз посчитать 10^5 цепочек быстрее чем один раз 10^6 цепочек, что конечно какой-то глюк в интерпретаторе/трансляторе gp2c

Но видно что с ростом количества 10^4->10^5 скорость растёт, что обнадёживает.

[wrest]

Да, ещё с матрицами или векторами векторов, есть такая засада в pari/gp
Они не могут быть короткими, увы. Ну или я не смог.

[Dmitriy40]

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

Видно что 10 раз посчитать 10^5 цепочек быстрее чем один раз 10^6 цепочек, что конечно какой-то глюк в интерпретаторе/трансляторе gp2c

Не обязательно: это типичное поведение при переполнении кэша, т.е. 10^5 цепочек ещё (почти) влезают в кэш, а 10^6 уже ну совсем не влезают.

В таких случаях рекомендуют поменять адресацию, переставить колонки и строки местами в chunk[][] (т.е. чтобы обращение было chuck[col][row]), тогда каждая итерация внешнего цикла по колонке будет обрабатывать лишь часть всего массива и будет выше локальность и может эта обрабатываемая часть снова станет помещаться в кэш. Хотя при таком размазанном доступе я не думаю что это чем-то поможет, кэши достаточно умные чтобы это нивелировать. Но опробовать можно.
И внутренний вектор по row тогда можно объявить vectorsmall, должен быть заметный выигрыш в размере.

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

[wrest]

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

И кстати матрица должна работать быстрее вектора векторов

Есть у меня подозрение, что матрица и вектор векторов - одно и то же во внутренней реализации. Во всяком случае типа T_Matsmall по аналогии T_Vecsmall нет. А тип T_Vec - он может содержать элементы T_INT и другие типы.
Но вы наверное правы насчёт поменять индексы местами, и попробовать внутренний вектор (который будет из 10^6 элементов) сделать vectorsmall -- надо попробовать.

[wrest]

Dmitriy40
Я попробовал:

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

Но вы наверное правы насчёт поменять индексы местами, и попробовать внутренний вектор (который будет из 10^6 элементов) сделать vectorsmall -- надо попробовать.

Не получается. В интерпретаторе работает, после компиляции не работает, хотя и компилируется. Причем сам вектор из vecsmall-ов создаётся. Но дальше изменить конкретный элемент типа my_v[3][5]=8 нельзя - segmentation fault.

[Dmitriy40]

wrest
Поменять индексы можно и независимо от small, это два разных лайфхака.

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

Есть у меня подозрение, что матрица и вектор векторов - одно и то же во внутренней реализации.

Разумеется нет:

? vector(100,i,vector(100,j,0));getstack()
%1 = 82400
? matrix(100,100);getstack()
%2 = 81984

[wrest]

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

Разумеется нет:

М-да, ну я ориентировался на гайд по типам в gp2c https://pari.math.u-bordeaux.fr/pub/par ... html#sec11 где написано что тип vec в gp2c это "vector and matrices (excluding vecsmall)"

Но вот интересно, если chunk это вектор векторов, то вызов
chunk[3][5]=111
транслируется в Си как два вызова gel (get element, я думаю)
gel(gel(chunk, 3), 5) = stoi(111);
А если chunk_mat это матрица, то
chunk_mat[3,5]=222 транслируется в Си как один вызов gcoef (наверное get coeff)
gcoeff(chunk_mat_h, 3, 5) = stoi(222);

[Dmitriy40]

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

М-да, ну я ориентировался на гайд по типам в gp2c

Значит есть разница в типах между PARI/GP и gp2c и тогда говорите точнее про что имеете в виду.

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

Но вот интересно,

Что же тут интересного, всё как и должно быть, или два вызова преобразования адреса, или один.