#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define GRIDSZ 80
#define TESTCASE 0

typedef struct {
   double x;
   double y;
} point;

typedef struct {
   int e;
   int f;
} pair;

typedef point quadrilateral[4], polygon[8];

/* Intersection of the segments P1P2 and P3P4 (not lines).
   If the segments do not intersect, the point (-1,-1) is returned */
point segint ( point P1, point P2, point P3, point P4 )
{
   double a1, b1, c1, a2, b2, c2, d;
   point Q;

   if (P1.x == P2.x) {
      a1 = 1; b1 = 0; c1 = -P1.x;
   } else {
      a1 = (P2.y - P1.y) / (P2.x - P1.x); b1 = -1; c1 = -a1 * P1.x - b1 * P1.y;
   }
   if (P3.x == P4.x) {
      a2 = 1; b2 = 0; c2 = -P3.x;
   } else {
      a2 = (P4.y - P3.y) / (P4.x - P3.x); b2 = -1; c2 = -a2 * P3.x - b2 * P3.y;
   }
   d = a1 * b2 - a2 * b1;
   if (d == 0) return (point){-1,-1};
   Q.x = (b1 * c2 - b2 * c1) / d;
   Q.y = (a2 * c1 - a1 * c2) / d;
   if (P1.x <= P2.x) {
      if ((Q.x < P1.x) || (Q.x > P2.x)) return (point){-1,-1};
   } else {
      if ((Q.x < P2.x) || (Q.x > P1.x)) return (point){-1,-1};
   }
   if (P3.x <= P4.x) {
      if ((Q.x < P3.x) || (Q.x > P4.x)) return (point){-1,-1};
   } else {
      if ((Q.x < P4.x) || (Q.x > P3.x)) return (point){-1,-1};
   }
   return Q;
}

/* On which side of the directed line segment P1P2 does P lie.
   Return values are as follows:
      0: P lies on P1P2
      1: P lies to the left of P1P2
     -1: P lies to the right of P1P2
*/
int side ( point P1, point P2, point P )
{
   double x1, y1, x2, y2, x, y, d;

   x1 = P1.x; y1 = P1.y;
   x2 = P2.x; y2 = P2.y;
   x = P.x; y = P.y;
   d = x2 * y - y2 * x - x1 * y + y1 * x + x1 * y2 - y1 * x2;
   if (d == 0) return 0;
   if (d > 0) return 1;
   return -1;
}

/* Position of point P relative to quadrilateral R.
   Return values are as follows:
      0: P is outside R
      1: P is inside R
      2: P is on an edge
      3: P is a corner
*/
int ptpos ( quadrilateral R, point P )
{
   int p0, p1, p2, p3, p;

   p0 = side(R[0],R[1],P); if (p0 < 0) return 0;
   p1 = side(R[1],R[2],P); if (p1 < 0) return 0;
   p2 = side(R[2],R[3],P); if (p2 < 0) return 0;
   p3 = side(R[3],R[0],P); if (p3 < 0) return 0;
   p = p0 + p1 + p2 + p3;
   if (p == 4) return 1;
   return p;
}

/*** GRAPHICAL METHOD ***/

double aoi1 ( quadrilateral R1, quadrilateral R2, int N )
{
   double pxldim;
   int i, j, npxl;
   point Q;

   pxldim = 1 / (double)N; npxl = 0;
   for (j=0; j<N; ++j) {
      for (i=0; i<N; ++i) {
         Q.x = ((double)j + 0.5) * pxldim;
         Q.y = ((double)i + 0.5) * pxldim;
         if ((ptpos(R1,Q) == 1) && (ptpos(R2,Q) == 1)) ++npxl;
      }
   }
   return npxl * pxldim * pxldim;
}

/*** NUMERICAL METHOD ***/

/* Intersection of the vertical segment P1P2 with the sides of R. The
   array stores the y-coordinates of the intersection points. The number
   of intersection points is returned. */
int verticalint ( quadrilateral R, point P1, point P2, double y[] )
{
   int j, k, l;
   point Q;

   l = 0;
   for (j=0; j<4; ++j) {
      k = (j == 3) ? 0 : j+1;
      Q = segint(R[j],R[k],P1,P2);
      if ((Q.x >= 0) && (Q.x <= 1) && (Q.y >= 0) && (Q.y <= 1)) y[l++] = Q.y;
   }
   return l;
}

