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

[Dmitriy40]

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

А чем сложение по модулю отличается от короткого деления?

Скоростью.
Например код
d+=b; if(d>=p, d-=p);
должен транслироваться в (x64, справа время выполнения в потоке)

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

mov RAX,d       ;0.5
add RAX,b       ;0.5
mov RDX,RAX     ;0.25
sub RAX,p       ;0.5
cmovnc RDX,RAX  ;0.5

с суммарным временем выполнения 0.5+0.5+0.25+0.5+0.5=2.25 тактов.
Даже если компилятор не любит cmov, то вместо неё поставит две команды jc+mov, которые выполняются за (0.5+0.25)...2=0.75...2 тактов и суммарно тогда 2.5...3.75.
Я взял худший случай, на самом деле c cmov будет 0.25+0.25+0+0.25+0.5=1.25.

А код с коротким делением
d=(d+b)%p;
наверняка странcлируется в код

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

mov RAX,d       ;0.5
add RAX,b       ;0.5
xor RDX,RDX     ;0.25
div p           ;21...74

с суммарным временем выполнения 22.25...75.25 тактов.

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

Вот делить-то можно например вычитанием,

В общем случае это ещё дольше аппаратного деления потому что цикл (возможно развёрнутый, не суть) из 64 сравнений/вычитаний. Даже оптимизация для пропуска старших нулевых битов не сильно поможет (и очень сомневаюсь что компилятор до неё додумается). В частном случае - будет код как первый выше без деления, но до такой оптимизации ни один компилятор не додумается, у него просто нет информации что d и b в этом месте кода оба всегда строго меньше p. Хотя, если это ему прямо указать ... Не знаю, всё равно не думаю.

-- 09.12.2025, 01:53 --

(Векторизация для AVX2)

А так как эти операции выполняются в цикле, то компилятор в принципе может их векторизировать для AVX2 до:

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

vmovdqu ymm1,[ p ]      ;0.5
vmovdqu ymm0,[ d ]      ;0.5
vpaddd  ymm0,ymm0,[ b ] ;0.5
vpcmpgtd ymm2,ymm1,ymm0 ;0.5
vpandn  ymm1,ymm2,ymm1  ;0.33
vpsubd  ymm0,ymm0,ymm1  ;0.5

и за 2.83 такта получит 8 результатов или по 0.35 такта на число.

Я же применяю финт ушами и ограничиваю диапазон до на бит меньше, до 7/15/31 битов, и использую такой код:

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

vmovdqu ymm0,[ d ]      ;0.5
vpaddd  ymm0,ymm0,[ b ] ;0.5
vpsubd  ymm1,ymm0,[ p ] ;0.5
vpminud ymm0,ymm0,ymm1  ;0.5

и трачу всего 2 такта на 8 результатов или 0.25 такта на число. Или 0.125 такта на число если хватает 15 битов на каждое.

-- 09.12.2025, 02:13 --

wrest
Выше намекал что ценой памяти можно реже выполнять сложение всего вектора остатков по модулю. Идея вот в чём: делим цикл сложения и проверки на два отдельных, только проверки и только сложения, тогда цикл проверки память не меняет и его можно прервать в любой момент, по обнаружению причины отбросить n. А цикл сложения тогда не содержит условий и выполняется быстрее.
Так это выигрыша не даёт, но тут появляется основная идея: цикл проверки разворачиваем в память раз скажем 1000, т.е. храним в памяти большую таблицу размером например 1900х1000 элементов (для простых до 2^14) в которой указаны добавки к базовому значению остатка по модулю каждого простого для всех 1000 итераций (т.е. tab[p,i]=(m*i)%p). Тогда цикл сложения выполняем лишь 1/1000 итерацию, складываем сразу n+=(m*1000)%p, это тоже предвычисленная константа. А в цикле проверки выполняем операцию (n[p]+tab[p,i])%p<#v, где сложение по модулю делается конечно без деления.
Выигрыш в том что сложения всего вектора делаются в 1000 раз реже, а цикл проверки практически для каждого i обрывается задолго до конца вектора. Так что сумма получается не 1.001 сложений по модулю на каждый элемент как было бы в худшем случае, а сильно меньше, типа 0.101. Цена - таблица в памяти и скорость доступа к ней (потому желательно чтобы она помещалась хотя бы в кэш L3 так как она общая для всех потоков).
На PARI это выигрыша не даёт, насколько помню.

