HDU5017(2014年西安网络赛K题)

第三场网络赛又水掉了……

说下这题吧。题意很简单，就是求椭球面a x^2 + b y^2 + c z^2 +d yz + e zx + f xy = 1到原点距离的最小值。比赛开始时就看到了这题，结果开始的思路完全错浪费了好多时间，甚至还以为是精度问题。到了最后才终于想出做法，无奈没时间敲。赛后敲一下子就AC了，只能说略可惜……

思路其实就是几代里面学的非退化线性变换。首先可以写出这个二次型的系数矩阵A，然后因为系数矩阵是对称的，所以一定存在一个正交矩阵C使得C^T*A*C等于一个对角阵，并且对角线上的元素就是A所有的特征值。这里的C事实上就代表了一个R^3空间的等距变换。加上如果d,e,f=0,那么这题的答案就是a,b,c三数的倒数中的最小值再开方，这样这题就转化成求A的最大特征值了。因为没有直接求特征值的模板，所以我推了下那个特征方程，然后用高次方程求实根的模板写了这题。

#include<stdio.h>
#include<math.h>
#include<iostream>
#include<algorithm>
#include<stdlib.h>
#include<vector>
using namespace std;

const double eps=1e-12;
const double inf=1e+12;
inline int sign(double x)
{
    return x<-eps?-1:x>eps;
}
inline double get(const vector<double> &coef,double x)
{
    double e=1,s=0;
    for(int i=0;i<coef.size();++i) s+=coef[i]*e,e*=x;
    return s;
}
double find(const vector<double> &coef,int n,double lo,double hi)
{
    double sign_lo,sign_hi;
    if((sign_lo=sign(get(coef,lo)))==0) return lo;
    if((sign_hi=sign(get(coef,hi)))==0) return hi;
    if(sign_lo*sign_hi>0) return inf;
    for(int step=0;step<100&&hi-lo>eps;step++)
    {
        double m=(lo+hi)*.5;
        int sign_mid=sign(get(coef,m));
        if(sign_mid==0) return m;
        if(sign_lo*sign_mid<0) hi=m;
        else lo=m;
    }
    return (lo+hi)*.5;
}
vector<double> equation(vector<double> coef,int n)
{
    vector<double> ret;
    if(n==1)
    {
        if(sign(coef[1])) ret.push_back(-coef[0]/coef[1]);
        return ret;
    }
    vector<double> dcoef(n);
    for(int i=0;i<n;i++) dcoef[i]=coef[i+1]*(i+1);
    vector<double> droot=equation(dcoef,n-1);
    droot.insert(droot.begin(),-inf);
    droot.push_back(inf);
    for(int i=0;i+1<droot.size();i++)
    {
        double tmp=find(coef,n,droot[i],droot[i+1]);
        if(tmp<inf) ret.push_back(tmp);
    }
    return ret;
}

int main()
{
    double a,b,c,d,e,f;
    while(scanf("%lf %lf %lf %lf %lf %lf",&a,&b,&c,&d,&e,&f)!=EOF)
    {
        d/=2;e/=2;f/=2;
        vector<double> tet;
        tet.push_back(-a*b*c+e*e*b+f*f*c+d*d*a-2*d*e*f);
        tet.push_back(a*b+b*c+c*a-e*e-f*f-d*d);
        tet.push_back(-a-b-c);
        tet.push_back(1.0);
        vector<double> ret=equation(tet,3);
        double res=-1000;
        for(int i=0;i<ret.size();i++) res=max(ret[i],res);
        printf("%.8lf\n",sqrt(1.0/res));
    }
    return 0;
}