√m≈√x+√y

[juna]

Реализовал пока алгоритм $x^2\equiv R\mod p$ по простому модулю. Заявленная вычислительная сложность $O(\log^4(p))$. Считает мгновенно на числах большой разрядности.

/* solve x^2=a mod p, p is prime */
modulop(a,p):=block([h, N1, N2,jacob, N,a1,a2, k,p1,j,i, b, c, d, ji],
  if primep(p) then
    block(
             jacob: jacobi (p, a),
        if jacob = -1 then
          "No solutions"
        else
          block(
             h:p-1,
        k:0,
        while mod(h,2)=0 do block(h:h/2,k:k+1),
             N:5,
        while jacobi (p, N)=1 do N:N+1,
        a1:power_mod(a, (h+1)/2,p),
        N1:power_mod(N,h,p),
        N2:1, j:0, a2:inv_mod(a,p),
        for i:0 thru k-2 do
          block(
                    b:mod(a1*N2,p), c:mod(a2*b^2,p),
               d:power_mod(c,2^(k-2-i),p),
               if d = 1 then ji:0,
               if d=p-1 then ji:1,
               N2:mod(N2*power_mod(N1,(2^i*ji),p),p)
                  ),
      [mod(a1*N2,p), p-mod(a1*N2,p)]
                   )
        )
  else
    "Not prime"
    )$
modulop(87,101);

Воспользоваться можно через любой online сервис, например http://maxima.cesga.es/

Вставляем вышеуказанный код, в функцию modulop подставляем нужные a и p (правда в online версии символ Якоби как-то для длинных простых плохо считается, поэтому все-таки лучше оффлайн последняя версия maxima).

-- Ср янв 11, 2023 22:48:28 --

(Оффтоп)

Утундрий в сообщении #1576796 писал(а):

Надо будет на досуге и себе завести птичью тему. Что-нибудь вроде разложения кубического корня свободного от кубов целого числа на сумму трёх квадратных корней свободных от квадратов целых чисел.

И у Вас по заявленному вопросу есть какой-то задел? В противном случае содержательное общение можно продолжить только в пургатории...

[Andrey A]

juna
Я пока свои тупички не обойду, не уймусь, на будущее очень интересно. Но кое что по горячим следам. Положим, нам не известно простое $m$ или составное. Берем произвольный квадрат (или "неудобный" для процедуры, если такое возможно), делим его на $m$, получаем вычет $R$. Запускаем процедуру по $R$ и получаем другой, удобный ей квадрат. Программе ничего не известно о природе $m$, но мы получили пару сравнимых квадратов и можем делать выводы о делимости. Если процедура быстрее известных методов факторизации, нет ли в том противоречия? А что если $R$ — целый квадрат. Влияет это на ее поведение?
На alpertron.com помните как было — программа делает массу только ей нужных вещей, все по уму и очень солидно, а в самом корне сидит маленький такой переборчик.

[juna]

Andrey A в сообщении #1576830 писал(а):

Если процедура быстрее известных методов факторизации, нет ли в том противоречия?

Конечно, конкретно эта процедура не быстрее алгоритма проверки на простоту (для больших $p$ там, по-моему, вероятностный тест Миллера-Рабина запускается). В функции modulop первой же строчкой идет:

if primep(p) then

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

Andrey A в сообщении #1576830 писал(а):

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

Конкретно здесь переборчика всего 3, не считая всяких внутренних power_mod (модульное возведение в степень), inv_mod(мультипликативно обратное):

      h:p-1,
        k:0,
        while mod(h,2)=0 do block(h:h/2,k:k+1),

Находим степень двойки в p-1.

N:5,
        while jacobi (p, N)=1 do N:N+1,

Находим первый попавшийся квадратичный невычет.

for i:0 thru k-2 do

Цикл по количеству найденных ранее 2 в p-1. Отсюда, в частности, понятно, что худшие случаи для алгоритма - это простые числа вида $p=2^n+1$

[Andrey A]

