/* Problem 4--The Fence Builder
   This was not only the hardest program in this contest (not a single
   team in the entire region got it) but also the hardest problem I have
   ever seen in a regional competition.  It is based on the geometrical
   fact that the maximum area occurs when all the vertices lie on the
   perimeter of the same circle.  I don't know what solution they had in
   mind, but I set it up as a series of n trig equations in n unknowns
   and solved the problem with the multivariable Newton's method.  This
   sometimes gives the "wrong" root, so I have to try several different
   initial estimates to get the correct one.  The matrix algorithms here
   are very crude; the running time could be greatly improved, but it
   runs fast enough, and I'm just glad I finally finished the thing.  To
   be quite honest, this problem was why I delayed so long in putting up
   solutions to this contest! */

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

#define PI 3.1415926535897932384626433832795

int main (int argc, char **argv);
double FindArea (double *l, int sz);
void ge (double **M, int rows, int cols);
int inv (double **M, int sz);
double ** Subtract (double **A, double **B, int rows, int cols);
double ** Multiply (double **A, double **B, int r, int rc, int c);
double ** FindAngles (double *l, int sz, int random);
double ** MakeMat (int r, int c);
void FreeMat (double **M, int r);
double L2 (double **A, double ** B, int rows, int cols);

int main (int argc, char **argv) {

  FILE *in, *out;
  int r, sz, cs, i;
  double *l;

 srand (time (NULL)); /* I use a randomized Newton's method. */
 in = fopen ("prob4.in","r");
 out = fopen ("prob4.out","w");
 cs = 0;
 for (;;) {
  fscanf (in,"%d",&sz); /* get number of fences */
 if (sz==0) break;
  l = malloc (sz*sizeof (double));
  i = 0;
  for (;;) {
  if (i==sz) break;  /* read in array of fence lengths */
   fscanf (in,"%lf",&l[i]);
   i++;
  }
  fprintf (out,"Case %d: maximum area = %.3f\n",++cs,FindArea(l,sz));
  free (l);
 }
 fclose (in);
 fclose (out);
 r = EXIT_SUCCESS;
 return r;
}

/* FindArea takes l, an array of size sz of fence lengths and returns
   the maximum area that can be enclosed by these fence lengths. */
double FindArea (double *l, int sz) {

  double **M, area, hrsq;
  int i, nct;

 M = FindAngles (l,sz,0); /* Use a straightforward starting estimate */
 for (;;) { /* to find a possible solution.  If no solution is found */
  nct = 0;  /* with that estimate or if the solution is out of bounds */
  i = 0;    /* for our problem by containing a angle outside the */
  for (;;) {/* range of -pi to pi, or by containing more than one */
  if (i==sz || nct > 1) break; /* negative angle, start using random */
   nct += M[i][0] < 0; /* starting estimates until you get the answer.*/
   if (fabs (M[i][0]) > PI) {
    nct = 2;
   } else {
   }
   i++;
  }
 if (nct <= 1) break;
  FreeMat (M,sz);
  M = FindAngles (l,sz,1);
 } /* we now have the angles from the center of the circle to the */
 hrsq = pow(l[0],2)/(4*(1+cos(M[0][0]))); /* vertices of the fence */
 i = 0; /* we compute the area of the polygon by summing the areas */
 area = 0; /* of the component triangles */
 for (;;) {
 if (i==sz) break;
  area += sin (M[i][0]);
  i++;
 }
 area *= hrsq;
 FreeMat (M,sz);
 return area;
}

/* ge takes a Matrix M with rows rows and cols columns and puts it into
   row-reduced echelon form using Gaussian Elimination with Partial
   Pivoting.  The new matrix overwrites M. */
