интерактивный курс: введение в программирование на PARI/GP

[waxtep]

Dmitriy40, wrest, большое спасибо! Имплементировал ваши рекомендации частично: предвычисление остатков и присвоение B[x]=g в forstep-е. Заметно полегчало: для #B=84 - быстрее в 4,7 раза, для #B=124 - в 7,1 раза, 6,5 минут вместо 45. Следующая цель, #B=172 обещает считаться в районе двух суток, можно еще наверное чуть повозиться с оптимизацией и запустить. Код выглядит более устрашающе, зато работает! Типа (самый внутренний цикл, над ним еще два):

1. for(j=1,2*k-2, 
2. g=gi[j]+2; lr=lri[j]; 
3. i=(g+1)\2; s=(g+1)%2+1; d=6*i+2*s-3; 
4. r01=2*lr*i; r02=bL-bR-2; 
5. r1=r01%d; r2=r02%d; r3=(r02-r01)%d; r4=(r02+r01)%d; r5=(-r01)%d; 
6.  
7. if(j%2==1, \\one shot left and one shot right of the island; strong assumption for the optimum strategy, valid for k=2,3,4 
8. forstep(x=bL-1,1,-d,if(B[x]==0,B[x]=g)); 
9. forstep(x=bL-1-r1,1,-d,if(B[x]==0,B[x]=g)); 
10. forstep(x=bR+1+r2,sb,d,if(B[x]==0,B[x]=g)); 
11. forstep(x=bR+1+r3,sb,d,if(B[x]==0,B[x]=g)); 
12. , 
13. forstep(x=bR+1,sb,d,if(B[x]==0,B[x]=g)); 
14. forstep(x=bR+1+r5,sb,d,if(B[x]==0,B[x]=g)); 
15. forstep(x=bL-1-r2,1,-d,if(B[x]==0,B[x]=g)); 
16. forstep(x=bL-1-r4,1,-d,if(B[x]==0,B[x]=g)); 
17. ); \\if j 
18.  
19. for(x=bR+1,sb+1,if(B[x]==0,bR=x-1;break)); 
20. forstep(x=bL-1,0,-1,if(B[x]==0,bL=x+1;break)); 
21. ); \\for j 

Тут наверняка есть еще, что улучшить, но уже приятно. Спасибо!

[waxtep]

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

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