Ага. Тест на простоту как базовая операция. И квадраты не только по $R$, но и по всем нужным вычетам (двойку-то зря что-ли в $ (p-1)$ возводили?). Отлично. Осталось подставить всё это в известные нам формулы, упорядочить по величине погрешности, и... вот оно решение. Которое всё равно захочется проверить перебором, мало ли... ) Ведь захочется? По двойке тоже есть свои квазипростые. Знаете что, мне Ваши исследования интересны, по ссылкам буду еще разбираться, но если честно — много ли радости в таком решении? Вроде как древнюю перебор-пушку меняем на современные ПЗРК. Всё по тем же воробьям. Нет. Всё-таки хочется прямой алгоритм.

[juna]

Улучшил предыдущий код:

/* solve x^2=a mod p, p is prime */
modulop(a,p):=block([h, N1, N2, legandr, N,a1,a2, k,p1,i, b, c, d, j],
  if primep(p) then
    block(
             legandr: power_mod(a,(p-1)/2,p),
        if legandr = p-1 then
          "No solutions"
        else
          block(
             h:p-1, k:0,
        while mod(h,2)=0 do block(h:h/2,k:k+1),
             N:2,
        while power_mod(N,(p-1)/2,p)=1 do N:N+1,
        a1:power_mod(a, (h+1)/2,p), N1:power_mod(N,h,p),
        N2:1, j:0, a2:inv_mod(a,p),
        for i:0 thru k-2 do
          block(
                    b:mod(a1*N2,p), c:mod(a2*b^2,p),
               d:power_mod(c,2^(k-2-i),p),
               if d = 1 then j:0,
               if d=p-1 then j:1,
               N2:mod(N2*power_mod(N1,(2^i*j),p),p)),
           j:mod(a1*N2,p),
      [j, p-j]))
  else
    "Not prime")$

Реализовал нахождение решения сравнения $v^2\equiv R\mod m$ для $m=p_1\cdot p_2$:

/* solve x^2=a mod m, m = p1 * p2 */
modulop1p2(a,m):= block([fact:ifactors(m), ans1, ans2,res:[]],
  if length(fact) # 2 then
   "m is not equal p1*p2"
  else
    block(
      if (fact[1][2] = 1) and (fact[2][2] = 1) then
      block(
        ans1:modulop(a,fact[1][1]),
        ans2:modulop(a,fact[2][1]),
        res:cons(chinese([ans1[1],ans2[1]],[fact[1][1], fact[2][1]]), res),
        res:cons(chinese([ans1[1],ans2[2]],[fact[1][1], fact[2][1]]), res),
        res:cons(chinese([ans1[2],ans2[1]],[fact[1][1], fact[2][1]]), res),
        res:cons(chinese([ans1[2],ans2[2]],[fact[1][1], fact[2][1]]), res),
        res)
      else
        "m is equal p1^a*p2^b, a or b > 1" ))$

Для общего вида $m$ больно муторно перебирать все $2^n$ комбинаций, где $n$ - число простых делителей $m$.

Для эффективного перебора нужно еще разработать правила останова для $R$.

Пусть известно решение $v_1^2\equiv R_1\mod m$, найдем ограничение для $R_2, v_2$, для которых возможно лучшее приближение.

Имеем:

$$2m-\sqrt{(m-v_1)^2-R_1}-\sqrt{(m+v_1)^2-R_1}>2m-\sqrt{(m-v_2)^2-R_2}-\sqrt{(m+v_2)^2-R_2}$$

т.е. должно выполняться:

$$\sqrt{(m-v_1)^2-R_1}+\sqrt{(m+v_1)^2-R_1}<\sqrt{(m-v_2)^2-R_2}+\sqrt{(m+v_2)^2-R_2}$$

Обозначим $a=\sqrt{(m-v_1)^2-R_1}+\sqrt{(m+v_1)^2-R_1}, b=(m-v_2)^2, c=(m+v_2)^2$

Возводим в квадрат:

$$a^2<b-R_2+c-R_2+2\sqrt{(b-R_2)(c-R_2)}$$

$$a^2-b-c+2R_2<2\sqrt{(b-R_2)(c-R_2)}$$

