POJ 1408スラグ計算ジオメトリ

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;

struct Point2D{ double x,y; }mesh[33][33],p[4][33];
struct Vector2D{ double x,y; };
Vector2D operator-(Point2D a,Point2D b)
{
    Vector2D v;
    v.x=b.x-a.x;v.y=b.y-a.y;
    return v;
}

double operator^(Vector2D a,Vector2D b){
    return fabs(a.x*b.y-b.x*a.y);
}

Point2D getPoint( Point2D a,Point2D b,Point2D c,Point2D d )
{
    Point2D temp;
    temp.x=((a.x*((c-b)^(d-b)))+(b.x*((c-a)^(d-a))))/(((c-b)^(d-b))+((c-a)^(d-a)));
    temp.y=((a.y*((c-b)^(d-b)))+(b.y*((c-a)^(d-a))))/(((c-b)^(d-b))+((c-a)^(d-a)));
    return temp;
}

double Area( Point2D a,Point2D b,Point2D c,Point2D d ){
    return (fabs((c-a)^(d-a))+fabs((c-b)^(d-b))+fabs((a-c)^(b-c))+fabs((a-d)^(b-d)))/4;
}

int main()
{
    //freopen( "input.txt","r",stdin );
    //freopen( "output.txt","w",stdout );
    int N;
    while( scanf("%d",&N)!=EOF )
    {
        if( N==0 )
            break;
        memset(p,0,sizeof(p));
        p[0][0].x=0;     p[0][0].y=0;     p[0][N+1].x=1;   p[0][N+1].y=0;
        p[1][0].x=0;     p[1][0].y=1;     p[1][N+1].x=1;   p[1][N+1].y=1;
        for( int i=0;i<2;i++ ) for( int j=1;j<=N;j++ )
            scanf( "%lf",&p[i][j].x ),p[i][j].y=i%2;
        for( int i=2;i<4;i++ ) for( int j=1;j<=N;j++ )
            scanf( "%lf",&p[i][j].y ),p[i][j].x=i%2;

        for( int i=0;i<=N+1;i++ )
        {
            for( int j=0;j<=N+1;j++ )
            {
            //printf( "i=%d j=%d ",i,j );
            if( i==0 )   { mesh[i][j]=p[0][j];continue; }
            if( i==N+1 ) { mesh[i][j]=p[1][j];continue; }
            if( j==0 )   { mesh[i][j]=p[2][i];continue; }
            if( j==N+1 ) { mesh[i][j]=p[3][i];continue; }
            //mesh[i][j].x=0.1;mesh[i][j].y=0.1;printf( "(%lf,%lf) ",mesh[i][j].x,mesh[i][j].y );
            mesh[i][j]=getPoint( p[0][j],p[1][j],p[2][i],p[3][i] );
            //printf( "(%lf,%lf) ",mesh[i][j].x,mesh[i][j].y );
            }
            //printf( "
" );
        }
/*
        for( int i=0;i<=N+1;i++ ){
        for( int j=0;j<=N+1;j++ )
            printf( "(%lf,%lf) ",mesh[i][j].x,mesh[i][j].y );
        printf( "
" );
        }
        */
        double ans=0;
        for( int i=0;i<=N;i++ )
        for( int j=0;j<=N;j++ )
                ans=max( ans,Area(mesh[i][j],mesh[i+1][j+1],mesh[i][j+1],mesh[i+1][j]) );
        printf( "%.6f
",ans );
    }
    return 0;
}