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

[dmd]

Параллельной версии hyperellratpoints нету в текущей 2.17.1, а хотелось бы. Если вдруг есть идеи, подскажите хотя бы примерное направление, как самому реализовать параллельный аналог этой функции.

[wrest]

dmd в сообщении #1671197 писал(а):

Если вдруг есть идеи, подскажите хотя бы примерное направление, как самому реализовать параллельный аналог этой функции.

Как самому не знаю, но я бы предложил задать вопрос разработчикам, сп. https://pari.math.u-bordeaux.fr/lists.html
И вот: https://pari.math.u-bordeaux.fr/archive ... 00006.html

[gris]

Иногда хочется определить целочисленный отрезок длины 1, в котором задано действительное число, например, квадратный корень из неквадрата.
gp > [floor(sqrt(150)),floor(sqrt(150))+1]
%75 = [12, 13]
А вот нужно типа
75886^17+487^13 = 917...143
[floor(sqrt(75886^17+487^13)),floor(sqrt(75886^17+487^13))+1]
И тут
*** floor: precision too low in truncr (precision loss in truncation).
Хотя бы до Е300

[wrest]

gris
? prec=floor(log(75886^17+487^13)/log(10))+10
%1 = 92
? default(realprecision, prec)
? [floor(sqrt(75886^17+487^13)),floor(sqrt(75886^17+487^13))+1]
%2 = [302952208377199909263956604869077972350392, 302952208377199909263956604869077972350393]
?

-- 18.06.2025, 19:02 --

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

Хотя бы до Е300

default(realprecision,300)
или
\p 300

[Dmitriy40]

Достаточно было заменить
floor(sqrt(...))
на
sqrtint(...)
и вообще не связываться с плавающими числами, оставаться в целых, которые считаются точно независимо от realprecision, в том числе и sqrtint():

? sqrtint(75886^17+487^13)
%1 = 302952208377199909263956604869077972350392
? sqrtint(75886^27+487^23)
%2 = 762399196777595305347944476443436697974097553827078983457348436251

[gris]

wrest, спасибо! Посчиталось быстро.
Проблема дурацкая, но очень важная. Найти расстояние от праймориала до окаймляющих квадратов.

{pp=1;
forprime(p=2, 1000, pp=pp*p;
   s1=floor(sqrt(pp));
   print(p,"# = ...",pp%1000," (",s1,",", s1+1,")",
         " d= (",pp-s1^2,",", (s1+1)^2-pp,")");
)}
2# = ...2 (1,2) d= (1,2)
3# = ...6 (2,3) d= (2,3)
5# = ...30 (5,6) d= (5,6)
7# = ...210 (14,15) d= (14,15)
11# = ...310 (48,49) d= (6,91)
13# = ...30 (173,174) d= (101,246)
17# = ...510 (714,715) d= (714,715)
19# = ...690 (3114,3115) d= (2694,3535)
...
997# = ...910 (442609767....

Все проверенные лежат далеко от точных квадратов. Ну кроме самых первых. ВотЪ Я даже составил последовательность кратчайших расстояний.
[1, 2, 5, 14, 6, 101...]
Dmitriy40, спасибо, подумаю над этим.
+++ а вдруг какой-то праймик вплотную приблизится к точной нечётной степени? Кстати, я заметил, что некоторые пары расстояний в точности равны парам квадратных корней. Нет ли тут чего интересного? Больше такого не встречалось... Ну то, что сумма двух расстояний равна 2s1+1 это понятно из формулы разности квадратов. Видно, что праймориал не тяготеет к определённому концу отрезка последовательных квадратов. Надо это изучить.

[Dmitriy40]

? pp=1; forprime(p=2,67, pp*=p; x=sqrtint(pp); print1(min(pp-x^2,(x+1)^2-pp),", "); );
1, 2, 5, 14, 6, 101, 714, 2694, 8774, 64874, 175470, 980201, 3368769, 36419829, 212951314, 5138843466, 10999470170, 103728576730, 2555438890809,
? pp=1; forprime(p=2,67, pp*=p; x=sqrtint(pp); print(min(pp-x^2,(x+1)^2-pp)/pp*1.0); );
0.50000000000000000000000000000000000000
0.33333333333333333333333333333333333333
0.16666666666666666666666666666666666667
0.066666666666666666666666666666666666667
0.0025974025974025974025974025974025974026
0.0033633033633033633033633033633033633034
0.0013986013986013986013986013986013986014
0.00027774083501637681204244671736931798851
3.9328912663143380601988759210457958607 E-5
1.0027368793806008635126583892139194983 E-5
8.7489814113569049244125917314340579987 E-7
1.3208942051222199473251216928963594430 E-7
1.1072361814729278073480506101913681429 E-8
2.7838029049598939396557674819216990663 E-9
3.4632436581323949475747003560141913485 E-10
1.5768567542475216729386661477646483043 E-10
5.7206662125718376531535587903857047973 E-12
8.8438919122043606018501622491359578324 E-13
3.2518889360765711493415252504127058019 E-13
? pp=1; forprime(p=2,67, pp*=p; x=sqrtint(pp); print(min(pp-x^2,(x+1)^2-pp)/x*1.0); );
1.0000000000000000000000000000000000000
1.0000000000000000000000000000000000000
1.0000000000000000000000000000000000000
1.0000000000000000000000000000000000000
0.12500000000000000000000000000000000000
0.58381502890173410404624277456647398844
1.0000000000000000000000000000000000000
0.86512524084778420038535645472061657033
0.58743974290305302624531333690412426352
0.80654946912996991322077728323843150906
0.39181491562816101322126925077985615366
0.35982523421471214561270260338907890047
0.19313267370190206590455151879251295565
0.31841109599161300186145469493726340849
0.27156993377895701442516898135598274074
0.90017887277718249298595613281771172821
0.25084715935818523067565388633157643326
0.30288022728654018477801684788053120325
0.91159220355618851699971718159342535187

