Ellipsoid

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others) Total Submission(s): 648    Accepted Submission(s): 194 Special Judge
Problem Description Given a 3-dimension ellipsoid(橢球面)

your task is to find the minimal distance between the original point (0,0,0) and points on the ellipsoid. The distance between two points (x₁,y₁,z₁) and (x₂,y₂,z₂) is defined as
Input There are multiple test cases. Please process till EOF.

For each testcase, one line contains 6 real number a,b,c(0 < a,b,c,< 1),d,e,f(0 ≤ d,e,f < 1), as described above.It is guaranteed that the input data forms a ellipsoid. All numbers are fit in double.

Output For each test contains one line. Describes the minimal distance. Answer will be considered as correct if their absolute error is less than 10^-5.
Sample Input 1 0.04 0.01 0 0 0
Sample Output 1.0000000模擬退火演算法：（Simulated Annealing，簡稱SA）是一種通用概率演算法，用來在一個大的搜尋空間內找尋命題的最優解。設：Y（x）是功能函式；要求最大值，x[i+1]是x[i]的下一個點，若Y(x[i+1])>=Y(x[i])，移動後可以得到更優解，總是接收移動；若Y(x[i+1])<Y(x[i])(即移動後的解比當前解要差)，則以一定的概率接受移動，而且這個概率隨著時間推移逐漸降低（逐漸降低才能趨向穩定） 下面給出模擬退火的虛擬碼：

/*
* J(y)：在狀態y時的評價函式值
* Y(i)：表示當前狀態
* Y(i+1)：表示新的狀態
* r： 用於控制降溫的快慢
* T： 系統的溫度，系統初始應該要處於一個高溫的狀態
* T_min ：溫度的下限，若溫度T達到T_min，則停止搜尋
*/
while( T > T_min )
{
　　dE = J( Y(i+1) ) - J( Y(i) ) ; 

　　if ( dE >=0 ) //表達移動後得到更優解，則總是接受移動
Y(i+1) = Y(i) ; //接受從Y(i)到Y(i+1)的移動
　　else
　　{
// 函式exp( dE/T )的取值範圍是(0,1) ，dE/T越大，則exp( dE/T )也
if ( exp( dE/T ) > random( 0 , 1 ) )
Y(i+1) = Y(i) ; //接受從Y(i)到Y(i+1)的移動
　　}
　　T = r * T ; //降溫退火 ，0<r<1 。r越大，降溫越慢；r越小，降溫越快
　　/*
　　* 若r過大，則搜尋到全域性最優解的可能會較高，但搜尋的過程也就較長。若r過小，則搜尋的過程會很快，但最終可能會達到一個區域性最優值
　　*/
　　i ++ ;
}
 

程式：

#include"string.h"
#include"stdio.h"
#include"queue"
#include"stack"
#include"vector"
#include"algorithm"
#include"iostream"
#include"math.h"
#include"stdlib.h"
#define M 100009
#define inf 100000
#define eps 1e-10
#define PI acos(-1.0)
using namespace std;
int disx[9]={0,0,1,-1,1,-1,1,-1};
int disy[9]={1,-1,0,0,1,1,-1,-1};
double a,b,c,d,e,f;
double fun(double x,double y,double z)
{
    return sqrt(x*x+y*y+z*z);
}
double cal(double x,double y)
{
    double A=c;
    double B=d*y+e*x;
    double C=f*x*y+a*x*x+b*y*y-1;
    double ff=B*B-4*A*C;
    if(ff<0.0)return inf*100;
    ff=sqrt(ff);
    double z1=(-B+ff)/(2*A);
    double z2=(-B-ff)/(2*A);
    if(fun(x,y,z1)<fun(x,y,z2))
        return z1;
    return z2;
}
double solve()
{
    double x=0,y=0,z=sqrt(1/c);
    double step=1,rate=0.99;
    while(step>eps)
    {
        for(int i=0;i<8;i++)
        {
            double xx=x+step*disx[i];
            double yy=y+step*disy[i];
            double zz=cal(xx,yy);
            if(zz>=inf*99)continue;
            if(fun(xx,yy,zz)<fun(x,y,z))
            {
                x=xx;
                y=yy;
                z=zz;
            }
        }
        step*=rate;
    }
    return fun(x,y,z);
}
int main()
{
    while(scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f)!=-1)
    {
        double ans=solve();
        printf("%.7lf\n",ans);
    }
    return 0;
}