bzoj1336/1337 && hdu3932 最小圆覆盖 计算几何：随机增量法&模拟退火

题目链接

bzoj1336
bzoj1337
hdu3932

知识点：随机增量法

增量法。
如果当前点不在当前的圆上，重新构造的最小圆上必定有此点。
然后枚举另外两个点即可。
枚举方法就是从1~i-1枚举j，再从1~j-1枚举k。
期望复杂度是O(n)在随机环境下是o的k的。

代码：随机增量法

bzoj1336和1337代码类似，这里只贴一个。

/**************************************************************
    Problem: 1336
    User: 2014lvzelong
    Language: C++
    Result: Accepted
    Time:452 ms
    Memory:3012 kb
****************************************************************/

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<cmath>
#include<set>
using namespace std;
const int N = 110000;
const double eps = 1e-10;
int read() {
    char ch = getchar(); int x = 0, f = 1;
    while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0'; ch = getchar();}
    return x * f;
}
int dblcmp(double x) {return x > eps ? 1 : (x < -eps ? -1 : 0);}
int n;
double ans, r;
struct point {
    double x, y;
    point(double a = 0, double b = 0) : x(a), y(b) {}
    double operator * (point a) const {return x * a.y - y * a.x;}
    double operator ^ (point a) const {return x * a.x + y * a.y;}
    point operator - (point a) const {return point(x - a.x, y - a.y);}
    point operator + (point a) const {return point(x + a.x, y + a.y);}
    point operator / (double a) const {return point(x / a, y / a);}
}o, p[N];
double sqr(double x) {return x * x;}
double sqr(point a) {return a ^ a;}
double len(point a) {return sqrt(sqr(a));}

void geto(point a, point b, point c) {
    double a1 = 2 * (b.x - a.x), b1 = 2 * (b.y - a.y), c1 = sqr(b.x) + sqr(b.y) - sqr(a.x) - sqr(a.y);
    double a2 = 2 * (c.x - b.x), b2 = 2 * (c.y - b.y), c2 = sqr(c.x) + sqr(c.y) - sqr(b.x) - sqr(b.y); //ax+by=c
    o.x = (c1 * b2 - c2 * b1) / (b2 * a1 - b1 * a2);
    o.y = (c1 * a2 - c2 * a1) / (b1 * a2 - a1 * b2);
}

int main() {
    srand(20010927);
    scanf("%d", &n);
    for(int i = 1;i <= n; ++i) scanf("%lf%lf", &p[i].x, &p[i].y);
    for(int i = 1;i <= n; ++i) swap(p[rand() % n + 1], p[rand() % n + 1]);
    o = p[1];
    for(int i = 2;i <= n; ++i) {
        if(dblcmp(len(o - p[i]) - r) <= 0) continue;
        o = (p[i] + p[1]) / 2; r = len(p[i] - o);
        for(int j = 2;j < i; ++j) {
            if(dblcmp(len(o - p[j]) - r) <= 0) continue;
            o = (p[i] + p[j]) / 2; r = len(p[j] - o);
            for(int k = 1;k < j; ++k) {
                if(dblcmp(len(o - p[k]) - r) <= 0) continue;
                geto(p[i], p[j], p[k]); r = len(p[i] - o);
            }
        }
    }
    printf("%.10lf\n%.10lf %.10lf\n", r, o.x, o.y);
    return 0;
}

知识点2：模拟退火

一个有技巧的退火。
首先扔十个随机的点。
然后每次枚举多个方向跑一定步长就可以了。
注意到最小圆覆盖只有一个极小值点，所以不需要接受答案更差的点。

代码2：模拟退火

bzoj1336小数点后10位什么的模拟退火是不存在的。
hdu3932才一位而且有限定范围嘻嘻，所以当然是挑软柿子捏咯。
不过仍然嫌弃跑得慢，考试当然是选择随机增量咯，不过思路是可以借鉴的。

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<ctime>
#include<cmath>
using namespace std;
const double eps = 1e-8, pi = acos(-1);
const int N = 1010, L = 40;
int n, x, y;
struct point {
    double x, y;
    point(double a = 0, double b = 0) : x(a), y(b) {}
    double operator ^ (point a) {return x * a.x + y * a.y;}
    point operator - (point a) {return point(x - a.x, y - a.y);}
    point Move(double step, double ang) {return point(x + step * cos(ang), y + step * sin(ang));}
}p[N], rp[N];
double len[N];
double Rand() {return (rand() % 10000 + 1) / 10000.0;}
double sqr(point p) {return p ^ p;}
double vlen(point p) {return sqrt(sqr(p));}
double Dis(point o) {
    double ret = 0;
    for(int i = 1; i <= n; ++i)
        ret = max(ret, vlen(p[i] - o));
    return ret;
}
bool check(point p) {return p.x > eps && p.y > eps && p.x + eps < x && p.y + eps < y;}
int main() {
    srand(20010927);
    while(~scanf("%d%d%d", &x, &y, &n)) {
        for(int i = 1;i <= n; ++i) scanf("%lf%lf", &p[i].x, &p[i].y);
        for(int i = 1;i <= 10; ++i) {
            rp[i] = point(Rand() * x, Rand() * y);
            len[i] = Dis(rp[i]);
        }
        for(double step = max(x, y) / sqrt(n); step > 0.001; step *= 0.8) {
            for(int i = 1;i <= 10; ++i) {
                point cur = rp[i];
                for(int j = 0;j < L; ++j) {
                    point nxt = cur.Move(step, Rand() * 2 * pi);
                    if(!check(nxt)) continue;
                    double dis = Dis(nxt);
                    if(dis + eps < len[i]) {
                        rp[i] = nxt;
                        len[i] = dis;
                    }
                }
            }
        }
        int k = 1;
        for(int i = 2;i <= 10; ++i) 
        if(len[i] + eps < len[k])
            k = i;
        printf("(%.1lf,%.1lf).\n", rp[k].x, rp[k].y);
        printf("%.1lf\n", len[k]);
    }
    return 0;
}