[wrest]

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

но тут появляется основная идея: цикл проверки разворачиваем в память раз скажем 1000, т.е. храним в памяти большую таблицу размером например 1900х1000 элементов (для простых до 2^14) в которой указаны добавки к базовому значению остатка по модулю каждого простого для всех 1000 итераций (т.е. tab[p,i]=(m*i)%p). Тогда цикл сложения выполняем лишь 1/1000 итерацию, складываем сразу n+=(m*1000)%p, это тоже предвычисленная константа.

Идея-то понятна - вместо деления с остатком делать таблицу умножения и по ней смотреть. Но заранее мы не знаем вся ли таблица понадобится, соответственно, надо смотреть на отношение деления с остатком

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

с суммарным временем выполнения 22.25...75.25 тактов.

к сложению и вычитанию

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

с суммарным временем выполнения 0.5+0.5+0.25+0.5+0.5=2.25 тактов.

Т.е. выгодно если на одно вычисление делением не больше 35 вычислен й сложением-сраунением-вычитанием.

[Dmitriy40]

Ну как бы да.
Но смотрите, на 100 тысяч kpod (т.е. n) выполняется 7.82млн делений для вычисления d, т.е .в среднем по 78 делений на каждое n (из которых реально разделились лишь 4.5 на n). Значит их всё ещё выгодно заменять на 78*35=2730 сложений по модулю. А это простые почти до 25000.
Но что реально накомпилит gcc никому неизвестно.

-- 09.12.2025, 16:59 --

Интерпретатор проигрывает уже при размере массива простых до 2^9, хотя при 2^8 ещё вровень. Компиляцию проверить не могу.

[Yadryara]

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

Вот, кстати, вопрос.

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]));

Почему здесь уважаемый wrest объявил смолом вектор, где куча чисел явно больше 255 ?

И вот всё равно при переходе к счёту по одному из своих паттернов получилось 4-кратное замедление. Пока не разобрался в чём дело:

