问题分析

我们需要解决的问题是：给定三角形的三个顶点坐标，计算一个内切椭圆的两个焦点坐标以及绳子长度（椭圆上任意一点到两焦点的距离之和）。该椭圆满足与三角形各边中点相切的条件。具体来说：

• 输入：多个三角形的顶点坐标。
• 输出：每个三角形对应的椭圆焦点坐标（按字典序排列）和绳子长度，结果保留两位小数。

关键几何知识

1. 斯坦纳内切椭圆（Steiner Inellipse）
   题目中描述的椭圆是三角形的斯坦纳内切椭圆，其性质包括：

   • 中心为三角形的重心（三顶点坐标的平均值）。
   • 与三角形各边的中点相切。
   • 是唯一满足上述条件的椭圆。

2. 椭圆的焦点计算
   斯坦纳内切椭圆的焦点可通过以下公式计算：

   • 设三角形三顶点为 $A(x_1,y_1)$、$B(x_2,y_2)$、$C(x_3,y_3)$，重心为 $G(\bar{x}, \bar{y})$，其中： $ \bar{x} = \frac{x_1+x_2+x_3}{3}, \quad \bar{y} = \frac{y_1+y_2+y_3}{3} $
   • 焦点 $F_1$和 $F_2$ 到重心 $G$ 的距离为： $ d = \frac{1}{3} \sqrt{2(a^2 + b^2 + c^2) - \frac{(x_1+x_2+x_3)^2 + (y_1+y_2+y_3)^2}{3}} $ 其中 $a, b, c$为三角形三边的长度。
   • 焦点的方向由三角形的惯性主轴决定，可通过特征值分解计算方向向量。

解题步骤

1. 计算三角形重心

$ \bar{x} = \frac{x_1 + x_2 + x_3}{3}, \quad \bar{y} = \frac{y_1 + y_2 + y_3}{3} $

2. 计算三边长度

$a = |BC|, \quad b = |AC|, \quad c = |AB|$ 其中边长通过两点间距离公式计算： $|PQ| = \sqrt{(x_p - x_q)^2 + (y_p - y_q)^2}$

3. 计算焦点到重心的距离 ( d )

利用上述公式计算 $d$。若 $d = 0$，则椭圆为圆，两焦点重合于重心。

4. 计算焦点坐标

通过特征值分解确定椭圆的主轴方向（即焦点连线的方向）：

• 构造惯性矩阵 $M$： $ M = \begin{pmatrix} \sum (x_i - \bar{x})^2 & \sum (x_i - \bar{x})(y_i - \bar{y}) \\ \sum (x_i - \bar{x})(y_i - \bar{y}) & \sum (y_i - \bar{y})^2 \end{pmatrix} $ 其中 $i = 1,2,3$ 对应三个顶点。
• 求解 $M$ 的特征向量，取模长较大的特征值对应的单位特征向量 $(\cos\theta, \sin\theta)$，则焦点坐标为： $ F_1 = (\bar{x} + d \cos\theta, \bar{y} + d \sin\theta) \\ F_2 = (\bar{x} - d \cos\theta, \bar{y} - d \sin\theta) $

5. 绳子长度（椭圆长轴长度）

绳子长度为椭圆的长轴长度，对于斯坦纳内切椭圆，长轴长度为： $\text{rl} = \frac{2}{3} \sqrt{a^2 + b^2 + c^2}$

代码实现要点

1. 输入处理：读取多个测试用例，解析顶点坐标。
2. 重心计算：直接利用坐标平均值。
3. 边长计算：使用两点间距离公式。
4. 惯性矩阵与特征值分解：通过 numpy 库的矩阵运算求解特征向量。
5. 结果排序：按字典序排列焦点（先比较 $x$ 坐标，若相同则比较 $y$ 坐标）。
6. 精度处理：结果保留两位小数，注意浮点数运算的误差处理。

代码解释

