Ellipsoid 退火模拟

Problem Description
Given a 3-dimension ellipsoid(椭球面)

your task is to find the minimal distance between the original point (0,0,0) and points on the ellipsoid. The distance between two points (x1,y1,z1) and (x2,y2,z2) is defined as

Input
There are multiple test cases. Please process till EOF.

For each testcase, one line contains 6 real number a,b,c(0 < a,b,c,< 1),d,e,f(0 ≤ d,e,f < 1), as described above. It is guaranteed that the input data forms a ellipsoid. All numbers are fit in double.

Output
For each test contains one line. Describes the minimal distance. Answer will be considered as correct if their absolute error is less than 10-5.

Sample Input

1 0.04 0.01 0 0 0

Sample Output

1.0000000

Source
2014 ACM/ICPC Asia Regional Xi’an Online

---

直接求解思路很简单：
拉格朗日乘数法即可
不过这个计算量极大；
因为是一般式，而且四个未知数；

---

所以不妨不断逼近最优解；

----> 退火模拟；

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<string>
#include<cmath>
#include<map>
#include<set>
#include<vector>
#include<queue>
#include<bitset>
#include<ctime>
#include<deque>
#include<stack>
#include<functional>
#include<sstream>
#include<cctype>
//#pragma GCC optimize("O3")
using namespace std;
#define maxn 200005
#define inf 0x3f3f3f3f
#define INF 0x7fffffff
#define rdint(x) scanf("%d",&x)
#define rdllt(x) scanf("%lld",&x)
typedef long long  ll;
typedef unsigned long long ull;
typedef unsigned int U;
#define ms(x) memset((x),0,sizeof(x))
const int mod = 1e9 + 7;
#define Mod 20100403
#define sq(x) (x)*(x)
#define eps 1e-10
const int N = 1505;

inline int rd() {
	int x = 0;
	char c = getchar();
	bool f = false;
	while (!isdigit(c)) {
		if (c == '-') f = true;
		c = getchar();
	}
	while (isdigit(c)) {
		x = (x << 1) + (x << 3) + (c ^ 48);
		c = getchar();
	}
	return f ? -x : x;
}

ll gcd(ll a, ll b) {
	return b == 0 ? a : gcd(b, a%b);
}
ll sqr(ll x) { return x * x; }

double a, b, c, d, e, f;
double Mindist;
double dx[] = { 0,1,0,-1,1,1,-1,-1 };
double dy[] = { 1,0,-1,0,1,-1,1,-1 };

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

double calZ(double x, double y) {
	double B = 1.0*d*y + 1.0*e * x;
	double A = c;
	double C = a * x*x + b * y*y + f * x*y - 1.0;
	double delt = B * B - 4.0 * A*C;
	if (delt < 0.0)return INF;
	double tmp1 = (-B + sqrt(delt)) / A / 2.0;
	double tmp2 = (-B - sqrt(delt)) / A / 2.0;
//	if (dist(x, y, tmp1) < dist(x, y, tmp2))return tmp1;
	if (tmp1 - tmp2 < 0.0)return tmp1;
	else return tmp2;
}

void Test() {
	double x = 0.0, y = 0.0, z = 1.0*sqrt(1.0 / c);
	double tmp = 100.0;
	while (tmp > eps) {
		double nx = x + (rand() * 2 - RAND_MAX)*tmp;
		double ny = y + (rand() * 2 - RAND_MAX)*tmp;
		double nz = calZ(nx, ny);
		double dt = dist(nx, ny, nz) - dist(x, y, z);
		if (dt < 0.0)x = nx, y = ny, z = nz;
		else if (exp(-dt / tmp)*RAND_MAX > rand())x = nx, y = ny, z = nz;
		tmp *= 0.998;
	}
	Mindist = dist(x, y, z);
}

int main()
{
	//ios::sync_with_stdio(false);
	while (cin >> a >> b >> c >> d >> e >> f) {
	    srand((int)time(NULL));
		Test();
		printf("%.7lf\n", 1.0*Mindist);
	//	cout << 1.0*sqrt(1.0 / c) << endl;
	}
	return 0;
}

或者：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<string>
#include<cmath>
#include<map>
#include<set>
#include<vector>
#include<queue>
#include<bitset>
#include<ctime>
#include<deque>
#include<stack>
#include<functional>
#include<sstream>
#include<cctype>
//#pragma GCC optimize("O3")
using namespace std;
#define maxn 200005
#define inf 0x3f3f3f3f
#define INF 0x7fffffff
#define rdint(x) scanf("%d",&x)
#define rdllt(x) scanf("%lld",&x)
typedef long long  ll;
typedef unsigned long long ull;
typedef unsigned int U;
#define ms(x) memset((x),0,sizeof(x))
const int mod = 1e9 + 7;
#define Mod 20100403
#define sq(x) (x)*(x)
#define eps 1e-9
const int N = 1505;