Возводим еще раз в квадрат:

$$(a^2-b-c+2R_2)^2<4(b-R_2)(c-R_2)$$

$$(a^2-b-c)^2+4R_2(a^2-b-c)+4R_2^2<4bc-4bR_2-4cR_2+4R_2^2$$

$$(a^2-b-c)^2+4R_2a^2<4bc$$

$$R_2<\frac{4bc-(a^2-b-c)^2}{4a^2}$$

После всех подстановок и сокращений окончательно получаем (если не ошибся):

$$R_2<m^2-\frac{a^2}{4}+v_2^2-\frac{4m^2v_2^2}{a^2}$$

Поскольку $4m^2>a^2$, то функция $f(v_2)=v_2^2(1-\frac{4m^2}{a^2})$ достигает максимума в точке $v_2=0$.

Поэтому имеем следующее ограничение:

$$R_2<m^2-\frac{a^2}{4}$$

Но это не особенно полезно.

$$v_2<\sqrt{\frac{m^2-\frac{a^2}{4}-R_2}{\frac{4m^2}{a^2}-1}}$$

А вот на это можно опираться. Для этого добавим еще порцию кода:

delta(m,R,L):=block([dd,v,vv,i:1,j],
  while not(integerp(xy(m,L[i],R)[1])) and (i<length(L)) do i:i+1,
  dd:bfloat(2*m-sqrt((m-L[i])^2-R)-sqrt((m+L[i])^2-R)),vv:L[i],
  for j:i+1 thru length(L) do
    block(
    if (dd>bfloat(2*m-sqrt((m-L[j])^2-R)-sqrt((m+L[j])^2-R))) and (integerp(xy(m,L[j],R)[1])) then
      block(dd:bfloat(2*m-sqrt((m-L[j])^2-R)-sqrt((m+L[j])^2-R)),vv:L[j])),
  if integerp(xy(m,vv,R)[1]) then
  [vv,dd]
  else
  false)$
xy(m,v,R):=([((m-v)^2-R)/(4*m),((m+v)^2-R)/(4*m)])$
fpprec:30$
v2(m,v1,R1,R2):=block([a], a:bfloat(sqrt((m-v1)^2-R1)+sqrt((m+v1)^2-R1)),bfloat(sqrt((m^2-a^2/4-R2)/(4*m^2/a^2-1))))$

Функция delta(m,R,L) выбирает из списка L решений сравнения $v^2\equiv R\mod m$ целочисленное решение с минимальной ошибкой $2m-a$, если целочисленных решений нет, то возвращается false.

Функция xy(m,v,R) просто получает список решений $x,y$ для заданных $m,v,R$.

Функция v2(m,v1,R1,R2) считает собственно верхнее ограничение для $v_2$ по полученной формуле для известных $v_1,R_1$ и нового $R_2$.

Продемонстрируем на примере треугольного числа (максимального в указанной ранее таблице) $m=714195944281=597577\cdot 1195153$

(%i7) delta(714195944281,1, modulop1p2(1, 714195944281));
(%o7)         [714193553974, 2.09178501637662800360928372356b-7]
(%i8) v2(714195944281, 714193553974, 1, 2);
(%o8)                 7.14191163658994066022843741247b11
(%i9) delta(714195944281,2, modulop1p2(2, 714195944281));
(%o9)                                false
(%i10) v2(714195944281, 714193553974, 1, 3);
(%o10)                7.14188773335993878063157787733b11
(%i11) delta(714195944281,3, modulop1p2(3, 714195944281));
(%o11)                               false
(%i12) v2(714195944281, 714193553974, 1, 4);
(%o12)                7.14186383004993475900066029253b11
(%i13) delta(714195944281,4, modulop1p2(4, 714195944281));
(%o13)        [714191163667, 4.18357703366352418061779872005b-7]
(%i14) v2(714195944281, 714191163667, 4, 5);
(%o14)                7.14189968508510919964289242979b11
(%i15) delta(714195944281,5, modulop1p2(5, 714195944281));
(%o15)                               false
(%i16) v2(714195944281, 714191163667, 4, 6);
(%o16)                7.14188773348010840036651324299b11
(%i17) delta(714195944281,6, modulop1p2(6, 714195944281));
(%o17)                               false
(%i18) v2(714195944281, 714191163667, 4, 7);
(%o18)                7.14187578185510713252166124388b11
(%i19) delta(714195944281,7, modulop1p2(7, 714195944281));
(%o19)                               false
(%i20) v2(714195944281, 714191163667, 4, 8);
(%o20)                7.14186383021010529569837096844b11
(%i21) delta(714195944281,8, modulop1p2(8, 714195944281));
(%o21)        [214899754753, 1.23165381973183718500308714283b-11]
(%i22) v2(714195944281, 214899754753, 8, 9);
(%o22)               1.08650074855601945444046486567b11 %i

