найти корешок по двум элементам цепной дроби

[Andrey A]

P.S. Во второй версии оплошность. Ограничив процесс пятью итерациями, урезал заодно и конечную проверку, а этого делать не следовало.
На всякий случай: если на заключительном этапе (наименьшее $m=...$) в столбце "J" под заголовком $a_2$ все пять строк заполнены числами, лучше проверить результат в первой версии. Бывает редко, не хочется переписывать таблицу.

[Andrey A]

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

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

Вроде бы и не редко. Переписал. Ссылка прежняя https://yadi.sk/i/D_CzGy4qt8reiw

В случае одного знака $ задача решается устно:
Для четного $u\ \ m=(u/2)^2+1.$
Для нечетного $u\ \ m=u^2+2.$

Для $3$-х знаков постараюсь всё же написать, хотя тема не для Excel.

[Andrey A]

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

Для $3$-х знаков постараюсь всё же написать

Готово. Ссылка та же.

[Andrey A]

P.S. Исправлена ошибка 27.10.2020.

[Dmitriy40]

(Тупой перебор)

Сделал втупую перебор $a_0$ на PARI/GP в соответствии с наблюдением mihaild и Null хотя бы просто для проверки:

? sq(b,c)=my(a,m,x,y); for(a=0,oo, x=(a+c/(b*c+1))^2; y=(a+(c+1)/(b*c+b+1))^2; m=ceil(min(x,y)); if(m<max(x,y), printf("sqrt(%d)=[%d,%d,%d,...]\n", m,a,b,c); break) );
? sq(5,4)
sqrt(174)=[13,5,4,...]
? sq(7,11)
sqrt(2129)=[46,7,11,...]
? sq(997,999)
sqrt(249000003002)=[498999,997,999,...]
time = 1,295 ms.
? sq(999,997)
sqrt(994014984005)=[997003,999,997,...]
time = 2,606 ms.

Не быстро конечно. Вроде линейно по $a_0$, которое очевидно квадратично по $a_1, a_2$ (точнее скорее линейно по их произведению $a_1 a_2$), так что давать $a_1,a_2 \approx 10^9$ не стоит.

-- 27.10.2020, 07:24 --

Более умный перебор (для совсем малых аргументов может ошибаться UPD: проверил, не ошибается) на PARI/GP:

? sq(b,c)=my(a,k,m,x,y,bc1); bc1=b*c+1; for(k=1,oo, a=floor((k*bc1/c-c/bc1)/2); x=(a+c/(b*c+1))^2; y=(a+(c+1)/(b*c+b+1))^2; m=ceil(min(x,y)); if(m<max(x,y), printf("sqrt(%d)=[%d,%d,%d,...]\n", m,a,b,c); break) );
? sq(999,997)
sqrt(994014984005)=[997003,999,997,...]
? sq(9999,9997)
sqrt(9994001499840005)=[99970003,9999,9997,...]
time = 62 ms.
? sq(99999,99997)
sqrt(99994000149998400005)=[9999700003,99999,99997,...]
time = 765 ms.
? sq(999999,999997)
sqrt(999994000014999984000005)=[999997000003,999999,999997,...]
time = 8,019 ms.

Кажется уже линейно по коэффициентам ... С коэффициентами под $10^9$ (каждый) должно справиться за пару-тройку часов ...
В принципе можно даже и к операциям в целых перейти, тогда вопрос точности вычислений и возможной ошибки отпадёт, но анализировать что-то лень ... С другой стороны PARI похоже и сам всё в целых (точнее конечно рациональных) считает, поленился проверить типы данных, но раз не ошибается и время не зависит от установленной точности вычислений с плавающей точкой, то очень похоже.

-- 27.10.2020, 08:06 --

