Не более пяти решений

[Dmitriy40]

Я немного пописал на бумажке, прямо по исходному уравнению, оно фактически квадратное относительно $y$. Причём из медленных операций надо одно извлечение корня и одно деление, на каждую пару $H,x$. Но, что и хочу заметить, эти медленные операции иногда вполне можно и не выполнять. Ведь при некоторых значениях под корнем заведомо отрицательная величина, а при некоторых (других) деление гарантированно не даст целого числа, и оба этих условиях можно проверить быстрее вычисления самих операций. А если вывести точную конечную формулу (я пока поленился), то там вроде бы можно использовать чётность/нечётность (или другой признак делимости), т.е. сразу исключать некоторые комбинации $H,x$ из рассмотрения, не вычисляя медленных операций.
В общем есть над чем подумать даже кроме математики.

Кстати корень у меня вычислялся поразрядным взвешиванием, через 32/64 умножений, плюс оценка начального приближения (самого старшего бита). Когда будет быстрее метод Герона я не разбирался, но тут на форуме была тема примерно об этом и там он кажется был быстрее (на языке С).

Если Вы можете из Maple вызвать произвольную функцию из внешней dll, то могу сделать такую dll сразу с решением квадратного уравнения, правда параллелить тогда придётся Вам. Ну или сделать сразу перебор всех $x$ в задаваемом диапазоне для заданного $H$, уже многопоточно. Прикидочно это делов на день-два-три (отладить, особенно если с многопоточным перебором). Правда не знаю как вернуть массив результатов, это надо договариваться.

-- 13.04.2020, 21:09 --

nnosipov
Сразу добавлю вопрос: у меня под корнем получается кубическая парабола $8x^3-4(H+1)x+1$, она всегда имеет положительные значения, при любых $H$, причём после извлечения корня оно ещё и растёт быстрее линейного, значит корни $y$ потенциально могут расти в бесконечность, значит на $x$ не может быть ограничения на перебор сверху, почему у Вас какие-то ограничения на циклы? Или я что-то не так понял?

[nnosipov]

Dmitriy40
Пожалуйста, посмотрите файл, который я приложил --- там описание алгоритма. Все ли там понятно? Насчет математики (т.е. почему этот алгоритм корректно решает задачу) объясню обязательно, но позже. Сейчас важно Ваше профессиональное (как программиста) мнение о реализации этого алгоритма. Я предварительно прикинул размер дискриминанта обоих квадратных уравнений --- это по порядку $|H|^{3/2}$. Тотально в алгоритме решается примерно $4|H|^{1/4}$ квадратных уравнений. Можно ли теперь обработать все значения $H$ до $10^9$?

-- Вт апр 14, 2020 01:43:59 --

Сейчас сам прикинул: если исходить из этого:

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

например до $10^{19}$ корень извлекается 15 млн.раз/с

то где-то сутки счета. В принципе, реалистично

Расчеты таковы. Суммарное число уравнений равно

$$\sum_{H=1}^{M} H^{1/4} \sim M^{5/4}$$

При $M=10^9$ получим $10^{12}$ уравнений. Поделив на 15 млн., получим что-то около суток.

-- Вт апр 14, 2020 02:04:14 --

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

Или я что-то не так понял?

Ну, немного математики здесь все-таки есть. Но она, поверьте, школьная. Прилагаю еще один файл со слайдами доклада моих школьников, который они сделали вот здесь http://www.school.ioffe.ru/readings/2018/ Там есть теорема, на которой основан алгоритм (см. слайд 8). Доказательство, правда, не приведено, но оно написано в самой работе. Если нужно, прикреплю и файл самой работы.

[Dmitriy40]

nnosipov
Ну так этот алгоритм прямо виден в Вашем коде выше на Maple, в нём всё понятно (кроме смысла пределов перебора, но тут уж проще поверить Вам на слово чем самому разбираться откуда они взялись). Фактически просто несколько циклов решений квадратных уравнений. Ничего принципиально сложного. Скорость оценить пока не могу, более-менее с потолка можно взять несколько миллионов решений уравнений в секунду (корень, деление, умножения/сложения/вычитания, проверки).