Итак, на $R=8$ уперлись. Имеем в порядке убывания точности:

$$\sqrt{714195944281}\approx\sqrt{87264806974}+\sqrt{302164561727}, R=8$$

$$\sqrt{714195944281}\approx\sqrt{2}+\sqrt{714193553976}, R=1$$

$$\sqrt{714195944281}\approx\sqrt{8}+\sqrt{714191163675}, R=4$$

(%i33) bfloat(sqrt(714195944281));
(%o33)                 8.45101144408762471059270771886b5
(%i34) bfloat(sqrt(87264806974)+sqrt(302164561727));
(%o34)                 8.45101144408762471059263484866b5
(%i35) bfloat(sqrt(2)+sqrt(714193553976));
(%o35)                 8.45101144408762470935511324262b5
(%i36) bfloat(sqrt(8)+sqrt(714191163675));
(%o36)                 8.45101144408762470811751462432b5

-- Вс янв 15, 2023 17:12:19 --

Andrey A в сообщении #1576869 писал(а):

много ли радости в таком решении

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

На алгоритм, который я реализовал в виде функции modulop, оказывается указывал еще waxtep:

waxtep в сообщении #1570810 писал(а):

Тогда, для простого $m$, можно попробовать перебирать значения $R$ по возрастанию (причем, не все $R$ будут разрешены по соображениям взаимности) и в лоб извлекать из них корни по модулю $m$, - здесь же есть быстрые алгоритмы (не имел реального опыта реализации, но вот пишут: Алгоритм Тонелли—Шенкса
; есть некоторые оговорки, но, на первый взгляд, алгоритм не выглядит устрашающим/безнадежным).

https://ru.wikipedia.org/wiki/%D0%90%D0 ... 1%81%D0%B0

[Andrey A]

juna
Я насчет "против лома нет приема" был неправ, конечно. Есть методы, и у Давенпорта кое-что описано, но они либо не элементарны, либо предполагают знание квадрата по единице и т.д. Боюсь, мы тут ставим телегу впереди лошади, и задача в таком виде теряет привлекательность. У меня тут еще одна химера рухнула, не стал и выкладывать. Знаете, я, пожалуй, возьму паузу. Надо успокоиться, заодно и в Ваших выкладках поразбераюсь. От такого беспрерывного бодания добра не жди, иногда нужно отвлечься или просто отдохнуть.

[juna]

Andrey A в сообщении #1577225 писал(а):

От такого беспрерывного бодания добра не жди, иногда нужно отвлечься или просто отдохнуть.

На мой взгляд имело место обычное критическое осмысление фактов...
Кстати, думаю, для простых $m$, аналог подхода Тонелли—Шенкса может существовать и в интерпретации цепных дробей, но уровень требуемого инсайта там наверное потребует 10 дана.

[juna]

И все-таки оценка

$$R_2(m,v_1, R_1)\leq\left\lfloor m^2-\frac{\left(\sqrt{(m-v_1)^2-R_1}+\sqrt{(m+v_1)^2-R_1}\right)^2}{4}\right\rfloor$$

довольна удобна для нахождения новых ограничений на $R_2$ и можно пользоваться только ей.
Продемонстрирую на паре примеров треугольных чисел (хотя она справедлива для любых чисел).

