P4518 [JSOI2018]绝地反击

模拟退火+网络流瞎跑 $\texttt{40pts}$

求退火优化（$or$ 调参 $or$ 卡常技巧）

---

主要思路：

退火圆上在 $x$ 轴上方第一个点。（退火角度）

然后二分求最小值，将问题转化成判定性问题。

然后求圆外到圆上的点的最大距离可以跑二分图匹配。

复杂度大概是一个 $O(Time \times \log{n} \times m \times \sqrt{n})$ 的样子。 $Times$ 指模拟退火的次数。

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <queue>

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 310, M = N * N * 2;
const double eps = 1e-8, INF = 2e9;;
const double pi = acos(-1);
double delta;
PDD p[N], c[N];
int h[N], e[M], ne[M], idx = 1;
int d[N], S, T;
int cur[N], n, r;
double ans = INF;
int f[M];

double dist(PDD a, PDD b) {
	return (double)sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

void add(int a, int b, int c) {
	e[ ++ idx] = b, ne[idx] = h[a], f[idx] = c, h[a] = idx;
}

namespace Max_Flow {
	bool bfs() {
	    queue<int> q;
	    memset(d, -1, sizeof d);
	    q.push(S), d[S] = 0, cur[S] = h[S];
	    
	    while (q.size()) {
	        int t = q.front();
	        q.pop();
	        
	        for (int i = h[t]; i; i = ne[i]) {
	            int ver = e[i];
	            if ((d[ver] == -1) && f[i]) {
	                d[ver] = d[t] + 1;
	                cur[ver] = h[ver];
	                if (ver == T) return true;
	                q.push(ver);
	            }
	        }
	    }
	    
	    return false;
	}
	
	int find(int u, int limit) {
	    if (u == T) return limit;
	    
	    int flow = 0;
	    for (int i = cur[u]; i && flow < limit; i = ne[i]) {
	        int ver = e[i];
	        cur[u] = i;
	        if ((d[ver] == d[u] + 1) && f[i]) {
	            int t = find(ver, min(limit - flow, f[i]));
	            if (!t) d[ver] = -1;
	            f[i] -= t, f[i ^ 1] += t, flow += t;
	        }
	    }
	    
	    return flow;
	}
	
	int dinic() {
	    int res = 0, flow;
	    while (bfs()) while (flow = find(S, n + 10)) res += flow;
		return res; 
	}
}
using namespace Max_Flow;

double rand(double l, double r) {
	return (double)rand() / RAND_MAX * (r - l) + l;
}

namespace Simulate_Anneal {
	double calc(double x) { // 初始角为x度时的值 
		for (int i = 1; i <= n; i ++ ) { // 处理每一个圆上点的值 
			double angle = (double)(i - 1) * delta + x;
			c[i + n] = {(double)cos(angle) * r, (double)sin(angle) * r};
		}
		
		double l = 0, r = 100.00;
		while (r - l > eps) {
			double mid = (l + r) / 2.0;
			memset(h, 0, sizeof h);
			idx = 1;
			
			for (int i = 1; i <= n; i ++ ) {
				for (int j = n + 1; j <= n + n; j ++ ) {
					if (dist(p[i], c[j]) < mid) 
						add(i, j, 1),
						add(j, i, 0);
				}
			}
			
			for (int i = 1; i <= n; i ++ )
				add(S, i, 1), add(i, S, 0);
			for (int i = n + 1; i <= n + n; i ++ )
				add(i, T, 1), add(T, i, 0);
			
			if (dinic() == n) r = mid;
			else l = mid;
		}
		ans = min(ans, l);
		return l;
	}
	
	void simulate_anneal() {
		double x = rand(0, delta);
		for (double t = 2.0 * pi; t >= 1e-7; t *= 0.99) {
			double y = rand(0, delta);
			double delta = calc(y) - calc(x);
			if (exp( - delta / t) < rand() / RAND_MAX)
				x = y;
			if ((double)clock() / CLOCKS_PER_SEC > 0.98) break;
		}
	}
}
using namespace Simulate_Anneal;

int main() {
//	freopen("anneal.in", "r", stdin);
//	freopen("anneal.out", "w", stdout);
	
	srand(1000000009);
	scanf("%d%d", &n, &r);
	S = 0, T = n + n + 1;
	delta = 2 * pi / n;
	for (int i = 1; i <= n; i ++ )
		scanf("%lf%lf", &p[i].x, &p[i].y);
	while ((double)clock() / CLOCKS_PER_SEC <= 0.98)
		simulate_anneal(); 
	
	printf("%.8lf\n", ans);
	
	return 0;
}