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

[wrest]

Улучшил функцию,

Фунция sqm(a,m) находит все решения n<m сравнения n^2 = a (mod m) включая составные модули

в части поиска базовых решений для нечетных простых делителей модуля.
DeepSeek там изначально устроил перебор решений, что для больших простых делителей модуля будет неэффективно. Сейчас перебор устранён использованием функции issquare() которая не только определяет есть ли решение, но если есть сразу и возвращает его.
можно было бы воспользоваться функцией kronecker() для определения есть ли решение, но не знаю будет ли какой-то выигрыш по производительности.
Думаю что функция теперь годная.
Структурированный вид (отступы, комментарии)

(Оффтоп)

1. \\ v2 Написано wrest с помощью DeepSeek  
2. \\ Переделан поиск базовых решений для нечетных простых делителей. 
3. \\ Теперь вместо перебора используется встроенная функция 
4. \\ Это должно существенно ускорить работу при очень больших простых делителях модуля 
5. \\ Фунция sqm(a,m) находит все решения n<m сравнения n^2 = a mod m и работает для составных модулей  
6. \\ модуль факторизуется, затем отдельно обрабатываются степени двойки и степени нечетных простых.  
7. \\ Найденные корни для нечетных простых делителей поднимаются до соответстующих степеней   
8. \\ с использованием леммы Гензеля. Затем корни по простым делителям модуля объединяются с использованием КТО  
9. \\ китайской теоремы об остатках  
10. sqm(a, m) = { 
11.     my(factors, solutions_list, mod_list, p, e, sols, solutions, new_solutions, crt_val, i, k, j, x0, x1, x2); 
12.  
13.     \\ Специальные случаи 
14.     if(m == 1, return([0])); 
15.     if(a == 0, return([0])); 
16.      
17.     factors = factor(m); 
18.     if(#factors[,1] == 0, return([])); 
19.      
20.     solutions_list = []; 
21.     mod_list = []; 
22.      
23.     for(i = 1, #factors[,1], 
24.         p = factors[i,1]; 
25.         e = factors[i,2]; 
26.         sols = []; 
27.          
28.         if(p == 2, 
29.             \\ Решение для 2^1 
30.             if(e == 1, 
31.                 if(a%2 == 0, sols = [0]) 
32.             , 
33.             \\ Решение для 2^2 
34.             if(e == 2, 
35.                 if(a%4 == 0, sols = [0, 2]); 
36.                 if(a%4 == 1, sols = [1, 3]) 
37.             , 
38.             \\ Решение для 2^k, k>=3 
39.                 \\ Рекурсивный поиск решений 
40.                 my(x0 = [], x0_temp = []); 
41.                 if(a%4 == 0, x0 = [0, 2]); 
42.                 if(a%4 == 1, x0 = [1, 3]); 
43.                  
44.                 for(k_temp = 3, e, 
45.                     x0_temp = []; 
46.                     for(idx = 1, #x0, 
47.                         x = x0[idx]; 
48.                         if((x^2 - a) % (2^k_temp) == 0, 
49.                             x0_temp = concat(x0_temp, [x]); 
50.                             x0_temp = concat(x0_temp, [x + 2^(k_temp-1)]); 
51.                         ) 
52.                     ); 
53.                     x0 = x0_temp; 
54.                     if(#x0 == 0, break); 
55.                 ); 
56.                 sols = x0; 
57.             )) 
58.         , 
59.             \\ Оптимизированная обработка нечетных простых 
60.             if(issquare(Mod(a,p),&x1),x1=lift(x1),return([])); 
61.             x2 = p - x1; 
62.             x0 = [x1, x2]; 
63.              
64.             \\ Подъем Гензеля для p^k 
65.             for(k_step = 2, e, 
66.                 my(temp = []); 
67.                 for(j = 1, #x0, 
68.                     my(x = x0[j], f = x^2 - a, df = 2*x); 
69.                     if(Mod(df, p) != 0, 
70.                         my(inv_df = lift(Mod(df, p^k_step)^(-1))); 
71.                         x = (x - f * inv_df) % (p^k_step); 
72.                         temp = concat(temp, [x]); 
73.                     , 
74.                         if(Mod(f, p^k_step) == 0, 
75.                             for(t = 0, p-1, 
76.                                 my(new_x = (x + t*p^(k_step-1)) % (p^k_step)); 
77.                                 temp = concat(temp, [new_x]); 
78.                             ) 
79.                         ) 
80.                     ) 
81.                 ); 
82.                 x0 = temp; 
83.                 if(#x0 == 0, break); 
84.             ); 
85.             sols = x0; 
86.         ); 
87.          
88.         if(#sols == 0, return([])); 
89.         solutions_list = concat(solutions_list, [sols]); 
90.         mod_list = concat(mod_list, [p^e]); 
91.     ); 
92.      
93.     \\ Объединение решений 
94.     if(#solutions_list == 0, return([])); 
95.     if(#solutions_list == 1, return(solutions_list[1])); 
96.      
97.     solutions = solutions_list[1]; 
98.     for(k = 2, #solutions_list, 
99.         new_solutions = []; 
100.         for(i = 1, #solutions, 
101.             for(j = 1, #solutions_list[k], 
102.                 crt_val = chinese(Mod(solutions[i], mod_list[1]), Mod(solutions_list[k][j], mod_list[k])); 
103.                 if(crt_val != -1, 
104.                     new_solutions = concat(new_solutions, [lift(crt_val)]); 
105.                 ) 
106.             ) 
107.         ); 
108.         solutions = new_solutions; 
109.         if(#solutions == 0, return([])); 
110.     ); 
111.      
112.     vecsort(solutions,,8); 
113. } 

Вид в одну строку, без пробелов и комментариев (может быть удобней для копипаста в pari)

(Оффтоп)

1. sqm(a,m)=my(factors,solutions_list,mod_list,p,e,sols,solutions,new_solutions,crt_val,i,k,j,x0,x1,x2);if(m==1,return([0]));if(a==0,return([0]));factors=factor(m);if(#factors[,1]==0,return([]));solutions_list=[];mod_list=[];for(i=1,#factors[,1],p=factors[i,1];e=factors[i,2];sols=[];if(p==2,if(e==1,if(a%2==0,sols=[0]),if(e==2,if(a%4==0,sols=[0,2]);if(a%4==1,sols=[1,3]),my(x0=[],x0_temp=[]);if(a%4==0,x0=[0,2]);if(a%4==1,x0=[1,3]);for(k_temp=3,e,x0_temp=[];for(idx=1,#x0,x=x0[idx];if((x^2-a)%(2^k_temp)==0,x0_temp=concat(x0_temp,[x]);x0_temp=concat(x0_temp,[x+2^(k_temp-1)]);));x0=x0_temp;if(#x0==0,break););sols=x0;)),if(issquare(Mod(a,p),&x1),x1=lift(x1),return([]));x2=p-x1;x0=[x1,x2];for(k_step=2,e,my(temp=[]);for(j=1,#x0,my(x=x0[j],f=x^2-a,df=2*x);if(Mod(df,p)!=0,my(inv_df=lift(Mod(df,p^k_step)^(-1)));x=(x-f*inv_df)%(p^k_step);temp=concat(temp,[x]);,if(Mod(f,p^k_step)==0,for(t=0,p-1,my(new_x=(x+t*p^(k_step-1))%(p^k_step));temp=concat(temp,[new_x]);))));x0=temp;if(#x0==0,break););sols=x0;);if(#sols==0,return([]));solutions_list=concat(solutions_list,[sols]);mod_list=concat(mod_list,[p^e]););if(#solutions_list==0,return([]));if(#solutions_list==1,return(solutions_list[1]));solutions=solutions_list[1];for(k=2,#solutions_list,new_solutions=[];for(i=1,#solutions,for(j=1,#solutions_list[k],crt_val=chinese(Mod(solutions[i],mod_list[1]),Mod(solutions_list[k][j],mod_list[k]));if(crt_val!=-1,new_solutions=concat(new_solutions,[lift(crt_val)]);)));solutions=new_solutions;if(#solutions==0,return([])););vecsort(solutions,,8) 

[wrest]

Натолкнулся тут на такую штуку.
В принципе, вместо скажем n=sqrt(Mod(a,m)) или n=0;issquare(Mod(a,m),&n) которые возвращают одно решение n сравнения $n^2 = a \pmod{m}$ можно было использовать n=polrootsmod(x^2-a,m) на случай простого нечётного m
Но тут у pari/gp какая-то крайне непрозрачная шняга с "виртуальными" (символьными?) переменными фигурирующими в многочленах. В частности, конструкция с polrootsmod() может вообще не работать если x использовался где-то раньше. Тогда может помочь kill(x) перед запуском программы. Но я натокнулся на другое. При каждом использовании polrootsmod() как бы создавался новый x. Всё работало, 65535 раз. А на 65536-й система сказала что-то типа "закончилось место под переменные".
Иллюстрация

? prm(a,m)=kill(x);polrootsmod(x^2-a,m)
%1 = (a,m)->kill(x);polrootsmod(x^2-a,m)
? variables
%2 = [x, y]
? for(a=1,10,n=prm(a,5))
? variables
%4 = [x, y, x, x, x, x, x, x, x, x, x, x]
?

Решение, видимо, в том что надо делать объявление my(x='x) в функциях где используются "символьные" вычисления. Иначе могут быть сюрпризы.

? prm(a,m)=my(x);polrootsmod(x^2-a,m)
%18 = (a,m)->my(x);polrootsmod(x^2-a,m)
? prm(2,5)
  ***   at top-level: prm(2,5)
  ***                 ^--------
  ***   in function prm: my(x);polrootsmod(x^2-a,m)
  ***                          ^--------------------
  *** polrootsmod: incorrect type in factormod (t_INT).

Но:

? prm(a,m)=my(x='x);polrootsmod(x^2-a,m)
%19 = (a,m)->my(x='x);polrootsmod(x^2-a,m)
? prm(2,5)
%20 = []~
?

[wrest]

В общем, что-то я погрузился не на шутку в решение n^2 = a mod m
Предыдущая функция с ошибками
Вот новая, всё работает на встртенных функциях, текст магический, как это случается с pari/gp
Но зато короткий.

1. sqrtmod(n, m) = { 
2.     my(p, r, t, l, q, roots, x='x, s); 
3.     /* функция sqrtmpd(n,m) вычисляет все корни сравнения x^2 = n mod m */ 
4.     /* в том числе для составных m */ 
5.     /* если корней нет, возвращает пустой вектор */ 
6.  
7.     /* Проверяем, является ли n квадратичным вычетом по модулю m */ 
8.     /* На будущее: тут по сути делается факторизация, а потом она же */ 
9.     /* делается ещё раз явно если корни есть -- надо бы оптимизировать */ 
10.     if(!issquare(Mod(n, m)), return([]));  
11.  
12.     /* Факторизуем модуль m */ 
13.     p = factorint(m); 
14.     l = matsize(p)[1]; 
15.  
16.     /* Находим квадратные корни для каждой степени простого делителя */ 
17.     t = vector(l,i, sqrt(n + O(p[i,1]^p[i,2])) ); 
18.     /* Понижаем модули до минимальной p-адической точности */ 
19.     t = apply(x->Mod(lift(x), x.p^padicprec(x, x.p)), t); 
20.  
21.     /* Комбинируем корни с помощью КТО */ 
22.     roots = []; 
23.     forvec(s = vector(l, i, [0,1]), 
24.         x = chinese(vector(l, i, t[i] * (-1)^s[i])); 
25.         roots = concat(roots, [x]); 
26.     ); 
27.  
28.     /* Удаляем дубликаты корней */ 
29.     roots = Set(roots); 
30.     if(#roots == 0, return([])); 
31.  
32.     /* Генерируем все решения, добавляя кратные корни (если нужно) */ 
33.     /* и преобразуем в обычные числа */ 
34.     r = []; 
35.     for(i = 1, #roots, 
36.         x = roots[i]; 
37.         step = x.mod; 
38.         for(j = 0, m \ step - 1, 
39.             r = concat(r, lift(x) + j * step); 
40.         ) 
41.     ); 
42.  
43.     /* Сортируем и возвращаем результат */ 
44.     return(vecsort(r)); 
45. }; 

[wrest]

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

Вот новая, всё работает на встртенных функциях, текст магический, как это случается с pari/gp

Кажется, это код ув. maxal отсюда post430305.html#p430305

[gris]

Супернубовский вопрос: как очистить историю? В смысле вызова предыдущих команд(?) стрелочкой вверх после gp>

[wrest]

gris в сообщении #1693163 писал(а):

Супернубовский вопрос: как очистить историю? В смысле вызова предыдущих команд(?) стрелочкой вверх после gp>

История хранится в файле (зависит от операционной системы).
По команде
? default(histfile)
можно посмотреть как он называется

В файле можно очистить. Чтобы не накапливалось, можно поставить histsize = 10 например
? default(histsize,10)

[gris]

wrest, спасибо, но не смейтесь. Называется от gp_history.txt, но где он находится? В C:\Program Files\Pari32-2-13-4\data не могу найти
Вероятно, он скрытый. А я не вижу скрытые файлы и не могу изменить доступ
Но воспользовался вашим советом в другую сторону. Пусть будет не 5000, а 6000 А дальше видно будет
При перезапуске не меняется

[wrest]

gris в сообщении #1693166 писал(а):

Называется от gp_history.txt, но где он находится? В C:\Program Files\

У меня в виндовсе нет pari/gp, так что не знаю. Воспользуйтесь поиском в проводнике виндовс...

-- 03.07.2025, 16:50 --

gris в сообщении #1693166 писал(а):

При перезапуске не меняется

Надо править файл gprc.txt

[gris]

Получилось! Там проблема в доступе к скрытым файлам. Надо новый комп и всё переустановить:facepalm: Но теперь файл истории новенький! Спасибо за дельный совет!

[maxal]

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

Кажется, это код ув. maxal отсюда post430305.html#p430305

Вот новая версия, которая использует библиотечную функцию Zn_quad_roots:

install(Zn_quad_roots, GGG);
isNULL(x=[1,[]])=x;

\\ all square roots of q modulo N
{ asqrtintmod2(N,q) = my(r); 
    r = isNULL(Zn_quad_roots(N, 0, -q));
    if(#r!=2, error("asqrtintmod2: check input ",N," ",q));
    r[2] * Mod(1,r[1]);
}

[wrest]

maxal
Непонятно... Как при помощи новой версии найти все числа, меньшие 10000, квадраты которых дают остаток 2121 при делении на 10000?
И зачем вы поменяли местами модуль и остаток?

Выглядит конечно загадочно это. Что за библиотека и что это значит, откуда берется?

[maxal]

wrest
Функция из библиотеки PARI, которой нет в GP, но она подгружается через функцию install. Ваша задача решается например так:

? foreach(asqrtintmod2(10000,2121),r, forstep(s=0,10000,r, print(s) ))
1011
6011
2261
7261
2739
7739
3989
8989

[wrest]

Посмотрел документацию на libpari ("Users' Guide to the PARI library")
Оказывается, у них там несколько вариантов функций поднятия корней по лемме Гензеля.
Например

GEN Zp_sqrtlift(GEN b, GEN a, GEN p, long e) let $a,b,p$ be t_INTs, with $p > 1$ odd, such that $a^2 \equiv b \pmod{p}$. Returns a t_INT $A$ such that $A^2 \equiv b \pmod{p^e}$. Special case of Zp_sqrtnlift.

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

А вот использованная ув. maxal функция Zn_quad_roots в документации не описана.

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

[maxal]

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

А вот использованная ув. maxal функция Zn_quad_roots в документации не описана.

Описана на стр. 104 в https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.17.2/libpari.pdf

[wrest]

maxal в сообщении #1693359 писал(а):

Описана на стр. 104

Хм.. действительно.
A... Ясно. Я почему-то смотрел в User’s Guide to the PARI library (version 2.7.7).