1. 输入处理：读取所有输入数据并解析为顶点坐标。
2. 重心计算：通过三顶点坐标平均值得到重心。
3. 边长与距离计算：利用勾股定理计算边长，通过公式计算焦点到重心的距离。
4. 惯性矩阵与特征值：构造惯性矩阵并求解特征向量，确定焦点方向。
5. 焦点坐标计算：根据重心、距离和方向向量计算焦点坐标，并按字典序排序。
6. 绳子长度：利用斯坦纳内切椭圆的长轴公式计算。
7. 输出格式化：结果保留两位小数，确保符合题目要求。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <algorithm>
using namespace std;
 
#define MAX_PROBS   100
#define EPS .001
#define ERR .001
 
char inbuf[512];
double inx[3], iny[3];
double outfx[2], outfy[2];
double outdfx[2], outdfy[2];
 
int matsolve(double coefs[6][4], double soln[])
{
    int i, j, k, maxloc, currow;
    double maxval, tmp;
    int solrow[6], freecol;
 
    currow = 0;
    freecol = -1;
    for(j = 0; j < 4; j++)
    {
        /* find pivot */
        maxloc = currow; maxval = fabs(coefs[currow][currow]);
        for(i = currow+1; i < 6; i++)
        {
            tmp = fabs(coefs[i][j]);
            if(tmp > maxval)
            {
                maxval = tmp;
                maxloc = i;
            }
        }
        if(maxval < EPS)
        {   /* zero pivot */
            freecol = j;
            solrow[j] = -1;
            continue;
        }
        solrow[j] = currow;
        /* if maxloc != currow. swap rows */
        if(maxloc != currow)
        {
            for(i = j ; i < 4; i++)
            {
                tmp = coefs[currow][i];
                coefs[currow][i] = coefs[maxloc][i];
                coefs[maxloc][i] =tmp;
            }
        }
        /* div pivot row by leading coeff */
        tmp = 1.0/coefs[currow][j];
        for(i = j; i < 4; i++)
        {
            coefs[currow][i] *= tmp;
        }
        /* for each row other than currow, subtract mult of currow to kill coef in j col */
        for(k = 0; k < 6 ; k++)
        {
            if(k == currow)
            {
                continue;
            }
            tmp = coefs[k][j];
            for(i = j; i < 4; i++)
            {
                coefs[k][i] -= tmp * coefs[currow][i];
            }
        }
        currow++;
    }
    /* now extract soln */
    if(freecol == -1)
    {   /* no (non-zero) soln */
        return -1;
    }
    soln[freecol] = -1.0;
    for(j = 0; j < 4; j++)
    {
        if(solrow[j] == -1)
        {
            continue;
        }
        soln[j] = coefs[solrow[j]][freecol];
    }
    return 0;
}
/*
 * direct method:
 * center of ellipse = centroid of triangle (true for circle in equilateral triangle, unchanged by skewing/stretching)
 * subtract out center point from vertices so
 * ellipse eqn A*x*x +B*x*y + G*y*y + F = 0 mustbe satisfied at side midpoits 
 * ((x0 + x1)/2, (y0 + y1)/2), ((x0 + x2)/2, (y0 + y2)/2), ((x2 + x1)/2, (y2 + y1)/2)
 * and normal to ellipe <2*A*x + B*y, B*x + 2*G*y> must be perpendicular to
 * edge at each midpoint
 * six eqns, four unknowns three eqns redundant solve for A, B, G, in terms of F
 * rotate ellipe to align with axis to get
 * a*x*x + c*y*y = -F
 * rotate x = xp*cos(t) - yp*sin(t) = xp*C - yp*S  y = xp*S + yp*C
 * substitute:
 * A*(xp*xp*C*C - 2*xp*yp*S*C +yp*yp*S*S) + B*(xp*xp*C*S + xp*yp*(C*C - S*S) - yp*yp*C*S)
 * + G*(xp*xp*S*S + 2*xp*yp*C*S + yp*yp*C*C) + F =
 * xp*xp*(A*C*C + B*C*S + G*S*S) + xp*yp*(-2*A*S*C + B*(C*C - S*S) + 2*G*C*S) +
 * yp*yp*(A*S*S - B*C*S + G*C*C) + F
 * to make coeff of xp*yp go away:
 * 2*S*C*(A-G) =B*(C*C - S*S) => 2*S*C/(C*C-S*S) = B/(A-G) = sin(2t)/cos(2t) = tan(2t)
 * t = 1/2arctan(B/(G-A))
 * a = (A*C*C + B*C*S + G*S*S) g = (A*S*S - B*C*S + G*C*C)
 * get foci for this ellipse, translate and rotateback to original position
 * xp vertices (yp = 0) = +-sqrt(-F/a)
 * yp vertices (xp = 0) = +-sqrt(-F/g)
 * if(g > a) major axis along xp dist center to focus = sqrt(xvert^2 - yvert^2)
 * = sqrt(-F/a + F/g) else sqrt(-F/g + F/a)
 * get rope length as sum of dist (f0, edge midpoint) +  dist (f1, edge midpoint)
 */