Не знаю надо ли пояснять идею, она проста: перебирать только те $a_0$, которые после возведения в квадрат дают максимально близкие снизу к целым числам, для чего оценивается слагаемое $k=\frac{2 a_0 a_2}{a_1 a_2 + 1}+(\frac{a_2}{a_1 a_2 + 1})^2$ из раскрытого $(a_0+\frac{a_2}{a_1 a_2 +1})^2$.

[Dendr]

Все не получается выделить достаточно времени, поэтому весьма коряво и не для слишком больших значений, зато для произвольной (тоже ограниченной, впрочем, если нужна точность) длины цепи.
Код тут:
http://cpp.sh/85mur
Экзешник тут:
https://yadi.sk/d/pX5RloDX8ja6aw

[Dmitriy40]

По подсказкам Andrey A и небольшой своей доработкой реализовал его алгоритм на PARI/GP:

? {sq3(b,c)=my(x,pq,r,m,p,q2,v,i,pf);
   pf=default(realprecision);default(realprecision,300);\\Для больших чисел может понадобиться и большая точность
   v=contfracpnqn([b,c,b]); if((v[2,2]*v[2,1])%2==0, v[2,2]=-v[2,2], v=contfracpnqn([b,c,c,b]));
   m=round(((v[2,2]*v[2,1]/2)%v[1,1]+v[2,1]/v[1,1])^2);\\Инициализация m каким-то большим точно подходящим числом
   while(1,\\Цикл максимального уменьшения найденного числа сохраняя его "подходящесть"
      p=0; pq=contfracpnqn(contfrac(sqrt(4*m),12),10+2);\\Разложение sqrt(4m) в 10 (+2 запас) подходящих дробей, пока 10-ти хватало
      forstep(i=10,1,-1,\\Перебираем подходящие дроби в обратном порядке, от больших
         r=m-pq[1,i]+pq[2,i]^2;\\Это будет будущее m=m-P+Q^2
         if(r<2 || r>=m, next);\\Такие не интересуют
         v=contfrac(sqrt(r),5); if(#v>3 && v[2..3]==[b,c], m=r; p=pq[1,i]; q2=pq[2,i]^2; break);\\Если новое подходит, то его и используем
      );
      if(p==0, break);\\Нашли ли новое m
      r=p%(2*q2);\\Попытаемся уменьшить p до минимума
      while(r<p,\\Пока оно меньше уже найденного проверяем подходит ли такое p
         x=round(r^2/q2/4); if(x>=m, break);\\Это новое m, если больше уже найденного, то не интересно
         v=contfrac(sqrt(x),5); if(#v>3 && v[2..3]==[b,c], m=x; break);\\Проверяем новое m и если подходит, то используем и выходим
         r+=2*q2;\\Продолжаем перебор в сторону увеличения до нахождения наименьшего подходящего, пока хватало 3-4 циклов
      );
   );
   default(realprecision,pf);\\Вернём точность вычислений обратно
   m
};
? m=sq33(7,11);printf("sqrt(%d)=%d\n", m,contfrac(sqrt(m),5))
sqrt(2129)=[46,7,11,3]
? m=sq33(22,45);printf("sqrt(%d)=%d\n", m,contfrac(sqrt(m),5))
sqrt(1048669)=[1024,22,45,2,7]
? m=sq33(222,11);printf("sqrt(%d)=%d\n", m,contfrac(sqrt(m),5))
sqrt(6522939)=[2554,222,11,2,9]

Для всех тестовых векторов оказалось достаточно 10-ти элементов в цепной дроби разложения корня, 12 элементов берётся исключительно из-за возможного усечения результата contfrac если последнее значение окажется 1.
Опять же для всех тестовых векторов достаточно 5-ти итераций цикла "спуска" по $m$, для больших чисел обычно хватает 2-3 итерации.
Использование встроенной функции contfrac потребовало увеличения точности операций с плавающей точкой, поставил с запасом, на скорость это влияет мало.
Отдельно с помощью Andrey A было проверено что возможна реализация только в длинных рациональных числах, уже не зависящая от точности вычислений с плавающей точкой. Но она медленнее (на не слишком больших числах типа как вон ниже где-то вдвое) и код вдвое объёмнее, потому не привожу, на алгоритм оно не влияет (за одним исключением: перебирать подходящие дроби можно в прямом порядке, причём тогда и разложение в цепную дробь тоже становится итерационным и можно оборвать процесс ранее 10-ти членов, но на скорость это влияет слабо).

Отличия этой реализации от алгоритма Andrey A:
подходящие дроби перебираются от конца к началу;
они формируются сразу все 10 (+2 для обхода оптимизации PARI/GP), а не итерационно;
вместо предложенного им двоичного поиска (деления пополам) наименьшего подходящего $P$ берётся сразу минимальное положительное $P$ и перебирается подряд (с шагом арифметической прогрессии) в сторону увеличения, обычно хватает 1-2 итераций, больше 4 на тестовых векторах не понадобилось ни разу.

Замеры скорости.
Сравнивал со своей реализацией полного перебора выше (оптимизированной примерно вдвое).
98901 комбинацию $[?,1,1]-[?,999,99]$ мой перебирает за 11с, этот за 10.5с.
98901 комбинацию $[?,1,1]-[?,99,999]$ мой перебирает за 105с, этот за 9.7с.
9801 комбинацию $[?,9901,9901]-[?,9999,9999]$ мой перебирает за 224с, этот за 1.0с.
9801 комбинацию $[?,999999901,999999901]-[?,999999999,999999999]$ этот перебирает за 1.2с (тут не хватило точности в 200 значащих цифр, поставил 300, PARI/GP позволяет хоть миллион значащих цифр использовать).
Т.е. зависимость от величины чисел практически отсутствует (в пределах десятка процентов) и в среднем решение занимает порядка 100мкс (3.7ГГц CPU).

[Andrey A]

Dmitriy40 Спасибо! Очень интересно и поучительно.
Фразу ... зависимость от величины чисел практически отсутствует я для себя перевожу так: программа написана толково.

[Dmitriy40]

(Оффтоп)

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

Фразу ... зависимость от величины чисел практически отсутствует я для себя перевожу так: программа написана толково.

Нет, правильным переводом будет "Алгоритм выполняет много сложных операций (корень, цепная и подходящие дроби) и мало простых (арифметика), потому накладные расходы на операции с длинными числами теряются на общем фоне.". Это достоинство не реализации программы, а наличия встроенных функций в PARI/GP. И разумеется Вашего алгоритма, что он не $O(a^n)$, не $O(n^2)$, не $O(n \log n)$, и даже не $O(n)$, а всего лишь $O(\log n)$ или менее (навскидку, просто столько занимают операции с длинными числами, а тут они необходимы в том или ином виде). Т.е. зависимость времени работы от длины чисел присутствует лишь в реализации чисел, но не в алгоритме, он сам по себе занимает близкое к константе время, хотя это чисто эмпирический факт, без доказательства (независимости количества итераций спуска по $m$ от величины чисел). Хорошей аналогией будет операция взятия остатка (дробной части числа), её можно вычислить как последовательность вычитаний с линейной зависимостью от величины чисел, а можно одной операцией деления (с округлением и вычитанием), вообще не зависящей от величины чисел (в частном случае их реализации).
Я программу почти не оптимизировал, больше боролся за простоту кода.

Да, хотелось бы обратить внимание всех, я детально в алгоритме не разбирался, лишь по минимуму для написания работающего кода, потому не понял есть ли нормальное доказательство что он выдаёт именно минимальное значение. По-моему это лишь эмпирический факт, но возможно это гарантируется формулой спуска $m=m-P+Q^2$, не знаю (в остальных местах гарантия есть).

[Andrey A]

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

... есть ли нормальное доказательство

Доказательства нет и контрпримера нет. Тонкое место.