hdu5017 Ellipsoid(旋转)

比赛的时候跳进这个大坑里，最后代码是写出来了。看到好像很多都是模拟退火做的，下面提供一个奇怪的思路吧。

ax^2+by^2+cz^2+dyz+exz+fxy=1(*)

通过一些奇特的YY我们可以知道这是由一个标准的椭球ax^2+by^2+cz^2=1旋转得到的，之所以有交叉项是因为绕了X,Y,Z轴旋转，所以我们的想法是要是能将这个椭球旋回去那么就可以知道最短的距离了。

首先我们是旋转z轴，希望使得xy项消掉，这个时候令z=0可以得到  ax^2+by^2+fxy=1

通过一些椭圆的旋转知识我们可以求出当旋转角为phi,将x,y做一个旋转代换就可以使得xy被消掉，具体的phi值求法可以百度一下椭圆，公式这里就是tan(2p)=f/(a-b)（对于上面那个式子），将这个代换弄到＊里其实就会得到一个新的方程:

a'x^2+b'y^2+c'z^2+d'yz+e'xz+f'xy=1

这个时候f'=0，然后我们希望通过旋转x轴将d'=0，方法同样是令x=0得到

b'y^2+c'z^2+d'yz=1

利用相同的代换可以得到新的方程

a''x^2+b''y^2+c''z^2+d''yz+e''xz+f''xy=1

这个时候注意,d''=0，但原先的f''却不一定等于0了。

同样我们再令旋转y轴，使得e''=0。

经过一轮这样的操作后，我们可以发现，最后e＝0是一定的，但是d和f不一定等于0，但是直觉上这样的操作只要重复的做最后一定会收敛，因此我们就重复地对这个椭球做这样的操作，做个100次（实测其实只要转4，5次）精度就已经非常足够了。

#pragma warning(disable:4996)
#include <iostream>
#include <cstring>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <cmath>
#include <queue>
#include <string>
#include <map>
using namespace std;

#define eps 1e-6
int dcmp(double x){
	return (x > eps) - (x < -eps);
}

double a, b, c, d, e, f;
double p, w, k;

double pi = acos(-1.0);

double get(double A, double B, double C)
{
	if (dcmp(B) == 0) return 0;

	if (dcmp(A - C) == 0) {
		return pi / 4;
	}
	double ret = atan(B / (A - C));
	return ret / 2;
}

double a1, a2, a3, b1, b2, b3, c1, c2, c3;

double A, B, C, D, E, F;