$m=12403$

Мы начинаем с $R_1=1$ находим $v_1=12088, \delta=1.6e-3$, тогда $R_2(12403,12088,1)\leq 19$.

-- Чт янв 19, 2023 08:45:50 --

Берем $R_1=2$, нет решений
Берем $R_1=3$, нет решений
Берем $R_1=4$, находим $v_1=11773, \delta=3.26e-3$, поскольку погрешность больше, оставляем старое ограничение.
Берем $R_1=5$, нет решений
Берем $R_1=6$, нет решений
Берем $R_1=7$, нет решений
Берем $R_1=8$, нет решений
Берем $R_1=9$, находим $v_1=11458, \delta=4.95e-3$, поскольку погрешность больше, оставляем старое ограничение.
Берем $R_1=10$, нет решений
Берем $R_1=11$, нет решений
Берем $R_1=12$, нет решений

-- Чт янв 19, 2023 08:52:29 --

Берем $R_1=13$, находим $v_1=4060, \delta=1.17e-3$, поскольку погрешность оказалась меньше, обновляем ограничение $R_2(12403, 4060, 13)\leq 14$
Берем $R_1=14$, нет решений
Таким образом, процесс закончен на $R=13$.

$m=1891$

Опять начинаем с $R_1=1$ находим $v_1=1768, \delta=4.2e-3$, тогда $R_2(1891,1768,1)\leq 7$.
Берем $R_1=2$, нет решений
Берем $R_1=3$, нет решений
Берем $R_1=4$, находим $v_1=1645, \delta=8.69e-3$, поскольку погрешность больше, оставляем старое ограничение
Берем $R_1=5$, находим $v_1=1246, \delta=4.67e-3$, поскольку погрешность больше, оставляем старое ограничение
Берем $R_1=6$, нет решений
Берем $R_1=7$, нет решений
Таким образом, остается $R=1$.

[Andrey A]

juna
Хорошо, да больно хлопотно. И зачем проверять всё подряд? Вы же сами открыли чудесную последовательность A003658.

[juna]

Andrey A в сообщении #1577885 писал(а):

И зачем проверять всё подряд? Вы же сами открыли чудесную последовательность A003658
.

Программе она неизвестна (а заморачиваться, как-то ее генерить - ленно), да и если ставить вопрос не только о лучшем приближении, нужно рассматривать все $R$.

[Andrey A]

juna в сообщении #1577889 писал(а):

нужно рассматривать все $R$.

Лучшее — всегда приведенное. Последовательность свободных от квадратов программе известна? Просто домножить ее на $4$ кроме $4k+1$-ых членов.

[juna]

Andrey A в сообщении #1577895 писал(а):

Последовательность свободных от квадратов программе известна?

Нет, конечно. Для реализации нужно факторизацией $R$ заниматься, маленьких правда, но кто знает.

Я вот думаю о попеременном уточнении $R$ через $v$ и обратно.

[juna]

Проверил все треугольные, представимые произведением двух простых в пределах доступных нам таблиц
Максимальное $m$ с $R=1$ - это $m=25651$, дальше уже действительно нет ни одного.
Максимальное $R=92$ для $m=25239989503$.

Вот код программки с учетом разработанного ранее:

sqrtmm(m):=block([Rm,dm,vm,R1,v1,d,RR2, dd,xx],
   R1:1,dd:delta(m,R1,modulop1p2(R1,m)),
   d:dd[2],v1:dd[1],
   Rm:R1, dm:d,
   vm:v1,RR2:R2(m,v1,R1),
   while R1<RR2 do 
     block(
       R1:R1+1,
       dd:delta(m,R1,modulop1p2(R1,m)),
       if dd # false then
         block(
         d:dd[2],v1:dd[1],
         if d<dm then
           block(
             Rm:R1, dm:d, vm:v1,RR2:R2(m,v1,R1)))),
   xx:xy(m,vm,Rm),
   [m,dm,Rm,vm,xx[1],xx[2]])$