По подсказкам
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-ти итераций цикла "спуска" по

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

берётся сразу минимальное положительное

и перебирается подряд (с шагом арифметической прогрессии) в сторону увеличения, обычно хватает 1-2 итераций, больше 4 на тестовых векторах не понадобилось ни разу.
Замеры скорости.
Сравнивал со своей реализацией полного перебора выше (оптимизированной примерно вдвое).
98901 комбинацию
![$[?,1,1]-[?,999,99]$ $[?,1,1]-[?,999,99]$](https://dxdy-04.korotkov.co.uk/f/b/3/0/b3033ec1978825929579e1ffddb14e0f82.png)
мой перебирает за 11с, этот за 10.5с.
98901 комбинацию
![$[?,1,1]-[?,99,999]$ $[?,1,1]-[?,99,999]$](https://dxdy-01.korotkov.co.uk/f/0/7/d/07da52106684194c182cf314b4eda84f82.png)
мой перебирает за 105с, этот за 9.7с.
9801 комбинацию
![$[?,9901,9901]-[?,9999,9999]$ $[?,9901,9901]-[?,9999,9999]$](https://dxdy-01.korotkov.co.uk/f/8/4/4/844e85ce2e086ab03330d9a1bd0282a382.png)
мой перебирает за 224с, этот за 1.0с.
9801 комбинацию
![$[?,999999901,999999901]-[?,999999999,999999999]$ $[?,999999901,999999901]-[?,999999999,999999999]$](https://dxdy-02.korotkov.co.uk/f/9/e/7/9e7c8b2ff1e105b4469ff044235863b782.png)
этот перебирает за 1.2с (тут не хватило точности в 200 значащих цифр, поставил 300, PARI/GP позволяет хоть миллион значащих цифр использовать).
Т.е. зависимость от величины чисел практически отсутствует (в пределах десятка процентов) и в среднем решение занимает порядка 100мкс (3.7ГГц CPU).