Ellipsoid HDU - 5017 [模拟退火] 有个不知名BUG

Ellipsoid HDU - 5017 [模拟退火] 

题意:求一三维空间中,椭球面上到(0,0,0) 最近的一点

思路:模拟退火做,这题对Rate要求设置为0.99,精度要求高

但是发现有个BUG,不能加区域限制,不知道什么情况.

#include<cstdio>
#include<vector>
#include<cmath>
#include<math.h>
#include<time.h>
#include<string>
#include<string.h>
#include<iostream>
#include<algorithm>
#include<map>
#define PI acos(-1.0)
#define pb push_back
#define F first
#define S second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=1005;
const int MOD=1e9+7;
const double eps=1e-10;
int sign(double x) {  //三态函数，减少精度问题
    return abs(x) < eps ? 0 : x < 0 ? -1 : 1;
}

double a,b,c,d,e,f;
int n;
int dir[8][2]={0,-1,0,1,1,0,-1,0,1,1,1,-1,-1,1,-1,-1};
double cul(double x,double y,double z){
    return sqrt(x*x+y*y+z*z);
}
double myrand(){
    return rand()%10000/10000.0;
}
double get_nz(double x,double y){
    double A=c;
    double B=d*y+e*x;
    double C=-1.0*(1-f*x*y-a*x*x-b*y*y);
    double delta=B*B-4*A*C;
    if(sign(delta)<0)   return (double)(MOD+10);
    double z1=(-B-sqrt(delta))/2/A;
    double z2=(-B+sqrt(delta))/2/A;
    return sign(cul(x,y,z1)-cul(x,y,z2))<0?z1:z2;
}
void SA(double &res){
    double T=1;
    double x=0,y=0,z=sqrt(1.0/c);
    double E=cul(x,y,z);
//    printf("stz=%f\n",z);
    while(T>eps){
        double mn=1e20;
        double tx,ty,tz;
        for(int i=0;i<8;i++){
            double nx=x+dir[i][0]*T;
            double ny=y+dir[i][1]*T;
            ///这里不能加
//            if(sign(nx-sqrt(1.0/a))>0)  nx=sqrt(1.0/a);
//            if(sign(nx+sqrt(1.0/a))<0)  nx=-sqrt(1.0/a);
//            if(sign(ny-sqrt(1.0/b))>0)  ny=sqrt(1.0/b);
//            if(sign(ny+sqrt(1.0/b))>0)  ny=-sqrt(1.0/b);
            double nz=get_nz(nx,ny);
            if(nz>=MOD)    continue;
            double tE=cul(nx,ny,nz);
//            printf("tE=%f\n",nz);
            if(tE < mn){
                mn=tE;
                tx=nx,ty=ny,tz=nz;
            }
        }
        if(sign(mn-E)<0 || (mn-E)/T<myrand()){
            E=mn;
            x=tx,y=ty,z=tz;
        }
        T*=0.99;//要成0.99
    }
    res=E;
}
void mian(){
    double res;
    SA(res);
    printf("%.10f\n",res);
}
int main(void){
    srand(time(0));
    while(scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f)!=EOF){
        mian();
    }

    return 0;
}