/* This function does not handle degeneracy, so the output may be poorer
   than the actual approximation, in particular when N is small. */
double aoi2 ( quadrilateral R1, quadrilateral R2, int N )
{
   double stepsize, area, y1[4], y2[4], t;
   int i, l1, l2;
   point P1, P2;

   stepsize = 1 / (double)N; area = 0;
   P1.y = 0; P2.y = 1;
   for (i=0; i<N; ++i) {
      P1.x = P2.x = ((double)i + 0.5) * stepsize;
      l1 = verticalint(R1,P1,P2,y1);
      l2 = verticalint(R2,P1,P2,y2);
      if ((l1 == 2) && (l2 == 2)) {  /* Degenerate cases are not considered */
         if (y1[0] > y1[1]) { t = y1[0]; y1[0] = y1[1]; y1[1] = t; }
         if (y2[0] > y2[1]) { t = y2[0]; y2[0] = y2[1]; y2[1] = t; }
         if ((y2[0] <= y1[0]) && (y1[0] <= y2[1]) && (y2[1] <= y1[1]))
            area += y2[1] - y1[0];
         else if ((y1[0] <= y2[0]) && (y2[0] <= y1[1]) && (y1[1] <= y2[1]))
            area += y1[1] - y2[0];
         else if ((y2[0] <= y1[0]) && (y1[0] <= y1[1]) && (y1[1] <= y2[1]))
            area += y1[1] - y1[0];
         else if ((y1[0] <= y2[0]) && (y2[0] <= y2[1]) && (y2[1] <= y1[1]))
            area += y2[1] - y2[0];
      }
   }
   return area * stepsize;
}

/*** GEOMETRIC METHOD ***/

/* Find all edge intersections between R1 and R2. The intersection
   points are stored in I[]. The dge numbers are stored in J[]. The
   number of intersections is returned. */
int edgeint ( quadrilateral R1, quadrilateral R2, polygon I, pair J[] )
{
   int l, j1, k1, j2, k2;
   point Q;
   pair q;
   double angle[8], t;

   l = 0;
   for (j1=0; j1<4; ++j1) {
      k1 = (j1 == 3) ? 0 : j1+1;
      for (j2=0; j2<4; ++j2) {
         k2 = (j2 == 3) ? 0 : j2+1;
         Q = segint(R1[j1],R1[k1],R2[j2],R2[k2]);
         if ((Q.x >= 0) && (Q.x <= 1) && (Q.y >= 0) && (Q.y <= 1)) {
            I[l] = Q;
            J[l] = (pair){j1,j2};
            ++l;
         }
      }
   }
   if (l >= 2) {
      Q.x = Q.y = 0;
      for (j1=0; j1<l; ++j1) {
         Q.x += I[j1].x;
         Q.y += I[j1].y;
      }
      Q.x /= (double)l;
      Q.y /= (double)l;
      for (j1=0; j1<l; ++j1) angle[j1] = atan2(I[j1].y-Q.y,I[j1].x-Q.x);
      for (k1=l-2; k1>=0; --k1) {
         for (j1=0; j1<=k1; ++j1) {
            if (angle[j1] > angle[j1+1]) {
               t = angle[j1]; angle[j1] = angle[j1+1]; angle[j1+1] = t;
               Q = I[j1]; I[j1] = I[j1+1]; I[j1+1] = Q;
               q = J[j1]; J[j1] = J[j1+1]; J[j1+1] = q;
            }
         }
      }
   }
   return l;
}

/* Area of a triangle with vertices A, B, C */
double tarea ( point A, point B, point C )
{
   double area;

   area = A.x * (B.y - C.y) + B.x * (C.y - A.y) + C.x * (A.y - B.y);
   if (area < 0) area = -area;
   area /= 2;
   return area;
}

/* Area of a polygon P with n sides */
double parea ( polygon P, int n )
{
   int i;
   double area;


   printf("\n--- Number of corners on the intersection polygon = %d\n", n);
   for (i=0; i<n; ++i)
      printf("    P[%d] = (%18.16lf, %18.16lf)\n", i, P[i].x, P[i].y);
   if (n <= 2) return 0;
   area = 0;
   for (i=2; i<n; ++i) area += tarea(P[0],P[i-1],P[i]);
   return area;
}

