интерактивный курс: введение в программирование на 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 на спасает. Ладно, что-нибудь еще придумаем