void ge (double **M, int rows, int cols) {

  int r, c, tr, tc, maxr;
  double max, td;

 r = c = 0; /* We are computing column c, from row r on down. */
 for (;;) {
 if (r==rows || c==cols) break;
  max = fabs (M[r][c]);
  maxr = r;
  tr = r+1;
  for (;;) { /* Find the element in the column with largest absolute */
  if (tr==rows) break; /* value.  This will be our pivot. */
   if (fabs (M[tr][c]) > max) {
    max = fabs (M[tr][c]);
    maxr = tr;
   } else {
   }
   tr++;
  }
  if (max > 1e-8) {
   tc = c;
   for (;;) {
   if (tc==cols) break;
    td = M[maxr][tc];
    M[maxr][tc] = M[r][tc];
    M[r][tc] = td;
    tc++;
   }
   tc = c+1;
   for (;;) { /* divide this row by the leading element */
   if (tc==cols) break;
    M[r][tc] /= M[r][c];
    tc++;
   }
   M[r][c] = 1;
   tr = 0;
   for (;;) { /* subtract this row from the other rows to zero out */
   if (tr==rows) break; /* column */
    if (tr != r) {
     tc = c+1;
     for (;;) {
     if (tc==cols) break;
      M[tr][tc] -= M[tr][c] * M[r][tc];
      tc++;
     }
     M[tr][c] = 0;
    } else {
    }
    tr++;
   }
   r++;
  } else {
  }
  c++;
 }
}

/* inv computes the inverse of sz x sz matrix M, by pasting an identity
   matrix to the right of it, row reducing the larger matrix, and
   pulling the inverse off the right.  The inverse overwrites the
   original matrix.  inv returns 1 if there was an inverse, 0 if M was
   not invertible. */
int inv (double **M, int sz) {

  double **TM;
  int tr, tc, valid;

 TM = MakeMat (sz,2*sz);
 tr = 0;
 for (;;) {
 if (tr==sz) break;
  tc = 0;
  for (;;) {
  if (tc==sz) break;
   TM[tr][tc] = M[tr][tc];
   tc++;
  }
  tr++;
 }
 tr = 0;
 for (;;) { /* Tack identity matrix to the right */
 if (tr==sz) break;
  tc = sz;
  for (;;) {
  if (tc==2*sz) break;
   if (tr == tc-sz) {
    TM[tr][tc] = 1;
   } else {
    TM[tr][tc] = 0;
   }
   tc++;
  }
  tr++;
 }
 ge (TM,sz,2*sz); /* row reduce it */
 tr = 0;
 for (;;) { /* put inverse back into M */
 if (tr==sz) break;
  tc = 0;
  for (;;) {
  if (tc==sz) break;
   M[tr][tc] = TM[tr][tc+sz];
   tc++;
  }
  tr++;
 }
 tr = 0; /* Make sure there really is an identity matrix on the left */
 valid = 1;
 for (;;) {
 if (tr==sz || !valid) break;
  valid = TM[tr][tr]==1;
  tr++;
 }
 FreeMat (TM,sz);
 return valid;
}

/* Subtract computes and returns A-B, where A and B each have rows rows
   and cols columns. */
double ** Subtract (double **A, double **B, int rows, int cols) {

  int tr, tc;
  double **diff;

 diff = MakeMat (rows,cols);
 tr = 0;
 for (;;) {
 if (tr==rows) break;
  tc = 0;
  for (;;) {
  if (tc==cols) break;
   diff[tr][tc] = A[tr][tc] - B[tr][tc];
   tc++;
  }
  tr++;
 }
 return diff;
}
/* Multiply computes and returns A*B where A has r rows and rc columns
   and B has rc rows and c columns. */
double ** Multiply (double **A, double **B, int r, int rc, int c) {

  double **prod;
  int tr, tc, td;

 prod = MakeMat (r,c);
 tr = 0;
 for (;;) {
 if (tr==r) break;
  tc = 0;
  for (;;) {
  if (tc==c) break;
   prod[tr][tc] = 0;
   td = 0;
   for (;;) {
   if (td==rc) break;
    prod[tr][tc] += A[tr][td] * B[td][tc];
    td++;
   }
   tc++;
  }
  tr++;
 }
 return prod;
}

/* FindAngles takes an array of fence sizes of size sz and computes and
   returns the angles the center of the circle makes with the vertices
   of the polygon, assuming all vertices lie on the perimeter of the
   circle.  It does this by setting up n trig equations in n unknowns,
   and solving them with multivariable Newton's method.  There is no
   guarantee that a solution will be found or that the solution will
   be the one we are looking for.  FindArea calls this method repeatedly
   until a good solution is found.  If random is true, a random initial
   estimate is made, otherwise the initial estimate is that all angles
   are equal.  The solution is returned in the form of a column
   vector. */
