[THOJ 1589] 椭球面 三分套三分

题意

　　现在给出一个椭球面:　$ax ^ 2 + by ^ 2 + cz ^ 2 + dyz + exz + fxy = 1$ .

　　求椭球面到 $(0, 0, 0)$ 的距离.

　　$T \le 200, 0 < a, b, c < 1, 0 \le d, e, f < 1$ .

　　假装 $(0, 0, 0)$ 在椭球内部.

分析

　　二次的式子通常都是单峰的.

　　猜测椭球到 $(0, 0, 0)$ 的距离也是单峰的.

　　三分 x , 三分 y , 利用一元二次方程解出 z 并计算距离.

实现

　　实现小结:

　　　　1.　对于三分, 模拟退火之类的, 最好记录访问到的最大的答案, 而不是最后再求一次.

　　　　2.　浮点数上, 如果运算出错, 会返回 nan , 判断一个数是不是 nan , 就判断 !(x < 0) && !(x >= 0) .

　　最终实现:

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <cstdlib>
 4 #include <cctype>
 5 #include <cmath>
 6 #include <algorithm>
 7 using namespace std;
 8 #define F(i, a, b) for (register int i = (a); i <= (b); i++)
 9 #define db double
10
11 const db INF = 1e30;
12
13 db a, b, c, d, e, f, ans;
14
15 inline bool nan(db x) { return !(x < 0) && !(x >= 0); }
16 inline db Dist(db x, db y, db z) { return sqrt(x * x + y * y + z * z); }
17 inline db Calc(db x, db y) {
18     db A = c;
19     db B = d * y + e * x;
20     db C = a * x * x + b * y * y + f * x * y - 1;
21     db D = B * B - 4 * A * C;
22     db z = (-B - sqrt(D)) / (2 * A);
23     db tmp = Dist(x, y, z);
24     return ans = min(ans, tmp), tmp;
25 }
26 inline db Y(db x) {
27     db l = 0, r = 1000 / sqrt(b);
28     while (r-l > 1e-10) {
29         db _l = (l+l+r)/3, ansL = Calc(x, _l);
30         db _r = (l+r+r)/3, ansR = Calc(x, _r);
31         nan(ansL) || nan(ansR) || ansL < ansR ? r = _r : l = _l;
32     }
33     db ans = Calc(x, l);
34
35     l = -1000 / sqrt(b), r = 0;
36     while (r-l > 1e-10) {
37         db _l = (l+l+r)/3, ansL = Calc(x, _l);
38         db _r = (l+r+r)/3, ansR = Calc(x, _r);
39         nan(ansL) || nan(ansR) || ansR < ansL ? l = _l : r = _r;
40     }
41     db ans2 = Calc(x, l);
42
43     return nan(ans) ? ans2 : nan(ans2) ? ans : min(ans, ans2);
44 }
45 inline db X(void) {
46     db l = 0, r = 1000 / sqrt(a);
47     while (r-l > 1e-10) {
48         db _l = (l+l+r)/3, ansL = Y(_l);
49         db _r = (l+r+r)/3, ansR = Y(_r);
50         nan(ansL) || nan(ansR) || ansL < ansR ? r = _r : l = _l;
51     }
52     db ans = Y(l);
53
54     l = -1000 / sqrt(a), r = 0;
55     while (r-l > 1e-10) {
56         db _l = (l+l+r)/3, ansL = Y(_l);
57         db _r = (l+r+r)/3, ansR = Y(_r);
58         nan(ansL) || nan(ansR) || ansR < ansL ? l = _l : r = _r;
59     }
60     db ans2 = Y(l);
61
62     return nan(ans) ? ans2 : nan(ans2) ? ans : min(ans, ans2);
63 }
64
65 int main(void) {
66     #ifndef ONLINE_JUDGE
67         freopen("e.in", "r", stdin);
68     #endif
69
70     while (~scanf("%lf %lf %lf %lf %lf %lf", &a, &b, &c, &d, &e, &f))
71         ans = INF, X(), printf("%0.8lf\n", ans);
72
73     return 0;
74 }