По диапазону. В первом уравнении под корнем максимум будет $8x^3$, хотелось бы его иметь до $10^{19}$ (это примерно $2^{64}$, чтобы не заморачиваться с длинной арифметикой), соответственно хотелось бы $x<1.3 \cdot 10^6$, с небольшими ухищрениями можно расширить до $x<10^7$, но примерно с полутора-двухкратным замедлением вычислений можно и до $x<3.5 \cdot 10^{12}$ (это даст подкоренное выражение до $2^{128}$). Это под x64 винду. Под x32/x86 будет ещё вдвое-втрое медленнее, потому не интересно. Второе уравнение под корнем будет максимум примерно $4k^2$, надеюсь это меньше чем $8x^3$ и оценку по $x$ не ухудшит.
И это оценка под те алгоритмы которые мне уже известны и опробованы, а ведь можно и с плавающей точкой считать и потом лишь проверить что результат целочисленный (это быстро, только умножения, ни делений ни корней). Но тут вопрос потери точности, его надо обдумать. Зато в этом варианте можно использовать векторизацию AVX даже на одном ядре.

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

-- 13.04.2020, 22:19 --

Упс, с $8x^3$ я ошибся, там 4-я степень. Значит придётся сразу делать 128 бит под корнем, это несложно, процедура давно написана, просто она дольше (примерно втрое).

[nnosipov]

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

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

Да, само собой, примеры для тестирования будут. Буду очень рад, если возьметесь довести все это до ума.

В приложенном файле примеры всех отрицательных значений $H$ (до какой-то границы), когда уравнение имеет 5 или 6 решений (в каждой строке: значение $H$ и соответствующее ему множество решений).

-- Вт апр 14, 2020 02:49:51 --

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

Упс, с $8x^3$ я ошибся, там 4-я степень. Значит придётся сразу делать 128 бит под корнем, это несложно, процедура давно написана, просто она дольше (примерно втрое).

Выше я оценивал максимальное значение дискриминанта обоих уравнений как $H^{3/2}$. При $H=10^9$ получим $10^{14}$. Что укладывается в диапазон до $10^{19}$. Или Вы о чем-то другом?

-- Вт апр 14, 2020 02:54:14 --

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

соответственно хотелось бы $x<1.3 \cdot 10^6$

Без проблем: ведь $x<Q(m)$, а $Q(m)$ при $m=H^{1/4}$ тоже порядка $H^{1/4}$ --- в этом вся фишка алгоритма.

-- Вт апр 14, 2020 03:03:41 --

Еще один момент, по поводу оценки общего времени работы: можно прогнать до, например, границы $10^6$, замерить время, а потом оценивать, исходя из того, что при увеличении границы на порядок время возрастет в $10^{5/4} \approx 18$ раз (конечно, это нижняя оценка).

[Dmitriy40]

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

Буду очень рад, если возьметесь довести все это до ума.

Ну совсем до ума лень, это много работы, а вот сварганить утилитку которая за недельку куда-то там досчитает, это могу, пока всё более-менее понятно. Или сразу в виде dll с многопоточным перебором по аргументу квадратного уравнения в заданных пределах с выдачей списка решений плюс оболочка с оставшимися циклами. Я ещё подумаю насколько второе сложнее. Хотя, его же тоже можно потом прикрутить к решению уравнения. Подумаю в общем.
"До ума" обычно подразумевает добавку нормального интерфейса к библиотеке, обработку ошибок, документацию, демопримеры, может быть тесты, иногда контроль правильности вычислений, балансировку нагрузки, оптимизацию под разные процессоры, да много всего ... Обычно это всё лень, пишу только саму вычислительную часть и всё, часто даже с необходимостью постпроцессинга в чём-то ещё (эксель, PARI/GP, а то и отдельные свои утилиты), так для себя получается удобнее чем делать всё "по уму".

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

Выше я оценивал максимальное значение дискриминанта обоих уравнений как $H^{3/2}$. При $H=10^9$ получим $10^{14}$. Что укладывается в диапазон до $10^{19}$. Или Вы о чем-то другом?

Нет, похоже именно об этом, именно про подкоренное выражение. Просто у меня есть две процедурки, извлечение целочисленного корня из 64 и 128 битного числа, вторая почти вчетверо медленнее (как оказалось, только что померил, думал лишь вдвое), остальное то не сильно важно какой разрядности (64 или 128).

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

Еще один момент, по поводу оценки общего времени работы: можно прогнать до, например, границы $10^6$, замерить время, а потом оценивать, исходя из того, что при увеличении границы на порядок время возрастет в $10^{5/4} \approx 18$ раз (конечно, это нижняя оценка).

Да я обычно даже проще делаю: запускаю счёт скажем тысячи $H$ подряд (на десятки-сотни секунд в сумме) и замеряю время при старте с разных начальных точек. Если для теста скорости.

[nnosipov]

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

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

Ну нет, это уже перебор

[Dmitriy40]

Так, на PARI/GP код написал (для разбирательства как там и что и для будущей проверки), скорость около $H=10^9$ всего 130 $H$ в секунду. Теперь переношу на асм, пока без особой оптимизации.

