[hdu3644 A Chocolate Manufacturer's Problem]模拟退火，简单多边形内最大圆

题意：判断简单多边形内是否可以放一个半径为R的圆

思路：如果这个多边形是正多边形，令r(x,y)为圆心在(x,y)处多边形内最大圆的半径，不难发现，f(x,y)越靠近正多边形的中心，r越大，所以可以利用模拟退火法来逼近最优点。对于一般的多边形，由于可能存在多个这样的"局部最优点"，所以可以选不同的点作为起点进行多若干次模拟退火即可。

模拟退火的过程：每次由原状态S生成一个新状态T，如果T比S优，那么接受这一次转移，否则以一定概率P接受这次转移，因为这样可能会跳过局部最优解而得到全局最优解。

PS:步长每次改变的系数一般设为0.8~0.9，eps不能设太高。

#pragma comment(linker, "/STACK:10240000")
#include <bits/stdc++.h>
using namespace std;

#define X                   first
#define Y                   second
#define pb                  push_back
#define mp                  make_pair
#define all(a)              (a).begin(), (a).end()
#define fillchar(a, x)      memset(a, x, sizeof(a))

typedef long long ll;
typedef pair<int, int> pii;

namespace Debug {
void print(){cout<<endl;}template<typename T>
void print(const T t){cout<<t<<endl;}template<typename F,typename...R>
void print(const F f,const R...r){cout<<f<<" ";print(r...);}template<typename T>
void print(T*p, T*q){int d=p<q?1:-1;while(p!=q){cout<<*p<<", ";p+=d;}cout<<endl;}
}
template<typename T>bool umax(T&a, const T&b){return b<=a?false:(a=b,true);}
template<typename T>bool umin(T&a, const T&b){return b>=a?false:(a=b,true);}
/* -------------------------------------------------------------------------------- */

const double eps = 1e-4;/** 设置比较精度 **/
struct Real {
    double x;
    double get() { return x; }
    int read() { return scanf("%lf", &x); }
    Real(const double &x) { this->x = x; }
    Real() {}
    Real abs() { return x > 0? x : -x; }

    Real operator + (const Real &that) const { return Real(x + that.x);}
    Real operator - (const Real &that) const { return Real(x - that.x);}
    Real operator * (const Real &that) const { return Real(x * that.x);}
    Real operator / (const Real &that) const { return Real(x / that.x);}
    Real operator - () const { return Real(-x); }

    Real operator += (const Real &that) { return Real(x += that.x); }
    Real operator -= (const Real &that) { return Real(x -= that.x); }
    Real operator *= (const Real &that) { return Real(x *= that.x); }
    Real operator /= (const Real &that) { return Real(x /= that.x); }

    bool operator < (const Real &that) const { return x - that.x <= -eps; }
    bool operator > (const Real &that) const { return x - that.x >= eps; }
    bool operator == (const Real &that) const { return x - that.x > -eps && x - that.x < eps; }
    bool operator <= (const Real &that) const { return x - that.x < eps; }
    bool operator >= (const Real &that) const { return x - that.x > -eps; }

    friend ostream& operator << (ostream &out, const Real &val) {
        out << val.x;
        return out;
    }
    friend istream& operator >> (istream &in, Real &val) {
        in >> val.x;
        return in;
    }
};

struct Point {
    Real x, y;
    int read() { return scanf("%lf%lf", &x.x, &y.x); }
    Point(const Real &x, const Real &y) { this->x = x; this->y = y; }
    Point() {}
    Point operator + (const Point &that) const { return Point(this->x + that.x, this->y + that.y); }
    Point operator - (const Point &that) const { return Point(this->x - that.x, this->y - that.y); }
    Real operator * (const Point &that) const { return x * that.x + y * that.y; }
    Point operator * (const Real &that) const { return Point(x * that, y * that); }
    Point operator += (const Point &that)  { return Point(this->x += that.x, this->y += that.y); }
    Point operator -= (const Point &that)  { return Point(this->x -= that.x, this->y -= that.y); }
    Point operator *= (const Real &that)  { return Point(x *= that, y *= that); }

    bool operator == (const Point &that) const { return x == that.x && y == that.y; }

    Real cross(const Point &that) const { return x * that.y - y * that.x; }
    Real dist() { return sqrt((x * x + y * y).get()); }
};
typedef Point Vector;

struct Segment {
    Point a, b;
    Segment(const Point &a, const Point &b) { this->a = a; this->b = b; }
    Segment() {}
    bool intersect(const Segment &that) const {
        Point c = that.a, d = that.b;
        Vector ab = b - a, cd = d - c, ac = c - a, ad = d - a, ca = a - c, cb = b - c;
        return ab.cross(ac) * ab.cross(ad) < 0 && cd.cross(ca) * cd.cross(cb) < 0;
    }
    Point getLineIntersection(const Segment &that) const {
        Vector u = a - that.a, v = b - a, w = that.b - that.a;
        Real t = w.cross(u) / v.cross(w);
        return a + v * t;
    }
    Real Distance(Point P) {
        Point A = a, B = b;
        if (A == B) return (P - A).dist();
        Vector v1 = B - A, v2 = P - A, v3 = P - B;
        if (v1 * v2 < 0) return v2.dist();
        if (v1 * v3 > 0) return v3.dist();
        return v1.cross(v2).abs() / v1.dist();
    }
};

const int maxn = 55;
double PI = acos(-1.0);

Point p[maxn];
int n;

Real getAngel(Point o, Point a, Point b) {
    a -= o;
    b -= o;
    Real ans = acos((a * b / a.dist() / b.dist()).get());
    return a.cross(b) <= 0? ans : -ans;
}

bool inPolygon(Point o) {
    Real total = 0;
    for (int i = 0; i < n; i ++) {
        total += getAngel(o, p[i], p[(i + 1) % n]);
    }
    return total.abs() > PI;
}

Real getR(Point o) {
    Real ans = 1e9;
    for (int i = 0; i < n; i ++) {
        Segment seg(p[i], p[(i + 1) % n]);
        umin(ans, seg.Distance(o));
    }
    return ans;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
#endif // ONLINE_JUDGE
    while (cin >> n, n) {
        p[0].read();
        Real maxx = p[0].x, minx = p[0].x, maxy = p[0].y, miny = p[0].y;
        for (int i = 1; i < n; i ++) {
            p[i].read();
            umax(maxx, p[i].x);
            umin(minx, p[i].x);
            umax(maxy, p[i].y);
            umin(miny, p[i].y);
        }
        Real R;
        R.read();
        Point a(minx, miny), b(maxx, maxy);
        bool ok = false;
        for (int i = 0; !ok && i < n; i ++) {
            Real deta = (b - a).dist() / 4;
            Point O = (p[i] + p[(i + 1) % n]) * 0.5;
            int cnt = 0;
            while (!ok && deta > 0 && cnt < 100) {
                for (int j = 0; ; j ++) {
                    double randnum = rand();
                    Point newp(O.x + deta * sin(randnum), O.y + deta * cos(randnum));
                    if (!inPolygon(newp)) continue;
                    Real buf = getR(newp);
                    if (buf > getR(O) || j > 4) { /** 这里考虑了概率因素 **/
                        if (buf >= R) ok = true;
                        O = newp;
                        break;
                    }
                }
                deta *= 0.8;
                cnt ++;
            }
        }
        puts(ok? "Yes" : "No");
    }

}

[hdu3644 A Chocolate Manufacturer's Problem]模拟退火，简单多边形内最大圆