/* BHS is a program that combines ham sandwich partitioning with grid
   coarsening to partition a graph quickly and efficiently.  It is
   invoked as follows:
   bhs <# of points> <# of pieces> <input filename> {<coarsening
   factor, def: 650>}
   bhs with no arguments prompts the user for the arguments.
   The input file consists of a series of lines each with four
   floating-point numbers: the x- and y-coordinates of the point and the
   two weights of that point.
   The output is sent to a file ham.out.dualcr.<# of pieces> .  Each
   line consists of a single integer between (inclusive) 0 and <# of
   pieces>-1 indicating the piece in which the corresponding input point
   should be placed.

   The program is written in ANSI C with NO non-ANSI feature of any kind
   and NO C++ feature of any kind.  It should compile and run on system
   supporting ANSI C. */

/**********************************************************************/

#include <stdio.h>                   /* Standard include files */
#include <stdlib.h>
#include <math.h>
#include <string.h>

/**********************************************************************/

#define TOLERANCE 1e-10  /* FP comparision tolerance */
#define SHRINKSIZE 25    /* How small does a set have to be to do
                            brute-force division */
#define RANDOMGUESS 100  /* How many random points are chosen during the
                            selection phase */
#define MAX(a,b) ((a) > (b) ? (a) : (b))  /* Standard macros */
#define MIN(a,b) ((a) < (b) ? (a) : (b))

/**********************************************************************/

typedef struct lntp {double wt,m,b; int index;} lntp; /* weighted line*/
typedef lntp *lna;     /* runtime array of weighted lines */
typedef struct lnst {int ct; double wt; lna ln;} lnst;
              /* weighted line structure */
typedef struct pttp *pta;  /* linked list of biweighted points */
typedef struct pttp {double wt1, wt2; pta n; double x,y;} pttp;
          /* biweighted points */
typedef struct ptst {int ct; pta pt;} ptst; /* biweighted point struct*/
typedef struct wr {double wt, val; int index;} wr;
    /* weight type for median */
typedef wr *dwarr;  /* array of weight type */
typedef struct inftype {int inf; double val;} inftype;
    /* type fp arithmetic with infinities */
typedef inftype *infarr; /* array of infinities */
typedef int (*cfunc) (const void*,const void*); /*general cmp function*/
typedef double (*wfunc) (const void*); /* general weight function */
typedef struct trapcoord {inftype x; double y, dwt; int index;}
               trapcoord;   /* Coordinates of a trap corner*/
typedef trapcoord trap[4]; /* the four corners of trapezoid */
typedef struct llt *lltptr;  /*general linked list */
typedef struct llt {lltptr n;} llt; /* general linked list node */
typedef struct bound {double minx, miny, maxx, maxy;} bound;
                  /* boundary region */
typedef struct ptlin *ptptr; /* linked list of point arrays */
typedef struct ptlin {pta pt; ptptr n;} ptlin; /* individual point arr*/
typedef void *(*nfunc)(const void *); /* general next function */
typedef void (*snfunc)(void *,void *);/* general setnext function */

/**********************************************************************/

int main (int argc, char **argv);
static void rham (ptst *Q, int p);
static void rhamo (ptptr partition, int size1, int p, ptptr pt,
                   int *Pct, ptst *P, ptst *Q, int *PP, lltptr *L,
                   int *Lsize, lltptr bigll);
static pta sb (pta *P1pt, int *P1ct, int *P2ct, ptst *P, ptst *Q,
               int *PP, lltptr *L, int *Lsize, lltptr bigll);
static pta sbfrac (pta *P1pt, int *P1ct, int *P2ct, double frac,
                   ptst *P, ptst *Q, int *PP, lltptr *L, int *Lsize,
                   lltptr bigll);
static int ham (pta P, int ct, double *m, double *b, pta bigp);
static int hamfrac (pta *P, int ct, double *m, double *b, double frac,
                    int *above, int *dln, pta bigp);
static int tcmp (double a, double b);
static void VertCheck (lnst *L, pta P, int first,
                       double *ma, double *mb, double frac, pta bigp);
static pttp IntersPt (lnst *L1, lnst *L2, double frac);
static int IsEqual (double a, double b);
static int median (void *A, int size, int n, double k, cfunc cmp,
                   wfunc wt, int *d, double *dwt);
static int mcmp (const void *x, const void *y);
static double Lwt (const void *x);
static int Shrink (lnst *L1, lnst *L2, double *k1, double *k2,
                   inftype *l, inftype *r, pttp *pt);
static pttp FindPt (lnst *L1, lnst *L2, double k1, double k2,
                    inftype *l, inftype *r);
static int partition (void *A, int size, int f, int l, double *prevwt,
                      cfunc cmp, wfunc wt);
static void BuildArray (lntp *ln1, lntp *ln2, infarr V, int *vct,
                        inftype *l, inftype *r);
static int SearchPt (lnst *L1, lnst *L2, infarr V, int *vct, double k1,
                     double k2, pttp *pt, inftype *l, inftype *r,
                     int *f, int *e);
static void GetCorner (lnst *L1, inftype *i, double k1, trapcoord *t1,
                       trapcoord *t2);
static int DCheck (lnst *L, trap T);
static double IntersTrap (lntp *ln, trap TV, double *bottom);
static int infcmp (inftype *a, inftype *b);
static void sort (infarr A, int *n);
static void computeintersections (lnst *L, inftype *x, dwarr A);
static int Acmp (const void *x, const void *y);
static double Awt (const void *x);
static double Sint (lntp *ln, trapcoord t1, trapcoord t2);
static double wtbelow (lntp *ln, trapcoord t);
static int infcmpc (inftype *a, inftype *b, int exact);
static int infexcmp (const void *a, const void *b);
static void linsplit (pta *A, double frac, pta *l, pta *r, int *wrap);
static int dcmp (const void *a, const void *b);
static void totalwt (pta A, double *twt1, double *twt2);
static int wtOK (double iwt, double owt, double twt, double frac);
static void Advance (double *iwt1, double *iwt2, double *owt1,
                     double *owt2, pta l, pta *r, pta A, int *wrap);
static void *emalloc (size_t size);
static void *ecalloc (size_t size);
static ptst *coarsen (pta pt, int ct, bound *grid, lltptr **L,
                      int *Lsize, lltptr *bigll, int **PP);
static void msort (void *A,cfunc cmp,nfunc NEXT,snfunc SETNEXT);
static void *msortsplit (void *A,nfunc NEXT,snfunc SETNEXT);
static void *msortmerge (void *A,void *B,cfunc cmp,nfunc NEXT,
                         snfunc SETNEXT);
static void *ptanext (const void *p);
static void ptasetnext (void *p, void *n);
static double distpointline (double x, double y, int vert, double m,
                             double b);
static int llcmp (const void *a, const void *b);
static void *llnext (const void *p);
static void llsetnext (void *p, void *n);
static void sreadln (char *line, int max);
static int nreadln (void);
static int nextchar(void);

/**********************************************************************/

static int CSIZE;            /* Coarsening size */
static lltptr gbll;          /* global linked list */
static ptst *gP;             /* global point list */

