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

[genk]

В сообщении https://dxdy.ru/topic14229-345.html wrest писал:
Для разложения числа в сумму двух квадратов натурпльных чисел, насколько мне известно, простых путей нет.
Я тут навскидку составил программку на PARI(наверняка не оптимальную,но Dmitriy40 сразу увидит)

sq(n)={my(v=apply(x->x^2,vector(sqrtint(n-1),i,i)));u=[];f=0;
for(i=1,n-1,v=concat(v,v[i]));
forsubset([#v,2],x,u=concat(u,[vecsort(vecextract(v,x))]));u=Set(u);
for(i=1,#u,if(vecsum(u[i])==n,f++;print(n,"=",sqrtint(u[i][1]),"^",2,"+",sqrtint(u[i][2]),"^",2)));
if(f==0,print("Not squares")}

Например

sq(25)
25=3^2+4^2
sq(32)
32=4^2+4^2
sq(82)
82=1^2+9^2
sq(185)
185=4^2+13^2
185=8^2+11^2
sq(91)
Not squares

Главная головная боль-в создании массива из сочетаний.А как иначе?

[Dmitriy40]

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

А как иначе?

Ну ... Я бы сделал так:

? sq(n)=my(x,y); for(x=1,sqrt(n/2), if(issquare(n-x^2,&y), print(n,"=",x,"^2+",y,"^2")));
? sq(91)
? sq(185)
185=4^2+13^2
185=8^2+11^2
? sq(1234567890123456)
time = 6,491 ms.
? sq(98765432^2+12345678^2)
9907026323426308=12345678^2+98765432^2
9907026323426308=23606018^2+96694272^2
9907026323426308=30889448^2+94619598^2
9907026323426308=41551008^2+90446338^2
time = 18,706 ms.

Время $O(\sqrt{n})$, памяти не требуется вообще.

[genk]

Dmitriy40
Класс! Ex ungue leonem pingere.

[wrest]

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

Ну ... Я бы сделал так:

? sq(n)=my(x,y); for(x=1,sqrt(n/2), if(issquare(n-x^2,&y), print(n,"=",x,"^2+",y,"^2")));

Я бы добавил проверку на наличие решений

sq(n)=my(x,y); if(setsearch(Set(factor(n)[, 1]%4), 1), print("Solutions:");for(x=1,sqrt(n/2), if(issquare(n-x^2,&y), print(n,"=",x,"^2+",y,"^2"))),print("No solutions"));

Тогда перебор будет только если решения есть. А если решений нет, будет быстрее.

[wrest]

Ой, что-то я накосячил. Предыдущее сообщение недействительно.

-- 27.05.2023, 15:24 --

Попытка номер 2.
Функция sqn() вычисляет количество представлений входного числа суммой двух натуральных квадратов

sqn(c)=my(v=factor(c),a0=if(v[1,1]==2,v[1,2],0),v1=v,b=1;b4=0);for(i=1,#v[,1],b4=v[i,1]%4;if(b4==1,b=b*(v[i,2]+1));if(b4==3,if(v[i,2]%2!=0,return(0))));if(b%2==0,return(b/2));return((b-(-1)^a0)/2)

Функция sqw() обращается к функции sqn(), затем перебирает и печатает решения пока не найдёт их все.

sqw(n)=my(x,y,c=sqn(n),c0=0);if(c==0,print("No solutions");return());print(c," solutions:");for(x=1,sqrt(n/2), if(issquare(n-x^2,&y), print(n,"=",x,"^2+",y,"^2");c0++;if(c0==c,return())));

Надеюсь, косяков нет
Ускорение за счет:
1. Ускоряется если решений нет.
2. Перебор заканчивается когда найдены все решения.

P.S. Я не силён в while и until, так что заколхозил просто выход из for
P.P.S. Идея как считать количество решений взята отсюда: https://mathworld.wolfram.com/SumofSquaresFunction.html

[wrest]

Вольфрам разлагает число 489886592769624989886788 в сумму двух квадратов влёт.
https://www.wolframalpha.com/input?i=Po ... 2C2%2C2%29
Получается 2 решения (хотя именно это число в доках по функции, так может просто из базы достают).
489886592769624989886788=179075995118^2+676622775808^2
489886592769624989886788=180837148448^2+676154212078^2

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

Функция sqn() из поста выше, конечно, вычисляет количество решений мгновенно. Так что вопрос именно в поиске их.

[zykov]

Как представить произвольное число суммой квадратов
Теорема Ферма — Эйлера
Представление чисел суммой двух квадратов и эллиптические кривые

[wrest]

zykov
Ссылки не дают ответа на вопрос как разложить число на сумму двух квадратов натуральных чисел...
Например, берём "пробное" число
$489886592769624989886788=2^2\cdot 6361^1 \cdot 19253521174721937977^1$
Остатки от деления на 4 простых множителей
$2;1;1$
И что дальше? Знаем, что искомых представлений два штуки.

[zykov]

Если вместо перебора "x=1,sqrt(n/2)" делать перебор "x=sqrt(n/2),sqrt(n)", то будет несколько короче (где-то в $2.4$ раза).

[wrest]

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

Если вместо перебора "x=1,sqrt(n/2)" делать перебор "x=sqrt(n/2),sqrt(n)", то будет несколько короче (где-то в $2.4$ раза).

Хорошее замечание, действительно.
Переопубликую с учетом этого (а то у меня там помарка нашлась)

sqn(c)=my(v=factor(c),a0=if(v[1,1]==2,v[1,2],0),v1=v,b=1,b4=0);for(i=1+(a0>1),#v[,1],b4=v[i,1]%4;if(b4==1,b=b*(v[i,2]+1));if(b4==3,if(v[i,2]%2!=0,return(0))));if(b%2==0,return(b/2));return((b-(-1)^a0)/2)

sqw(n)=my(x,y,c=sqn(n),c0=0);if(c==0,print("No solutions");return());print(c," solutions:");for(x=floor(sqrt(n/2))+1, sqrt(n),if(issquare(n-x^2,&y), print(n,"=",x,"^2+",y,"^2");c0++;if(c0==c,return())));

[zykov]

Если $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$$

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

[wrest]

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

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

Ну вот в "пробном числе" есть очень большой простой множитель:

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

Например, берём "пробное" число
$489886592769624989886788=2^2\cdot 6361^1 \cdot 19253521174721937977^1$

При этом $2^2\cdot 6361^1=138^2+80^2$
Вольфрам раскладывает
19253521174721937977=3098660171^2+3106738856^2
Но как это сделать нам? Хотя тут нам конечно повезло, ибо оба числа практически равны корню из половины раскладываемого и перебор быстро отрабатывает.

[zykov]

Да, $6361$ легко раскладывается. Единственный вариант $6361=40^2+69^2$.

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

у вот в "пробном числе" есть очень большой простой множитель

Его и надо раскладывать.
Всё же меньше, чем исходное число.

[wrest]

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

Единственный вариант

Ессно, простые если и раскладываются, то единственным образом.

[zykov]

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

Вольфрам раскладывает
19253521174721937977=3098660171^2+3106738856^2

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

-- 27.05.2023, 19:18 --

Вот здесь приводят какой-то алгоритм: https://math.stackexchange.com/question ... to-a-prime