double aoi3 ( quadrilateral R1, quadrilateral R2 )
{
   polygon I, P;
   pair J[8];
   int i, l, e, f, n, type;

   l = edgeint(R1,R2,I,J);
   printf("\n--- Number of edge intersections = %d\n", l);
   for (i=0; i<l; ++i)
      printf("    I(%d,%d) = (%18.16lf, %18.16lf)\n", J[i].e, J[i].f, I[i].x, I[i].y);

   if (l == 0) {
      if (ptpos(R2,R1[0]) > 0) /* R1 is inside R2 */
         return parea(R1,4);
      if (ptpos(R1,R2[0]) > 0) /* R2 is inside R1 */
         return parea(R2,4);
      return 0;
   }

   n = 0; type = 0; i = 0;
   while (i < l) {
      if (type == 0) {          /* Side-Side intersection */
         P[n] = I[i];
         e = J[i].e + 1; if (e == 4) e = 0;
         f = J[i].f + 1; if (f == 4) f = 0;
         if (ptpos(R2,R1[e]) > 0) type = 1;
         else if (ptpos(R1,R2[f]) > 0) type = 2;
         else ++i;
      } else if (type == 1) {   /* Corner of R1 */
         P[n] = R1[e];
         ++e; if (e == 4) e = 0;
         if (ptpos(R2,R1[e]) == 0) { ++i; type = 0; }
      } else if (type == 2) {   /* Corner of R2 */
         P[n] = R2[f];
         ++f; if (f == 4) f = 0;
         if (ptpos(R1,R2[f]) == 0) { ++i; type = 0; }
      }
      ++n;
   }

   return parea(P,n);
}

/* Generate a random convex quadrilateral. This is not part of the
   assignment. But how will the teachers generate the samples?
   This function does it. The generated quadrilateral has its four
   corners in counterclockwise order.
*/

int badangle ( point A, point B, point C )
{
   double area, a, c, theta;

   area = tarea(A,B,C);
   a = sqrt((B.x-C.x)*(B.x-C.x) + (B.y-C.y)*(B.y-C.y));
   c = sqrt((A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y));
   theta = 180 * asin(2*area/(a*c)) / M_PI;
   return (theta < 60);
}

