интерактивный курс: введение в программирование на 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)