Много точных квадратов

[nnosipov]

dmd
Для полноты картины сообщите, пожалуйста, сколько времени (и в какой CAS --- видимо, pari/gp?) работал Ваш i5, пока не решил этот миллион уравнений.

Как я понял, Ваш алгоритм основан на оценке $|y| \leqslant H$, но здесь, действительно, возникает вопрос о корректности:

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

Выходит что $|y|$ не может быть больше $H$, но почему?

Это можно обосновать буквально голыми руками. А именно, можно рассуждать так: если $|y|>H$, то число $y^4+2Hy+1$ не может быть точным квадратом, потому что его можно поместить строго между двумя последовательными квадратами. И в самом деле: если, например, $y>H$, то имеем $(y^2)^2<y^4+2Hy+1<(y^2+1)^2$ (левое неравенство очевидно, а правое станет очевидным, если раскрыть скобки). Если же $y<-H$, то $(y^2-1)^2<y^4+2Hy+1<(y^2)^2$, что проверяется столь же просто. Таким образом, для любого решения a priori имеем $|y| \leqslant H$, и можно просто перебрать все такие $y$. Этим рассуждениям можно придать общий характер, что сделано, например, в статье Poulakis D. A simple method for solving the diophantine equation $y^2=x^4+ax^3+bx^2+cx+d$ // Elem. Math. 1999. Vol. 54. P. 32-36. Эта статья абсолютно элементарна (доступна школьникам). Алгоритм, который в ней описывается, я привык называть стандартным алгоритмом. Он известен очень давно. Я думал, что он существует на уровне фольклора, пока не наткнулся на статью Пулакиса.

Интрига, однако, в том, что этот алгоритм иногда может быть существенно улучшен. По-настоящему мы это поняли год назад, вот в этой теме: topic136587.html Применительно к нашему уравнению $y^4+2Hy+1=z^2$ этот улучшенный алгоритм выглядит так. Пусть $(y,z)$ --- некоторое решение, при этом $z \geqslant 0$. Рассмотрим целое число

$$ x=-y^2+z=-y^2+\sqrt{y^4+2Hy+1}=\frac{H}{y}+O\left(\frac{1}{y^2}\right), \quad y \to \infty. $$

Тогда

$$ x^2+2y^2x-2Hy-1=0. \eqno(*) $$

Кроме того, если $|y|>c_1H^{1/2}$, то $|x|<c_2H^{1/2}$ (здесь $c_i$ --- абсолютные константы, близкие к единице). Более точно, нетрудно доказывается, что:

а) если $y>H^{1/2}$, то $0<x<H^{1/2}$;

б) если $y<-H^{1/2}-1$, то $-H^{1/2}<x<0$.

Как следствие, можно предложить следующий алгоритм решения нашего уравнения.

1. Для каждого целого $y \in [-H^{1/2}-1,H^{1/2}]$ определить, будет ли число $y^4+2Hy+1$ точным квадратом.

2. Для каждого целого $x \in (-H^{1/2},H^{1/2})$, $x \neq 0$, решить уравнение $(*)$ относительно $y$, а именно, определить, является ли число

$$ y=\frac{H \pm \sqrt{H^2-2x^3+2x}}{2x} $$

целым.

Предложенный алгоритм решения требует всего $O(H^{1/2})$ извлечений квадратного корня. Используя Maple и процессор i7, можно решить миллион таких диофантовых уравнений за (примерно) 70 минут. Подозреваю, что если написать это на pari/gp, то будет быстрее.

Вот такая история

P.S. Еще один побочный результат нынешнего обсуждения: пример из статьи Stroeker R., de Weger B. Solving elliptic diophantine equations: the general cubic case // Acta Arith. 1999. Vol. 87. P. 339—365 о котором я рассказываю вот здесь http://dxdy.ru/post1421247.html#p1421247 оказывается совсем неудачным, ибо само уравнение можно свести к элементарной "разности квадратов" и, тем самым, решать уравнение Морделла в данном случае как-то неприлично.

Upd. На самом деле все еще проще: все решения уравнения $xy^2-Hx^2-y=0$ могут быть легко описаны следующим образом. Пусть $S$ --- множество всех целых делителей числа $H$. Тогда любое решение $(x,y)$ находится по формулам

$$x=\sqrt[3]{\frac{H+t}{t^2}}, \quad y=tx^2 \quad (t \in S).$$

Таким образом, все сводится к перебору делителей $H$.

[dmd]

nnosipov в сообщении #1487129 писал(а):

Для полноты картины сообщите, пожалуйста, сколько времени (и в какой CAS --- видимо, pari/gp?) работал Ваш i5, пока не решил этот миллион уравнений.

В данном случае больше суток, при этом на 3-х ядрах, и прошел только первые три четверти от миллиона $H$, дальше не стал ждать, выключил. Чуть позже проверю Ваш алгоритм и затем сообщу время выполнения в pari/gp (теперь уже весь миллион конечно проверю).

Возможно ли алгоритм для оценки границ из статьи Пулакиса модифицировать так, чтоб учитывался отличный от единицы старший коэффициент при $x^4$?

[nnosipov]

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

Возможно ли алгоритм для оценки границ из статьи Пулакиса модифицировать так, чтоб учитывался отличный от единицы старший коэффициент при $x^4$?

Только если этот старший коэффициент будет точным квадратом. Иначе задача станет, вообще говоря, существенно не элементарной. Конечно, бывают случаи, когда такие уравнения все же решаются с помощью элементарных рассуждений (например, уравнение $y^2=2x^4+1$), но это специально сконструированные примеры. В любом случае решать такие уравнения миллионами будет затруднительно.

-- Чт окт 15, 2020 13:46:19 --

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

Чуть позже проверю Ваш алгоритм и затем сообщу время выполнения в pari/gp

