Трёхмерный треугольник

#include <iostream>
#include <cmath>
 
struct point {
    double x;
    double y;
    double z;
};
 
point operator* (const double &d, const point &p) {
    return {p.x*d, p.y*d, p.z*d};
}
 
point operator* (const point &p, const double &d) {
    return {p.x*d, p.y*d, p.z*d};
}
 
point operator/ (const point &p, const double &d) {
    return {p.x/d, p.y/d, p.z/d};
}
 
point operator+ (const point &a, const point &b) {
    return {a.x+b.x, a.y+b.y, a.z+b.z};
}
 
point operator- (const point &a, const point &b) {
    return {a.x-b.x, a.y-b.y, a.z-b.z};
}
 
double f(const point A, const point B, const point C, const point P) {
    return hypot(hypot(A.x-P.x, A.y-P.y), A.z-P.z) + 
           hypot(hypot(B.x-P.x, B.y-P.y), B.z-P.z) +
           hypot(hypot(C.x-P.x, C.y-P.y), C.z-P.z);
}
 
point getGradient(const point A, const point B, const point C, const point P, double eps = 0.0001) {
    point result;
    
    result.x = (f(A, B, C, {P.x+eps,P.y,P.z}) - f(A, B, C, {P.x-eps,P.y,P.z})) / (2*eps);
    result.y = (f(A, B, C, {P.x,P.y+eps,P.z}) - f(A, B, C, {P.x,P.y-eps,P.z})) / (2*eps);
    result.z = (f(A, B, C, {P.x,P.y,P.z+eps}) - f(A, B, C, {P.x,P.y,P.z-eps})) / (2*eps);
    
    return result;
}
 
int main()
{
    // 38.7644 // 14.6261 // 9.1231
    point A = {38, 38, 0}; // 2, 7, 3 // 1, 3, 0 
    point B = {54, 34, 0}; // 1, -2, 0 // 5, 6, 0
    point C = {66, 53, 0}; // 3, 1, -5 //6, 10, 0
    point P = A;
    double eps = 0.0001, step = 2000, len = 0.0, cnt = 0.0;
    do {
        len = f(A,B,C,P);
        P = P - step*getGradient(A,B,C,P);
        step *= 0.9;
        cnt++;
    } while (fabs(len - f(A,B,C,P)) > eps);
    std::cout << P.x << ' ' << P.y << ' ' << P.z << ' ' << f(A,B,C,P) << ' ' << cnt;
}

#include <iostream>
#include <cmath>
 
struct Point {
    double x, y, z;
};
 
struct Triangle {
    Point A, B, C;
};
 
Point operator* (const double &d, const Point &p) {
    return {p.x*d, p.y*d, p.z*d};
}
 
Point operator* (const Point &p, const double &d) {
    return {p.x*d, p.y*d, p.z*d};
}
 
Point operator/ (const Point &p, const double &d) {
    return {p.x/d, p.y/d, p.z/d};
}
 
Point operator+ (const Point &a, const Point &b) {
    return {a.x+b.x, a.y+b.y, a.z+b.z};
}
 
Point operator- (const Point &a, const Point &b) {
    return {a.x-b.x, a.y-b.y, a.z-b.z};
}
 
double f(const Triangle T, const Point P) {
    return std::hypot(T.A.x-P.x, T.A.y-P.y, T.A.z-P.z) + 
           std::hypot(T.B.x-P.x, T.B.y-P.y, T.B.z-P.z) +
           std::hypot(T.C.x-P.x, T.C.y-P.y, T.C.z-P.z);
}
 
Point getGradient(const Triangle T, const Point P, double eps = 0.0001) {
    Point result;
    result.x = f(T, {P.x+eps,P.y,P.z}) - f(T, {P.x-eps,P.y,P.z});
    result.y = f(T, {P.x,P.y+eps,P.z}) - f(T, {P.x,P.y-eps,P.z});
    result.z = f(T, {P.x,P.y,P.z+eps}) - f(T, {P.x,P.y,P.z-eps});
    return result / (2*eps);
}
 
double df(const Triangle T, const Point P, const Point grad, double x, double eps = 0.0001) {
    return (f(T, P - (x+eps)*grad) - f(T, P - (x-eps)*grad)) / (2*eps);
}
 