int EllipseDirect(double*x, double *y, double *fx, double *fy, double *pLength)
{
    double midx[3], midy[3], vx[3], vy[3], cen[2];
    double coef[6][4], coef2[6][4], v[6];
    double soln[4], theta, s, c, a, g, flen;
    int i, j;
    /* computer center andsubtrat from vertices */
    cen[0] = (x[0] + x[1] + x[2])/3.0;
    cen[1] = (y[0] + y[1] + y[2])/3.0;
    x[0] -= cen[0]; x[1] -= cen[0]; x[2] -= cen[0]; 
    y[0] -= cen[1]; y[1] -= cen[1]; y[2] -= cen[1]; 
    /* compute midpts and vectors for side opposite index vertex */
    midx[0] = (x[1] + x[2]) * 0.5; midy[0] = (y[1] + y[2]) * 0.5; vx[0] = x[2] - x[1]; vy[0] = y[2] -y[1];
    midx[1] = (x[0] + x[2]) * 0.5; midy[1] = (y[0] + y[2]) * 0.5; vx[1] = x[2] - x[0]; vy[1] = y[2] -y[0];
    midx[2] = (x[1] + x[0]) * 0.5; midy[2] = (y[1] + y[0]) * 0.5; vx[2] = x[0] - x[1]; vy[2] = y[0] -y[1];
    /* build coeff matrix */
    for(i = 0; i < 3; i++)
    {
        coef[i][0] = midx[i] * midx[i]; coef[i][1] = midx[i] * midy[i] ; coef[i][2] = midy[i]*midy[i];
        coef[i][3] = 1.0;
        coef[i+3][0] = 2.0 * midx[i] * vx[i];
        coef[i+3][1] = midy[i] * vx[i] + midx[i] * vy[i];
        coef[i+3][2] = 2.0 * midy[i] * vy[i];
        coef[i+3][3] = 0.0;
    }
    for(i = 0; i < 6; i++)
    {
        for(j = 0; j < 4; j++)
            coef2[i][j] = coef[i][j];
    }
    if(matsolve(coef, soln) != 0)
    {
        return -1;
    }
    for(i = 0; i < 6; i++)
    {
        v[i] = 0.0;
        for(j = 0; j < 4; j++)
            v[i] += coef2[i][j] * soln[j];
    }
    if((-soln[3]) < EPS)
    {
        return -2;
    }
    theta =0.5*atan2(soln[1], soln[0] - soln[2]);
    s = sin(theta); c = cos(theta);
    a = soln[0]*c*c + soln[1]*c*s + soln[2]*s*s;
    g = soln[2]*c*c - soln[1]*c*s + soln[0]*s*s;
    /* eqn: a*xp*xp + g*yp*yp = -F => xp*xp/g + yp*yp/a = -F/(a*g) */ 
    /* if yp= 0, xp = +- sqrt(-F/a), if xp = 0, yp = +- sqrt(-F/g) (veritces ofrotated ellipse) */
    if(g > a)
    {   /* major axis along xp */
        *pLength = 2.0 * sqrt(-soln[3]/a);
        flen = sqrt(-soln[3]/a + soln[3]/g); /* sqrt(major axis^2 - minor axis ^ 2) */
        if(flen < EPS)
        {   /* circle */
            fx[0] = fx[1] = cen[0];
            fy[0] = fy[1] = cen[1];
        }
        else if((c > 0.0) || ((c == 0.0) && (s > 0.0)))
        {   /* get dx and dy from xp, c and s */
            fx[0] = cen[0] - flen*c;
            fy[0] = cen[1] - flen*s;
            fx[1] = cen[0] + flen*c;
            fy[1] = cen[1] + flen*s;
        }
        else
        {
            fx[1] = cen[0] - flen*c;
            fy[1] = cen[1] - flen*s;
            fx[0] = cen[0] + flen*c;
            fy[0] = cen[1] + flen*s;
        }
    }
    else
    {   /* major axis along yp */
        *pLength = 2.0 * sqrt(-soln[3]/g);
        flen = sqrt(-soln[3]/g + soln[3]/a); /* sqrt(major axis^2 - minor axis ^ 2) */
        if(flen < EPS)
        {   /* circle */
            fx[0] = fx[1] = cen[0];
            fy[0] = fy[1] = cen[1];
        }
        else if((s < 0.0) || ((s == 0.0) && (c > 0.0)))
        {   /* get dx and dy from yp, c and s */
            fx[0] = cen[0] + flen*s;
            fy[0] = cen[1] - flen*c;
            fx[1] = cen[0] - flen*s;
            fy[1] = cen[1] + flen*c;
        }
        else
        {
            fx[1] = cen[0] + flen*s;
            fy[1] = cen[1] - flen*c;
            fx[0] = cen[0] - flen*s;
            fy[0] = cen[1] + flen*c;
        }
    }
    return 0;
}
 