/**********************************************************************/

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

  int i,SIZE=0,PARTS=0;
  ptst P;
  FILE *rfile;
  char fn[256], ifn[256];

 if (argc != 1 && argc != 4 && argc != 5) {
  fprintf (stderr,"ERROR:  Improper invocation\n"); return EXIT_FAILURE;
 }
 CSIZE = 650;  /* Default coarsening size */
 switch (argc) {              /* Input verification */
  case 1:
   printf ("Enter size:  "); SIZE = nreadln ();
   printf ("Enter parts:  "); PARTS = nreadln ();
   printf ("Enter filename:  ");sreadln (ifn,sizeof ifn);
   printf ("Enter coarsening size:  "); CSIZE = nreadln ();break;
  case 5: CSIZE = atoi (argv[4]);
  case 4:
   SIZE = atoi (argv[1]);
   PARTS = atoi (argv[2]);
   strncpy (ifn,argv[3],255); ifn[255]=0; break;
 }
 if (NULL==(rfile = fopen (ifn,"r"))) {
  printf ("Input file cannot be opened.\n"); exit (EXIT_FAILURE);}
 P.pt = emalloc (SIZE*sizeof(pttp));
 for (i=0;i<SIZE;i++) {     /* reading information from the pointfile */
  fscanf (rfile,"%lf %lf %lf %lf\n",&(P.pt[i].x),&(P.pt[i].y),
         &(P.pt[i].wt1),&(P.pt[i].wt2));
  P.pt[i].n = &P.pt[i+1];
 }
 P.pt[SIZE-1].n= NULL;
 fclose (rfile);
 P.ct = SIZE;
 rham (&P,PARTS);   /* time the program run */
 sprintf (fn,"ham.out.dualcr.%d",PARTS);
 if (NULL==(rfile = fopen (fn,"w"))) {
  fprintf (stderr," ERROR:  Output file cannot be opened.\n");
  exit (EXIT_FAILURE);
 }
 for (i=0;i<P.ct;i++) fprintf (rfile,"%d\n",(int)P.pt[i].n);
 fclose (rfile);                 /* write to the output file */
 free (P.pt);
 return EXIT_SUCCESS;
}

/**********************************************************************/

/* rham divides the point set Q into p nearly equal pieces. */

static void rham (ptst *Q, int p) {

  int i, Lsize, *PP, pct;
  pta j;
  ptst *P;
  lltptr *L, Pctr, bigll;
  bound grid;
  ptptr partition;

 P = coarsen (Q->pt,Q->ct,&grid,&L,&Lsize,&bigll,&PP); /* shrink the */
 partition = ecalloc (p*sizeof (ptlin));  /* number of points */
 pct = 1;
 partition[0].pt = P->pt;
 rhamo (partition,P->ct,p,&partition[0],&pct,P,Q,PP,L,&Lsize,bigll);
 i = 0;
 for (i=0;i<p;i++)  /* Give the points their correct partition number */
  for (j=partition[i].pt;j!=NULL;j=j->n)
   for (Pctr = L[PP[j-P->pt]];Pctr!=NULL;Pctr=Pctr->n)
    Q->pt[Pctr-bigll].n=(pta)i;
 free (PP);
 free (bigll);
 free (L);
 free (P->pt);
 free (P);
 free (partition);
}

/**********************************************************************/

/* rhamo divides a point set in half and recursively divides these
   sets in half until the proper division is obtained.
   partition will contain the array of a linked list of points
   corresponding to the partition.  size1 is the number of points to
   divide.  p is the number of pieces in the partition, pt is the
   list of points to divide.  Pct is the index of partition to put the
   split points.  P is the list of coarsened points; Q is the original
   list of points.
   PP is the index map to expand coarsened points.
   L is the list of expandable points.  Lsize is the size of PP.
   bigll is a contiguous list of linked list nodes corresponding to
   points.*/

static void rhamo (ptptr partition, int size1, int p, ptptr pt,
                   int *Pct, ptst *P, ptst *Q, int *PP, lltptr *L,
                   int *Lsize, lltptr bigll) {

  int size2;

 if (p<2) return;
 partition[*Pct].pt = p%2==0
                 ? sb (&pt->pt,&size1,&size2,P,Q,PP,L,Lsize,bigll)
  : sbfrac (&pt->pt,&size1,&size2,p/2/(double)p,P,Q,PP,L,Lsize,bigll);
 partition[*Pct].n = pt->n;
 pt->n = &partition[*Pct];
 (*Pct)++;
 rhamo (partition,size2,(p+1)/2,pt->n,Pct,P,Q,PP,L,Lsize,bigll);
 rhamo (partition,size1,p/2,pt,Pct,P,Q,PP,L,Lsize,bigll);
}

/**********************************************************************/

/* sb divides a point set in half.
   P1pt is the array of points to divide.  This is shrunk to contain
   only one half.
   P1ct is the size of this array.
   P2ct is the size of the newly-formed split array.  The array itself
   is returned.
   P is the original array of points.  Q is the coarsened array.
   PP is the index map to expand coarsened points.
   L is the list of expandable points.  Lsize is the size of PP.
   bigll is a contiguous list of linked list nodes corresponding to
   points.*/

static pta sb (pta *P1pt, int *P1ct, int *P2ct, ptst *P, ptst *Q,
               int *PP, lltptr *L, int *Lsize, lltptr bigll) {

  int vert, tot;
  double P1wt1=0, P1wt2=0, P2wt1=0, P2wt2=0, totx, toty, totwt1, totwt2;
  double m,b;
  pta P2pt = NULL, P3pt=NULL, i, j, thispt;
  lltptr ii, jj;

 *P2ct = 0;
 if (*P1ct==0) return NULL;
 vert = ham (*P1pt,*P1ct,&m,&b,P->pt); /* Find HS cut */
 for (i = *P1pt,j=NULL;i!=NULL;i = i==NULL ? *P1pt : i->n) {
  if (tcmp (distpointline (i->x,i->y,vert,m,b),0) <= 0) {
   if (j==NULL) *P1pt = i->n; else j->n = i->n;
   i->n = P3pt; /* If the point is on the dividing line */
   P3pt = i;    /* push onto the P3 list */
   i = j;
   (*P1ct)--;
  } else if (!vert ? i->y > m*i->x+b : i->x > b) {
   P1wt1 += i->wt1; /* One side of the line */
   P1wt2 += i->wt2;
  } else {
   P2wt1 += i->wt1;
   P2wt2 += i->wt2;
   if (j==NULL) *P1pt = i->n; else j->n = i->n;
   i->n = P2pt; /* Other side of the line */
   P2pt = i;    /* Push onto the P2 list */
   i = j;
   (*P2ct)++;
   (*P1ct)--;
  }
  j = i;
 }
 for (;(i=P3pt)!=NULL;) {
  gbll = bigll; /* Dividing the points on the line */
  gP = Q;        /* Sort points by weight */
  msort (&L[PP[i-P->pt]],llcmp,llnext,llsetnext);
  for (ii = L[PP[i-P->pt]],jj = NULL;ii!=NULL;ii = ii==NULL ?
                                        L[PP[i-P->pt]] : ii->n) {
   thispt = &Q->pt[ii-bigll];
   if (P1wt1 + P1wt2 < P2wt1 + P2wt2) {
    P1wt1 += thispt->wt1;
    P1wt2 += thispt->wt2; /* Divide points on the line into two sets */
   } else {
    P2wt1 += thispt->wt1;
    P2wt2 += thispt->wt2; /* Put point into other set */
    if (jj == NULL) L[PP[i-P->pt]] = ii->n; else jj->n = ii->n;
    ii->n = L[*Lsize];
    L[*Lsize] = ii;
    ii = jj;
   }
   jj = ii;
  }
  totx = toty = totwt1 = totwt2 = tot = 0;
  if ((ii = L[PP[i-P->pt]])==NULL) { /* None of the points remained */
   P3pt = i->n;                      /* in this set. */
   i->n = NULL;
  } else {
   for (;ii!=NULL;ii=ii->n) {
    thispt = &Q->pt[ii-bigll];
    totx += thispt->x;
    toty += thispt->y;
    totwt1 += thispt->wt1;
    totwt2 += thispt->wt2;
    tot++;
   }
   i->x = totx/tot; /* Finding average location and total weight */
   i->y = toty/tot;
   i->wt1 = totwt1;
   i->wt2 = totwt2;
   P3pt = i->n;
   i->n = *P1pt;
   *P1pt = i;
   (*P1ct)++;
  }
  totx = toty = totwt1 = totwt2 = tot = 0;
  if ((ii = L[*Lsize])!=NULL) {
   for (;ii!=NULL;ii=ii->n) {
    thispt = &Q->pt[ii-bigll];
    totx += thispt->x;
    toty += thispt->y;
    totwt1 += thispt->wt1;
    totwt2 += thispt->wt2;
    tot++;
   }
   P->pt[P->ct].x = totx/tot;
   P->pt[P->ct].y = toty/tot;
   P->pt[P->ct].wt1 = totwt1;
   P->pt[P->ct].wt2 = totwt2;
   P->pt[P->ct].n = P2pt;
   P2pt = &P->pt[P->ct];
   PP[P->ct] = *Lsize;
   P->ct++;
   (*Lsize)++;
   (*P2ct)++;
  }
 }
 return P2pt;
}