void genquadrilateral ( quadrilateral R, int testcase )
{
   double x1, y1, x2, y2, x3, y3, x4, y4, dy;

   switch (testcase) {
      case 1:
         x1 = 0.10; x2 = 0.60; x3 = 0.50; x4 = 0.20;
         y1 = 0.80; y2 = 0.50; y3 = 0.85; y4 = 0.90;
         break;
      case 101:
         x1 = 0.30; x2 = 0.80; x3 = 0.90; x4 = 0.40;
         y1 = 0.30; y2 = 0.05; y3 = 0.40; y4 = 0.50;
         break;
      case 2:
         x1 = 0.20; x2 = 0.90; x3 = 0.50; x4 = 0.10;
         y1 = 0.10; y2 = 0.20; y3 = 0.90; y4 = 0.70;
         break;
      case 102:
         x1 = 0.25; x2 = 0.45; x3 = 0.60; x4 = 0.40;
         y1 = 0.50; y2 = 0.25; y3 = 0.60; y4 = 0.75;
         break;
      case 3:
         x1 = 0.20; x2 = 0.70; x3 = 0.30; x4 = 0.10;
         y1 = 0.30; y2 = 0.40; y3 = 0.90; y4 = 0.70;
         break;
      case 103:
         x1 = 0.15; x2 = 0.80; x3 = 0.90; x4 = 0.60;
         y1 = 0.10; y2 = 0.20; y3 = 0.60; y4 = 0.80;
         break;
      case 4:
         x1 = 0.10; x2 = 0.50; x3 = 0.40; x4 = 0.15;
         y1 = 0.20; y2 = 0.30; y3 = 0.80; y4 = 0.75;
         break;
      case 104:
         x1 = 0.30; x2 = 0.90; x3 = 0.60; x4 = 0.20;
         y1 = 0.10; y2 = 0.40; y3 = 0.70; y4 = 0.50;
         break;
      case 5:
         x1 = 0.20; x2 = 0.50; x3 = 0.80; x4 = 0.60;
         y1 = 0.40; y2 = 0.10; y3 = 0.20; y4 = 0.80;
         break;
      case 105:
         x1 = 0.40; x2 = 0.90; x3 = 0.30; x4 = 0.10;
         y1 = 0.30; y2 = 0.50; y3 = 0.90; y4 = 0.70;
         break;
      case 6:
         x1 = 0.20; x2 = 0.80; x3 = 0.90; x4 = 0.05;
         y1 = 0.20; y2 = 0.10; y3 = 0.90; y4 = 0.80;
         break;
      case 106:
         x1 = 0.50; x2 = 0.95; x3 = 0.40; x4 = 0.25;
         y1 = 0.05; y2 = 0.30; y3 = 0.70; y4 = 0.50;
         break;
      case 7:
         x1 = 0.70; x2 = 0.15; x3 = 0.45; x4 = 0.90;
         y1 = 0.90; y2 = 0.70; y3 = 0.30; y4 = 0.40;
         break;
      case 107:
         x1 = 0.80; x2 = 0.85; x3 = 0.35; x4 = 0.10;
         y1 = 0.10; y2 = 0.75; y3 = 0.95; y4 = 0.20;
         break;
      case 8:
         x1 = 0.95; x2 = 0.40; x3 = 0.05; x4 = 0.50;
         y1 = 0.60; y2 = 0.95; y3 = 0.50; y4 = 0.05;
         break;
      case 108:
         x1 = 0.10; x2 = 0.70; x3 = 0.80; x4 = 0.20;
         y1 = 0.30; y2 = 0.10; y3 = 0.90; y4 = 0.80;
         break;
      default:
         while (1) {
            x1 = ((double)rand() / (double)RAND_MAX) * 0.5;
            x2 = 0.5 + ((double)rand() / (double)RAND_MAX) * 0.5;
            x3 = 0.5 + ((double)rand() / (double)RAND_MAX) * 0.5;
            x4 = ((double)rand() / (double)RAND_MAX) * 0.5;
            y1 = ((double)rand() / (double)RAND_MAX) * 0.75;
            y2 = ((double)rand() / (double)RAND_MAX) * 0.75;

            dy = 1-y1;
            dy *= 0.5 * (1 + ((double)rand() / (double)RAND_MAX));
            y4 = y1 + dy;

            dy = 1-y2;
            dy *= 0.5 * (1 + ((double)rand() / (double)RAND_MAX));
            y3 = y2 + dy;

            if (badangle((point){x1,y1}, (point){x2,y2}, (point){x3,y3})) continue;
            if (badangle((point){x2,y2}, (point){x3,y3}, (point){x4,y4})) continue;
            if (badangle((point){x3,y3}, (point){x4,y4}, (point){x1,y1})) continue;
            if (badangle((point){x4,y4}, (point){x1,y1}, (point){x2,y2})) continue;

            break;
         }
   }

   R[0].x = x1; R[0].y = y1;
   R[1].x = x2; R[1].y = y2;
   R[2].x = x3; R[2].y = y3;
   R[3].x = x4; R[3].y = y4;

   printf("    %18.16lf %18.16lf\n", x1, y1);
   printf("    %18.16lf %18.16lf\n", x2, y2);
   printf("    %18.16lf %18.16lf\n", x3, y3);
   printf("    %18.16lf %18.16lf\n", x4, y4);
}

/*** The following functions before main() are for text-based visualization
     of the intersection region. These are not part of the assignment.
 ***/

void printgrid ( char (*G)[GRIDSZ+2] )
{
   int x, y;

   printf("\n");
   for (y=GRIDSZ+1; y>=0; --y) {
      for (x=0; x<=GRIDSZ+1; ++x) printf("%c%c", G[y][x], G[y][x]);
      printf("\n");
   }
}

void initgrid ( char (*G)[GRIDSZ+2] )
{
   int x, y;

   for (y=1; y<=GRIDSZ; ++y) for (x=1; x<=GRIDSZ; ++x) G[y][x] = ' ';
   G[0][0] = G[0][GRIDSZ+1] = G[GRIDSZ+1][0] = G[GRIDSZ+1][GRIDSZ+1] = '+';
   for (x=1; x<=GRIDSZ; ++x) G[0][x] = G[GRIDSZ+1][x] = '=';
   for (y=1; y<=GRIDSZ; ++y) G[y][0] = G[y][GRIDSZ+1] = '|';
}