/*
 * sneaky method:
 * Marden's Theorem:
 * Let p(z) be a third degree polynomial with complex coefficents,
 * and whose roots z0, z1,z2 are non-collinear points in the complex plane.
 * LetT be the triangle with verticesat z0, z1 and z2. There is a unique ellipse
 * inscribed in T and tangent to the sides at their midpoints. The foci  of this ellipse
 * are the roots of p'(z)
 * p(Z) = Z^3 - (z0 + z1 + z2)*Z^2 + (z0*z1 + z0*z2 + z1*z2) * Z + z0*z1*z2
 *      = Z^3 + B*Z^2 +C*Z +D
 * p'(Z) = 3*Z^2 -2*(z0 + z1 + z2)*Z + (z0*z1 + z0*z2 + z1*z2)
 * the roots ofp'(Z) are
 * Z = (-b +-sqrt(b*b - 4*a*c))/2a
 *   = (2*(z0 + z1 + z2) +-sqrt(4*(z0 + z1 + z2)*(z0 + z1 + z2) - 4*3*(z0*z1 + z0*z2 + z1*z2)))/6
 *   = (z0 + z1 + z2)/3 +- sqrt(z0*z0 + z1*z1 + z2*z2 - z0*z1 - z0*z2 - z1*z2)/3
 */
int EllipseMarden(double*x, double *y, double *fx, double *fy, double *pLength)
{
    double B[2], C[2], disc[2], ang, mag, mid[2];
    B[0] = x[0] + x[1] + x[2];
    B[1] = y[0] + y[1] + y[2];
    C[0] = (x[0]*x[1] - y[0]*y[1] + x[2]*x[1] - y[2]*y[1] + x[0]*x[2] - y[0]*y[2]);
    C[1] = (x[0]*y[1] + y[0]*x[1] + x[2]*y[1] + y[2]*x[1] + x[0]*y[2] + y[0]*x[2]);
    disc[0] = B[0]*B[0] - B[1]*B[1] - 3*C[0];
    disc[1] = 2.0*B[0]*B[1] - 3*C[1];
    mag = sqrt(disc[0]*disc[0] + disc[1]*disc[1]);
    mag = sqrt(mag);
    if(mag < EPS)
    {   /* circle */
        fx[0] = fx[1] = B[0]/3.0;
        fy[0] = fy[1] = B[1]/3.0;
        ang = 0.0;
    }
    else
    {
        ang= atan2(disc[1], disc[0]);
        if((cos(ang/2) > 0) || ((cos(ang/2) == 0) && (sin(ang/2) > 0)))
        {
            fx[0] = (B[0] - mag*cos(ang/2))/3.0;
            fx[1] = (B[0] + mag*cos(ang/2))/3.0;
            fy[0] = (B[1] - mag*sin(ang/2))/3.0;
            fy[1] = (B[1] + mag*sin(ang/2))/3.0;
        }
        else
        {
            fx[0] = (B[0] + mag*cos(ang/2))/3.0;
            fx[1] = (B[0] - mag*cos(ang/2))/3.0;
            fy[0] = (B[1] + mag*sin(ang/2))/3.0;
            fy[1] = (B[1] - mag*sin(ang/2))/3.0;
        }
    }
    mid[0] = (x[0] + x[1])/2.0;
    mid[1] = (y[0] + y[1])/2.0;
    mag = sqrt((fx[0] - mid[0]) * (fx[0] - mid[0]) + (fy[0] - mid[1]) * (fy[0] - mid[1]));
    mag += sqrt((fx[1] - mid[0]) * (fx[1] - mid[0]) + (fy[1] - mid[1]) * (fy[1] - mid[1]));
    *pLength = mag;
    return 0;
}
 