/**********************************************************************/
/* sb divides a point set in half.
   P1pt is the array of points to divide.  This is shrunk to contain
   only one half.
   P1ct is the size of this array.
   P2ct is the size of the newly-formed split array.  The array itself
   is returned.
   frac is the fraction to divide.
   P is the original array of points.  Q is the coarsened array.
   PP is the index map to expand coarsened points.
   L is the list of expandable points.  Lsize is the size of PP.
   bigll is a contiguous list of linked list nodes corresponding to
   points.*/

static pta sbfrac (pta *P1pt, int *P1ct, int *P2ct, double frac,
                   ptst *P, ptst *Q, int *PP, lltptr *L, int *Lsize,
                   lltptr bigll) {

  int vert, tot, dln, above;
  double P1wt1=0, P1wt2=0, P2wt1=0, P2wt2=0, totx, toty, totwt1, totwt2;
  double m,b;
  pta P2pt=NULL, P3pt=NULL, i, j, thispt;
  lltptr ii, jj;

 *P2ct = 0;
 if (*P1ct==0) return NULL;
 vert = hamfrac (P1pt,*P1ct,&m,&b,frac,&above,&dln,P->pt);
 for (i = *P1pt,j=NULL;i!=NULL;i = i==NULL ? *P1pt : i->n) {
  if (!vert ? tcmp(i->y,m*i->x+b)==0 :
       dln ? tcmp(i->x,m)==0 || tcmp(i->x,b)==0 :
       tcmp(i->x,b)==0) {
   if (j==NULL) *P1pt = i->n; else j->n = i->n;
   i->n = P3pt;
   P3pt = i;
   i = j;
   (*P1ct)--;
  } else if (!vert ? (above ? i->y > m*i->x+b : i->y < m*i->x+b) :
        dln ? (above ? i->x > m&&i->x < b : i->x < m||i->x > b):
        above ? i->x > b : i->x < b) {
   P1wt1 += i->wt1;
   P1wt2 += i->wt2;
  } else {
   P2wt1 += i->wt1;                      /* This is similar to sb */
   P2wt2 += i->wt2;
   if (j==NULL) *P1pt = i->n; else j->n = i->n;
   i->n = P2pt;
   P2pt = i;
   i = j;
   (*P2ct)++;
   (*P1ct)--;
  }
  j = i;
 }
 for (;(i=P3pt)!=NULL;) {
  gbll = bigll;
  gP = Q;
  msort (&L[PP[i-P->pt]],llcmp,llnext,llsetnext);
  for (ii = L[PP[i-P->pt]],jj = NULL;ii!=NULL;ii = ii==NULL ?
                                            L[PP[i-P->pt]] : ii->n) {
   thispt = &Q->pt[ii-bigll];
   if (frac*(P1wt1+P1wt2+P2wt1+P2wt2)-(P1wt1+P1wt2) >
        (1-frac)*(P1wt1+P1wt2+P2wt1+P2wt2)-(P2wt1+P2wt2)) {
    P1wt1 += thispt->wt1;
    P1wt2 += thispt->wt2;
   } else {
    P2wt1 += thispt->wt1;
    P2wt2 += thispt->wt2;
    if (jj == NULL) L[PP[i-P->pt]] = ii->n; else jj->n = ii->n;
    ii->n = L[*Lsize];
    L[*Lsize] = ii;
    ii = jj;
   }
   jj = ii;
  }
  totx = toty = totwt1 = totwt2 = tot = 0;
  if ((ii = L[PP[i-P->pt]])==NULL) {
   P3pt = i->n;
   i->n = NULL;
  } else {
   for (;ii!=NULL;ii=ii->n) {
    thispt = &Q->pt[ii-bigll];
    totx += thispt->x;
    toty += thispt->y;
    totwt1 += thispt->wt1;
    totwt2 += thispt->wt2;
    tot++;
   }
   i->x = totx/tot;
   i->y = toty/tot;
   i->wt1 = totwt1;
   i->wt2 = totwt2;
   P3pt = i->n;
   i->n = *P1pt;
   *P1pt = i;
   (*P1ct)++;
  }
  totx = toty = totwt1 = totwt2 = tot = 0;
  if ((ii = L[*Lsize])!=NULL) {
   for (;ii!=NULL;ii=ii->n) {
    thispt = &Q->pt[ii-bigll];
    totx += thispt->x;
    toty += thispt->y;
    totwt1 += thispt->wt1;
    totwt2 += thispt->wt2;
    tot++;
   }
   P->pt[P->ct].x = totx/tot;
   P->pt[P->ct].y = toty/tot;
   P->pt[P->ct].wt1 = totwt1;
   P->pt[P->ct].wt2 = totwt2;
   P->pt[P->ct].n = P2pt;
   P2pt = &P->pt[P->ct];
   PP[P->ct] = *Lsize;
   P->ct++;
   (*Lsize)++;
   (*P2ct)++;
  }
 }
 return P2pt;
}

/**********************************************************************/
/* ham applies the ham sandwich theorem to biweighted points.
   P is the set of points to divide.  ct is the size of the array.
   m and b will be the slope and y-intercept of the line.  bigp is the
   master array of points. It returns whether the line is vertical.*/

