poj1408 计算几何

题意很简单，算法也很简单，就是枚举所有的点然后求得最大的四边形的面积，但是我的数组一开始开小了，结果WA了好多次，而且没有找到错误，无语了都！

#include<iostream>
#include<cmath>
#include<cstdio>
using namespace std;
double Max;
#define exp 0.00000000001
int n;
typedef struct point
{
	double x,y;
}rr;
point a[4][40];	 
point s[5];
void fun(point p1,point p2,point p3,point p4,int j)
{
	double k,b,x,y;
	k=(p2.y-p1.y)/(p2.x-p1.x);
	b=p2.y-p2.x*k;
	if(fabs(p3.x-p4.x)<exp)
	{
         x=p3.x;
		 y=k*x+b;
	}
	else
	{
		double k1,b1;
		k1=(p3.y-p4.y)/(p3.x-p4.x);//k于k1肯定是不相等的
		b1=p3.y-k1*p3.x;
        x=(b-b1)/(k1-k);
		y=k*x+b;
	}
	s[j].x=x;
	s[j].y=y;
}
void area(point p1,point p2,point p4,point p3)
{
	 double len=sqrt(pow(p3.y-p1.y,2)+pow(p3.x-p1.x,2));
	 double a;
	 if(fabs(p3.x-p1.x)>exp)
	 {
       double k=(p3.y-p1.y)/(p3.x-p1.x);
	   double b=p3.y-p3.x*k;
	   double h1=fabs(k*p2.x-p2.y+b)/(sqrt(k*k+1));
	   double h2=fabs(k*p4.x-p4.y+b)/(sqrt(k*k+1));
	   a=(h1*len+h2*len)/2;
	 }
	 else
	 {
		 double h1=fabs(p2.x-p1.x);
		 double h2=fabs(p4.x-p1.x);
         a=(h1*len+h2*len)/2;
	 }
	 if(Max<a)
		 Max=a; 
}
int main()
{
     int i,j; 
	 while(scanf("%d",&n))
	 {
		 if(n==0)
			 break;
		  for(i=1; i<=n; i++)
		  {
			  scanf("%lf",&a[0][i].x);
			  a[0][i].y=0;
		  }
	      a[0][0].x=0;a[0][0].y=0;
		  a[0][n+1].x=1;
		  a[0][n+1].y=0;
		  Max=0;
		  for(i=1; i<=n; i++)
		  {
			  scanf("%lf",&a[1][i].x);
			  a[1][i].y=1;
		  }
		  a[1][0].x=0;a[1][0].y=1;
		  a[1][n+1].x=1;a[1][n+1].y=1;
		  for(i=1; i<=n; i++)
		  {
			  scanf("%lf",&a[2][i].y);
			  a[2][i].x=0;
		  }
		  a[2][0].x=0;a[2][0].y=0;
		  a[2][n+1].x=0;a[2][n+1].y=1;
		  for(i=1; i<=n; i++)
		  {
			  scanf("%lf",&a[3][i].y);
			  a[3][i].x=1;
		  }
	//	  n=30;
		  a[3][0].x=1;a[3][0].y=0;
		  a[3][n+1].x=1;a[3][n+1].y=1;
		  for(i=0; i<=n; i++)
		  {
			  for(j=0; j<=n; j++)//求交点的
			  {
				  //分别求得是四个交点的
                   fun(a[2][i],a[3][i],a[0][j],a[1][j],0);
				   fun(a[2][i],a[3][i],a[0][j+1],a[1][j+1],1);
				   fun(a[2][i+1],a[3][i+1],a[0][j],a[1][j],2);
				   fun(a[2][i+1],a[3][i+1],a[0][j+1],a[1][j+1],3);
				   area(s[0],s[1],s[2],s[3]);
			  }
		  }
		  printf("%0.6lf\n",Max);
	 }
	 return 0;
}