double getStep(const Triangle T, const Point P, double eps = 0.0001) {
    Point grad = getGradient(T, P, eps);
    double lx = 0.01, rx = lx, cx, df1;
    if (df(T, P, grad, lx, eps) >= 0) return lx;
    while (df(T, P, grad, rx, eps) < 0) rx *= 2;
    do {
        cx = (lx + rx) / 2;
        df1 = df(T, P, grad, cx, eps);
        if (df1 > 0) {
            rx = cx;
        } else {
            lx = cx;
        }
    } while (fabs(df1) > eps);
    return cx;
}
 
int main()
{
    Triangle Ts[] = {
                        { {1, 3, 0}, {5, 6, 0}, {6, 10, 0} },       // 9.1231
                        { {38, 38, 0}, {54, 34, 0}, {66, 53, 0} },  // 38.7644
                        { {2, 7, 3}, {1, -2, 0}, {3, 1, -5} },      // 14.6261
                        { {20, 73, 0}, {17, -42, 0}, {-12, 51, 0} } // 135.5567
                    };
    Triangle T = Ts[3];
    double crutch = 10.0;
    T.A = T.A*crutch; T.B = T.B*crutch; T.C = T.C*crutch; // Первый костыль
    Point P = T.A;
    double eps = 0.0001, step = 0.0, len = 0.0, cnt = 0.0;
    do {
        len = f(T,P);
        step = getStep(T, P, eps);
        P = P - step * getGradient(T, P, eps);
        cnt++;
    } while (fabs(len - f(T,P)) > eps);
    T.A = T.A/crutch; T.B = T.B/crutch; T.C = T.C/crutch; P = P/crutch;  // Второй костыль
    std::cout << P.x << ' ' << P.y << ' ' << P.z << ' ' << f(T,P) << ' ' << cnt << std::endl;
}

#include <iostream>
#include <cmath>
 
#define D60 1.0471975511965977461542144610932 // 60° в радианах
 
struct Vector3D {
    double x, y, z;
};
 
struct Matrix {
    Vector3D A, B, C;
};
 
// квадрат числа
double sqr(const double x) {
    return x*x;
}
 
// формирование матрицы поворота
Matrix getRotationMatrix(const Vector3D &V, double angle) {
    double sina = sin(angle);
    double cosa = cos(angle);
    return {
                { cosa + (1 - cosa)*sqr(V.x), (1 - cosa)*V.x*V.y - sina*V.z, (1 - cosa)*V.x*V.z + sina*V.y },
                { (1 - cosa)*V.x*V.y + sina*V.z, cosa + (1 - cosa)*sqr(V.y), (1 - cosa)*V.y*V.z - sina*V.x },
                { (1 - cosa)*V.x*V.z - sina*V.y, (1 - cosa)*V.y*V.z + sina*V.x, cosa + (1 - cosa)*sqr(V.z) }
           };
}
 
// транспонирование матрицы. Решил поразвлекаться и запихал это в оператор
Matrix operator! (const Matrix &M) {
    return { {M.A.x, M.B.x, M.C.x},
             {M.A.y, M.B.y, M.C.y},
             {M.A.z, M.B.z, M.C.z} };
}
 
// векторное произведение
Vector3D operator* (const Vector3D &A, const Vector3D &B) {
    return {A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x};
}
 
// произведение матрицы 3х3 и вектора
Vector3D operator* (const Matrix &M, const Vector3D &V) {
    return { M.A.x * V.x + M.A.y * V.y + M.A.z * V.z,
             M.B.x * V.x + M.B.y * V.y + M.B.z * V.z,
             M.C.x * V.x + M.C.y * V.y + M.C.z * V.z };
}
 
// частное вектора и числа
Vector3D operator/ (const Vector3D &A, const double x) {
    return {A.x / x, A.y / x, A.z / x};
}
 
// сумма двух векторов
Vector3D operator+ (const Vector3D &A, const Vector3D &B) {
    return {A.x + B.x, A.y + B.y, A.z + B.z};
}
 
// разность двух векторов
Vector3D operator- (const Vector3D &A, const Vector3D &B) {
    return {A.x - B.x, A.y - B.y, A.z - B.z};
}
 
// вывод вектора
std::ostream& operator << (std::ostream& out, const Vector3D V) {
    return out << V.x << "," << V.y << "," << V.z;
}
 
// расстояние между двумя точками
double getDistance(const Vector3D &A, const Vector3D &B) {
    return std::hypot(A.x - B.x, A.y - B.y, A.z - B.z);
}
 