fun ( potok:small, starti:small, koli:small, maxpredp:small ) = {

my(t0=getwalltime());

my(zama:vecsmall=Vecsmall([13]));
\\my(zamb:vecsmall=Vecsmall([61]));
my(kdel = 48, maxlen = 5):small;
my(bolv=[13,2,3,28,5,18,1,32,3,50,7,12,1,338,45,88,1,6,17,20,3,2]);
my(mkv:vecsmall=Vecsmall(Set(select(x->(numdiv(x)%3<>0),bolv,1))));
my(lenmkv = #mkv):small;
my(rkp:vecsmall = vectorsmall(lenmkv,x,0));
rkp[1] = nextprime(maxlen+1); for(j=2, lenmkv,
if(nextprime(rkp[j-1]+1)<>zama[1], rkp[j] = nextprime(rkp[j-1]+1),
rkp[j] = nextprime(nextprime(rkp[j-1]+1)+1) ));
\\rkp = setminus(rkp,zam1); rkp = concat(rkp,zam2);
my(rassprob:vecsmall=Vecsmall([1..lenmkv]));

my(v:vecsmall=Vecsmall([3757,242,1587,28,4205,
18,1369,32,5547,50,11767,12,40931,338,45,7688,361,4998,2809,20,6627,6962]));
my(ip = 17):small;
my(aa = 4998):small;
my(m=132115864604174642234399073835381356548570400):int;
my(n0=68182721761635739290794149646849536294407271):int;
my(maa=132115864604174642234399073835381356548570400*aa):int;
my(n0aa=68182721761635739290794149646849536294407271*aa):int;

my(nu:vecsmall=Vecsmall([3, 3, 3, 3, 3, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 4, 1, 4, 3, 3, 3]));

my(lenv=#v):small;

my(predp:vecsmall = Vecsmall(primes( [67, maxpredp] )));
my(kolpredp = #predp):small; 

my(ostmaa:vecsmall  = vectorsmall(kolpredp, nompr, maa  % predp[nompr] ));
my(ostn0aa:vecsmall = vectorsmall(kolpredp, nompr, n0aa % predp[nompr] ));

my(n):int;
my(mes):small;
my(kolfac:vecsmall=vectorsmall(lenv,mes,1)); 
my(pro:vec=vector(lenv,mes,1));
my(q=q1=best=valids=0):small;

print(starti,"   ", koli,"   ", maxpredp, "    ",strtime(getwalltime()-t0), "        ");

t0=getwalltime();

for(i:small=starti, starti + koli - 1, 


kanp = n0 + m * i;

\\if(!ispseudoprime(kanp,1), next);

n = aa * kanp - ip;

\\n = aa * (n0 + m * i) - ip;


kolfac=vectorsmall(lenv,mes,1); pro=vector(lenv,mes,1);

for(nompr=1, kolpredp,

\\mes = lenv - (ostn0[nompr]+2*i*ostm[nompr]+lenv-1) % predp[nompr];
\\mes = lenv - (aa * (ostn0[nompr] + ostm[nompr]*i) - ip + lenv - 1) % predp[nompr];

mes = lenv - (ostn0aa[nompr] + ostmaa[nompr]*i - ip + lenv - 1) % predp[nompr];

if(mes<1, next); kolfac[mes]++; pro[mes] *= predp[nompr]; 

if( kolfac[mes] == nu[mes], q1++;
if(!ispseudoprime((n+mes-1)/(v[mes]*pro[mes]),1), next(2))));

for(mes=1,lenv,
if(kolfac[mes]<nu[mes],q1++;
if(ispseudoprime((n+mes-1)/(v[mes]*pro[mes]),1),next(2))));

q++; 
print;


valids=0; for(mes=1,lenv,if(numdiv(n+mes-1)==48,valids++));

if(valids>=lenv-3,best++;print(n, "   ", valids);print));

print("Potok ", potok, "      ",best,"        ", q,"      ",q1,"    ", 1000*koli\(getwalltime()-t0),"    ",strtime(getwalltime()-t0));
return(best);
}


/*
GP;default(threadsize,32M);
GP;default(debugmem,0);
GP;init_RabFun_66_7();
*/

{t0=getwalltime();print;

\\print1("      ");

parfor(i=0,5,fun(i,i*10^5,10^5,2^20));


print;
print("                                                ",
strtime(getwalltime()-t0));print;
}quit

Кстати, специально так заточено, чтоб он находил известное приближение с valids=19.

[wrest]

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

Почему здесь уважаемый wrest объявил смолом вектор, где куча чисел явно больше 255 ?

Потому что small это всё, что меньше 2^64 на 64-разрядных версиях pari.
"Маленькие" числа это такие, которые помещаются в регистр процессора - тогда вычисления с ними максимально быстрые.
В Си- нотации это long (на 64битных системах соответственно 64 бита) или ulongc Там есть некоторая путаница, но для вас безопасно считать что pari-шный small это от -2^63 до 2^63

-- 10.12.2025, 10:45 --

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

Пока не разобрался в чём дело:

Добавьте пож-ста комментарии и отступы.

[Yadryara]

Отступы ещё больше не люблю, потому что код норовит за правый край вылезти. У меня ведь шрифт крупный. Но это всё равно периодически происходит. Недавно опять было. Я не заметил и возникла ошибка. Это обсуждалось на форуме.

Как видите, я строки добавляю для лучшей разбивки кода. Так что вертикальные отступы ставлю.

А горизонтальные для чего нужны? Для понимания иерархии? Я и так пока что понимаю, что куда вложено. Здесь пока сохраняется то, что было у вас. Чуть по-другому считается $n$. На скорость ещё может влиять, что паттерн v длиной 22. Другие числа больше, но для того ведь и забрались на арифмост, чтобы это не играло.

Ставить горизонтальные отступы не вижу смыла, извините.

Коменты дам. Скорее всего, они нужны буквально для 3-5 строк. Что именно непонятно?

[Yadryara]

Зато, кстати, есть очень маленькое падение по скорости. Один поток:

7838280910875426566029915787276966511735589626244441   19

Potok 0   100000   1048576    9 ms

Potok 0      1        1      6510    41511    2,409 ms

                                              2,431 ms

Два потока:

7838280910875426566029915787276966511735589626244441   19

Potok 0   100000   1048576    12 ms
Potok 1   100000   1048576    10 ms

Potok 1      0        0      6325    120336     831 ms
Potok 0      1        1      6510    41528    2,408 ms

                                              2,433 ms

6 потоков:

Potok 4   100000   1048576    32 ms
Potok 0   100000   1048576    34 ms
Potok 1   100000   1048576    33 ms
Potok 2   100000   1048576    33 ms
Potok 3   100000   1048576    32 ms
Potok 5   100000   1048576    32 ms

Potok 5      0        0      6364    113378     882 ms
Potok 2      0        0      6394    112107     892 ms
Potok 4      0        0      6299    112107     892 ms
Potok 1      0        0      6325    109529     913 ms
Potok 3      0        0      6257    101214     988 ms

7838280910875426566029915787276966511735589626244441   19

Potok 0      1        1      6510    39572    2,527 ms

                                              2,579 ms

12 потоков:

Potok 0   100000   1048576    69 ms
Potok 2   100000   1048576    70 ms
Potok 1   100000   1048576    71 ms
Potok 3   100000   1048576    72 ms
Potok 10  100000   1048576    69 ms
Potok 6   100000   1048576    69 ms
Potok 4   100000   1048576    72 ms
Potok 5   100000   1048576    71 ms
Potok 9   100000   1048576    70 ms
Potok 8   100000   1048576    70 ms
Potok 7   100000   1048576    71 ms
Potok 11  100000   1048576    71 ms

Potok 8      0        0      6174    74128    1,349 ms
Potok 1      0        0      6325    74019    1,351 ms
Potok 3      0        0      6257    73637    1,358 ms
Potok 4      0        0      6299    73152    1,367 ms
Potok 7      0        0      6240    73099    1,368 ms
Potok 9      0        0      6281    72992    1,370 ms
Potok 11     0        0      6156    73046    1,369 ms
Potok 2      0        0      6394    72727    1,375 ms
Potok 5      0        0      6364    72886    1,372 ms
Potok 6      0        0      6189    72202    1,385 ms
Potok 10     0        0      6300    71994    1,389 ms

7838280910875426566029915787276966511735589626244441   19

Potok 0      1        1      6510    34916    2,864 ms

                                              2,957 ms

[wrest]

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

Я и так пока что понимаю, что куда вложено.

Ок, тогда и разбирайтесь в вашем коде сами

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

А горизонтальные для чего нужны?

Для того что такой код:

rkp[1] = nextprime(maxlen+1);
for(j = 2, lenmkv,
    if(nextprime(rkp[j-1]+1)<>zama[1],
        rkp[j] = nextprime(rkp[j-1]+1)
        ,
        rkp[j] = nextprime(nextprime(rkp[j-1]+1)+1)
      )
   );
\\rkp = setminus(rkp,zam1); rkp = concat(rkp,zam2);
my(rassprob: vecsmall = Vecsmall([1..lenmkv]));

понятней и легче читать чем такой:

if(nextprime(rkp[j-1]+1)<>zama[1], rkp[j] = nextprime(rkp[j-1]+1),
rkp[j] = nextprime(nextprime(rkp[j-1]+1)+1) ));
\\rkp = setminus(rkp,zam1); rkp = concat(rkp,zam2);
my(rassprob:vecsmall=Vecsmall([1..lenmkv]));

-- 10.12.2025, 13:40 --

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

Коменты дам. Скорее всего, они нужны буквально для 3-5 строк. Что именно непонятно?

Непонятно для чего переменные. Вам имена rkp, mkv, n0aa и прочее о чём-то говорят, мне - нет.
Непонятен и ход вычисления - что для чего и почему делается.
Так что если хотите, чтобы ваш код был понятен (и можно было в него вникнуть) -- потрудитесь сделать его таковым. Ну а на нет и суда нет.

[Yadryara]

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

Для того что такой код:
понятней и легче читать чем такой:

Согласен. Но только эти 4 строчки можно пока игнорить. Они пока не используются. Никак.

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

Непонятно для чего переменные. Вам имена rkp, mkv, n0aa и прочее о чём-то говорят, мне - нет.

Я уже писал, но мне не трудно повторить:

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

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

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

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

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

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

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

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

n0aa это вспомогательная переменная: n0aa = n0*aa

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

Непонятен и ход вычисления - что для чего и почему делается.

С какого места непонятно? Делается почти ровно то же самое что у вас в функции em3. Только m и n по-другому вычисляются.

[wrest]

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

Я уже писал, но мне не трудно повторить:

Желательно, чтобы это было комментариями прямо в коде.

-- 10.12.2025, 16:58 --

Yadryara
Я (и надеюсь впоследствии и вы тоже) смотрю ещё и на сгенерированный Си-шный код.
Комментарии gp2c переносит из pari-шного в Си-шный код, так что там можно понять о каком участке pari-шного кода идёт речь. Это ещё один плюс комментирования в этом конкретном случае.

[Yadryara]

Согласен про коменты. Но сейчас важнее другое.

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

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

Нужно 12 разных паттернов считать в 12 потоках. Собственно, в работающей сейчас программе так и считается. Только потоков сейчас 9.

Имеется 162 комплекта по 10! паттернов в каждом. 10! и на 9 отлично делится, и на 12, и на 24.

[wrest]

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

Согласен про коменты. Но сейчас важнее другое.

Ок, хорошо, тогда вернёмся как только более важную проблему решите.

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

Нужно 12 разных паттернов считать в 12 потоках. Собственно, в работающей сейчас программе так и считается. Только потоков сейчас 9.

Я этого не понимаю. Задача была сформулирована (мной) вот тут: post1711508.html#p1711508
Да, сейчас в (моей) функции em3() специфичные для "паттерна" константы заданы явно в коде функции, так что она подходит только для конкретной пары (n0;m) , и в коде задан конкретный шаг 2*m для прогрессии, и конкретная длина цепочки для проверки, равная 21.
Если это не подходит и условия меняются - ну значит надо передавать в функцию паттерно-специфичные константы через параметры.
При этом, для эффективной параллельной работы надо, чтобы функция работала заметное время (единицы-десятки секунд), а не вызывалась десятки-тысячи или более раз в секунду.

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

Я чувствую, что вы внутренне против этого подхода, так же как и против отступов, комментариев и т.п.
Когда созреете, можем поговорить о правильном подходе, буду ждать.

[Yadryara]

Здесь довольно сложная ситуация.

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

длина цепочки для проверки, равная 21.

Которая в настоящий момент давно уже неактуальна. Ищут D(48,22) и D(48,23). В дальнейшем, как и в моём коде выше, речь о длине 22.

Есть серия паттернов 1-1-17-3-10! в которой идёт счёт. В ней 4 группы комплектов паттернов. Итого 162 комплекта по 10! паттернов, то есть всего более полумиллиарда паттернов только в этой серии.

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

Да, сейчас в (моей) функции em3() специфичные для "паттерна" константы заданы явно в коде функции, так что она подходит только для конкретной пары (n0;m)

То бишь для конкретного паттерна. Потому что при смене паттерна они перевычисляются.

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

в коде задан конкретный шаг 2*m для прогрессии,

Вот шаг как раз другой, выше это чётко показано. Рядом с прежним шагом — новый.

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

Если это не подходит и условия меняются - ну значит надо передавать в функцию паттерно-специфичные константы через параметры.

Конечно. Но я пока даже этого не сделал. Потому и считал что пока вам будет легко разобраться в коде, раз он не сильно изменён. И на входе по сути те же 4 параметра.

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

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

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

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

То есть надо спроектировать какой-то универсальный интерфейс (API) - что передавать в функцию.

Не против.

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

Я чувствую, что вы внутренне против этого подхода, так же как и против отступов, комментариев и т.п.

Неправильно чувствуете. Не надо всё в одну кучу валить. Быстроту и удобство программы и совершенно необязательные коменты и/или горизонтальные отступы.

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

Когда созреете, можем поговорить о правильном подходе, буду ждать.

Смотря какой подход считаете правильным. Собственный? Выше я с вами в главном согласился.

[wrest]

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

Смотря какой подход считаете правильным. Собственный?

Ессно. И я про подход к обмену кодом и как его писать, в рамках локальных задач ускорения, про моё удобство. Чтобы я участвовал, мне должно быть удобно. Если мне неудобно/непонятно и т.п. -- я пас :) Но я и не требую многого, вроде. Более того, это что-то навроде правил "хорошего тона" (как например старательная письменная речь без проблем с орфографией и пунктуацией), которые вы можете подвергать сомнению, но тогда кооперация будет затруднена.

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

[Dedekind]

Yadryara, wrest, переходите на Python. Там и проблема с отступами автоматически решится:)