Note: This article is available in Chinese only. 本文暂无英文版本。
View original
题意：给出椭球方程的 6 个参数 a,b,c,d,e,f
问椭球上的点到原点 (0,0,0) 的最小距离是多少。
思路：感觉难点在于，如何保证搜到的点一直在椭球上。
一开始我考虑用椭球的参数方程，然后发现不记得是什么了 2333。
然后看了题解，发现比较巧妙的做法是，只搜索 x,y，然后从椭球方程中解出 z。
x,y 确定以后，椭球方程就变成了一个关于 z 的一元二次方程，可解。
由于是要求距离原点的最小距离，而现在可能得到的两个解是关于 xoy 平面对称的，只有 z 坐标不同，因此我们取距离原点近的那个 z。
以及，感觉在平面上搜 4 个方向就好，没必要 8 个方向？
WA 到死是因为计算距离时忘记开根号了。呵呵呵呵，我是傻逼。
代码实现
1/* ***********************************************
2Author :111qqz
3Created Time :2016年08月30日 星期二 20时46分05秒
4File Name :code/poj/5017.cpp
5 ************************************************ */
6#include <cstdio>
7#include <cstring>
8#include <iostream>
9#include <algorithm>
10#include <vector>
11#include <queue>
12#include <set>
13#include <map>
14#include <string>
15#include <cmath>
16#include <cstdlib>
17#include <ctime>
18#define fst first
19#define sec second
20#define lson l,m,rt<<1
21#define rson m+1,r,rt<<1|1
22#define ms(a,x) memset(a,x,sizeof(a))
23typedef long long LL;
24#define pi pair < int ,int >
25#define MP make_pair
26#define INF 1E20
27#define MAX 1
28using namespace std;
29const double eps = 1E-8;
30const int dx4[4]={1,0,0,-1};
31const int dy4[4]={0,-1,1,0};
32const int inf = 0x3f3f3f3f;
33double a,b,c,d,e,f;
34int dblcmp(double d)
35{
36 return d<-eps?-1:d>eps;
37}
38struct point
39{
40 double x,y,z;
41 point(){}
42 point (double _x,double _y,double _z):x(_x),y(_y),z(_z){}
43 double dis(point b){
44 double res = 0 ;
45 res = (x-b.x)*(x-b.x)+(y-b.y)*(y-b.y)+(z-b.z)*(z-b.z);
46 return sqrt(res);
47 }
48 void look()
49 {
50 printf("x:%.7f y:%.7f z:%.7f ",x,y,z);
51 }
52};
53double getZ(double x,double y)
54{
55 double A,B,C,delta;
56 A = c;
57 B = d*y+e*x;
58 C = a*x*x+b*y*y+f*x*y-1;
59 delta = B*B-4*A*C;
60 // cout<<"A:"<<A<<" B:"<<B<<" C:"<<C<<" delta:"<<delta<<endl;
61 if (dblcmp(delta)<0) return 1E31; //无解
62 delta = sqrt(delta);
63 double z1 = (-B+delta)/2.0/A;
64 double z2 = (-B-delta)/2.0/A;
65 if (dblcmp(z1*z1-z2*z2)<0) return z1; //x,y的坐标已经确定，我们选距离(0,0,0）近的解。
66 else return z2;
67}
68double SA(point &ret)
69{
70 double step,cur,ans=INF;
71 const double delta=0.99;
72 bool tag;
73 point nxt;
74 point target = point(0,0,0);
75 for ( ret = point(0,0,0),step = 1 ; step>eps ; step*=delta){
76 for (tag = true; tag ;){
77 tag = false;
78 for ( int i = 0 ; i < 4 ; i++){
79 nxt.x = ret.x + dx4[i]*step; //只搜x,y，通过解方程得到z啊巧妙地保证了点一直在椭球面上。
80 nxt.y = ret.y + dy4[i]*step;
81 nxt.z = getZ(nxt.x,nxt.y);
82 if (nxt.z>=1E30) continue; //无解。
83 cur = nxt.dis(target);
84 if (dblcmp(cur-ans)<0)
85 {
86 ans = cur;
87 ret = nxt;
88 tag = true;
89 break;
90 }
91 }
92 }
93 }
94 return ans;
95}
96int main()
97{
98#ifndef ONLINE_JUDGE
99 freopen("code/in.txt","r",stdin);
100#endif
101 while (~scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f))
102 {
103 double ans;
104 point ret;
105 ans = SA(ret);
106 printf("%.7f\n",ans);
107 }
108#ifndef ONLINE_JUDGE
109 fclose(stdin);
110#endif
111 return 0;
112}