Это нормальный баг я беру размер полоски с запасом, чтобы почти наверняка было B[1]=B[#B]=0. Сделал его на всякий случай еще неправильнее, чтобы в таком случае программа крэшилась на попытке обращения к B[0] или B[#B+1]. Это не со зла, я просто не знаю максимально возможный размер массива, есть только какие-то эвристические соображения.

PS: у меня оказывается выше по потоку был страшный баг, в самом внешнем цикле, - из-за него пропускал в итоге значительную часть возможных комбинаций. Чудом это не влияло на результат (искомых максимальных комбинаций несколько; какая-то находилась), и тем не менее. Я пытался получить все возможные комбинации $\pm1$ так:
for(lrc=0,2^(2*k-2)-1,
lri=vectorsmall(2*k-2,i,2*(iferr(binary(lrc)[i],EE,0)-1/2));
...
Но так это не работает; рабочий непричесанный вариант (двоичного представления с лидирующими нулями) такой:
lri=vectorsmall(2*k-2,i,lrb=binary(lrc);slrb=#lrb;if(i<=2*k-2-slrb,-1,2*(lrb[i-2*k+2+slrb]-1/2)));

Может быть, эту забавную багу можно будет как-то использовать для сокращения перебираемых вариантов, надо будет подумать

[Dmitriy40]

waxtep в сообщении #1717364 писал(а):

Я пытался получить все возможные комбинации $\pm1$ так:

Проще так:
lri=Vecrev([bittest(lrc,t)*2-1|t<-[0..2*k-3]]);
Или так:
lri=vectorsmall(2*k-2,t, bittest(lrc,2*k-2-t)*2-1);

[waxtep]

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

Проще так:

Точно, биттест же, спасибо!

-- 06.02.2026, 04:01 --

waxtep в сообщении #1717364 писал(а):

Сделал его на всякий случай еще неправильнее, чтобы в таком случае программа крэшилась на попытке обращения к B[0] или B[#B+1].

Кстати, для #B=172 это произошло, спасибо, что и это отметили; расширил в полтора раза, время счета вроде растет заметно слабее, чем линейно по #B

[wrest]

Как говорится, хозяйке на заметку.
Обнаружил следующее, в механизме трансляции-компиляции при помощи скрипта командной строки gp2c-run который есть в составе транслятора gp2c
Этот скрипт (это обычный shell скрипт), среди прочей активности, ищет строки в файле .c начинающиейся с "GP;" и вставляет то, что идёт после "GP;" в файл *.gp.run
В доках https://pari.math.u-bordeaux.fr/pub/par ... /gp2c.html это описано так:

Цитата:

More generally, gp2c-run automatically passes any line in the C file starting with ’GP;’ to GP at start-up.

Казалось бы, есть инструмент, при помощи которого можно управлять созданием файла *.gp.run просто вставляя строки в исходный файл *.gp такого содержания

/*
GP;Here goes my command for gp interpreter
*/

Но не тут-то было. При создании файла .c транслятор gp2c действительно переносит комментарии из *.gp в *.c но из-за расстановки отступов, перед комментариями могут оказаться пробелы (а может и табуляции, тут не проверял) и в результат трансляции, файл *.c вышеописанное попадёт как

... some already indented .С code
   /*
   GP;Here goes my command for gp interpreter
   */

И тогда комбинация "GP;" будет не в начале строки.
В скрипте gp2c-run обработка это устроена так (уже отработал транслятор в *.c и создал файл $name.c):

grep "^GP;" $name.c | sed 's/^GP;//' >$name.run

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

grep "^[ \t]*GP;" "$name.c" | sed 's/^[ \t]*GP;//' >"$name.run"

Тогда в файл *.gp.run попадёт то, что нам надо и наши команды будут выполняться при инициализации интерпретатора.
Понадобилось это для того, чтобы после команды gp2c-run *.gp мы видели не промт "?" интерпретатора pari/gp, а могли сразу запустить исполнение нашего кода/функции, ну и если надо, перед этим настроив например размер стека и сделав ещё что-то нам нужное, а после исполнения, могли бы исполнить команду quit ну и т.п.
Разработчики, конечно, думали о другом: что люди накомпилируют себе функций, и будут их вызывать в коде обычного *.gp скрипта, не забыв "установить" эти функции перед вызовом и т.п. Но нам-то виднее как надо пользоваться

[waxtep]

waxtep в сообщении #1717341 писал(а):

Код выглядит более устрашающе, зато работает!

Кажется, в человеческом виде этот код может выглядеть так, с предрасчетом всего, что можно, и отказом от лишнего функционала:

1. for(j=1,2*k-2, 
2. g=gi[j]+2; lr=lri[j]; 
3. i=ii[g]; d=dd[g]; 
4. jc=1-jc; \\one shot left and one shot right of the island; strong assumption for the optimum strategy, valid for k=2,3,4 
5.  
6. \\do period d shots with fork i, knowing position of the lr-th tine bb 
7. bb=if(jc,bL-1,bR+1); 
8. x0L=(bb-(1+lr)*i)%d; if(!x0L,x0L=x0L+d); 
9. x0R=sb-(sb-bb-(1-lr)*i)%d; 
10. forstep(x=x0L,x0R,d,B[x]=g;B[x+2*i]=g); 
11. if((xf=x0L-d+2*i)>0,B[xf]=g); if((xf=x0R+d-2*i)<=sb,B[xf]=g);  
12.  
13. \\calcualte new island bounds 
14. for(x=bR+1,sb+1,if(B[x]==0,bR=x-1;break)); 
15. forstep(x=bL-1,0,-1,if(B[x]==0,bL=x+1;break)); 
16. ); \\for j 

Некорректное поведение в последних двух строках (выход индекса за границы диапазона) оставил сознательно, пусть крэшится, если уткнулись в края. Старый код еще сколько-то часов посчитает, и попробую

[waxtep]

waxtep в сообщении #1717647 писал(а):

Старый код еще сколько-то часов посчитает, и попробую

Стало быстрее... в 1,2 раза (270 с вместо 329 с для k=5, #B=6(k^2+k+1)=186) Похоже, насильно мил не будешь отсюда уже ничего не выжать (нет оптимизации в разы без тотальной переделки алгоритма). Ну хоть код теперь выглядит попритянее.

При этом, старый код считал для k=6, #B=258 52 часа, ускорение на 20% здесь и для больших k на спасает. Ладно, что-нибудь еще придумаем

[wrest]

maxal
Помогите советом, пож-ста.
Хочу задействовать встроенную факторизацию методом MPQS
Для этого устанавливаю встроенную функцию
install(mpqs,G)
Она возвращает вектор
? type(mpqs(nextprime(10^8)*nextprime(10^10)))
time = 1 ms.
%3 = "t_VEC"
?
но присвоить это переменной не получается
? v=mpqs(nextprime(10^8)*nextprime(10^10))
*** at top-level: v=mpqs(nextprime(10^8)*nextprime(10^10))
*** ^----------------------------------------
*** bug in PARI/GP (Segmentation Fault), please report.
?
Однако, покомпонентно доступ есть:
? component(mpqs(nextprime(10^8)*nextprime(10^10)),1)
time = 1 ms.
%4 = 10000000019
?
Однако, так получается получить доступ только однократно, вектор результатов может меняться (множители меняются местами). Поэтому, я в обёртке делаю так:
my_mpqs(n:int)=my(my_div:int=component(mpqs(n),1));return([my_div,n\my_div]);

Насколько я понял (но это не точно), вектор 6-компоненьный, компоненты:
1 - первый множитель
2 - степень первого множителя
3 - вызывает ошибку при отображении, что-то типа NULL или фиг знает
4 - второй множитель
5 - степень второго множителя
6 - вызывает ошибку при отображении, что-то типа NULL или фиг знает

Как ещё можно установить библиотечную mpqs() чтобы её использовать, в какую обёртку обернуть?

Суть в том, что mpqs() не делает trial division и не пробует другие методы, ищет один множитель и останавливается. А именно это и надо.

[Dmitriy40]

wrest
Замените покомпонентно NULL на скажем 0:

? print(mpqs(nextprime(10^8)*nextprime(10^10)))
[10000000019, 1, NULL, 100000007, 1, NULL]
? isnull(x=0)=x
%2 = (x=0)->x
? [isnull(t)|t<-mpqs(nextprime(10^8)*nextprime(10^10))]
%3 = [10000000019, 1, 0, 100000007, 1, 0]
? apply(t->isnull(t),mpqs(nextprime(10^8)*nextprime(10^10)))
%4 = [10000000019, 1, 0, 100000007, 1, 0]

-- 22.03.2026, 21:58 --

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

ищет один множитель и останавливается.

Не обязательно останавливается и может вернуть длинный вектор:
? print(n=vecprod([nextprime(10^t)|t<-[5..11]]),": ",mpqs(n))
100003497908545547874191502297671839116672752648391577603: [100003, 1, NULL, 10000000019, 1, NULL, 1000003, 1, NULL, 100000007, 1, NULL, 10000019, 1, NULL, 100000000003, 1, NULL, 1000000007, 1, NULL]
time = 1,998 ms.
А может одно число:
? print(n=vecprod([nextprime(10^t)|t<-[4..10]]),": ",mpqs(n))
10007350035407932474708105660686777770584839024407: 10007

[wrest]

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

Замените покомпонентно NULL на скажем 0:

Точно, apply() помогает, спасибо. Каждый треий компонент это флаг простоты. NULL - неизвестно. 0 - составной.
Заодно выяснилось, что mpqs может вернуть и больше множителей, чем два.

-- 22.03.2026, 22:03 --

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

А может одно число:

Засада, однако. Я разочарован

[Dmitriy40]

Попытался сейчас зайти на PARI/GP online, 32-битные версии работают, обе 64-битные нет, и

stable писал(а):

wasm streaming compile failed: CompileError: WebAssembly.instantiateStreaming(): invalid table elements limits flags @+9493
falling back to ArrayBuffer instantiation
failed to asynchronously prepare wasm: CompileError: WebAssembly.instantiate(): invalid table elements limits flags @+9493
Aborted(CompileError: WebAssembly.instantiate(): invalid table elements limits flags @+9493)

и

development писал(а):

wasm streaming compile failed: CompileError: WebAssembly.instantiateStreaming(): invalid table elements limits flags @+9869
falling back to ArrayBuffer instantiation
failed to asynchronously prepare wasm: CompileError: WebAssembly.instantiate(): invalid table elements limits flags @+9869
Aborted(CompileError: WebAssembly.instantiate(): invalid table elements limits flags @+9869)

Беда однако! 32-битные версии заметно медленнее при работе с числами больше 2^32 (например isprime).

[Dmitriy40]

Забавный глюк функции gettime(), итерации цикла:

4296.066s
4282.009s
-4290691.042s
4282.212s
4279.873s
4273.694s
4281.651s

Очевидно произошло переполнение 32 битов результата.
Странно что за много лет впервые наткнулся. Да и другие три потока, работающие столько же времени как и этот, не наткнулись.
При этом нигде в доке не сказано что она возвращает лишь 32 бита. Единственный намёк - "The library syntax is long gettime(void)." в описании функции. Как и про getabstime(). А вот getwalltime() возвращает нормальное длинное число (тип результата GEN).

PS. Кстати новая версия вышла в начале недели, 28.09.2026, уже 2.19.0.
Никаких интересных/полезных изменений не заметил.

[wrest]

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

Никаких интересных/полезных изменений не заметил.

А они есть :) bitclear, bitflip, bitset как раз полезны по теме простых чисел (для реализации разных решёт). Неясно пока (не пробовал) быстро ли работают. Но заявлено, что все три работают по месту (по ссылке), без копирования битового вектора, в отличие от bitand и bitor
Я в это упёрся при реализации продвинутого primepi (по мотивам primecount). Битовые операции которые очень быстро работают в Си (или тем более аккуратно оптимизированные на асме), в pari/gp тормозят из-за перекидывания битового вектора туда-сюда.
И я ещё так понял что там наконец-то начались подвижки пока под капотом) по оптимизации параллельности. Но это не точно, не уверен.

[Dmitriy40]

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

А они есть :) bitclear, bitflip, bitset

Хм, в changelog 2.19.0 их нету, я смотрел туда. А, их добавили в тестовой 2.18.1, которая stable не стала. Эти гады скрывают changelogs кроме самого последнего, фиг найдёшь.
Да, в 2.18.1 много чего поправили и добавили.
И опять исправляли проблемы с forprime. Когда они наконец их порешают уже, уже страшно его использовать. Хотя ни одной ошибки и не вижу.
Опять поменяли поведение Mod(a,b)==Mod(A,B).

[wrest]

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

Эти гады скрывают changelogs кроме самого последнего, фиг найдёшь.

Да, бардачно немного. Вот тут сводка 2.17.4 ->2.19.0
http://pari.math.u-bordeaux.fr/archives ... 00005.html