static int ham (pta P, int ct, double *m, double *b, pta bigp) {

  int vert;
  double m1a, m1b, m2a, m2b;
  pttp p;
  lnst L1, L2;

 L1.ln = emalloc (ct*sizeof(lntp));  /* Convert points to dual lines. */
 L2.ln = emalloc (ct*sizeof(lntp));
 VertCheck (&L1,P,1,&m1a,&m1b,0.5,bigp);/* Check for vertical lines*/
 VertCheck (&L2,P,0,&m2a,&m2b,0.5,bigp);
 if (tcmp(m1a,m2b)<=0 && tcmp(m2a,m1b)<=0) {
  vert = 1;
  *b = (MAX(m1a,m2a) + MIN(m1b,m2b))/2;
  *m = 0;
 } else {
  p = IntersPt (&L1,&L2,0.5);
  vert = 0;
  *m = p.x;
  *b = -p.y;
 }
 free (L1.ln);
 free (L2.ln);
 return vert;
}

/**********************************************************************/
/* hamfrac applies the ham sandwich theorem to biweighted points.
   P is the set of points to divide.  ct is the size of the array.
   m and b will be the slope and y-intercept of the line.  bigp is the
   master array of points. frac is the fraction in which the points
   should be divided.  above indicates the side of the line on which
   the smaller piece is located.  dln indeicates whether the submarine
   sandwich theorem had to be invoked.  It returns whether the line is
   vertical.*/

static int hamfrac (pta *P, int ct, double *m, double *b, double frac,
                    int *above, int *dln, pta bigp) {

  int vert;
  pta l, r;
  double m1a, m1b, m2a, m2b, m3a, m3b, m4a, m4b;
  pttp p;
  lnst L1, L2;

 L1.ln = emalloc (ct*sizeof(lntp));
 L2.ln = emalloc (ct*sizeof(lntp));
 VertCheck (&L1,*P,1,&m1a,&m1b,frac,bigp);
 VertCheck (&L2,*P,0,&m2a,&m2b,frac,bigp);
 if (tcmp(m1a,m2b)<=0 && tcmp(m2a,m1b)<=0) {
  vert = 1;
  *above = 0;
  *dln = 0;
  *m = 0;
  *b = (MAX(m1a,m2a) + MIN(m1b,m2b))/2;
 } else {
  VertCheck (&L1,*P,1,&m3a,&m3b,1-frac,bigp);
  VertCheck (&L2,*P,0,&m4a,&m4b,1-frac,bigp);
  if (tcmp(m3a,m4b)<=0 && tcmp(m4a,m3b)<=0) {
   vert = 1;
   *above = 1;
   *dln = 0;
   *b = (MAX(m3a,m4a) + MIN(m3b,m4b))/2;
  } else if ((m1a > m2a && m3a > m4a) || (m1a < m2a && m3a < m4a)) {
   p = IntersPt (&L1,&L2,frac);
   vert = 0;
   *above = 1;
   *dln = 0;
   *m = p.x;
   *b = -p.y;
  } else {
   vert = 1;
   linsplit (P,frac,&l,&r,above);
   *m = MIN (l->x,r->x);
   *b = MAX (l->x,r->x);
   *dln = 1;
  }
 }
 free (L1.ln);
 free (L2.ln);
 return vert;
}

/**********************************************************************/
/* tcmp compares a and b with respect to a comparison tolerance. */

static int tcmp (double a, double b) {

 if (IsEqual (a,b)) return 0;
 if (a < b) return -1;
 return 1;
}

/**********************************************************************/
/* VertCheck checks whether the ham sandwich cut is vertical while
   converting points to dual lines.  L is the list of dual lines.  P
   is the list of points.  first indicates which weight is
   being used. frac is the fraction to which they are to be divided.
   bigp is the master array of points. */

static void VertCheck (lnst *L, pta P, int first,
                       double *ma, double *mb, double frac, pta bigp) {

  int pos, d, ii;
  double dwt;
  pta i;

 L->wt = 0;
 ii=0;
 for (i=P;i!=NULL;i=i->n) {
  L->ln[ii].m = i->x;
  L->ln[ii].b = -i->y;
  L->ln[ii].wt = first ? i->wt1 : i->wt2;
  L->wt += L->ln[ii].wt;
  L->ln[ii].index = i-bigp;
  ii++;
 }
 L->ct = ii;
 pos = median (L->ln,sizeof(lntp),L->ct,L->wt*frac,mcmp,Lwt,&d,&dwt);
 *ma = L->ln[pos].m;
 *mb = d?L->ln[pos+1].m:*ma;
}

/**********************************************************************/
/* IntersPt locates the median interesection point of two sets of
   weighted lines.  L1 and L2 are the lines; frac is the fraction in
   which they are to be divided.  This returns the intersection point.*/

static pttp IntersPt (lnst *L1, lnst *L2, double frac) {

  double k1, k2;
  int abort;
  inftype l,r;
  pttp p;

 abort = 0;
 k1 = L1->wt*frac;
 k2 = L2->wt*frac;
 l.inf = -1;
 l.val = 0;
 r.inf = 1;
 r.val = 0;
 while (1) {   /* Shrink the lines until they are small enough.*/
  if (L1->ct < SHRINKSIZE && L2->ct < SHRINKSIZE) break;
  if ((abort = L1->ct > L2->ct ?  Shrink (L1,L2,&k1,&k2,&l,&r,&p) :
                          Shrink (L2,L1,&k2,&k1,&l,&r,&p)) > 0) break;
 }              /* Then use brute force to find the intersection. */
 if (abort != 1) p = FindPt (L1,L2,k1,k2,&l,&r);
 return p;
}

/**********************************************************************/
/* This checks whether a and b are equal with respect to a comparision
   tolerance. */

static int IsEqual (double a, double b) {

 if (fabs (a) > fabs (b)) {
   double t = b;
  b = a;
  a = t;
 }
 if (b==0) return 1;
 if (a==0) {
  if (fabs (b) < TOLERANCE) return 1;
  return 0;
 }
 if (fabs (b-a) / fabs (a) < TOLERANCE) return 1;
 return 0;
}

/**********************************************************************/
/* median is a generic weighted median function accepting an array A
   containing n elements of size bytes each. k is the weight at which
   the median should be found.  cmp and wt are the comparision and
   weight functions for the elements.  d indicates whether the median
   falls between two elements. dwt represents how much weight overlaps
   the median. */

static int median (void *A, int size, int n, double k, cfunc cmp,
                   wfunc wt, int *d, double *dwt) {

  int f,l,q,min,i;
  double prevwt, twt, w;
  void *t;

 prevwt = 0;
 f = 0;
 l = n-1;
 while (1) {
  twt = prevwt;
  q = partition (A,size,f,l,&twt,cmp,wt); /* locate a partition for */
  w = twt + wt ((char*)A+q*size);         /* this array */
  if ((q==0 || tcmp(twt,k)<0) && tcmp(w,k)>=0) break;
  if (w < k) {
   f = q+1;
   prevwt = w;
  } else l = q-1;
 }
 if (tcmp(w,k)==0 && q < n-1) {
  if (q==l) min = l+1;
  else {
   min = q+1;
   i = q+2;
   for (i=q+2;i<=l;i++)/* median falls between two elements, find both*/
    if (cmp((char*)A+i*size,(char*)A+min*size) < 0) min = i;
  }
  t = emalloc (size);
  memcpy (t,(char*)A+min*size,size);
  memcpy ((char*)A+min*size,(char*)A+(q+1)*size,size);
  memcpy ((char*)A+(q+1)*size,t,size);
  free (t);
  *d = 1;
 } else {
  *d = 0;
 }
 *dwt = k - twt;
 return q;
}