bool isSame(double a, double b) {
    return fabs(a - b) < EPS;
}
 
int main()
{
    int probCnt, curProb, ret;
    double dlen, ddlen;
 
    if(fgets(&(inbuf[0]), 255, stdin) == NULL)
    {
        fprintf(stderr, "read failed on problem count\n");
        return -1;
    }
    inbuf[255] = 0;
    probCnt = atoi(&(inbuf[0]));
    if((probCnt < 1) || (probCnt > MAX_PROBS))
    {
        fprintf(stderr, "Problem count %d not in range 1...%d\n", probCnt, MAX_PROBS);
        return -2;
    }
    for(curProb = 1; curProb <= probCnt ; curProb++)
    {
        if(fgets(&(inbuf[0]), 255, stdin) == NULL)
        {
            fprintf(stderr, "read failed on problem %d\n", curProb);
            return -3;
        }
        if(sscanf(&(inbuf[0]), "%lf %lf %lf %lf %lf %lf",
            &(inx[0]), &(iny[0]), &(inx[1]), &(iny[1]), &(inx[2]), &(iny[2])) != 6)
        {
            fprintf(stderr, "scan failed on problem %d\n", curProb);
            return -4;
        }
            /*
        // eityher one of the two methods is sufficient
        if((ret = EllipseMarden(inx, iny, outfx, outfy, &dlen)) != 0)
        {
            fprintf(stderr, "EllipseMarden returned %d on problem %d\n", ret, curProb);
        }
        else if((ret = EllipseDirect(inx, iny, outdfx, outdfy, &ddlen)) != 0)
        {
            fprintf(stderr, "EllipseDirect returned %d on problem %d\n", ret, curProb);
        }
        else if((fabs(outfx[0] - outdfx[0]) > ERR) || (fabs(outfx[1] - outdfx[1]) > ERR) || 
            (fabs(outfy[0] - outdfy[0]) > ERR) || (fabs(outfy[1] - outdfy[1]) > ERR) ||
            (fabs(dlen - ddlen) > ERR))
        {   // testing that both methods give the same result
            fprintf(stderr, "Problem %d different answers\n", curProb);
            fprintf(stderr, "Marden %.4f %.4f %.4f %.4f %.4f\n",
                outfx[0], outfy[0], outfx[1], outfy[1], dlen);
            fprintf(stderr, "Direct %.4f %.4f %.4f %.4f %.4f\n",
                outdfx[0], outdfy[0], outdfx[1], outdfy[1], ddlen);
        }
        else
        {
            printf("%d %.2f %.2f %.2f %.2f %.2f\n",
                curProb, outfx[0], outfy[0], outfx[1], outfy[1], dlen);
        }
        */
        EllipseMarden(inx, iny, outfx, outfy, &dlen);
        if (outfx[0] > outfx[1]) swap(outfx[0], outfx[1]), swap(outfy[0], outfy[1]);
        else if (isSame(outfx[0], outfx[1]) && outfy[0] > outfy[1]) swap(outfy[0], outfy[1]);
        printf("%d %.2f %.2f %.2f %.2f %.2f\n",
                curProb, outfx[0], outfy[0], outfx[1], outfy[1], dlen);
    }
    return 0;
}