模擬退火+圓與多邊形的面積交

#include"bits/stdc++.h"
using namespace std;
//圓在多邊形內，或者多邊形在圓內都可以算出
const double eps = 1e-15;
const double PI = acos( -1.0 ) ;
inline double sqr( double x ){ return x * x ; }
inline int dcmp( double x ){
    if ( fabs(x) < eps ) return 0 ;
    return x > 0? 1 : -1 ;
}

struct Point{
    double x , y ;
    Point(){}
    Point( double _x , double _y ): x(_x) , y(_y) {}
    void input() { scanf( "%lf%lf" ,&x ,&y ); }
    void print(){
        printf("%.6f %.6f\n",x,y);
    }
//    double norm() { return sqrt( sqr(x) + sqr(y) ); }

    friend Point operator + ( const Point &a , const Point &b ) { return Point( a.x + b.x , a.y + b.y ) ; }
    friend Point operator - ( const Point &a , const Point &b ) { return Point( a.x - b.x , a.y - b.y ) ; }
    friend Point operator * ( const Point &a , const double &b ) { return Point( a.x * b , a.y * b ) ; }
    friend Point operator * ( const double &a , const Point &b ) { return Point( b.x * a , b.y * a ) ; }
    friend Point operator / ( const Point &a , const double &b ) { return Point( a.x / b , a.y / b ) ; }
    friend bool operator == ( const Point &a , const Point &b ) { return dcmp( a.x - b.x ) == 0 && dcmp( a.y - b.y ) == 0 ; }

    bool operator < ( const Point &a )const{
        return ( dcmp( x - a.x ) < 0 ) || ( dcmp( x - a.x ) == 0 && dcmp( y - a.y ) < 0 ) ;
    }
};
typedef Point Vector;
double Dot( Point a , Point b ) { return a.x * b.x + a.y * b.y ; }
double Cross( Point a , Point b ) { return a.x * b.y - a.y * b.x ; }
double Length(Vector A) { return sqrt(Dot(A,A)) ; }

int n,m;
Point p[205];

int CircleInterLine( Point a, Point b, Point o, double r, Point *p )
{
    Point p1 = a - o ;
    Point d = b - a ;
    double A = Dot( d, d ) ;
    double B = 2 * Dot( d, p1 ) ;
    double C = Dot( p1, p1 ) - sqr(r) ;

    double delta = sqr(B) - 4*A*C ;
    if ( dcmp(delta) < 0 ) return 0 ;//相離
    if ( dcmp(delta) == 0 ) { //相切
        double t = -B / (2*A) ; // 0 <= t <= 1說明交點線上段上
        if ( dcmp( t - 1 ) <= 0 && dcmp( t ) >= 0 ) {
            p[0] = a + t * d ;
            return 1 ;
        }
    }
    if ( dcmp(delta) > 0 ) { //相交
        double t1 = ( -B - sqrt(delta) ) / (2*A) ;
        double t2 = ( -B + sqrt(delta) ) / (2*A) ; //0 <= t1, t2 <= 1說明交點線上段上
        int k = 0 ;
        if ( dcmp( t1 - 1 ) <= 0 && dcmp( t1 ) >= 0 )
            p[k++] = a + t1 * d ;
        if ( dcmp( t2 - 1 ) <= 0 && dcmp( t2 ) >= 0 )
            p[k++] = a + t2 * d ;
        return k ;
    }
    return 0;
}
double Triangle_area( Point a, Point b )
{
    return fabs( Cross( a , b ) ) / 2.0  ;
}
double Sector_area( Point a, Point b, double r )
{
    double ang = atan2( a.y , a.x ) - atan2( b.y, b.x  ) ;
    while ( ang <= 0 ) ang += 2.0 * PI ;
    while ( ang > 2.0 * PI ) ang -= 2.0 * PI ;
    ang = min( ang, 2.0*PI - ang ) ;
    return sqr(r) * ang/2.0 ;
}
double calc( Point a , Point b , double r )
{
    Point pi[2] ;
    if ( dcmp( Length(a) - r ) < 0 ) {
        if ( dcmp( Length(b) - r ) < 0 ) {
            return Triangle_area( a, b ) ;
        }
        else {
            CircleInterLine( a, b, Point(0,0), r, pi) ;
            return Sector_area( b, pi[0], r ) + Triangle_area( a, pi[0] ) ;
        }
    }
    else {
        int cnt = CircleInterLine( a, b, Point(0,0), r, pi ) ;
        if ( dcmp( Length(b) - r ) < 0 ) {
            return Sector_area( a, pi[0],r ) + Triangle_area( b, pi[0] ) ;
        }
        else {
            if ( cnt == 2 )
                return Sector_area( a, pi[0],r ) + Sector_area( b, pi[1], r ) + Triangle_area( pi[0], pi[1]) ;
            else
                return Sector_area( a, b, r ) ;
        }
    }
}
double area_CircleAndPolygon( Point *p , Point o , double r )
{
    double res = 0 ;
    p[n] = p[0] ;
    for ( int i = 0 ; i < n ; i++ ) {
        int tmp = dcmp( Cross( p[i] - o , p[i+1] - o ) ) ;
        if ( tmp )
            res += tmp * calc( p[i] - o , p[i+1] - o , r ) ;
    }
    return fabs( res ) ;
}

double PolygonArea(Point *p,int n) {
    double area=0;
    for(int i=1; i < n-1; i++)
        area+=Cross(p[i]-p[0],p[i+1]-p[0]);
    return fabs(area/2);
}

double get_max(Point *p, Point o)
{
    double ret = 0;
    for(int i = 0; i < n; i++)
    ret = max(ret,Length(o-p[i]));
    return ret;
}
void solve(int n, double r)
{
    double delta = 0.99;
    double ret = 0;
    for(int i = 0; i < n; i++){
        for(int k = 0; k < 20; k++){
            double step = 100;
            Point o = Point(p[i].x,p[i].y);
            double ans = area_CircleAndPolygon(p,o,r);
            while(step > eps){
                for(int j = 0; j < 10; j++){
                    double rad = double(rand()%360+1)*PI/180;
                    double nx = o.x + step*cos(rad);
                    double ny = o.y + step*sin(rad);
                    if(nx < 0 && ny < 0 || (nx > 100 && ny > 100)) continue;
                    double area = area_CircleAndPolygon(p,Point(nx,ny),r);
                    if(area > ans) {
                        ans = area;
                        o = Point(nx,ny);
                    }
                    step *= delta;
                }
            }
            ret = max(ret,ans);
        }
    }
    printf("%7f\n",ret);
}
int main()
{
#ifdef LOCAL
    freopen("in.txt","r",stdin);
#endif // LOCAL
    srand(time(0));
    double r;
    scanf("%d",&n);
    scanf("%lf",&r);
    for(int i = 0; i < n; i++){
        p[i].input();
    }
    solve(n,r);
}