hdu 5017 Ellipsoid

最新推荐文章于 2021-04-27 21:58:58 发布

原创最新推荐文章于 2021-04-27 21:58:58 发布 · 1.1k 阅读

本内容遵循CC 4.0 BY-SA版权协议

本文探讨了解决特定数学问题的方法，指出模拟退火算法虽然常用但效果不稳定，依赖于参数选择和运气。正确的解决方案应采用正交变换，通过矩阵变换将复杂问题简化为易于解决的形式。

这题大部分人是用模拟退火做的，但是参数的选择不能保证算法适合任何一种情况，比如下面的代码如果把step改成100，样例通不过却能ac。。。因为样例的那种情况收敛的点特别不 “平坦",所以从不同方向去逼近得到的结果相差很大，两种方法可以改进，一是增加r，二是增加方向的数量，不过最关键的ac或wa还是看脸，模拟退火对于这题不是一个很好的方法。

这题的正确做法应该是正交变换，正交变换不会改变一个向量的内积，把二次型方程变为标准型，解一个一元三次方程求出特征值，然后求特征向量，这是三维的情况，如果维数很高就用QR分解求特征值，然后高斯消元求特征向量。

#pragma comment(linker, "/STACK:102400000,102400000")
#include<iostream>
#include<vector>
#include<algorithm>
#include<cstdio>
#include<queue>
#include<stack>
#include<string>
#include<map>
#include<set>
#include<cmath>
#include<cassert>
#include<cstring>
#include<iomanip>
using namespace std;
#ifdef _WIN32
#define i64 __int64
#define out64 "%I64d\len"
#define in64 "%I64d"
#else
#define i64 long long
#define out64 "%lld\len"
#define in64 "%lld"
#endif
/************ for topcoder by zz1215 *******************/
#define foreach(c,itr)  for(__typeof((c).begin()) itr=(c).begin();itr!=(c).end();itr++)
#define FOR(i,a,b)      for( int i = (a) ; i <= (b) ; i ++)
#define FF(i,a)         for( int i = 0 ; i < (a) ; i ++)
#define FFD(i,a,b)      for( int i = (a) ; i >= (b) ; i --)
#define S64(a)          scanf(in64,&a)
#define SS(a)           scanf("%d",&a)
#define LL(a)           ((a)<<1)
#define RR(a)           (((a)<<1)+1)
#define pb              push_back
#define pf              push_front
#define X               first
#define Y               second
#define CL(Q)           while(!Q.empty())Q.pop()
#define MM(name,what)   memset(name,what,sizeof(name))
#define MC(a,b)			memcpy(a,b,sizeof(b))
#define MAX(a,b)        ((a)>(b)?(a):(b))
#define MIN(a,b)        ((a)<(b)?(a):(b))
#define read            freopen("out.txt","r",stdin)
#define write           freopen("out2.txt","w",stdout)

const int inf = 0x3f3f3f3f;
const i64 inf64 = 0x3f3f3f3f3f3f3f3fLL;
const double oo = 10e9;
const double eps = 1e-13;
const double pi = acos(-1.0);

double dir[360][2];
double a, b, c, d, e, f;

double dis(double x, double y, double z){
	return sqrt(x*x + y*y + z*z);
}

double getz(double x, double y) {
	double A = c, B = e * x + d * y,
		C = a * x * x + b * y * y + f * x * y - 1;
	double delta = B * B - 4 * A * C;
	if (delta < 0) return 1e60;
	double z1 = (-B + sqrt(delta)) / 2 / A,
		z2 = (-B - sqrt(delta)) / 2 / A;
	if (z1 * z1 < z2 * z2) return z1;
	else return z2;
}

double start(){
	double r = 0.99;
	double step = 1;
	double nowx=0,nowy=0,nowz=getz(nowx,nowy);
	double x, y, z;
	double temp;
	while (step > eps){
		temp = oo;
		for (int i = 0; i < 12; i++){
			x = nowx + dir[i][0] * step;
			y = nowy + dir[i][1] * step;
			z = getz(x, y);
			if (z >= oo){
				continue;
			} 
			if (dis(x, y, z) < dis(nowx,nowy,nowz)){
				nowx = x;
				nowy = y;
				nowz = z;
			}
		}
		step *= r;
	}
	return dis(nowx, nowy, nowz);
}

int main(){
	double sita = 2.0*pi / 12;
	for (int i = 0; i < 12; i++){
		dir[i][0] = cos(sita*i);
		dir[i][1] = sin(sita*i);
	}
	while (cin >> a >> b >> c >> d >> e >> f){
		printf("%.10lf\n", start());
	}
	return 0;
}