Автокомментарий: текст был переведен пользователем polosatic (оригинальная версия, переведенная версия, сравнить).

Если из одного большого вещественного числа вычесть другое большое вещественное число, относительная погрешность может сильно поменяться. Например, если окажется, что $$$a - b = 1$$$, при этом $$$a = 265265266$$$, а $$$b = 265265265$$$, и при вычислениях мы получили значения a = 265265266.1 и b = 265265264.9 (в обоих случаях ошибка меньше $$$10^{-9}$$$ относительно истинного значения), то a - b = 1.2, где ошибка уже $$$20\%$$$.

Ваша функция find возвращает выражение (lx * ry - ly * rx) / (k * (lx - rx) + ry - ly), где делимое вычисляется в целых числах, и там всё в порядке, а вот делитель -- в (лонг) даблах, и там может возникнуть этот эффект. Попробуем построить такой тест, где разность в делителе очень маленькая.

Заметим, что эта разность -- это $$$[P, R - L] / P_x$$$, поэтому просто попробуем сделать первое векторное произведение поменьше (меньше $$$\pm 1$$$ не получится, но мы сделаем $$$\pm 2$$$, что отличается всего в 2 раза), а $$$P_x$$$ побольше. Для этого можно взять конструкцию $$$R - L = (F_k, F_{k+1})$$$, а $$$P = (F_{k+1}, F_{k+2})$$$, где $$${F_i}$$$ есть последовательность чисел Фибоначчи. Однако чтобы прямая $$$\left<0, P\right>$$$ вообще пересекала отрезок через $$$L$$$ и $$$R$$$, мы можем взять $$$P = (F_{k+1}, F_{k+2})$$$, $$$L = P - (F_k, F_{k+1})$$$, $$$R = P + (F_k, F_{k+1})$$$ (в зависимости от чётности $$$k$$$ точки $$$L$$$ и $$$R$$$ могут быть в другом порядке, нам без разницы).

Итак, рассмотрим такой тест:

2
102334155 165580141
433494437 701408733
1
267914296 433494437

Значение py * (lx - rx) + px * (ry - ly) (которое явным образом не вычисляется в функции find) должно быть равно $$$2$$$, а k * (lx - rx) + ry - ly, соответственно, должно быть равно $$$2/267914296$$$. Но если мы выведем (k * (lx - rx) + ry - ly) * px, мы получим 1.99611705541610717773, то есть мы накопили около $$$0.2\%$$$ ошибки. В самом деле, ry - ly вычислено точно, а k * (lx - rx) вычислилось как 535828592.00000000745058059692 (должно быть $$$535828592 + 2/267914296$$$, что примерно 535828592.00000000746507383092390112694844772299870). Дело в том, что мантиссы в лонг даблах хватает только на 53 бита информации, а нам надо помнить 9 цифр целой части, потом ещё 8 нулей и ещё потом какое-то значимое количество цифр (в данном случае хватило только на две цифры). В 128-битном флоате мантисса занимает больше памяти, поэтому там хватает сохранить ещё сколько-то цифр, поэтому при вычитании всё равно будет нормальная погрешность.

Можно вылечить это решение, например, тем, что функция find будет возвращать не (lx * ry - ly * rx) / (k * (lx - rx) + ry - ly), а (ld)(lx * ry - ly * rx) / (py * (lx - rx) + px * (ry - ly)) (это значение в px раз меньше предыдущего, зато и делимое, и делитель вычисляются в интах, поэтому без ошибок). Тогда выводить надо будет не P.x / X, а просто (ld)1 / X.