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

[wrest]

Yadryara
Чтобы вектор был vectorsmall надо чтобы элементы были small.
Проверьте что здесь
vectorsmall( #predp, i, m % predp[i] )
m%predp[i] - имеет тип small, если нет то не получится.
У меня для того чтобы работало, вектор с простыми задекларирован как vectorsmall:
my(pr:vecsmall=Vecsmall(primes([59,lim]))); \\ вектор с простыми
а m задекларировано как длинное целое:
my(m=229403502054304075118105309394590566946400):int;
Поскольку остаток от деления длинного на короткое числа - короткий, то транслятор подставляет туда функцию smodis которая возвращает именно короткое число.

884 smodis(GEN x, long y)
885 {
886   pari_sp av = avma;
887   long r; (void)divis_rem(x,y, &r);
888   return gc_long(av, (r >= 0)? r: labs(y) + r);
889 }

А divis_rem в итоге вызывает функцию деления с остатком длинноготчисла на короткое в библиотеке gmp

Там всё ускорение построено именно на том, что 64-битные числа 64-битный процессор вычитает, складывает и т.п. быстро, а для больших чисел всё делается "руками". Вот представьте что у вас есть калькулятор (физический прибор). Вы на нём можете быстро посчитать произведение чисел которые по длине в сумме меньше чем ёмкость экрана. Ну допустим, экран там на 10 цифр. А вам надо умножить два 10-значных числа и получить точный ответ. Ну что тут делать: берёте листок бумаги и умножаете "в столбик", используя ваш аппаратный калькулятор для промежуточных (коротких) вычислений. Вот тут скорость и падает

[Yadryara]

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

Там всё ускорение построено именно на том, что 64-битные числа 64-битный процессор вычитает, складывает и т.п. быстро,

Да, я это понимаю. И я кстати, эти строчки уже преодолел, вспомнив про forperm и сделав вот так:

ostm:vecsmall = vectorsmall( #predp, nompr, m % Vec(predp)[nompr] );

[wrest]

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

ostm:vecsmall = vectorsmall( #predp, nompr, m % Vec(predp)[nompr] );

Это нехорошо, т.к. транслятор будет думать, что справа от знака оператора взятия остатка от деления % длинное число (вот это: Vec(predp)[nompr]).
Можете заглянуть в Си код, чтобы убедиться.

[Yadryara]

В Си код заглядывать не хочется. Я же ведь её ранее объявил:

maxprim = 2^10-2; predp = Vecsmall(primes( [67, maxprim] ));

[wrest]

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

В Си код заглядывать не хочется. Я же ведь её ранее объявил:

Хорошо, то есть всё работает уже, как надо?

[Yadryara]

Нет пока не всё работает. В других местах та же ошибка возникает.

Например, вот здесь я правильно объявил?

aa = v[ip+1]; m:int = 6*lcm(v)/aa;

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

m задекларировано как длинное целое:
my(m=229403502054304075118105309394590566946400):int;

Здесь что-то не догоняю. Число ведь намного больше чем $2^{64}$

[Dmitriy40]

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

Этот синтаксический сахар хорошо выглядит и даже имеет ускорительный смысл когда вы программируете интерпретатор.

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

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

Ну AVX это вы уже немного через край хватанули, имхо. У меня вот например под рукой две архитектуры aarch64 (она же arm64) и amd64 (она же x86-64), при том для развлечений с pari основная aarch64

Зато и у меня, и у Yadryara и у Демиса AVX2 вполне доступен. Вы же к нашему счёту помнится не присоединялись ...

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

А pcoul этот кстати как без gmp оходится - там длинная арифметика своя?

Да вроде и не обходится, использует.

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

Вот именно, а это 80 тысяч на каждую цепочку.
При том что в среднем нам нужны сотня-другая остатков.

Ну так переводим на остатки лишь пару сотен (или пару тысяч, до 2^14) первых простых, если цепочка не отброшена, запускаем обычный forprime. 80тыс конечно многовато без векторизации (SSE/AVX). Хотя способ решения как считать побыстрее тоже есть (и в программе на асме у меня использован), но требует много памяти (глобальной, read only, на асме у меня общей для всех потоков).

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

Я же ведь её ранее объявил:

А потом конструкцией Vec() превратили обратно в длинное число.

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

ostm:vecsmall = vectorsmall( #predp, nompr, m % Vec(predp)[nompr] );

И кстати это страшно медленно должно быть - Vec(predp) будет перевычисляться для каждого nompr, а не один раз и потом лишь индексироваться. В vector() инициирующее выражение вычисляется для каждого индекса заново целиком. Сможет ли gcc это оптимизировать и вычислить один раз - неизвестно.

-- 07.12.2025, 16:19 --

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

Здесь что-то не догоняю. Число ведь намного больше чем $2^{64}$

Вероятно в PARI int это указание на целое число произвольной длины, а не int64.

[wrest]

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

Здесь что-то не догоняю. Число ведь намного больше чем $2^{64}$

Да, поэтому и int. О типах в gp2c см. https://pari.math.u-bordeaux.fr/pub/par ... .html#sec4

Цитата:

int
Multi-precision integers represented by t_INT GENs.

-- 07.12.2025, 17:01 --

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

Зато и у меня, и у Yadryara и у Демиса AVX2 вполне доступен. Вы же к нашему счёту помнится не присоединялись ...

Да, ну я в этих ваших AVX не понимаю -- это ваше поле :)
Ассемблерные вставки для архитектуры amd64 (она же x86-64) определены в файле (даю путь с windows-слэшами)
...\pari\src\kernel\x86_64\asm0.h
А для x86 (32 бит) тут:
...\pari\src\kernel\ix86\asm0.h
Но как я понимаю, за арифметикой pari идёт в библиотеку gmp,
https://my.eng.utah.edu/~pajensen/ACM/D ... index.html
где есть и длинные и смешанные операции (это моё предположение). Вдруг там уже есть использование SSE/AVX и т.п. для чего-нибудь? Я не знаю...
Соответственно, если вы захотите остаться в pari, но при этом максимально использовать особенности вашей архитектуры, то вам надо будет разбираться как там это всё устроено, и собирать pari под себя, с вашим ассемблерным кодом для некоторых операций под архитектуру amd64

Ну например. Вот как сделано деление с остатком длинного на короткое в pari:
Последние две буквы в divis это типы аргументов i - длинное, s - короткое (32 или 64 бита в зависимости от архитектуры).

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

GEN
divis_rem(GEN y, long x, long *rem)
{
  long sy=signe(y),ly,s;
  GEN z;

  if (!x) pari_err_INV("divis_rem",gen_0);
  if (!sy) { *rem = 0; return gen_0; }
  if (x<0) { s = -sy; x = -x; } else s = sy;

  ly = lgefint(y);
  if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }

  z = cgeti(ly);
  *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
  if (sy<0) *rem = -  *rem;
  if (z[ly - 1] == 0) ly--;
  z[1] = evallgefint(ly) | evalsigne(s);
  return z;
}

mpn_divrem_1 - это функция из библиотеки GMP, вот она для x86-64, https://github.com/WinBuilds/gmplib/blo ... vrem_1.asm
Вот она же под SSE2
https://github.com/WinBuilds/gmplib/blo ... vrem_1.asm
Вот она же под core2
https://github.com/WinBuilds/gmplib/blo ... vrem_1.asm
Ну когда всё собирается в кучу, то собирается под ваш конкретный проц.

-- 07.12.2025, 17:13 --

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

Да вроде и не обходится, использует.

Тогда pcoul видимо тоже обёрнут в mingw/cygwin или какой-то подобный эмулятор, как и pari. Ну что же, все дороги ведут в линукс

[wrest]

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

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

Мне кажется так оно не работает, надо рабоче-крестьянским циклом менять каждый элемент последовательно.
В вашем (нашем) коде с модульной арифметикой например есть инициализация аккумуляторов для накопления омеги и остатка (частного) по числам цепочки:
nf=vectorsmall(len,d,1); nn=vector(len,d,1); \\ инициализируем аккумуляторы для цепочки
По-хорошему, её надо превратить в цикл присвоения единиц последовательно каждому элементу (векторы я кстати создаю до цикла), потому что gp2c сперва создаёт два новых вектора p9 и p10, заполняет их единицами (причем в случае nn не просто единицами а длинными единицами - gen_1) и потом присваивает нашим nf и nn:

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

      {
        GEN p9 = gen_0;   /* vecsmall */
        GEN p10 = gen_0;          /* vec */
        {
          long d;
          p9 = cgetg(len+1, t_VECSMALL);
          for (d = 1; d <= len; ++d)
            p9[d] = 1;
        }
        /* идём по цепочкам */
        /*Здесь начинается проверка конкретного n */
        nf = p9;
        {
          long d;
          p10 = cgetg(len+1, t_VEC);
          for (d = 1; d <= len; ++d)
            gel(p10, d) = gen_1;
        }
        nn = p10;

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

Ага. Ну вот как я и сказал.

Вот такое:

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

\\ делаем инициализацию аккумулятора nf через цикл
   for(j:small=1,len,nf[j]=1);
\\  делаем инициализацию аккумулятора nn через цикл
   for(j:small=1,len,nn[j]=1);
   \\ идём по простым для каждой цепочки

Транслируется в:

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

      /* делаем инициализацию аккумулятора nf через цикл */
      {
        pari_sp btop = avma;
        long j;
        for (j = 1; j <= len; ++j)
        {
          nf[j] = 1;
          if (gc_needed(btop, 1))
            nf = gerepilecopy(btop, nf);
        }
      }
      /*  делаем инициализацию аккумулятора nn через цикл */
      {
        pari_sp btop = avma;
        long j;
        for (j = 1; j <= len; ++j)
        {
          gel(nn, j) = gen_1;
          if (gc_needed(btop, 1))
            nn = gerepilecopy(btop, nn);
        }
      }
      /* идём по простым для каждой цепочки */

Но поскольку там всего по 21 элементу, на общую скорость не влияет.

-- 07.12.2025, 18:05 --

Dmitriy40
Кстати, я заглянул сейчас ещё раз в Си код функции em3(), и там оказывается при подсчёте остатков для цепочки используется таки деление короткого на короткое, вот это:
d=len-(pmods_n0[i]+2*k*pmods_m[i]+len-1)%pr[i];
транслируется в
d = len - smodss(((pmods_n0[i] + ((2*k)*pmods_m[i])) + len) - 1, pr[i]);
А smodss это по буквам в начале и конце -- все small (64 бит)
Так что тут всё вычисление короткое - и сложение-умножения в скобках, и деление результата с остатком на короткое простое из таблицы простых. Ну а деление 64-битных это как мне кажется слёзы по сравнению с isdseudoprime.

[Dmitriy40]

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

Так что тут всё вычисление короткое - и сложение-умножения в скобках, и деление результата с остатком на короткое простое из таблицы простых. Ну а деление 64-битных это как мне кажется слёзы по сравнению с isdseudoprime.

Ну да, я же об этом и сказал:

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

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

И не сказать чтобы слёзы, эти слёзы выполняются до $\pi(2^{20})$ раз, а ispseudoprime лишь один раз.

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

Да, ну я в этих ваших AVX не понимаю -- это ваше поле :)

Про AVX забудем, ОК.
Функции в gmp глянул, там под x64 делений вообще нет, только умножения, молодцы. А под SSE есть, но кажется не в основном цикле.

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

Мне кажется так оно не работает, надо рабоче-крестьянским циклом менять каждый элемент последовательно.

Как раз этого и хотелось избежать и использовать быстрые встроенные функции работы с векторами. Так то пройтись циклом по вектору несложно:
for(j=1,#nmodp, nmodp[j]+=s[j]; if(nmodp[j]>=pr[j], nmodp[j]-=pr[j]); if(nmodp[j]>=#v, next); d=#v-nmodp[j]; p=pr[j]; ... );
На месте ... вся проверка с d и p, только next там не должен обрывать цикл по j, а продолжать его.
Но тут явно не будет векторных операций, на что была надежда при использовании встроенных функций.
С другой стороны, тут сложение, сравнение, условное вычитание (все короткие, по модулю p же) - и это заменяет деление, что должно быть быстрее (на порядок). Фактически этот цикл по j заменяет собой цикл forprime, только обрывать это уже нельзя, ну так не будем делать его до 2^20, сделаем до 2^14 или ещё меньше, а что останется допроверим старым forprime.

[wrest]

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

эти слёзы выполняются до $\pi(2^{20})$ раз,

Сейчас в среднем остаток вычисляется 91 раз на цепочку (короткое деление). А pseudoprime в среднем 1,1 раз на цепочку.

[Dmitriy40]

Выходит можно подобрать такую длину вектора чтобы суммирование по модулю было быстрее 90 коротких делений.

[wrest]

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

Выходит можно подобрать такую длину вектора чтобы суммирование по модулю было быстрее 90 коротких делений.

Это платформозависимо, я думаю. Но экспериментом - можно.

[Dmitriy40]

Всё же есть преимущество у встроенных типов:

? pr=primes([59,2^10]); m=vectorsmall(#pr,j,j*3); s=vectorsmall(#pr,j,j);
? for(k=1,100000, for(j=1,#pr, m[j]+=s[j]; if(m[j]>=pr[j], m[j]-=pr[j]); ); );
time = 3,089 ms.
? pr=primes([59,2^10]); m=vector(#pr,j,j*3); s=vector(#pr,j,j);
? for(k=1,100000, for(j=1,#pr, m[j]+=s[j]; if(m[j]>=pr[j], m[j]-=pr[j]); ); );
time = 6,474 ms.
? pr=primes([59,2^10]); m=vector(#pr,j,j*3); s=vector(#pr,j,j);
? for(k=1,100000, m+=s; for(j=1,#pr, if(m[j]>=pr[j], m[j]-=pr[j]); ); );
time = 2,419 ms.
? pr=primes([59,2^10]); m=vector(#pr,j, Mod(j*3, pr[j])); s=vector(#pr,j, Mod(j, pr[j]));
? for(k=1,100000, m+=s; );
time = 718 ms.
? pr=primes([59,2^10]); n0=3294839204893274832748932748937498327948734; m=(n0*7)\3;
? for(k=1,100000, for(j=1,#pr, d=(n0+m*k+21-1)%pr[j]; ); );
time = 5,382 ms.
? pr=primes([59,2^12]); m=vectorsmall(#pr,j,j*3); s=vectorsmall(#pr,j,j);
? for(k=1,100000, for(j=1,#pr, m[j]+=s[j]; if(m[j]>=pr[j], m[j]-=pr[j]); ); );
time = 10,670 ms.
? pr=primes([59,2^12]); m=vector(#pr,j,j*3); s=vector(#pr,j,j);
? for(k=1,100000, for(j=1,#pr, m[j]+=s[j]; if(m[j]>=pr[j], m[j]-=pr[j]); ); );
time = 24,057 ms.
? pr=primes([59,2^12]); m=vector(#pr,j,j*3); s=vector(#pr,j,j);
? for(k=1,100000, m+=s; for(j=1,#pr, if(m[j]>=pr[j], m[j]-=pr[j]); ); );
time = 8,486 ms.
? pr=primes([59,2^12]); m=vector(#pr,j, Mod(j*3, pr[j])); s=vector(#pr,j, Mod(j, pr[j]));
? for(k=1,100000, m+=s; );
time = 2,528 ms.
? pr=primes([59,2^12]); n0=3294839204893274832748932748937498327948734; m=(n0*7)\3;
? for(k=1,100000, for(j=1,#pr, d=(n0+m*k+21-1)%pr[j]; ); );
time = 19,032 ms.

Видите, сложение векторов почти втрое быстрее.
К сожалению vectorsmall+vectorsmall не работает.
А сложение сразу модулярных типов ещё втрое быстрее.
Длинные деления почти на порядок медленнее.

[wrest]

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

Видите, сложение векторов почти втрое быстрее.
К сожалению vectorsmall+vectorsmall не работает.
А сложение сразу модулярных типов ещё втрое быстрее.
Длинные деления почти на порядок медленнее.

Да, всё так. Но всё немного меняется после компиляции. Уже быстрые функции почти не ускоряются. В итоге модулярные функции перестают явно выигрывать. Я сперва написал с Mod-ами, потом переписал по рабоче-крестьянски. После компиляции (с типизацией всего что можно как small для gp2c) второе оказалось быстрей, просто потому что Mod-ы это GEN-структуры, а не короткие 64-битные числа. И вам после вычислений ещё нужно достать из Mod-а результат через lift.
То есть если вы пишете под интерпретатор, то надо использовать встроенные функции и по максимуму запихнуть вычисления в них. Типа сложения векторов для интерпретатора надо писать просто как c=a+b а не покомпонентно for(i=1,#a, c[i]=a[i]+b[i]). А если под gp2c, то надо уже проверять.