/**********************************************************************/
/* Compare the slopes of two lines */

static int mcmp (const void *x, const void *y) {

 if (((lntp *)x)->m < ((lntp *)y)->m) return -1;
 if (((lntp *)x)->m > ((lntp *)y)->m) return 1;
 return ((lntp *)x)->index - ((lntp *)y)->index;
}

/**********************************************************************/
/* General weight function */

static double Lwt (const void *x) {

 return ((lntp *)x)->wt;
}

/**********************************************************************/
/* Shrink prunes the set of lines containing the ham sandwich
   intersection.  L1 and L2 are the sets of lines; k1 and k2 are the
   weights weights to lie below the median lines.  l and r are the
   endpoints of the region.  pt is the point found, if one is.  It
   returns true if an intersection point happened to be found during the
   shrinking process. */

static int Shrink (lnst *L1, lnst *L2, double *k1, double *k2,
                   inftype *l, inftype *r, pttp *pt) {

  int r1,r2,vct,i,abort,lct,f,e;
  trap TV;
  lna temp;
  lnst newL1;
  infarr V;
  double inner, bct, b;

 lct = 0;
 V = emalloc ((RANDOMGUESS+2)*sizeof(inftype));
 newL1.ln = emalloc (L1->ct*sizeof(lntp));
 while (1) {
  vct = 0;
  for (i=0;i < RANDOMGUESS;i++) {
   r1 = rand()%L1->ct;
   r2 = rand()%L1->ct;      /* Find random intersections */
   BuildArray (&L1->ln[r1],&L1->ln[r2],V,&vct,l,r);
  }
  if ((abort = SearchPt (L1,L2,V,&vct,*k1,*k2,pt,l,r,&f,&e))) break;
  newL1.ct = 0;
  newL1.wt = 0;
  GetCorner (L1,&V[f],*k1,&TV[0],&TV[3]);
  GetCorner (L1,&V[e],*k1,&TV[1],&TV[2]);
  if (DCheck (L1,TV)) {
   bct = 0;
   for (i=0;i < L1->ct;i++) {
    inner = IntersTrap (&(L1->ln[i]),TV,&b);
    if (tcmp(inner,0)>0) {
     newL1.ln[newL1.ct] = L1->ln[i];
     newL1.ln[newL1.ct].wt = inner;
     newL1.ct++;
     newL1.wt += inner;
    }
    bct += b;
   }
   L1->ct = newL1.ct;
   L1->wt = newL1.wt;
   temp = newL1.ln;
   newL1.ln = L1->ln;
   L1->ln = temp;
   *k1 -= bct;
   *l = V[f];
   *r = V[e];
   break;
  } else if (lct==L1->ct) {
   abort = 2; break;
  } else lct++;
 }
 free (newL1.ln);
 free (V);
 return abort;
}

/**********************************************************************/
/* FindPt uses a brute force search technique to the find the
   intersection point of the medians of two sets of lines.  L1 and L2
   are the lines; k1 and k2 represent the weight that falls below these
   medians.  l and r are the left and right endpoints and the fucntion
   returns an intersection point. */

static pttp FindPt (lnst *L1, lnst *L2, double k1, double k2,
                    inftype *l, inftype *r) {

  int i,j,vct,f,e;
  infarr V;
  pttp pt;

 V = emalloc ((L1->ct*L2->ct+2)*sizeof(inftype));
 vct = 0;
 for (i=0;i<L1->ct;i++)
  for (j=0;j<L2->ct;j++)
   BuildArray (&L1->ln[i],&L2->ln[j],V,&vct,l,r);
 if (!SearchPt (L1,L2,V,&vct,k1,k2,&pt,l,r,&f,&e)) {
  fprintf (stderr,"ERROR:  Can't find an Intersection Point\n");
  exit (EXIT_FAILURE);/* IN THEORY, this should never get executed. */
 }                    /* Intersection points always exist. */
 free (V);
 return pt;
}

/**********************************************************************/
/* Partition is a general partition for the weighted median function.
   A is the array containing elements of size bytes.  f is the first
   and l is the last element in the subarray of A that we need to
   partition. prevwt is the amount of weight preceding the partition,
   cmp and wt are generic comparision and weight functions. */

static int partition (void *A, int size, int f, int l, double *prevwt,
                      cfunc cmp, wfunc wt) {

  int r,p,q;
  void *Af, *t;

 Af = emalloc (size);
 t = emalloc (size);
 r = rand()%(l-f+1) + f; /*Choose a random element */
 memcpy (Af,(char*)A+r*size,size);/*Make this the pivot */
 memcpy ((char*)A+r*size,(char*)A+f*size,size);
 p = f+1;
 q = l;
 while (p<=q) {
  for (;p!=l+1 && cmp ((char*)A+p*size,Af)<= 0;p++)
   *prevwt += wt ((char*)A+p*size);
  for (;q!=f && cmp (Af,(char*)A+q*size) < 0;q--);
  if (p<q) {
   memcpy (t,(char*)A+p*size,size); /*move elements around pivot */
   memcpy ((char*)A+p*size,(char*)A+q*size,size);
   memcpy ((char*)A+q*size,t,size);
   *prevwt += wt ((char*)A+p*size);
   p++;
   q--;
  }
 }
 memcpy ((char*)A+f*size,(char*)A+q*size,size);
 memcpy ((char*)A+q*size,Af,size);/*put pivot in proper place */
 free (Af);
 free (t);
 return q;
}

/**********************************************************************/
/* BuildArray adds the abcissa of the intersection of two lines to
   an array.  ln1 and ln2 are the lines, V is the array already
   containing vct elements.  l and r are the left and right endpoints
   of the region. */

static void BuildArray (lntp *ln1, lntp *ln2, infarr V, int *vct,
                        inftype *l, inftype *r) {

  inftype m;

 if (tcmp(ln1->m,ln2->m)!=0) {/*lines are not parallel */
  m.inf = 0;
  m.val = (ln2->b-ln1->b)/(ln1->m-ln2->m);
  if (infcmp (l,&m)<0 && infcmp (&m,r)<0) {/*intersection is in region*/
   V[*vct] = m;
   (*vct)++;
  }
 }
}

/**********************************************************************/
/*SearchPt narrows donw the region in which the median intersection
  can possibly be.  L1 and L2 are the sets of lines.  V is the array
  of intersections of size vct.  k1 and k2 are the weights that must
  fall below the medians.  pt is the intersection point (if found)
  l and r are the endpoints of the region.  f and e are the endpoints
  of the smaller region containing the intersection.  It returns true
  if and only if an intersection point happens to be found during
  this procedure. */