int main()
{
	while (~scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e, &f))
	{
		int T = 100;
		while (T--){
			p = 0;
			w = 0;
			k = get(a, f, b);

			a1 = cos(p)*cos(k) - sin(p)*sin(w)*sin(k);
			a2 = -cos(p)*sin(k) - sin(p)*sin(w)*cos(k);
			a3 = -sin(p)*cos(w);

			b1 = cos(w)*sin(k);
			b2 = cos(w)*cos(k);
			b3 = -sin(w);

			c1 = sin(p)*cos(k) + cos(p)*sin(w)*sin(k);
			c2 = -sin(p)*sin(k) + cos(p)*sin(w)*cos(k);
			c3 = cos(p)*cos(w);

			A = a*a1*a1 + b*b1*b1 + c*c1*c1 + d*b1*c1 + e*a1*c1 + f*a1*b1;
			B = a*a2*a2 + b*b2*b2 + c*c2*c2 + d*b2*c2 + e*a2*c2 + f*a2*b2;
			C = a*a3*a3 + b*b3*b3 + c*c3*c3 + d*b3*c3 + e*a3*c3 + f*a3*b3;
			D = a * 2 * a2*a3 + b * 2 * b2*b3 + c * 2 * c2*c3 + d*(b2*c3 + b3*c2) + e*(a2*c3 + a3*c2) + f*(a2*b3 + a3*b2);
			E = a * 2 * a1*a3 + b * 2 * b1*b3 + c * 2 * c1*c3 + d*(b1*c3 + b3*c1) + e*(a1*c3 + a3*c1) + f*(a1*b3 + a3*b1);
			F = a * 2 * a1*a2 + b * 2 * b1*b2 + c * 2 * c1*c2 + d*(b1*c2 + b2*c1) + e*(a1*c2 + a2*c1) + f*(a1*b2 + a2*b1);

			a = A; b = B; c = C; d = D; e = E; f = F;

			//haha
			//cout << a << " " << b << " " << c << " " << d << " " << e << " " << f << endl;

			p = 0;
			w = get(b, d, c);
			k = 0;

			a1 = cos(p)*cos(k) - sin(p)*sin(w)*sin(k);
			a2 = -cos(p)*sin(k) - sin(p)*sin(w)*cos(k);
			a3 = -sin(p)*cos(w);

			b1 = cos(w)*sin(k);
			b2 = cos(w)*cos(k);
			b3 = -sin(w);

			c1 = sin(p)*cos(k) + cos(p)*sin(w)*sin(k);
			c2 = -sin(p)*sin(k) + cos(p)*sin(w)*cos(k);
			c3 = cos(p)*cos(w);

			A = a*a1*a1 + b*b1*b1 + c*c1*c1 + d*b1*c1 + e*a1*c1 + f*a1*b1;
			B = a*a2*a2 + b*b2*b2 + c*c2*c2 + d*b2*c2 + e*a2*c2 + f*a2*b2;
			C = a*a3*a3 + b*b3*b3 + c*c3*c3 + d*b3*c3 + e*a3*c3 + f*a3*b3;
			D = a * 2 * a2*a3 + b * 2 * b2*b3 + c * 2 * c2*c3 + d*(b2*c3 + b3*c2) + e*(a2*c3 + a3*c2) + f*(a2*b3 + a3*b2);
			E = a * 2 * a1*a3 + b * 2 * b1*b3 + c * 2 * c1*c3 + d*(b1*c3 + b3*c1) + e*(a1*c3 + a3*c1) + f*(a1*b3 + a3*b1);
			F = a * 2 * a1*a2 + b * 2 * b1*b2 + c * 2 * c1*c2 + d*(b1*c2 + b2*c1) + e*(a1*c2 + a2*c1) + f*(a1*b2 + a2*b1);
			a = A; b = B; c = C; d = D; e = E; f = F;

			//cout << a << " " << b << " " << c << " " << d << " " << e << " " << f << endl;

			p = get(a, e, c);
			w = 0;
			k = 0;

			a1 = cos(p)*cos(k) - sin(p)*sin(w)*sin(k);
			a2 = -cos(p)*sin(k) - sin(p)*sin(w)*cos(k);
			a3 = -sin(p)*cos(w);

			b1 = cos(w)*sin(k);
			b2 = cos(w)*cos(k);
			b3 = -sin(w);

			c1 = sin(p)*cos(k) + cos(p)*sin(w)*sin(k);
			c2 = -sin(p)*sin(k) + cos(p)*sin(w)*cos(k);
			c3 = cos(p)*cos(w);

			A = a*a1*a1 + b*b1*b1 + c*c1*c1 + d*b1*c1 + e*a1*c1 + f*a1*b1;
			B = a*a2*a2 + b*b2*b2 + c*c2*c2 + d*b2*c2 + e*a2*c2 + f*a2*b2;
			C = a*a3*a3 + b*b3*b3 + c*c3*c3 + d*b3*c3 + e*a3*c3 + f*a3*b3;
			D = a * 2 * a2*a3 + b * 2 * b2*b3 + c * 2 * c2*c3 + d*(b2*c3 + b3*c2) + e*(a2*c3 + a3*c2) + f*(a2*b3 + a3*b2);
			E = a * 2 * a1*a3 + b * 2 * b1*b3 + c * 2 * c1*c3 + d*(b1*c3 + b3*c1) + e*(a1*c3 + a3*c1) + f*(a1*b3 + a3*b1);
			F = a * 2 * a1*a2 + b * 2 * b1*b2 + c * 2 * c1*c2 + d*(b1*c2 + b2*c1) + e*(a1*c2 + a2*c1) + f*(a1*b2 + a2*b1);

			a = A; b = B; c = C; d = D; e = E; f = F;

			//cout << a << " " << b << " " << c << " " << d << " " << e << " " << f << endl;
		}


		A = sqrt(1.0 / A);
		B = sqrt(1.0 / B);
		C = sqrt(1.0 / C);

		double ans = min(min(A, B), C);
		printf("%.10lf
", ans);
	}
	return 0;
}