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

[wrest]

В статье https://www.mathnet.ru/php/getFT.phtml? ... what=fullt

М. Н. Вялый, О представлении чисел в виде суммы двух квадратов, Матем
просв., 2006, выпуск 10, 190-194

Написано:

Цитата:

Как найти представление числа п в виде суммы двух квадратов? Ко-
нечно, можно перебрать все пары $0\le x ,y \le n$ и проверить для каждой
пары равенство $x^2+y^2=n$. Нельзя ли придумать алгоритм, который
решает эту задачу быстрее?
Эффективных алгоритмов для представления числа в виде суммы двух
квадратов пока нет. Но задача значительно упрощается, если известен
корень из $-1$ по модулю $n$. В этом случае можно обойтись без перебора.

[zykov]

Да, по той ссылке они начинают с поиска корня из -1 (по модулю).
И говорят, что это легко сделать.

-- 27.05.2023, 19:46 --

Например по их формуле $18572338448117274863^2=-1 \mod 19253521174721937977$
Дальше надо искать GCD через алгоритм Евклида для $19253521174721937977$ и $18572338448117274863+i$.
Это и даст два числа, сумма квадратов которых даёт $19253521174721937977$.

[zykov]

Там у них кстати есть готовый код на Питоне.
Попробовал его, он сразу мгновенно разложил $19253521174721937977$ в $3098660171$ и $3106738856$.

Скопирую сюда их код, чтобы не потерялся:

Используется синтаксис Python

def mods(a, n):
    if n <= 0:
        return "negative modulus"
    a = a % n
    if (2 * a > n):
        a -= n
    return a

def powmods(a, r, n):
    out = 1
    while r > 0:
        if (r % 2) == 1:
            r -= 1
            out = mods(out * a, n)
        r //= 2
        a = mods(a * a, n)
    return out

def quos(a, n):
    if n <= 0:
        return "negative modulus"
    return (a - mods(a, n))//n

def grem(w, z):
    # remainder in Gaussian integers when dividing w by z
    (w0, w1) = w
    (z0, z1) = z
    n = z0 * z0 + z1 * z1
    if n == 0:
        return "division by zero"
    u0 = quos(w0 * z0 + w1 * z1, n)
    u1 = quos(w1 * z0 - w0 * z1, n)
    return(w0 - z0 * u0 + z1 * u1,
           w1 - z0 * u1 - z1 * u0)

def ggcd(w, z):
    while z != (0,0):
        w, z = z, grem(w, z)
    return w

def root4(p):
    # 4th root of 1 modulo p
    if p <= 1:
        return "too small"
    if (p % 4) != 1:
        return "not congruent to 1"
    k = p//4
    j = 2
    while True:
        a = powmods(j, k, p)
        b = mods(a * a, p)
        if b == -1:
            return a
        if b != 1:
            return "not prime"
        j += 1

def sq2(p):
    a = root4(p)
    return ggcd((p,0),(a,1))

[iifat]

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

И говорят, что это легко сделать

Ну да, вполне себе легко. Берём случайное $a$, считаем $a^\frac{p-1}2\pmod p$, пока не получим -1 — и вот он, корень.

[wrest]

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

Там у них кстати есть готовый код на Питоне.

Так раскладываются любые $n$ или только простые вида $4k+1$?
И как получить все разложения?
Допустим, у нас $n$ разложилось на множители так:
$n=2^{a_0}p_1^{2a_1}\dots p_i^{2a_i}q_1^{b_1}\dots q_s^{b_s}$
($p_i$ и $q_s$ простые).
Если $n$ так не представимо, то оно не раскладывается в сумму квадратов, так что его сразу откидываем.
Допустим мы разложили все $q_s$ в суммы квадратов ($q_s=x_s^2+y_s^2$).
Кстати $x_s$ и $y_s$ же необязател но уникальные будут.
Ну и что дальше? Что-то ещё надо проверять или можно выписать все представления?

-- 27.05.2023, 21:52 --

Почему Вялый пишет что

Вялый в сообщении #1595569 писал(а):

Эффективных алгоритмов для представления числа в виде суммы двух
квадратов пока нет. Но задача значительно упрощается, если известен
корень из $-1$ по модулю $n$.

Как мы видим, корень этот ищется влёт, так ведь?

[mathematician123]

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

Как мы видим, корень этот ищется влёт, так ведь?

