Разложение числа в сумму двух квадратов

[wrest]

zykov в сообщении #1595570 писал(а):

Дальше надо искать GCD через алгоритм Евклида для $19253521174721937977$ и $18572338448117274863+i$.
Это и даст два числа, сумма квадратов которых даёт $19253521174721937977$.

Победа всё ближе. Действительно:

? x=gcd(19253521174721937977,18572338448117274863+I)
%1 = 3098660171 + 3106738856*I
? real(x)^2+imag(x)^2
%2 = 19253521174721937977

[wrest]

Пока не забылось...
Разложим простое число p вида 4k+1 в сумму двух ненулевых квадратов. Предполагается, что проверки на простоту и остаток от деления на 4 выполнены ранее. Возвращается вектор из двух чисел, сумма квадратов которых равна p.

? sqp(p)=my(x=lift(p-sqrt(Mod(p-1,p))),y=gcd(p,x+I));return([real(y),imag(y)]);
? print(sqp(19253521174721937977))
[3098660171, 3106738856]
?

[dmd]

В pari/gp имеется команда thue, которая может решить сумму двух квадратов:

? thue(x^2+1,19253521174721937977)
%1 = [[-3106738856, -3098660171], [-3106738856, 3098660171], [-3098660171, -3106738856], [-3098660171, 3106738856], [3098660171, -3106738856], [3098660171, 3106738856], [3106738856, -3098660171], [3106738856, 3098660171]]
?

[Andrey A]

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

Победа всё ближе.

Если простое $p \equiv 1 \pmod 4,$ период разложения $\sqrt{p}$ нечетен, и в конце первого полупериода встретится подпоследовательность знаков вида $u,v,v,u.$ Назовем "интересной" подходящую дробь $x/y$ , соответствующую первому знаку этого мини-палиндрома. Она интересна тем, что остаток $x^2-py^2$ есть основание одного из двух искомых квадратов. Гипотеза. Подробности выложу в параллельной теме.

$\sqrt{61}=7,1,4,3,\underline{1,2,2,1},...$

$7,1,4,3,1,=\dfrac{7}{1},\dfrac{8}{1},\dfrac{39}{5},\dfrac{125}{16},\dfrac{164}{21},...\ 164^2-61 \cdot 21^2=-5.\ \sqrt{61-5^2}=6.$

[wrest]

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

В pari/gp имеется команда thue, которая может решить сумму двух квадратов:

Ничосе!
Тогда функция sqt(n) возвращает вектор из двухкомпонентных векторов (все натуральные решения, отсортированные по возрастанию, если есть, и пустой если нет):