Вот это было бы любопытно. Вот код на Maple:

solve_int:=proc(a,b,c)
local S,d,d2;
if a=0 and b=0 then return {} end if;
if a=0 then if irem(c,b)=0 then return {-c/b} else return {} end if end if;
S:={};
d:=b^2-4*a*c;
d2:=isqrt(d);
if d2^2=d then
if irem(-b+d2,2*a)=0 then S:=S union {(-b+d2)/2/a} end if;
if irem(-b-d2,2*a)=0 then S:=S union {(-b-d2)/2/a} end if;
end if;
return S;
end proc:

runge:=proc(H)
local S,Q,res,x,y;
S:={0};
Q:=isqrt(H);
for y from -Q-1 to Q do
if isqrt(y^4+2*H*y+1)^2=y^4+2*H*y+1 then S:=S union {y} end if
end do;
for x from -Q to -1 do res:=solve_int(2*x,-2*H,x^2-1); if res<>{} then S:=S union res end if
end do;
for x from 1 to Q do res:=solve_int(2*x,-2*H,x^2-1); if res<>{} then S:=S union res end if
end do;
return S
end proc:

Разумеется, процедуру solve_int решения квадратных уравнений можно не переписывать, в pari/gp есть встроенные функции, которые это делают адекватно (например, nfroots).

[dmd]

15 минут на одном ядре:

? #
   timer = 1 (on)
?
? \r nno5.gp
? nno5()

H = 4657; #T = 10
[[-68, 4555], [-68, 4555], [-49, 2304], [-40, 1479], [-24, 329], [-4657, 21687648], [0, 1], [4657, 21687650], [7, 260], [33, 1222]]

H = 53684; #T = 10
[[-84, 6385], [-53684, 2881971855], [0, 1], [53684, 2881971857], [8, 929], [10, 1041], [20, 1519], [64, 4863], [66, 5105], [416, 173185]]

H = 95281; #T = 10
[[-1148, 1317821], [-80, 5071], [-95281, 9078468960], [0, 1], [95281, 9078468962], [7, 1156], [40, 3191], [60, 4939], [76, 6917], [657, 431794]]
time = 15min, 22,090 ms.

Случайно артефакт $y=\lfloor \sqrt{H}\rfloor$ зацепило для H = 4657, исправлять не стал.

gp-код:

nno5()=
{
 for(H=1, 10^6,
  Y= [];
  q= sqrtint(H);
  for(x=-q-1, q,
   y= x;
   if(issquare(y^4+2*H*y+1), Y= concat(Y, [y]));
   t2= H^2-2*x^3+2*x;
   if(issquare(t2), if(x,
    t= sqrtint(t2);
    y= (H+t)/2/x;
    if(y==floor(y), if(issquare(y^4+2*H*y+1), Y= concat(Y, [y])));
    y= (H-t)/2/x;
    if(y, if(y==floor(y), if(issquare(y^4+2*H*y+1), Y= concat(Y, [y]))));
   ))
  );
  if(#Y>9,
   T= apply(y -> [y, sqrtint(y^4+2*H*y+1)], Y);
   print("\nH = "H"; #T = "#T"\n"T)
  )
 )
};

[nnosipov]

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

15 минут на одном ядре

Как-то неожиданно быстро.

[nnosipov]

И еще одно упражнение: предлагается решить уравнение $xy^2-Hx^2-y-1=0$ для всех $H \in [1,10^6]$ и найти те $H$, для которых число решений максимально. (Здесь должно быть быстрее, чем в предыдущем случае; на pari/gp вычисления должны занять всего несколько минут.)

[dmd]

Неужели оно проще предыдущего случая?

Выглядит очень похоже

$y^4 - 4 H y - 4 H = (2 H x - y^2)^2$

или даже сложнее, если в форме Вейерштрасса

$\Big(4 H x\Big)^3 + 16 H \Big(4 H x\Big) + 16 H^2 = \Big(4 H (2 x y - 1)\Big)^2$

[nnosipov]

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

Неужели оно проще предыдущего случая?

Выглядит очень похоже

$y^4 - 4 H y - 4 H = (2 H x - y^2)^2$

Да, можно быстрее. В предыдущем случае нужно было решать порядка $H^{1/2}$ квадратных уравнений, а здесь можно обойтись всего лишь решением порядка $H^{1/3}$ квадратных уравнений.

[dmd]

Одно решение, 33 секунды.

H = 89; #T = 5
[[0, -1], [1, -9], [1, 10], [4, 19], [5, -21]]

[nnosipov]

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

Одно решение, 33 секунды.

Не понял: за 33 секунды было решено миллион уравнений? Это что за алгоритм такой?

Хотя при $H=89$ имеем 6 решений (а не 5), это значение $H$ действительно рекордное.

[dmd]

Хм. Значит не правильно повторил идею для этой задачи. Буду думать. Действительно, пропущено [91,90].

Не миллион, а гораздо больше. У нас миллион только значений $H$. В прошлой задаче для каждого $H$ в моём коде решалось $2H^{1/2}$ подзадач, но не квадратных уравнений. В каждой просто проверялось на квадратность подкоренного выражения (где if(issquare(t2),), и только если оно квадрат, то идём дальше.

[nnosipov]

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

Значит не правильно повторил идею для этой задачи.

Видимо, да. Все-таки при таком подходе (когда фактически решается уравнение $y^4-4Hy-4H=z^2$) у Вас время счета должно быть порядка 15 минут. Но уравнение $xy^2-Hx^2-y-1=0$ не обязательно сводить к уравнению выше. На Maple, i7 на миллион таких уравнений уходит 5 минут. Поскольку pari/gp быстрее Maple раз в 5, действительно можно уложиться в одну минуту.