Только для простых $n$. А для составных $n$ даже определить, разрешимо ли сравнение $x^2 = -1 \mod{n}$ выглядит сложной задачей.

Например, пусть $n = pq$ для больших простых $p, q$. Если $p \equiv q \equiv 3 \pmod{4}$, то сравнение не имеет решений, а если $p \equiv q \equiv 1 \pmod{4}$, то решения есть. Как различить эти два случая без факторизации $n$? Непонятно.

[wrest]

mathematician123 в сообщении #1595587 писал(а):

Как различить эти два случая без факторизации $n$? Непонятно.

Факторизация не слишком больших чисел (ну не знаю, 64-битных например) - весьма быстрый процесс в сравнении с перебором при поиске квадратов. То есть, считаем что $n$ факторизовали, и уже знаем количество решений нашей задачи.

[zykov]

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

Так раскладываются любые $n$ или только простые вида $4k+1$?

Этот Питон-скрипт для простых вида $4k+1$ (простые вида $4k+3$ вообще не раскладываются).

Если число удалось разложить на простые и оно содержит только степени двойки и простые вида $4k+1$, то каждое такое простое раскладываем на сумму квадратов и потом получаем разложение для произведения по формуле.

Если разложение содержит простые вида $4k+3$ в чётной степени, то так не сработает.
Но может легко разложить на сумму квадаратов число вида $(4k+3)^{2n}$?

[mathematician123]

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

Ну и что дальше? Что-то ещё надо проверять или можно выписать все представления?

Можно
http://kvant.mccme.ru/pdf/1999/03/kv0399senderov.pdf
Вот здесь на последних двух страницах написано.

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

Если разложение содержит простые вида $4k+3$ в чётной степени, то так не сработает.
Но может легко разложить на сумму квадаратов число вида $(4k+3)^{2n}$?

Если $np^2 = a^2 + b^2$ для простого $p \equiv 3 \pmod{4}$, то $p \mid a$ и $p \mid b$. И задача сводится к представлению $n$ в виде суммы квадратов.

[wrest]

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

Но может легко разложить на сумму квадаратов число вида $(4k+3)^{2n}$?

Дык такие простые в сумму квадратов не раскладываются же

[zykov]

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

Дык такие простые в сумму квадратов не раскладываются же

А их чётные степени?
(Нечётны степени не раскладываются по Теореме Ферма — Эйлера).

-- 27.05.2023, 23:21 --

mathematician123
Спасибо, пронятно.
Если в разложении есть чётные степени простых, то их и раскладывать не надо.
Просто разложим остальное, а потом домножим.

[Dmitriy40]

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

for(x=floor(sqrt(n/2))+1,

Это неправильно, "+1" лишний:

? sqw(50)
2 solutions:
50=7^2+1^2
? sq(50)
50=5^2+5^2
50=7^2+1^2

Вместо floor(sqrt()) можно использовать sqrtint().
Чтобы не заводить лишнюю c0 можно делать так: c--;if(c==0,return).

В sqn() объявленная v1 не используется.

-- 28.05.2023, 00:25 --

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

Например по их формуле $18572338448117274863^2=-1 \mod 19253521174721937977$

Результат может быть не единственным:

? x=Mod(-1,19253521174721937977)
Mod(19253521174721937976, 19253521174721937977)
? sqrt(x)
Mod(681182726604663114, 19253521174721937977)
? (681182726604663114^2)%19253521174721937977-19253521174721937977
-1

Т.е. число $681182726604663114$ тоже подходит: $681182726604663114^2=-1 \mod 19253521174721937977$.

[wrest]

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

Это неправильно, "+1" лишний:

Да, спасибо. Я заметил, у себя использую пока round(sqrt(n/2))

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

Чтобы не заводить лишнюю c0 можно делать так: c--;if(c==0,return).

В sqn() объявленная v1 не используется.

Согласен, спасибо.

P.S. Оказалось, что количество разбиений на два квадрата есть в A025426

[zykov]

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

Т.е. число $681182726604663114$ тоже подходит: $681182726604663114^2=-1 \mod 19253521174721937977$.

$681182726604663114=19253521174721937977 - 18572338448117274863$, т.е. $681182726604663114=-18572338448117274863 \mod 19253521174721937977$.
Очевидно, что в квадрате оно даст тоже самое.

[Dmitriy40]

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

Очевидно, что в квадрате оно даст тоже самое.

Согласен, недоглядел.