static int SearchPt (lnst *L1, lnst *L2, infarr V, int *vct, double k1,
                     double k2,
                     pttp *pt, inftype *l, inftype *r, int *f, int *e) {

  int abort, t, Apos, Bpos, Cpos, Dpos, dA, dB, dC, dD;
  dwarr A, B, C, D;
  double dwt, vAa, vAb, vBa,  vBb, vCa, vCb, vDa, vDb;

 A = emalloc (L1->ct*sizeof(wr));
 B = emalloc (L1->ct*sizeof(wr));
 C = emalloc (L2->ct*sizeof(wr));
 D = emalloc (L2->ct*sizeof(wr));
 V[*vct] = *l;
 V[*vct+1] = *r;
 *vct += 2;
 sort (V,vct);
 *f = 0;
 *e = *vct-1;
 abort = 0;
 while (*e!=*f+1) { /*Use binary search to narrow down region */
  t = (*f+*e)/2;
  computeintersections(L1,&V[*f],A);
  Apos = median (A,sizeof(wr),L1->ct,k1,Acmp,Awt,&dA,&dwt);
  vAa = A[Apos].val;
  vAb = dA?A[Apos+1].val:vAa;
  computeintersections(L1,&V[t],B);
  Bpos = median (B,sizeof(wr),L1->ct,k1,Acmp,Awt,&dB,&dwt);
  vBa = B[Bpos].val;
  vBb = dB?B[Bpos+1].val:vBa;
  computeintersections(L2,&V[*f],C);
  Cpos = median (C,sizeof(wr),L2->ct,k2,Acmp,Awt,&dC,&dwt);
  vCa = C[Cpos].val;
  vCb = dC?C[Cpos+1].val:vCa;
  computeintersections(L2,&V[t],D);
  Dpos = median (D,sizeof(wr),L2->ct,k2,Acmp,Awt,&dD,&dwt);
  vDa = D[Dpos].val;   /* Check whether lines have crossed */
  vDb = dD?D[Dpos+1].val:vDa;
  if (V[*f].inf == 0 && tcmp(vAa,vCb)<=0 && tcmp(vCa,vAb) <= 0) {
   abort = 1;
   pt->x = V[*f].val;
   pt->y = (MAX(vAa,vCa) + MIN (vAb,vCb))/2;
   break;
  }                /* We have by chance found an intersection point*/
  if (V[t].inf == 0 && tcmp(vBa,vDb)<=0 && tcmp(vDa,vBb) <= 0) {
   abort = 1;
   pt->x = V[t].val;
   pt->y = (MAX(vBa,vDa) + MIN (vBb,vDb))/2;
   break;
  }
  if ((tcmp(vAa,vCa) < 0 &&
       tcmp(vBa,vDa) > 0)
   || (tcmp(vAa,vCa) > 0 &&   /* We have split the region in half */
       tcmp(vBa,vDa) < 0)) *e = t; else *f = t;
 }
 free (A);
 free (B);
 free (C);
 free (D);
 return abort;
}

/**********************************************************************/
/*This locates the corners of a trapezoid for narrowing.  L1 is the set
  of lines, i is the abcissa, k1 is the amount of weight below the
  median.  t1 and t2 are the corners. */

static void GetCorner (lnst *L1, inftype *i, double k1, trapcoord *t1,
                       trapcoord *t2) {

  dwarr A;
  int Apos;
  int d;

 t1->x = t2->x = *i;
 A = emalloc (L1->ct*sizeof(wr));
 computeintersections(L1,i,A);
 Apos = median (A,sizeof(wr),L1->ct,MAX(k1-L1->wt/6,0),Acmp,Awt,&d,
                &t1->dwt);
 t1->y = A[Apos].val;
 t1->index = A[Apos].index;
 Apos = median (A,sizeof(wr),L1->ct,MIN(L1->wt,k1+L1->wt/6),Acmp,Awt,&d,
                &t2->dwt);
 t2->y = A[Apos].val;
 t2->index = A[Apos].index;
 free (A);
}

/**********************************************************************/
/* This checks both sides of the trapezoid to make sure the proper
   amount of weight falls outside of it. L is the set of lines; T
   is the trapezoid. */

static int DCheck (lnst *L, trap T) {

  int i;
  double twt;

 twt = 0;
 for (i = 0;i != L->ct && tcmp(twt,L->wt/3)<=0;i++)
  twt += Sint (&L->ln[i],T[0],T[1]);
 if (tcmp(twt,L->wt/3)>0) return 0;
 twt = 0;
 for (i=0;i != L->ct && tcmp(twt,L->wt/3)<=0;i++)
  twt += Sint (&L->ln[i],T[3],T[2]);
 return tcmp(twt,L->wt/3) <= 0;
}

/**********************************************************************/
/* IntersTrap computes how much weight intersects the trapezoid.
   ln is the set of lines, TV is the trapezoid, bottom is the amount
   of weight falling below. */

static double IntersTrap (lntp *ln, trap TV, double *bottom) {

  double top;

 *bottom = MIN (wtbelow(ln,TV[0]),wtbelow(ln,TV[1]));
 top = ln->wt - MAX (wtbelow(ln,TV[2]),wtbelow(ln,TV[3]));
 return ln->wt - (*bottom+top);
}

/**********************************************************************/
/* infcmp uses tolerance arithmetic to compare two extended real
   numbers. */

int infcmp (inftype *a, inftype *b) {

 return infcmpc (a,b,0);
}

/**********************************************************************/
/* sort is a duplicate-deleting sort of inftypes.  A is the array of
   inftypes of size n. */

static void sort (infarr A, int *n) {

  int i,ct;

 if (*n==0) return;
 qsort (A,*n,sizeof(inftype),infexcmp);
 for (ct=i=0;i<*n;i++)
  if (infcmp (&A[i],&A[ct])!=0) A[++ct] = A[i];
 *n = ct+1;
}

/**********************************************************************/
/* computeintersections computes the intersections that a set of lines
   has with a vertical line.  L is the set of lines, x is the
   x-intercept of the vertical line.  The intersections are stored in
   A. */

static void computeintersections (lnst *L, inftype *x, dwarr A) {

  int i;

 for (i=0;i<L->ct;i++) {
  A[i].val = x->inf==0 ? L->ln[i].m*x->val + L->ln[i].b :
             x->inf * L->ln[i].m;
  A[i].wt = L->ln[i].wt;
  A[i].index = L->ln[i].index;
 }
}

/**********************************************************************/
/* Acmp compares two indexed points; first by value, second by index. */

static int Acmp (const void *x, const void *y) {

 if (((wr *)x)->val < ((wr *)y)->val) return -1;
 if (((wr *)x)->val > ((wr *)y)->val) return 1;
 return ((wr *)x)->index - ((wr *)y)->index;
}

/**********************************************************************/
/* This retrieves the weight of a weighted point */

static double Awt (const void *x) {

 return ((wr *)x)->wt;
}

/**********************************************************************/
/* This computes how much weight has crossed a side of a trapezoid.
   ln is the line; t1 and t2 are the coordinates of the trapezoid. */

static double Sint (lntp *ln, trapcoord t1, trapcoord t2) {

 return fabs (wtbelow(ln,t1) - wtbelow(ln,t2));
}

/**********************************************************************/
/* This finds how much weight of an individual line falls below the
   coordinate of a trapezoid.  ln is the line; t is the coordinate. */