// модуль вектора
double getDistance(const Vector3D &A) {
    return std::hypot(A.x, A.y, A.z);
}
 
// возвращает единичный вектор, путём деления исходного на модуль. Возможно, нужно назвать getOrt
Vector3D getOrth(const Vector3D &V) {
    return V / getDistance(V);
}
 
// вычисляем точку пересечения прямых AB и CD. Если прямые скрещиваются, то всё равно что-то вернёт. До конца не въехал ещё
// единственная функция, реализацию которой слямзил где-то на форуме. Перекурочил нулевые знаменатели. Вроде работает
Vector3D getIntersection(const Vector3D &A, const Vector3D &B, const Vector3D &C, const Vector3D &D){
    Vector3D P = B - A;
    Vector3D Q = D - C;
    return {
               ((P.y * Q.x - Q.y * P.x)?(A.x * P.y * Q.x - C.x * Q.y * P.x - A.y * P.x * Q.x + C.y * P.x * Q.x) /
               (P.y * Q.x - Q.y * P.x):0),
               ((P.x * Q.y - Q.x * P.y)?(A.y * P.x * Q.y - C.y * Q.x * P.y - A.x * P.y * Q.y + C.x * P.y * Q.y) /
               (P.x * Q.y - Q.x * P.y):0),
               ((P.y * Q.z - Q.y * P.z)?(A.z * P.y * Q.z - C.z * Q.y * P.z - A.y * P.z * Q.z + C.y * P.z * Q.z) /
               (P.y * Q.z - Q.y * P.z):0)
            };
}
 
int main()
{
    Vector3D P[] = { {2, 7, 3}, {1, -2, 0}, {3, 1, -5} };                  // массив вершин трегольника
    //std::cout << P[0] << ' ' << P[1] << ' ' << P[2] << ' ' << P[0] << "  \n";    // вывод для AutoCAD'a
    double dist[] = {0, 0, 0};                                                   // массив длин сторон, лежащих напротив вершин
    for (unsigned i = 0; i < 3; i++) {
        dist[i] = getDistance(P[(i + 1) % 3], P[(i + 2) % 3]);                   // вычисление массива длин
    }
    Vector3D R = {0, 0, 0};                                                      // вектор-результат
    
    // сначала ищем углы больше 120°. Перебираем вершины, по теореме косинусов находим угол. Если угол при вершине больше 120°, то это ответ
    bool find = false;
    for (unsigned i = 0; i < 3; i++) {
        if ((sqr(dist[(i + 1) % 3]) + sqr(dist[(i + 2) % 3]) - sqr(dist[i])) / (dist[(i + 1) % 3] * dist[(i + 2) % 3]) <= -1.0) {
            find = true;
            R = P[i];
            break;
        }
    }
    
    // если ничего не нашли, то ищем точку Ферма внутри треугольника
    if (!find) {
        Vector3D T1, T2, A_, B_;
        Vector3D orth = getOrth((P[1] - P[0]) * (P[2] - P[0])); // единичный вектор-нормаль к плоскости треугольника
        Matrix MR = getRotationMatrix(orth, D60);               // матрица поворота относительно нормали на 60° D60 = M_PI/3
        Matrix ML = !MR; // транспонирование                    // она же, но в другую сторону
 
        T1 = MR * (P[1] - P[0]) + P[0];                                // крутим первую точку относительно нулевой вправо
        T2 = ML * (P[1] - P[0]) + P[0];                                // и влево. Крутим в плоскости треугольника
        // Выбираем ту из двух, которая дальше от второй. Та, что дальше — выкручена наружу треугольника
        A_ = ((getDistance( T1, P[2] ) > getDistance( T2, P[2] )) ? T1 : T2 );
        
        // аналогично крутим вторую относительно первой. Оставляем ту, что дальше от нулевой
        T1 = MR * (P[2] - P[1]) + P[1];
        T2 = ML * (P[2] - P[1]) + P[1];
        B_ = ((getDistance( T1, P[0] ) > getDistance( T2, P[0] )) ? T1 : T2 );
        R = getIntersection(P[2], A_, P[0], B_); // вычисляем результат как пересечение отрезков [P[2], A_] и [P[0], B_]
    }
    
    std::cout << '(' << R << ')' << '\t' << getDistance(P[0],R) + getDistance(P[1],R) + getDistance(P[2],R);
}