hdu5017 Ellipsoid (模拟退火)

时间:2023-03-08 23:36:17
hdu5017 Ellipsoid (模拟退火)

Ellipsoid

原题链接

题目描述

  • 给定。一个要满足的椭球的方程\(ax^2+by^2+cz^2+dyz+exz+fxy=1\)
  • 求球面上一个点到原点\((0,0,0)\)的距离最小。
  • 有多组输入数据

解题思路

  • 这题有一个麻烦的地方: 如何求满足椭球方程的点!
  • 那就只能用自己浅薄的知识:求根公式。
  • 我们先随机x坐标和y坐标,再带入到方程中,用求根公式算出z坐标。
  • 这样就能保证所扩展的新状态一定满足条件。
  • 之后的更新就和普通的模拟退火一样了。
  • 具体求根方法,按照程序来看吧。反正就是移个项之后带公式。

代码

#include<bits/stdc++.h>
using namespace std;
inline void read(int &x)
{
x=0;
static int p;p=1;
static char c;c=getchar();
while(!isdigit(c)){if(c=='-')p=-1;c=getchar();}
while(isdigit(c)) {x=(x<<1)+(x<<3)+(c-48);c=getchar();}
x*=p;
}
double a,b,c,d,e,f;
double ansx,ansy,ansz;
const double eps=1e-10;
double dis(double x,double y,double z)
{
return sqrt(x*x+y*y+z*z);
}
double calc(double x,double y)
{
double A=c;
double B=d*y+e*x;
double C=a*x*x+b*y*y+f*x*y-1.0;
double delta=B*B-4.0*A*C;
if(delta<0)return 210000000.0;
double x1=(-B+sqrt(delta))/(2.0*A);
double x2=(-B-sqrt(delta))/(2.0*A);
if(dis(x,y,x1)>dis(x,y,x2))return x2;
return x1;
}
void MNTH()
{
for(int times=1;times<=1;times++)
{
double T=10000;
while(T>eps)
{
double nowx=ansx+(rand()*2-RAND_MAX)*T;
double nowy=ansy+(rand()*2-RAND_MAX)*T;
double nowz=calc(nowx,nowy);
if(nowz==210000000.0){T*=0.99;continue;}
double delta=dis(nowx,nowy,nowz)-dis(ansx,ansy,ansz);
if(delta<0)ansx=nowx,ansy=nowy,ansz=nowz;
else if(exp(delta/T)*RAND_MAX<rand())ansx=nowx,ansy=nowy,ansz=nowz;
T*=0.99;
}
}
printf("%.6lf\n",dis(ansx,ansy,ansz));
}
int main()
{
srand(19890604);
while(scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f)!=EOF)
{
ansx=0;
ansy=0;
ansz=sqrt(1.0/c);
MNTH();
}
return 0;
}