void addquadrilateral ( char (*G)[GRIDSZ+2], quadrilateral R, char c )
{
   point P1, P2, P3, P4, Q;
   int i, j, k, x, y;
   double stepsize;

   for (i=0; i<4; ++i) {
      y = 1 + round(GRIDSZ * R[i].y);
      x = 1 + round(GRIDSZ * R[i].x);
      G[y][x] = 'A'+i;
   }

   stepsize = 1 / (double)GRIDSZ;

   P3.y = 0; P4.y = 1;
   for (i=0; i<GRIDSZ; ++i) {
      P3.x = P4.x = ((double)i + 0.5) * stepsize;
      for (j=0; j<4; ++j) {
         k = j + 1;
         if (k == 4) k = 0;
         P1 = R[j]; P2 = R[k];
         Q = segint(P1,P2,P3,P4);
         if ((Q.x >= 0) && (Q.x <= 1) && (Q.y >= 0) && (Q.y <= 1)) {
            y = 1 + round(GRIDSZ * Q.y);
            x = 1 + round(GRIDSZ * Q.x);
            if (G[y][x] == ' ') G[y][x] = c;
         }
      }
   }

   P3.x = 0; P4.x = 1;
   for (i=0; i<GRIDSZ; ++i) {
      P3.y = P4.y = ((double)i + 0.5) * stepsize;
      for (j=0; j<4; ++j) {
         k = j + 1;
         if (k == 4) k = 0;
         P1 = R[j]; P2 = R[k];
         Q = segint(P1,P2,P3,P4);
         if ((Q.x >= 0) && (Q.x <= 1) && (Q.y >= 0) && (Q.y <= 1)) {
            y = 1 + round(GRIDSZ * Q.y);
            x = 1 + round(GRIDSZ * Q.x);
            if (G[y][x] == ' ') G[y][x] = c;
         }
      }
   }
}

void addintersection ( char (*G)[GRIDSZ+2], quadrilateral R1, quadrilateral R2 )
{
   int y, x;
   double stepsize;
   point Q;

   stepsize = 1 / (double)GRIDSZ;
   for (x=0; x<GRIDSZ; ++x) for (y=0; y<GRIDSZ; ++y) {
      Q.x = (double)x * stepsize;
      Q.y = (double)y * stepsize;
      if ((ptpos(R1,Q) == 1) && (ptpos(R2,Q) == 1) && (G[y+1][x+1] == ' '))
         G[y+1][x+1] = '.';
   }
}

int main ()
{
   quadrilateral R1, R2;
   char G[GRIDSZ+2][GRIDSZ+2];

   srand((unsigned int)time(NULL));

   printf("\n*** Case (%c)\n", 'a'+TESTCASE-1);
   printf("\n+++ Rectangle 1\n"); genquadrilateral(R1,TESTCASE);
   printf("\n+++ Rectangle 2\n"); genquadrilateral(R2,TESTCASE+100);

   initgrid(G);
   addquadrilateral(G,R1,'1');
   addquadrilateral(G,R2,'2');
   addintersection(G,R1,R2);
   printgrid(G);

   printf("\n+++ Graphical method\n");
   printf("    N = %7d. Area of intersection = %18.16lf\n", 10, aoi1(R1,R2,10));
   printf("    N = %7d. Area of intersection = %18.16lf\n", 100, aoi1(R1,R2,100));
   printf("    N = %7d. Area of intersection = %18.16lf\n", 1000, aoi1(R1,R2,1000));

   printf("\n+++ Numerical method\n");
   printf("    N = %7d. Area of intersection = %18.16lf\n", 10, aoi2(R1,R2,10));
   printf("    N = %7d. Area of intersection = %18.16lf\n", 100, aoi2(R1,R2,100));
   printf("    N = %7d. Area of intersection = %18.16lf\n", 1000, aoi2(R1,R2,1000));
   printf("    N = %7d. Area of intersection = %18.16lf\n", 10000, aoi2(R1,R2,10000));
   printf("    N = %7d. Area of intersection = %18.16lf\n", 100000, aoi2(R1,R2,100000));
   printf("    N = %7d. Area of intersection = %18.16lf\n", 1000000, aoi2(R1,R2,1000000));

   printf("\n+++ Geometric method\n");
   printf("\n    Area of intersection = %18.16lf\n", aoi3(R1,R2));

   exit(0);
}