[Dmitriy40]

В одном потоке до миллиона H считает 2.7с, до 2 млн 6.6с, до 10млн 55с, до 100млн почти 20минут. Около $H=10^9$ миллион H считается за 31с, около $H=10^{12}$ тот же миллион H считается за 237с.

[nnosipov]

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

Около $H=10^9$ миллион H считается за 31с

Так это же хорошо: тысяча миллионов --- это 30000 сек, т.е. меньше 9 часов. Я правильно понял?

-- Вт апр 14, 2020 10:37:17 --

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

до 10млн уже 55с.

О, это что, все значения $H$ до $10^7$ обработаны? На python на это ушло около часа. А значения $H$ точно отрицательны? Для положительных все гораздо быстрее.

[Dmitriy40]

Да, за 20 минут просчитались все H до минус ста миллионов. Новых 6 и более не найдено, пятёрок 166 штук (из которых 99 штук до 10млн).
Замедление примерно 20 раз на порядок по H. До миллиарда надо примерно 400 минут или 7 часов. Примерно. Как раз поспать.

[nnosipov]

Dmitriy40
круто!

[Dmitriy40]

Наблюдение: почти все пятёрки идут парами, для близких $H$, разница меньше 0.1%. При этом $x$ у них отличается на 1, а сами $x$ равны примерно $2/3m$, и некоторые $y$ совпадают, причём все такие пары найдены циклом по $Q$, первый цикл ничего не обнаружил (я его решения игнорирую если они попадают в диапазон $-Q..Q$, они гарантированно будут найдены вторым циклом). Но иногда проскакивают одиночные решения с очень большими $x$, из первого цикла по $m$ (обычно кстати для близких к нулю $k$), вот пример:

H=-93286619, m=99, Q=97, n=5: (-69,9659) (-1,9659) (-1,-9658) (0,-1) (70,9659)
H=-93325259, m=99, Q=97, n=5: (-70,-9661) (0,-1) (1,-9661) (1,9660) (69,-9661)
H=-96049797, m=99, Q=99, n=5: (-6931,13861) (6929,-13859) (-1,9801) (-1,-9800) (0,-1)
H=-98773779, m=100, Q=99, n=5: (-70,9939) (-1,9939) (-1,-9938) (0,-1) (71,9939)
H=-98813539, m=100, Q=99, n=5: (-71,-9941) (0,-1) (1,-9941) (1,9940) (70,-9941)

Это кстати конец перед 100 миллионами.
Наверное это можно как-то использовать для ускорения перебора ... если идти на рекорд по $H$.

-- 14.04.2020, 07:01 --

Только что нашлось новое:

H=-175481375, m=116, Q=114, n=6: (-9368,18735) (9366,-18733) (-23,13247) (-9,13247) (0,-1) (32,13247)

[nnosipov]

Что касается величины решений. При $H<0$ верна следующая оценка:

$$|x| \leqslant 1+\sqrt{\frac{|H|+3}{2}}.$$

Эта оценка точна, она достигается для значений $H=-2t^2+3$, где $t \geqslant 2$, одним из решений будет $(x,y)=(-1-t,2t+1)$. Значения $H=-96049797$ и $H=-175481375$ как раз такие.

Вот с числом решений все очень загадочно. Будет ли еще одна семерка? Можно делать ставки.

-- Вт апр 14, 2020 11:22:50 --

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

Наблюдение: почти все пятёрки идут парами, для близких $H$, разница меньше 0.1%. При этом $x$ у них отличается на 1, а сами $x$ равны примерно $2/3m$, и некоторые $y$ совпадают

Это для меня новый момент, надо будет обдумать. Возможно, есть какая-то аналитическая закономерность.

Кстати, с пятерками в случае положительных $H$ все просто: они все явно выписываются и оказываются экспоненциально редки (до $10^9$ их всего пять штук). А в случае отрицательных $H$ стоит ожидать пару сотен пятерок до $10^9$.

-- Вт апр 14, 2020 12:12:51 --

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

если идти на рекорд

Если хотим поставить рекорд по числу решений, то имеет смысл пройтись по $H=-2t^2+3$. Две имеющиеся семерки (для $H=-1$ и $H=-1219919$ как раз такого вида).

Увы, при $t \leqslant 10^5$ ничего нового нет.

[dmd]

Есть еще
H = -854699 [[0, -1], [1, -925], [1, 924], [21, -925], [-22, -925], [349, -1048]]

[nnosipov]

Это третья и последняя шестерка до $10^7$.

-- Вт апр 14, 2020 14:07:47 --

Есть экспоненциально редкие семейства четверок. А вот с пятерками ничего аналитического не находится, увы.