static double wtbelow (lntp *ln, trapcoord t) {

  double w, v;
  int c;

 v = t.x.inf==0 ? ln->m*t.x.val+ln->b :
     t.x.inf==1 * ln->m;
 c = tcmp (t.y,v);
 if (c<0 || (c==0 && t.index < ln->index)) return 0;
 if (c>0 || (c==0 && t.index > ln->index)) return w = ln->wt;
 return t.dwt;
}

/**********************************************************************/
/* infcmpc compares two extended real numbers a and b.  Exact is 1 when
   exact comparison is needed; 0 if tolerance is needed. */

static int infcmpc (inftype *a, inftype *b, int exact) {

 if (abs(a->inf) == 1) {
  if (b->inf == a->inf) return 0;
  return a->inf;
 }
 if (abs(b->inf) == 1) return -b->inf;
 if (exact) {
  if (a->val < b->val) return -1;
  if (a->val > b->val) return 1;
  return 0;
 }
 return tcmp (a->val,b->val);
}

/**********************************************************************/
/* This compares two extended real numbers with exact comparison. */

static int infexcmp (const void *a, const void *b) {

 return infcmpc ((inftype *)a,(inftype *)b,1);
}

/**********************************************************************/
/* linsplit finds a submarine sandwich cut for an array of weighted
   points.  A is the array, f is the fraction, l and r are the endpoints
   of the cut.  wrap indicates whether l precedes r or r precedes l. */

static void linsplit (pta *A, double frac, pta *l, pta *r, int *wrap) {

  double iwt1, iwt2, owt1, owt2, twt1, twt2;

 msort (A,dcmp,ptanext,ptasetnext);
 totalwt (*A,&twt1,&twt2);
 *l = *r = *A;
 iwt1 = iwt2 = 0;
 owt1 = twt1 - (*A)->wt1;
 owt2 = twt2 - (*A)->wt2;
 *wrap = 1;
 while (!wtOK (iwt1,owt1,twt1,frac) || !wtOK (iwt2,owt2,twt2,frac)) {
  if (iwt1 <= frac*twt1)
   Advance (&iwt1,&iwt2,&owt1,&owt2,*l,r,*A,wrap);
  else
   Advance (&owt1,&owt2,&iwt1,&iwt2,*r,l,*A,wrap);
 }
}

/**********************************************************************/
/* dcmp compares the abcissae of two points. */

static int dcmp (const void *a, const void *b) {

 if (((pttp *)a)->x==((pttp *)b)->x) return 0;
 if (((pttp *)a)->x < ((pttp *)b)->x) return -1;
 return 1;
}

/**********************************************************************/
/* This computes the totals of both weights of an array A of weighted
   points.  twt1 and twt2 are the total weights. */

static void totalwt (pta A, double *twt1, double *twt2) {

  pttp *i;

 *twt1 = *twt2 = 0;
 for (i=A;i!=NULL;i=i->n) {
  *twt1 += i->wt1;
  *twt2 += i->wt2;
 }
}

/**********************************************************************/
/* wtOK decides whether the amount of weight within and without a
   trapezoid is acceptable.  iwt is the inner weight; owt is the
   outer weight; twt is the total weight and frac is the desired
   fraction. */

static int wtOK (double iwt, double owt, double twt, double frac) {

 return iwt <= frac*twt && owt <= (1-frac)*twt;
}

/**********************************************************************/
/* Advance recomputes weights based on the slight advancing of a
   sub cut.  iwt1 and iwt2 are the inner weights; owt1 and owt2 are
   the outer weights.  A is the linked list of points; l and r are
   the left and right endpoints of the cut.  wrap indicates whether
   we may wrap back to the beginning. */

static void Advance (double *iwt1, double *iwt2, double *owt1,
                     double *owt2, pta l, pta *r, pta A, int *wrap) {

 *iwt1 += (l==*r)?0:(*r)->wt1;
 *iwt2 += (l==*r)?0:(*r)->wt2;
 *r = (*r)->n==NULL ? (*wrap=0,A) : (*r)->n;
 *owt1 -= (l==*r)?0:(*r)->wt1;
 *owt2 -= (l==*r)?0:(*r)->wt2;
}

/**********************************************************************/
/* emalloc attempts to allocate space, halting in error if it can't. */

static void *emalloc (size_t size) {

  void *p;

 if ((p=malloc (size))==NULL) {
  fprintf (stderr,"Out of memory!\n");
  exit (EXIT_FAILURE);
 }
 return p;
}

/**********************************************************************/
/* ecalloc attempts to allocate and initialize space, halting in error
   if it can't. */

static void *ecalloc (size_t size) {

  void *p;

 if ((p=calloc (size,1))==NULL) {
  fprintf (stderr,"Out of memory!\n");
  exit (EXIT_FAILURE);
 }
 return p;
}

/**********************************************************************/

/* Coarsen makes a grid and puts points into the appropriate grid
   region.
   pt is the array of points, ct is the size of the array, grid
   will contain the boundaries of the gridded region.  L will be a grid
   of linked lists of points corresponding to the points collapsed into
   a region.  Lsize will be the number of non-empty grid regions.
   bigll will contain a contiguous region of linked list nodes
   corresponding to points.  PP will contain an identity array of
   integers whose size is equal to the number of non-empty grid regions.
   Coarsen returns a an array of points corresponding to the non-empty
   grid regions. */

static ptst *coarsen (pta pt, int ct, bound *grid, lltptr **L,
                      int *Lsize, lltptr *bigll, int **PP) {

  ptst *Q;
  int i, xloc, yloc;
  double xwid, ywid;
  pttp *cgraph;
  lltptr P, *L2;

 *PP = emalloc (2*CSIZE*CSIZE*sizeof(int));
 Q = emalloc (sizeof (ptst));
 Q->pt = emalloc (2*CSIZE*CSIZE*sizeof (pttp));
 *L = ecalloc (2*CSIZE*CSIZE*sizeof (lltptr));
 if (ct ==0) {
  Q->ct = 0; return Q;
 }
 L2 = ecalloc (CSIZE*CSIZE*sizeof (lltptr));
 grid->minx = grid->maxx = pt[0].x;  /* Compute boundaries of grid */
 grid->miny = grid->maxy = pt[0].y;
 for (i=1;i<ct;i++) {
  if (pt[i].x < grid->minx)
   grid->minx = pt[i].x;
  else if (pt[i].x > grid->maxx)
   grid->maxx = pt[i].x;
  if (pt[i].y < grid->miny)
   grid->miny = pt[i].y;
  else if (pt[i].y > grid->maxy)
   grid->maxy = pt[i].y;
 }
 xwid = (grid->maxx - grid->minx)/CSIZE; /* Width of each section */
 ywid = (grid->maxy - grid->miny)/CSIZE;
 cgraph = ecalloc (CSIZE*CSIZE*sizeof (pttp));
 *bigll = emalloc (ct*sizeof (llt));
 for (i=0;i<ct;i++) {
  xloc = ceil ((pt[i].x-grid->minx)/xwid);/* Put each point in its */
  xloc = MIN (CSIZE-1,MAX (xloc,0));      /* appropriate location */
  yloc = ceil ((pt[i].y-grid->miny)/ywid);
  yloc = MIN (CSIZE-1,MAX (yloc,0));
  cgraph[CSIZE*xloc + yloc].wt1 += pt[i].wt1;
  cgraph[CSIZE*xloc + yloc].wt2 += pt[i].wt2;
  cgraph[CSIZE*xloc + yloc].x += pt[i].x;
  cgraph[CSIZE*xloc + yloc].y += pt[i].y;
  cgraph[CSIZE*xloc + yloc].n = /* sice of list */
               (pta)(1+(int)cgraph[CSIZE*xloc + yloc].n);
  P = &(*bigll)[i]; /* Add this to head of linked list of nodes */
  P->n = L2[CSIZE*xloc+yloc]; /* in this region */
  L2[CSIZE*xloc+yloc]=P;
 }
 for (i=xloc=0;xloc<CSIZE;xloc++)
  for (yloc=0;yloc<CSIZE;yloc++)/* If this grid is non-empty */
   if ((int)cgraph[CSIZE*xloc + yloc].n > 0) {
    (*L)[i] = L2[CSIZE*xloc+yloc]; /* Make a point that represents */
    Q->pt[i].x = cgraph[CSIZE*xloc + yloc].x/ /* all points in the */
                 (int)cgraph[CSIZE*xloc + yloc].n; /* grid */
    Q->pt[i].y = cgraph[CSIZE*xloc + yloc].y/
                 (int)cgraph[CSIZE*xloc + yloc].n;
    Q->pt[i].wt1 = cgraph[CSIZE*xloc + yloc].wt1;
    Q->pt[i].wt2 = cgraph[CSIZE*xloc + yloc].wt2;
    Q->pt[i].n = &Q->pt[i+1];
    (*PP)[i] = i;
    i++;
   }
 Q->pt[i-1].n=NULL;
 Q->ct = i;
 *Lsize = i;
 free (cgraph);
 free (L2);
 return Q;
}

