Algorithm Description

设$T(cos\theta ,sin\theta )$,则有
$$PT+QT=\sqrt{(P_x-cos\theta {)}^2+si{n}^2\theta }+\sqrt{(Q_x-cos\theta {)}^2+(Q_y-sin\theta {)}^2}$$
$$PT+QT=\sqrt{P_x^2-2P_{x}cos\theta +1}+\sqrt{Q_x^2+Q_y^2+1-2Q_{x}cos\theta -2Q_{y}sin\theta }$$
由费马原理，光线沿$PT+QT$最短的路径传播，因此只需对上式求导求极小值点。关于$\theta$求导得
$$\frac{P_{x}sin\theta }{\sqrt{P_x^2-2P_{x}cos\theta +1}}+\frac{Q_{x}sin\theta -Q_{y}cos\theta }{\sqrt{Q_x^2+Q_y^2+1-2Q_{x}cos\theta -2Q_{y}sin\theta }}$$
故只需用二分法解非线性方程
$$P_{x}sin\theta \sqrt{Q_x^2+Q_y^2+1-2Q_{x}cos\theta -2Q_{y}sin\theta }+(Q_{x}sin\theta -Q_{y}cos\theta )\sqrt{P_x^2-2P_{x}cos\theta +1}=0$$
$$T_x=cos\theta$$
$$T_y=sin\theta$$
由对称性易知
$$Rx=\frac{2Q_y-Qxtan\theta -k{Q}_x+k(2Tx-Qx)-2Ty}{k-tan\theta }$$
$$Ry=Qy-(Qx-Rx)\theta$$

Code

#include<iostream>
#include<cmath>
#include <iomanip>
using namespace std;

int main(int argc, char *argv[])
{
    if(argc != 4) {
        cout << "Usage: " << argv[0] << " <Px> <Qx> <Qy>\n";
        return 1;
    }

    long double Px = stold(argv[1]);
    long double Qx = stold(argv[2]);
    long double Qy = stold(argv[3]);

    if(Px >= -1 || Qx >= 0 || Qy <= 0 || Qx * Qx + Qy * Qy <= 1) {
        cout << "输入错误，Px应该小于-1，Qx应该小于0，Qy应该大于0，Qx^2+Qy^2应该大于1\n";
        return 1;
    }
    
    long double Tx = 0;
    long double Ty = 0;
    long double theta = 0;
    long double low = 3.14159265358979323 , high = 1.570796326794896;
    long double k = 0;
    while(1)
    {
        theta = (low + high) / 2;
        long double eq1 = Px * sin(theta);
        long double eq2 = sqrt(Qx * Qx + Qy * Qy +1 - 2 * Qx * cos(theta) - 2 * Qy * sin(theta));
        long double eq3 = (Qx * sin(theta) - Qy * cos(theta));
        long double eq4 = sqrt(Px * Px -2 * Px * cos(theta) + 1);
        long double res = eq1 * eq2 + eq3 * eq4;
        long double absres = abs(res);
        if(absres <= 1e-7)
        {
            Tx = cos(theta);
            Ty = sin(theta);
            k = 1 / tan(3.14159265358979323 - theta);
            break;
        }
        else if(res > 0)
        {
            low = theta;
        }
        else
        {
            high = theta;
        }
    }
    long double eq5 = 2 * Qy - Qx * tan(theta) + k * (2 * Tx - Qx) - 2 * Ty;
    long double eq6 = k - tan(theta); 
    long double Rx = eq5 / eq6;
    long double Ry = Qy - tan(theta) * (Qx - Rx);
    cout << "T = (" << fixed << setprecision(6) << Tx << " , " << fixed << setprecision(6) << Ty << ") , R = (" << fixed << setprecision(6) << Rx << " , " << fixed << setprecision(6) << Ry << ")";
    return 0;
}

Results

$$P=(-2,0),Q=(-1,1):T=(-0.885670,0.464316),R=(-0.380057,0.674993)$$

$$P=(-10,0),Q=(-2,1):T=(-0.959312,0.282350),R=(0.304214,0.321811)$$

$$P=(-1.000001,0),Q=(-2,2):T=(-1.000000,0.000002),R=(0.000007,1.999996)$$

$$P=(-2,0),Q=(-1,0.000001):T=(-1.000000,0.000001),R=(-1.000000,0.000001)$$

$$P=(-2.33,0),Q=(-3,1):T=(-0.989279,0.146038),R=(1.182424,0.382590)$$

$$P=(-3,0),Q=(-1,0.5):T=(-0.922615,0.385721),R=(-0.786920,0.410917)$$

$$P=(-3,0),Q=(-2,10):T=(-0.827028,0.562160),R=(8.380296,2.944148)$$

$$P=(-3,0),Q=(-3,1):T=(-0.987408,0.158192),R=(1.187435,0.329136)$$

$$P=(-10,0),Q=(-2,1):T=(-0.959312,0.282350),R=(0.304214,0.321811)$$

$$P=(-1024,0),Q=(-8,4):T=(-0.970066,0.242842),R=(7.000894,0.244735)$$

Conclusion

本实验提高精度的主要措施：

1. 使用long double 类型；

2. 用较多位数来表示$\pi$；

3. 将含有除法的方程交叉相乘，通过乘法代替除法，以减少在求商时的误差；

4. 在用二分法解非线性方程时，限制条件是结果的绝对值$\le 10^{-7}$。

但由于本实验的方程较为复杂，含有三角函数、根式、平方等易造成误差放大的因素，暂时还未找到其他较好的减少误差的办法，但还尝试了另一种思路，叙述如下：

设OT的延长线交PQ于R，则由角平分线定理，有$$\frac{QT}{PT}=\frac{QR}{RP}=\frac{Q_y}{R_y}-1$$
设$T(x,y),k=\frac{Qy}{Qx-Px}$，带入上述方程化简得：
$$\frac{Q_y^2(kx-y{)}^2+k^2{P}_x^2{y}^2-2k{Q}_{y}Pxy(kx-y)}{k^{2}P{x}^2{y}^2}=\frac{Q_y^2+Q_x^2+1-2y{Q}_y-2x{Q}_x}{1+P_x^2-2x{P}_x}$$
故只需用for循环遍历或二分法解非线性方程$$(Q_y^2(kx-y{)}^2+k^2{P}_x^2{y}^2-2k{Q}_{y}Pxy(kx-y))(1+P_x^2-2x{P}_x)-(k^{2}P{x}^2{y}^2)(Q_y^2+Q_x^2+1-2y{Q}_y-2x{Q}_x)=0$$
$$y=\sqrt{1-x^2}$$

由对称性易知
$$Rx=\frac{\frac{QxTy}{Tx}+\frac{Tx(2Tx-Qx)}{Ty}+2Ty-2Qy}{\frac{Tx}{Ty}+\frac{Ty}{Tx}}$$
$$Ry=Qy-\frac{Ty}{Tx(Qx-Rx)}$$
但验证结果时发现，上述方法在遇到一些极端情况如$P=(-1.000001,0),Q=(-2,2)$时，中间的运算过程会出现极小的浮点数导致无法继续运算，所以这种思路在细节上仍有待改进。