? sqt(n)=my(s=Set(apply(vecsort,abs(thue(x^2+1,n)))),l=List(s));if(#l>0 && l[1][1]==0,listpop(l,1));return(Vec(l));

Запуск:

? print(sqt(9907026323426308))
[[12345678, 98765432], [23606018, 96694272], [30889448, 94619598], [41551008, 90446338]]
? print(sqt(489886592769624989886788))
[[179075995118, 676622775808], [180837148448, 676154212078]]
? ##
  ***   last result: cpu time 0 ms, real time 1 ms.

1(одна!) миллисеунда на разбиение...

P.S. Львиная часть функции (в смысле текста, не времени) занимается исключением перестановок, нулевых и отрицательных решений.
P.P.S. В некотором смысле конечно читерство, но раз уж есть это в PARI/GP то почему бы и не воспользоваться.

[dmd]

Еще есть bnf-команды:

? bnfisnorm(bnfinit(x^2 + 1, 1), 19253521174721937977)
%1 = [Mod(3106738856*x + 3098660171, x^2 + 1), 1]
?

Внутри thue думаю bnf-код.

[wrest]

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

Внутри thue думаю bnf-код.

Фкнция thue() сама перечисляет все решения и выдаёт пустой вектор если их нет. Как раз то, что надо.
Так-то оно уже понятно: факторизуем, затем проверяем нет ли в факторизации нечётных степеней вида $4k+3$, затем раскладываем простые множители вида $4k+1$ в суммы квадратов и комбинируем, потом ещё смотрим на степень двойки, и учитываем её. запрограммировать-то несложно, просто немножко муторно. Всё что нужно, после разложения простых $4k+1$ в сумму квадратов, это

zykov в сообщении #1595553 писал(а):

Если $x^2 + y^2 = n$, то $(x+y)^2 + (x-y)^2 = 2n$.
Значит степени двойки можно сразу выкинуть и раскладывать нечётное.

-- 27.05.2023, 18:09 --

Ещё теорема Фиббоначи: Если целые числа N и M могут быть представлены в виде суммы двух квадратов, то их произведение так же может быть представлено в виде суммы двух квадратов.

$$(a^2+b^2) (c^2+d^2) = (ac+bd)^2 + (ad-bc)^2$$

Если удаётся разложить число на множители, то можно каждый множитель разложить на сумму.
Тут правда возможно некоторые пары будут потеряны.

Ничего там потеряно не будет, всё сойдется. Ессно, мы для произведения сумм квадратов запишем оба варианта.

$$(a^2+b^2) (c^2+d^2) = (ac+bd)^2 + (ad-bc)^2$$

и

$$(a^2+b^2) (c^2+d^2) = (ad+bc)^2 + (ac-bd)^2$$

Выкидывать степени двойки надо канеш не огульно, а тоже посмотреть четная она или нет. Четную часть переносим к остальным квадратам $4k+3$, и у нас остаются или только $4k+1$ или они же умноженные на два.

Т.е. приводим факторизацию

$$n=2^{a_0}p_1^{2a_1}\dots p_i^{2a_i}q_1^{b_1}\dots q_s^{b_s}$$

где ($p_i$ простые вида $4k+3$, а $q_s$ простые вида $4k+1$) к виду

$$n=2^k\cdot (2^m\cdot {p_1^{a_1}\dots p_i^{a_i})^2 \cdot (x_1^2+y_1^2)^{b_1}\dots (x_s^2+y_s^2)^{b_s}$$

где $k$ остаток от деления $a0$ на $2$, а $m$ - целая часть от деления $a0$ на $2$

Я уже так и собирался сделать... Но подоспела функция thue()

[genk]

Вариация по мотивам постов о разложении в сумму квадратов(Dmitriy40).
Программка для разложения натурального n в сумму двух натуральных k-х степеней.

sq(n,k)=my(x,y,f);f=0;
for(x=1,(n/2)^(1/k),y=n-x^k;
if(sqrtnint(y,k)==y^(1/k),f++;print(n,"=",x,"^",k,"+",y^(1/k),"^",k)));
if(f==0,print("No solutions"))

Например

sq(243,3)
243=3^3+6^3
sq(41405,2)
41415=14^2+203^2
41415=91^2+182^2
41415=133^2+154^2
sq(145141577656337,4)
145141577656337=1234^4+3457^4
sq(1234567890,5)
No solutions

Может,кому-нибудь пригодится.Ввиду общности программы,для шибко больших чисел работает медленно.

[wrest]

genk в сообщении #1596151 писал(а):

Программка для разложения натурального n в сумму двух натуральных k-х степеней.

...

genk в сообщении #1596151 писал(а):

Ввиду общности программы,для шибко больших чисел работает медленно.

Так а зачем, если есть thue(), которая работает мгновенно
? spt(n,k)=my(s=Set(apply(vecsort,abs(thue(x^k+1,n)))),l=List(s));if(#l>0 && l[1][1]==0,listpop(l,1));return(Vec(l));
? spt(145141577656337,4)
%3 = [[1234, 3457]]
? spt(1234567890,5)
%4 = []

[wrest]

genk в сообщении #1596151 писал(а):

if(sqrtnint(y,k)==y^(1/k),f++;print(n,"=",x,"^",k,"+",y^(1/k),"^",k)));

На заметку: то же самое (но без потенциальных ошибок округления) будет if(ispower(y,k,&y),f++;print(n,"=",x,"^",k,"+",y,"^",k)));
Экономится по крайней мере одно возведение в нецелую степень, что для цикла будет заметно, плюс исключается ошибка округления.
А вот тут

genk в сообщении #1596151 писал(а):

for(x=1,(n/2)^(1/k),

Как раз лучше будет for(x=1,sqrtint(n/2,k), - особо не влияет, но выглядит более, так сказать, целевым образм (понятней текст).

[genk]

wrest
Про thue знаю,читал,но особенно не вникал и не понравился формат вывода,поэтому для повседневной работы с не очень большими числами решил составить свой вариант.Про ispower как-то выпало из памяти,thanks.Кстати,почему вы все время(я имею ввиду разные посты) употребляете оператор $\&\&$ -у меня прекрасно работает $\&$,т.е.унарно,без повтора.

[zykov]

Цитата:

Finally, the logical operators ! (not prefix operator), && (and operator), || (or operator) exist, giving as results 1 (true) or 0 (false). Note that & and | are also accepted as synonyms respectively for && and ||.

[wrest]

genk в сообщении #1596298 писал(а):

и не понравился формат вывода

Ну это самое лёгкое, вывести на печать вы можете в любом формате. Тут смысл в том что функциявозвращает как раз в удобном формате вектора векторов, пригодном для дальнейшей обработки, в т.ч. вывода на печать.

genk в сообщении #1596298 писал(а):

Кстати,почему вы все время(я имею ввиду разные посты) употребляете оператор $\&\&$ -у меня прекрасно работает $\&$

Чтобы не было конфуза

[zykov]

В Pari "&&" и "&" синонимы.
В других языках нет. Обычно первый - логическое "ИЛИ", второй - побитовое "ИЛИ".