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;im,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