[wrest]

В соседней теме возникла потребность определять является ли число квадратичным вычетом по составному модулю, причем в разложении модуля есть степени двойки.
На удивление, встроенной в pari/gp такой функции нет. Есть только kronecker() которая проверяет по простому нечётному модулю.
Пришлось делать. Ну вот чтобы не пропала функция, даю её сюда.

1. is_quadratic_residue(a, m) = { 
2.     if(m <= 0, error("Modulus must be positive")); 
3.     if(m == 1, return(1)); 
4.     my(F = factor(m), p, k, s, a1, k_prime); 
5.     for(i = 1, #F~, 
6.         [p, k] = F[i,]; 
7.         s = valuation(a, p); 
8.         if(s >= k, \\ a делится на p^k 
9.             next() 
10.         ); 
11.         if(s % 2 == 1, \\ нечётная кратность 
12.             return(0) 
13.         ); 
14.         a1 = a / p^s; \\ Теперь a1 не делится на p 
15.          
16.         if(p == 2, 
17.             k_prime = k - s; 
18.             if(k_prime == 0, \\ Оставшаяся часть тривиальна 
19.                 next() 
20.             ); 
21.             \\ Проверяем условия для степени двойки 
22.             if(k_prime == 1, \\ Модуль 2 
23.                 next() 
24.             ); 
25.             if(k_prime == 2, \\ Модуль 4 
26.                 if(a1 % 4 != 1, return(0)); 
27.                 next() 
28.             ); 
29.             if(k_prime >= 3, \\ Модуль >=8 
30.                 if(a1 % 8 != 1, return(0)); 
31.                 next() 
32.             ); 
33.         , 
34.             \\ Для нечётных простых 
35.             if(kronecker(a1, p) != 1, return(0)) 
36.         ) 
37.     ); 
38.     return(1); 
39. }; 

Принимает два параметра a - проверяемый остаток (может быть отрицательным), m - модуль (должен быть неотрицательным, может быть составным числом).
Возвращает 1 если a - квадратичный вычет по модулю m и возвращает 0 если невычет.
Производится факторизация модуля, так что наверное для очень больших простых или псевдопростых модулей будет небыстро.
Функция написана с помощью ИИ, который как оказалось даже знает, что в pari есть функция valuation() и функция kronecker()

Пример. Может ли запись квадрата целого числа оканчиваться на 161616? Да, может:
? is_quadratic_residue(161616,10^6)
%3 = 1
А на 16161616? Нет, не может:
? is_quadratic_residue(16161616,10^8)
%2 = 0

[Dmitriy40]

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

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

Нет, (псевдо)простота определяется быстро, а вот если модуль будет произведением нескольких (строго больше одного) больших простых, вот тогда беда.

[wrest]

Dmitriy40
Кстати, якобы вот это должно определять является ли a квадратичным вычетом или нет по модулю m
issquare(Mod(a,m))
но это неточно

[Dmitriy40]

Охренеть умный ИИ, такую функцию забабахал, а можно было обойтись одной командой ...

[wrest]

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

Охренеть умный ИИ, такую функцию забабахал, а можно было обойтись одной командой ...

ага ну может есть какие нюансы, я не знаю :) но функция написанная ИИ работает "как положено": факторизует модуль, разбирает тонкости связанные со степенями двойки.
В документации на issquare() написано что модуль не проверяется на простоту. Я это читаю как "на ваш страх и риск".
Но ведь работает. Более того, и собсно сравнение решает. Например, найти квадрат какого числа оканчивается на 161616:
? n=0;if(issquare(Mod(161616,10^6),&n),print(lift(n)),print("No solutions"))
609796
?
А если надо 16161616? Пож-ста:
? n=0;if(issquare(Mod(16161616,10^8),&n),print(lift(n)),print("No solutions"))
No solutions
?
ахренеть просто, да? знают ли об этом разрабы pari ? что вот так вот просто раз - и определяется квадратичная вычетность. раз - и решилось сравнение.
и не найдёте вы в интернете что оно так работает. я искал...

[Dmitriy40]

Вопрос лишь как из одного корня
? n=0;issquare(Mod(84848484,100^4),&n);n
%1 = Mod(82292022, 100000000)
получить три других: 7292022, 32292022, 57292022.

[wrest]

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

Вопрос лишь как из одного корня

Никак
Тут писать надо.

[wrest]

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

Тут писать надо.

Ну что ж... пришлось написать.

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

Пример (найти такие числа что запись их квадратов оканчивается на 84848484):

? print(sqm(84848484,100^4))
[7292022, 9895478, 15104522, 17707978, 32292022, 34895478, 40104522, 42707978, 57292022, 59895478, 65104522, 67707978, 82292022, 84895478, 90104522, 92707978]
? ##
  ***   last result: cpu time 0 ms, real time 0 ms.
?

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

получить три других:

ну вот получилось ещё 15 других