/**********************************************************************/
/* msort runs merge sort on a generic linked list.  A is the list;
   cmp is the general compare function; NEXT and SETNEXT are the
   general functions for obtaining and setting the next pointer. */

static void msort (void *A,cfunc cmp,nfunc NEXT,snfunc SETNEXT) {

  void *B;

 if (*(void **)A == NULL || NEXT(*(void **)A) == NULL) return;
 B = msortsplit (*(void **)A,NEXT,SETNEXT); /* Split them. */
 msort (A,cmp,NEXT,SETNEXT); /* Sort the pieces. */
 msort (&B,cmp,NEXT,SETNEXT); /* Merge the pieces. */
 *(void **)A = msortmerge (*(void **)A,B,cmp,NEXT,SETNEXT);
}

/**********************************************************************/
/* msortsplit splits a linked list in half.  A is the list NEXT and
   SETNEXT are generic functions; it returns the other list. */

static void *msortsplit (void *A,nfunc NEXT,snfunc SETNEXT) {

  void *B, *P, *Q;

 B = NEXT (A);
 P = A;
 for (Q=B;Q!=NULL;) {
  SETNEXT (P,NEXT(Q));
  P = NEXT (P);
  if (P!=NULL) {
   SETNEXT (Q,NEXT(P));
   Q = NEXT (Q);
  } else
   Q = NULL;
 }
 return B;
}

/**********************************************************************/
/* This merges two sorted linked lists A and B.  cmp is the general
   compare function.  NEXT and SETNEXT are general pointer manipulation
   functions.  The merged list is returned. */

static void *msortmerge (void *A,void *B,cfunc cmp,nfunc NEXT,
                         snfunc SETNEXT) {

  void *C, *P;

 if (cmp (A,B)<=0) {
  C = A;
  A = NEXT(A);
 } else {
  C = B;
  B = NEXT(B);
 }
 for (P=C;A!=NULL || B!=NULL;P=NEXT(P)) {
  if (A==NULL) {
   SETNEXT (P,B);
   B = NEXT (B);
  } else if (B==NULL) {
   SETNEXT (P,A);
   A = NEXT (A);
  } else if (cmp(A,B)<=0) {
   SETNEXT (P,A);
   A = NEXT (A);
  } else {
   SETNEXT (P,B);
   B = NEXT (B);
  }
 }
 return C;
}

/**********************************************************************/
/* returns the next element in a linked list. */

static void *ptanext (const void *p) {

 return ( (pta) p)->n;
}

/**********************************************************************/
/* sets the next element in a linked list. */

static void ptasetnext (void *p, void *n) {

 ((pta)p)->n = n;
}

/**********************************************************************/
/* returns the distance between a point and a line.  x and y are the
   coordinates of the point.  m and b are the slope and intercept of
   the line; vert indicates whether the line is vertical. */

static double distpointline (double x, double y, int vert, double m,
                             double b){

  double deltax, deltay;

 if (vert) return fabs (x-b);
 if (tcmp (m,0)==0) return fabs (y-b);
 deltax = fabs (x - (y-b)/m);
 deltay = fabs (y - (m*x+b));
 return tcmp(deltax,0)==0||tcmp(deltay,0)==0 ? 0 :
        deltax*deltay/sqrt(deltax*deltax + deltay*deltay);
}

/**********************************************************************/
/* llcmp compares two points by dual weight. */

static int llcmp (const void *a, const void *b) {

  double d;
  pta ta, tb;

 ta = &gP->pt[(lltptr)a - gbll];
 tb = &gP->pt[(lltptr)b - gbll];
 d = ta->wt1+ta->wt2-tb->wt1-tb->wt2;
 if (d > 0) return -1;
 if (d < 0) return 1;
 return 0;
}

/**********************************************************************/
/* General linked list next function. */

static void *llnext (const void *p) {

 return ( (lltptr) p)->n;
}

/**********************************************************************/
/* General linked list setnext function. */

static void llsetnext (void *p, void *n) {

 ((lltptr)p)->n = n;
}

/**********************************************************************/
/* Reads a string from standard input.  line is the string read.  max
   is the maximum number of characters in the string.  If the line
   length exceeds max, the additional characters on the line are read
   and discarded. */

static void sreadln (char *line, int max) {

  int i,ci;

 memset (line,0,max);
 while ((ci=nextchar())== ' ');
 i = 0;
 while (1) {
  if (i >= max-1) break;
  if (ci=='\n') break;
  if (ci==EOF) break;
  if (ci >= 32 && ci <= 126) line[i++] = ci;
  ci = nextchar();
 }
 for (;ci!='\n' && ci != EOF;ci=nextchar());
 while (--i >= 0 && line[i]==' ') line[i]=0;
 if (line[0]==0 && ci==EOF) {
  printf ("\nEnd Of File\n"); exit (EXIT_FAILURE);}
}

/**********************************************************************/
/* nreadln returns a nonnegative integer read from the keyboard.  It
   checks for validity and prompts the user to reenter if invalid. */

static int nreadln (void) {

  int n,i,bad=0;
  char line[256];

 while (1) {
  sreadln (line,sizeof line);
  for (i=0;line[i]!=0;i++)
   if (line[i]<'0' || line[i] > '9') {bad=1;break;}
  if (!bad) break;
  printf ("Invalid number.  Reenter:  ");
  bad = 0;
 }
 sscanf (line,"%d",&n);
 return n;
}

/**********************************************************************/
/* nextchar gets the next character from the keyboard after checking
   for end-of-file. */

static int nextchar(void) {

  static int last = 0;

 if (last==EOF) return EOF;
 return (last=getchar());
}