double ** FindAngles (double *l, int sz, int random) {

  double **x0, **x1, **M, **f, **p;
  int tr, tc, maxr, ict;

 x0 = MakeMat (sz,1);
 if (!random) {
  tr = 0;
  for (;;) { /* All angles are equal. */
  if (tr==sz) break;
   x0[tr][0] = PI*(sz-2)/sz;
   tr++;
  }
 } else {
  tr = 0;
  maxr = 0;
  for (;;) { /* Angles are randomly generated. */
  if (tr==sz) break;
   x0[tr][0] = PI*rand()/(1+(double)RAND_MAX);
   tr++;
  }
 }
 ict = 0;
 for (;;) {
  M = MakeMat (sz,sz);
  tr = 0;
  for (;;) {
  if (tr==sz) break;
   tc = 0;
   for (;;) {
   if (tc==sz) break;
    M[tr][tc] = 0;
    tc++;
   }
   tr++;
  }
  tr = 0;
  for (;;) { /* Program trig equations into matrix */
  if (tr==sz-1) break;
   M[tr][0] = pow(l[0],2)*sin(x0[0][0])/pow(1+cos(x0[0][0]),2);
   M[tr][tr+1] =
      -pow(l[tr+1],2)*sin(x0[tr+1][0])/pow(1+cos(x0[tr+1][0]),2);
   tr++;
  }
  tc=0;
  for (;;) {
  if (tc==sz) break;
   M[sz-1][tc] = 1;
   tc++;
  }
  f = MakeMat (sz,1);
  tr = 0;
  for (;;) {
  if (tr==sz-1) break;
   f[tr][0] = pow(l[0],2)/(1+cos(x0[0][0])) -
              pow(l[tr+1],2)/(1+cos(x0[tr+1][0]));
   tr++;
  }
  f[sz-1][0] = -PI*(sz-2);
  tr = 0;
  for (;;) {
  if (tr==sz) break;
   f[sz-1][0] += x0[tr][0];
   tr++;
  } /* Plug into Newton's Method */
  if (!inv (M,sz)) { /* If not invertible, bad estimate */
   ict = 100;
  } else {
  }
  p = Multiply (M,f,sz,sz,1);
  x1 = Subtract (x0,p,sz,1);
  FreeMat (M,sz);
  FreeMat (f,sz);
  FreeMat (p,sz); /* If solution is close to previous solution or I've*/
 if (L2 (x0,x1,sz,1) < 1e-5 || ict == 100) break;/*iterated too long*/
  FreeMat (x0,sz);
  x0 = x1;
  ict++;
 }
 if (ict==100) { /* If I've iterated too long make this unacceptable */
  x1[0][0] = 2*PI;
 } else {
 }
 FreeMat (x0,sz);
 return x1;
}

/* MakeMat allocates space and returns an r x c matrix. */
double ** MakeMat (int r, int c) {

  double **M;
  int tr;

 M = malloc (r * sizeof (double *));
 tr = 0;
 for (;;) {
 if (tr==r) break;
  M[tr] = malloc (c * sizeof (double));
  tr++;
 }
 return M;
}

/* FreeMat frees up the space used by an r-row matrix. */
void FreeMat (double **M, int r) {

  int tr;

 tr = 0;
 for (;;) {
 if (tr==r) break;
  free (M[tr]);
  tr++;
 }
 free (M);
}

/* L2 computes the L2-norm of the difference between two matrices A and
   B with rows rows and cols columns. */
double L2 (double **A, double **B, int rows, int cols) {

  int tr, tc;
  double norm;

 norm = 0;
 tr = 0;
 for (;;) {
 if (tr==rows) break;
  tc = 0;
  for (;;) {
  if (tc==cols) break;
   norm += pow (A[tr][tc]-B[tr][tc],2);
   tc++;
  }
  tr++;
 }
 norm = sqrt (norm);
 return norm;
}