inline int rd() {
    int x = 0;
    char c = getchar();
    bool f = false;
    while (!isdigit(c)) {
        if (c == '-') f = true;
        c = getchar();
    }
    while (isdigit(c)) {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    return f ? -x : x;
}

ll gcd(ll a, ll b) {
    return b == 0 ? a : gcd(b, a%b);
}
ll sqr(ll x) { return x * x; }

double a, b, c, d, e, f;
double Mindist;
double dx[] = { 0,1,0,-1,1,1,-1,-1 };
double dy[] = { 1,0,-1,0,1,-1,1,-1 };

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

double calZ(double x, double y) {
    double B = 1.0*d*y + e * x;
    double A = c;
    double C = a * x*x + b * y*y + f * x*y - 1;
    double delt = B * B - 4 * A*C;
    if (delt < 0)return INF;
    double tmp1 = (-B + sqrt(delt)) / A / 2.0;
    double tmp2 = (-B - sqrt(delt)) / A / 2.0;
    if (dist(x, y, tmp1) < dist(x, y, tmp2))return tmp1;
    else return tmp2;
}

void Test() {
    double x = 0, y = 0, z = sqrt(1.0 / c);
    double sp = 100.0;
    while (sp > eps) {
        for (int i = 0; i < 8; i++) {
            double nx = x + sp * dx[i], ny = y + sp * dy[i];
            double nz = calZ(nx, ny);
            if (dist(nx, ny, nz) < dist(x, y, z))x = nx, y = ny, z = nz;
            //sp *= 0.998;
        }
        sp *= 0.998;
    }
    Mindist = dist(x, y, z);
}

int main()
{
    //ios::sync_with_stdio(false);
    while (cin >> a >> b >> c >> d >> e >> f) {
    //    srand((int)time(NULL));
        Test();
        printf("%.7lf\n", 1.0*Mindist);
    }
    return 0;
}

---

之前一道省选题考过退火模拟；
[JSOI2004]平衡点 / 吊打XXX
当然也可以按照上述方法处理；

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<string>
#include<cmath>
#include<map>
#include<set>
#include<vector>
#include<queue>
#include<bitset>
#include<ctime>
#include<deque>
#include<stack>
#include<functional>
#include<sstream>
#include<cctype>
//#pragma GCC optimize("O3")
using namespace std;
#define maxn 200005
#define inf 0x3f3f3f3f
#define INF 0x7fffffff
#define rdint(x) scanf("%d",&x)
#define rdllt(x) scanf("%lld",&x)
typedef long long  ll;
typedef unsigned long long ull;
typedef unsigned int U;
#define ms(x) memset((x),0,sizeof(x))
const int mod = 1e9 + 7;
#define Mod 20100403
#define sq(x) (x)*(x)
#define eps 1e-7
const int N = 1505;

inline int rd() {
	int x = 0;
	char c = getchar();
	bool f = false;
	while (!isdigit(c)) {
		if (c == '-') f = true;
		c = getchar();
	}
	while (isdigit(c)) {
		x = (x << 1) + (x << 3) + (c ^ 48);
		c = getchar();
	}
	return f ? -x : x;
}

ll gcd(ll a, ll b) {
	return b == 0 ? a : gcd(b, a%b);
}
ll sqr(ll x) { return x * x; }

struct node {
	double x, y, w;
}tg[maxn];
double ansx, ansy;
int n;
double dx[] = { 0,1,0,-1,1,1,-1,-1 };
double dy[] = { 1,0,-1,0,-1,1,1,-1 };



double cal(double x, double y) {
	double ans = 0;
	for (int i = 1; i <= n; i++) {
		double dx = x - tg[i].x;
		double dy = y - tg[i].y;
		ans += sqrt(dx*dx + dy * dy)*tg[i].w;
	}
	return ans;
}

void Test() {
	double tp = 180.0;
	while (tp > eps) {
		for (int i = 1; i < 8; i++) {
			double nx = ansx + tp * dx[i];
			double ny = ansy + tp * dy[i];
			double dt = cal(nx, ny) - cal(ansx, ansy);
			if (dt < 0)ansx = nx, ansy = ny;
		}
		tp *= 0.998;
	}
}

int main()
{
	//ios::sync_with_stdio(false);
	rdint(n);
	double sumx = 0, sumy = 0;
	for (int i = 1; i <= n; i++) {
		scanf("%lf%lf%lf", &tg[i].x, &tg[i].y, &tg[i].w);
		sumx += tg[i].x; sumy += tg[i].y;
	}
	ansx = 1.0* sumx / n; ansy = 1.0*sumy / n;
	Test();
	printf("%.3lf %.3lf\n", ansx, ansy);
	return 0;
}