Colin Champion, 7 Jan 2022

procedure : results : numerical analysis : software : index

This page describes a voting method based on Bayesian decision theory which can be applied to a certain spatial model of voting . It is not fully Bayesian but nonetheless outperforms all possible Condorcet methods.

The spatial model is as follows: An infinite number of voters satisfy a distribution which is a mixture of 3 components, each of which is uniform over a unit circular disc. The centres of the discs are drawn from a standard bivariate Gaussian. 4 candidates are drawn randomly from the voter distribution.

The parameters for each election therefore comprise 3 disc centres and 4 candidate locations, each of which has 2 degrees of freedom; but the entire system is invariant under translation and rotation, so the model has only 11 degrees of freedom.

Given the parameters of an election we can compute the proportion of ballots which will rank the candidates in each of 24 possible orderings. These 24 numbers are passed to two voting functions, Condorcet/Minimax and Semi-bayes; and each voting method reports the candidate it thinks should be elected (the Condorcet/Minimax function also provides Condorcet bounds). The accuracies of each method are recorded over 1000 simulated elections.

The model was chosen under the constraints that the voter distibution must not be symmetric, since then the Median Voter Theorem makes the problem too easy, and that the solution of the equations must be numerically tractable. It is horribly artificial but this has no bearing on the conclusions, which would have been much the same if I had overcome the computational obstacles to using a mixture of Gaussians.

[Update: I have subsequently programmed a similar solver using a GMM for voters and candidates: see A Spatial Solver. I again find that the Bayesian solver outperforms all Condorcet methods, and that the results are consistent with the belief that 100% accuracy would be obtained if only I could afford exact solution of the equations.]

Condorcet LBMinimaxCondorcet UBSemi-Bayes
90.5%90.7%91.0%95.7%

Condorcet LB and UB are the lower and upper bounds on the performance of Condorcet methods, the former corresponding to a method which is invariably wrong when presented with a Condorcet cycle and the latter to one which is invariably right; neither bound is attainable. The difference between the two bounds is equal to the five Condorcet cycles in the elections generated.

The Semi-Bayes method tries to find the positions in space of the discs and candidates in order to obtain consistency with the supplied ballot frequencies. It knows the generative model but not the parameters.

The equations are not easy to solve, in part because the solution space is poorly connected (it is hard for a candidate to migrate from one disc to another). The results quoted are for 100 random starts for each election. The effect of the number of random starts on Semi-Bayes accuracy is shown in the following table.

1310100
64.9%90.8%94.8%95.7%

I would probably have done better to take the candidates from a simpler distribution, eg. a single Gaussian centred on the voter mean.

The equations have an interesting structure. Although there are 24 possible rankings, the number which occur in any given election varies from 4 to 18. I assume that every non-zero frequency reduces the dimension of the solution space until it becomes singular, while every zero frequency constrains the space without reducing its dimension. When there is more than one solution, it will often happen that all solutions agree in the best candidate to elect (sometimes as a result of additional symmetries such as reflection or relabelling).

The frequencies of the different numbers of non-zero equations, and the number of semi-Bayes errors for them (with 100 random starts), are as follows:

equations:45678 910111213 1415161718
elections:713237256 63593712839 8596140681141000
errors:12596 53–21 5211–43

Errors can arise in either of two ways: either no valid solution has been found for the equations, or solutions exist corresponding to different preferred candidates. Evidently with a single random start the former cause is the commoner, whereas the distribution of errors suggests that the latter factor dominates when there are 100 starts.

What my software does is to look for a solution from each random start. A putative solution with an associated error term is obtained for each starting position; the error term would be 0 if the equations were exactly satisfied. My code chooses the candidate for which the error term is least.

A truly Bayesian method would integrate over the solution space. A half-way house would be to use prior information to choose between alternative solutions when these exist. This would be an approximation to MAP, and could be implemented as follows: Set a threshold for the error term such that values below this threshold are assumed to correspond to true solutions. If no solution passes the threshold, accept the one which comes closest to it. If there is at least one solution passing the threshold, choose whichever has the greatest posterior likelihood. I have not attempted this.

What does this prove? I can hardly claim to have undermined the Condorcet criterion. So far as I can see its only possible merit is as a heuristic for devising simple and effective voting methods; and since Bayesian decision theory is anything other than simple (computationally speaking), the criterion’s merits are untouched by my results.

I wouldn’t dream of suggesting that any Bayesian method is worth considering in practice, and I strongly doubt than a reasonable approximation can be derived.

It is a mathematical truism with no need of experimental verification that for any evaluation condition, the optimal voting method is the Bayesian one whose prior is the evaluation’s generative model and whose loss function is its metric. (And what is real life besides an uncontrolled evaluation?) It would be possible to hypothesise that this optimum method would always satisfy the Condorcet criterion, and while this seems to me a far-fetched hypothesis, I can claim to have provided a counter-example.

Three numerical analysis problems arose in the context of my semi-Bayes method.

Equation solving: this is the big problem which in fact caused me no difficulty. I already had my own implementation of Powell’s NEWUOA algorithm for function minimisation, and I believe it to be of good quality. I applied it with a suitable objective function (validity), which is the sum of squared errors in the ballot frequencies plus the sum of squared distances of candidates from the circumference of the nearest voter disc whenever they fall outside the voter distribution.

Disc distances. The spatial model makes use of the average distance to locations on a disc from a given point. This is not numerically difficult if you know how to do it, but the knowledge was not on the internet when I started work (although a paper by Lews et al says that the full solution is present in an NBS internal report). I rederived the solution for myself and put it on Wikipedia for the benefit of anyone interested.

Disc integrals. It is also necessary to compute the area of the region of a unit disc satisfying a given set of linear inequalities. This is mathematically trivial but fussy to program (and not without numerical pitfalls).

bayestest.c : semibayes.c : discint.c : testint.c 

bayestest.c semibayes.c discint.c testint.c
utils.html: memory.h ranf.c random.h
condorcet.html: condorcet.h 
hillclimb.html: newuoa.c mdplib.c

bayestest does the evaluation; semibayes implements the semi-Bayesian method; discint computes disc integrals. The call is: 

bayestest ntrials nstarts

allowing you to specify the number of trials to perform (default 1000) and the number of random starts (default 100). Hence the following is an example call: 

bayestest 200 3

I compile with the following line: 

g++ -o bayestest -O -g bayestest.c semibayes.c discint.c ranf.c newuoa.c mdplib.c

testint was written to validate discint, which it does by comparing its result with numerical integration. I compile by 

g++ -o testint testint.c discint.c ranf.c

and execute testint without arguments; the result is 

3504:     2     0.139
4424:  3218  3438.722
5665:  2218  2048.312

This lists the most significant discrepancies from 10000 comparisons, each line showing the test number, then the number of samples satisfying the inequalities, and the expected number according to discint.

• euclid     • minimax     • discdist     • getvotes     • getlosses     • getwinner     • genelection     • main 

#include "memory.h"
#include <math.h>
#include "condorcet.h"
#include "random.h"

static double euclid(xy A,xy B) 
{ return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)) ; }
#define qscl 1                        // qscl is the radius of the voter discs

double radians(),euclid(xy A,xy B) ;  // radians is in discint
int factorial(int m) ;                // this is in condorlib
void decompress(int r,int *x,int m) ;
double discint(double **le,int neq) ; // this is a separate file
void semibayes(int **bal,double *freq,int nv,int nc,double *loss,int nstart) ;
static double pi=radians()/2 ; 

/* -------------------------------- minimax --------------------------------- */

int minimax(int **bal,double *freq,int nbal,int m)
{ int i,j,t,balno,*b,imax ;     
  double v,**r=matrix(m,m),*score=vector(m) ; 

  for(balno=0;balno<nbal;balno++) for(v=freq[balno],b=bal[balno],i=0;i<m-1;i++)
    for(j=i+1;j<m;j++) r[b[i]][b[j]] += v ; 
  for(i=0;i<m-1;i++) for(j=i+1;j<m;j++) 
  { v = r[i][j] - r[j][i] ; r[i][j] = v ; r[j][i] = -v ; }

  for(i=0;i<m;i++) for(score[i]=1<<24,j=0;j<m;j++) // arbitrary large number
    if(i!=j&&r[i][j]<score[i]) score[i] = r[i][j] ; 
  for(v=imax=t=0;t<m;t++) if(t==0||score[t]>v) { imax = t ; v = score[t] ; }
  freematrix(r) ; free(score) ; 
  if(v<0) return imax | (1<<24) ; else return imax ; 
}
/* -------------------------------- discdist -------------------------------- */

double discdist(double q)        // as seen on wikipedia
{ double x,a,b,c,oldc,adash,mult,csum,K,KminusE,qsq ; 
  if(q==0) return 2.0/3.0 ; else q = fabs(q) ;
  if(q==1) return 32/(9*pi) ;
  qsq = q * q ; 
  if(q<1) x = qsq ; else x = 1/qsq ; 
  
  // agm
  b = sqrt(1-x) ;
  c = sqrt(x) ; 
  for(a=1,oldc=c+1,csum=c*c/2,mult=1;c>1e-15&&c<oldc;csum+=mult*c*c,mult*=2)
  { oldc = c ; c = (a-b)/2 ; adash = (a+b)/2 ; b = sqrt(a*b) ; a = adash ; }
  K = pi / (2*a) ;
  KminusE = K * csum ; 
  
  if(q<1) return (4/(9*pi)) * ((5*qsq+3)*K   - (qsq+7)*KminusE) ;
  else    return (4/(9*pi)) * ((5*qsq+3)*K/q - (qsq+7)*KminusE*q) ;
}
/* -------------------------------- getvotes -------------------------------- */

void getvotes(xy *v,int nv,xy *c,int nc,int **b,double *freq)
{ int bno,cno,vno,eqno,N=factorial(nc) ;
  double q,x,y,**le=matrix(nc-1,3) ;
  xy C0,C1 ; 

  for(bno=0;bno<N;bno++) 
    for(decompress(bno,b[bno],nc),freq[bno]=vno=0;vno<nv;vno++)
  { for(x=v[vno].x,y=v[vno].y,eqno=0;eqno<nc-1;eqno++)
    { C0 = xy(c[b[bno][eqno]].x-x,c[b[bno][eqno]].y-y) ; 
      C1 = xy(c[b[bno][eqno+1]].x-x,c[b[bno][eqno+1]].y-y) ; 
      le[eqno][0] = 2 * (C1.x-C0.x) ; 
      le[eqno][1] = 2 * (C1.y-C0.y) ; 
      le[eqno][2] = C1.x*C1.x + C1.y*C1.y - C0.x*C0.x - C0.y*C0.y ; 
    }
    freq[bno] += discint(le,nc-1) ; 
  }
  freematrix(le) ;
  for(q=bno=0;bno<N;bno++) q += freq[bno] ; 
  for(bno=0;bno<N;bno++) freq[bno] /= q ; 
}
/* -------------------------------- getlosses ------------------------------- */

void getlosses(xy *v,int nv,xy *c,int nc,double *x)
{ int cno,vno ; 
  for(cno=0;cno<nc;cno++) for(x[cno]=vno=0;vno<nv;vno++) 
    x[cno] += qscl * discdist(euclid(c[cno],v[vno])/qscl) / nv ; 
}
int getwinner(double *loss,int nc)
{ int cno,winner ; 
  for(cno=0;cno<nc;cno++) if(cno==0||loss[cno]<loss[winner]) winner = cno ; 
  return winner ; 
}
/* ------------------------------- genelection ------------------------------ */

void genelection(xy *v,int nv,xy *c,int nc)
{ int vno,cno ; 
  double h,theta,x,y ; 
  for(y=x=vno=0;vno<nv;vno++) 
  { v[vno] = xy(gaussv(),gaussv()) ; x += v[vno].x ; y += v[vno].y ; }
  for(vno=0;vno<nv;vno++) v[vno] = xy(v[vno].x-x/nv,v[vno].y-y/nv) ;
  for(cno=0;cno<nc;cno++)
  { vno = rani(nv) ; 
    theta = 2*pi*ranf() ; 
    h = qscl * sqrt(ranf()) ; 
    c[cno] = xy(v[vno].x+h*cos(theta),v[vno].y+h*sin(theta)) ;
  }
}
/* ---------------------------------- main ---------------------------------- */

int main(int argc,char **argv)
{ int ntests=argc<2?1000:atoi(argv[1]),nstart=argc<3?100:atoi(argv[2]) ;
  int i,j,nv=3,nc=4,N=factorial(nc),testno,rightful,winner,ndi ; 
  int **b=imatrix(N,nc),tally[]={0,0,0,0} ; 
  double q,*freq=vector(N),*loss=vector(nc),*bayesloss=vector(nc) ; 
  xy *V=xyvector(nv),*C=xyvector(nc) ;
  char *algs[]={"condorcet","minimax","psychic","semibayes"} ;

  for(i=0;i<N;i++) decompress(i,b[i],nc) ; 

  for(testno=0;testno<ntests;testno++)
  { ranset(testno) ; 
    genelection(V,nv,C,nc) ; 
    getvotes(V,nv,C,nc,b,freq) ; 
    for(ndi=i=0;i<N;i++) if(freq[i]) ndi += 1 ; 
    getlosses(V,nv,C,nc,loss) ; 
    rightful = getwinner(loss,nc) ; 
    winner = minimax(b,freq,N,nc) ; 
    if(winner==rightful) tally[0] += 1 ;               // condorcet
    if((winner&0xffff)==rightful) tally[1] += 1 ;      // minimax
    if(winner==rightful||(winner>>24)) tally[2] += 1 ; // psychic
    semibayes(b,freq,nv,nc,bayesloss,nstart) ; 
    winner = getwinner(bayesloss,nc) ; 
    if(winner==rightful) tally[3] += 1 ; 

    printf("%4d: %2d | ",testno,ndi) ;
    for(q=i=0;i<nc;i++) 
    { q += (loss[i]-bayesloss[i])*(loss[i]-bayesloss[i]) ;
      printf("(%.3f:%.3f) ",loss[i],bayesloss[i]) ; 
    }
    printf(" [%.3f] | %d",sqrt(q/nc),rightful) ; 
    if(winner!=rightful) printf("(%d)\n",winner) ; else printf("\n") ; 
    if(testno%10==9) for(i=0;i<4;i++) printf("%s=%.1f%%%s",algs[i],
                              (100*tally[i])/(testno+1.0),(i==3?"\n":"; ")) ;
  }
}

• euclid     • flatten     • unflatten     • validity     • semibayes 

#include "memory.h"
#include <math.h>
#include "condorcet.h"
#define qscl 1

static double euclid(xy A,xy B) 
{ return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)) ; }
double ranf(),gaussv() ; 
int rani(int) ; 
void ranset(int) ;

void genelection(xy *v,int nv,xy *c,int nc) ; 
void getvotes(xy *v,int nv,xy *c,int nc,int **b,double *freq) ;
double uoamin(double (*func)(double*),double *x,int n,
              double xinc,double tol,int maxiter) ;
int uoacond(),factorial(int) ;
double discdist(double q),euclid(xy,xy) ;
void getlosses(xy *v,int nv,xy *c,int nc,double *x) ;

static int **truebal,NV,NC ; 
static double *truefreq ; 

static void flatten(xy *V,int nv,xy *C,int nc,double *x) 
{ int i,k ; 
  for(k=i=0;i<nv-1;i++) { x[k++] = V[i].x ; x[k++] = V[i].y ; }
  for(i=0;i<nc;i++) { x[k++] = C[i].x ; x[k++] = C[i].y ; }
}
static void unflatten(double *x,xy *V,int nv,xy *C,int nc) 
{ int i,k ; 
  for(k=i=0;i<nv-1;i++) { V[i].x = x[k++] ; V[i].y = x[k++] ; }
  for(V[nv-1]=xy(),i=0;i<nv-1;i++) 
  { V[nv-1].x -= V[i].x ; V[nv-1].y -= V[i].y ; }
  for(i=0;i<nc;i++) { C[i].x = x[k++] ; C[i].y =  x[k++] ; }
}
/* --------------------------------- validity ------------------------------- */

static double validity(double *x)
{ int i,j,N=factorial(NC),**b=imatrix(N,NC) ; 
  double q,qmin,qdist,r,*freq=vector(N) ;
  xy *V=xyvector(NV),*C=xyvector(NC) ;
  unflatten(x,V,NV,C,NC) ; 
  getvotes(V,NV,C,NC,b,freq) ; 
  for(q=i=0;i<N;i++) 
  { for(j=0;j<NC;j++) if(b[i][j]!=truebal[i][j])
      throw up("faulty ballot no %d",i) ; 
    q += (freq[i]-truefreq[i]) * (freq[i]-truefreq[i]) ;
  }
  for(j=0;j<NC;j++) 
  { for(qmin=-1,i=0;i<NV;i++) 
    { qdist = euclid(V[i],C[j]) ; if(qmin<0||qdist<qmin) qmin = qdist ; }
    if(qmin>qscl) q += (qmin-qscl)*(qmin-qscl) ;
  }
  freeimatrix(b) ; 
  free(freq,V,C) ; 
  return q ; 
}
/* --------------------------------- semibayes ------------------------------ */

void semibayes(int **bal,double *freq,int nv,int nc,double *loss,int nstart)
{ int i,tryno,nbal=factorial(nc),n=2*(nc+nv-1) ; 
  double *x=vector(n),q,qmin ;
  xy *V=xyvector(nv),*Vmin=xyvector(nv),*C=xyvector(nc),*Cmin=xyvector(nc) ;

  truefreq = freq ; 
  truebal = bal ; 
  NC = nc ; NV = nv ; 
  ranset(1<<24) ; // for reproducible results

  for(tryno=0;tryno<nstart;tryno++)
  { genelection(V,nv,C,nc) ; 
    flatten(V,nv,C,nc,x) ; 
    q = uoamin(validity,x,n,1e-2,1e-4,100000) ;
    if(tryno==0||q<qmin) { qmin = q ; unflatten(x,Vmin,nv,Cmin,nc) ; }
  }
  getlosses(Vmin,nv,Cmin,nc,loss) ;
  free(x,V,Vmin,C,Cmin) ;
}

• radians     • bearing     • linety     • swapends     • satisfiedby     • setpt     • setdtheta     • chordty     • revisept     • crossing     • chordsegment     • euclid     • discint 

#include "memory.h"
#include <math.h>
double radians() { return 2 * 3.1415926535897932384626433832795029 ; }
static double pi=radians()/2 ; 
double radians(double x,double qlo) // reduces x mod 2pi to betw qlo and qlo+2pi
{ while(x>=qlo+2*pi) x -= 2*pi*(int)((x-qlo)/(2*pi)) ; 
  while(x<qlo) x += 2*pi*(int)((qlo+2*pi-x)/(2*pi)) ; 
  return x ; 
}
double radians(double x) { return radians(x,0) ; } 
double bearing(xy A) { return atan2(A.y,A.x) ; }

struct linety 
{ double a,b,c,theta[2],dtheta ; xy pt[2] ; 
  linety() { ; }
  linety(double *x)
  { double A,B,C,q ; 
    a = x[0] ; b = x[1] ; c = x[2] ; 
    A = a*a+b*b ; 
    if(fabs(a)<fabs(b)) // careful with arithmetic precision
    { C = c*c-b*b ;
      B = a*c ;
      q = B*B-A*C ; 
      if(q<0) { pt[0] = xy() ; return ; }
      if(B<0) q = B-sqrt(q) ; else q = B+sqrt(q) ; 
      pt[0] = xy(q/A,(c-a*q/A)/b)  ; 
      pt[1] = xy(C/q,(c-a*C/q)/b) ;
    }
    else
    { C = c*c-a*a ; 
      B = b*c ; 
      q = B*B-A*C ; 
      if(q<0) { pt[0] = xy() ; return ; }
      if(B<0) q = B-sqrt(q) ; else q = B+sqrt(q) ; 
      pt[0] = xy((c-b*q/A)/a,q/A) ; 
      pt[1] = xy((c-b*C/q)/a,C/q) ; 
    }
    for(int i=0;i<2;i++) theta[i] = bearing(pt[i]) ; 
    dtheta = radians(theta[1]-theta[0],-pi) ; 
  }
  void swapends() 
  { swap(theta[0],theta[1]) ; swap(pt[0],pt[1]) ; dtheta = -dtheta ; }
  int satisfiedby(xy X) { return a*X.x+b*X.y<c ; }
} ;
genvector(linety,linvector) ; 
struct chordty 
{ xy pt[2] ; double theta[2],dtheta ; int lineno ; 
  void setpt(xy X,int i,int flag) 
  { pt[i] = X ; 
    if(lineno<0) 
    { theta[i] = bearing(X) ; if(flag) setdtheta(dtheta>=0) ; }
  }
  void setdtheta(int f) // f=0/1 => dtheta -ve/+ve
  { dtheta = radians(theta[1]-theta[0],f?0:-2*pi) ; }
  chordty(xy A,xy B,int f) // negative f encodes equation number
  { if(f<0) { dtheta = theta[1] = theta[0] = 0 ; lineno = -(1+f) ; }
    else lineno = -1 ; 
    setpt(A,0,0) ; 
    setpt(B,1,0) ; 
    if(f>=0) setdtheta(f) ;  
  }
  chordty(chordty *u,chordty *v,int i)
  { this[0] = chordty(u->pt[1],v->pt[0],-(1+i)) ; }
  void revisept(xy X,double phi,int i)
  { pt[i] = X ; theta[i] = phi ; setdtheta(dtheta>=0) ; }
} ;
genvector(chordty,chordvector) ; 
struct crossing 
{ xy pt ; int chordno ; 
  crossing() { chordno = -1 ; }
  crossing(xy P,int c) { pt = P ; chordno = c ; }
} ;
static int chordsegment(chordty *c,linety *l,int lineno)
{ c[0] = chordty(l->pt[0],l->pt[1],-(1+lineno)) ;
  c[1] = chordty(l->pt[1],l->pt[0],l->dtheta*l->c>=0) ;
  return 2 ;
}
/* -------------------------------------------------------------------------- */

static double euclid(xy A,xy B) 
{ return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)) ; }

double discint(double **le,int neq) 
{ static int visit=-1 ; visit += 1 ; 
  int i,j,t,eqno,nchord,cno,ndel,nX ;
  double q,dtheta,x,y,dx,dy,phi[2] ;
  linety l,c,*line=linvector(neq) ; 
  crossing X[2] ;

  // find intercepts with circle
  for(t=eqno=0;eqno<neq;eqno++,t++) 
  { line[t] = linety(le[eqno]) ;
    if(line[t].pt[0]==xy())
    { if(line[t].c>=0) { t -= 1 ; continue ; } free(line) ; return 0 ; } 
    if(line[t].a==0&&line[t].b==0) 
      throw up("degenerate eqn %d: 0x+0y<=%.1f",eqno,line[t].c) ; 
  }
  neq = t ; 
  if(neq==0) { free(line) ; return pi ; } 

  chordty *chord=chordvector(4*neq) ; // extra space for reordering
  nchord = chordsegment(chord,line,0) ; 

  for(eqno=1;eqno<neq;eqno++)
  { // find intersection of the new line with each chord of the rounded polygon
    for(l=line[eqno],nX=cno=0;cno<nchord;cno++) 
      if(chord[cno].lineno>=0) // the chord is linear
      { c = line[chord[cno].lineno] ; // other equation
        q = l.a*c.b - l.b*c.a ;
        x = ( c.b*l.c - l.b*c.c ) / q ; 
        y = ( l.a*c.c - c.a*l.c ) / q ; 
        dy = chord[cno].pt[1].y - chord[cno].pt[0].y ;
        dx = chord[cno].pt[1].x - chord[cno].pt[0].x ;
        if(fabs(dy)>fabs(dx)) q = (y-chord[cno].pt[0].y) / dy ; 
        else q = (x-chord[cno].pt[0].x) / dx ; 
        if(q>=0&&q<1)
        { if(nX==2) throw up("> %d:%d:%d",visit,eqno,cno) ; 
          X[nX++] = crossing(xy(x,y),cno) ; 
        }
      }
      else for(j=0;j<2;j++) // the chord is an arc, which the line crosses twice
      { dtheta = chord[cno].dtheta ; // the sign of dtheta is anticlockwiseness
        q = radians(l.theta[j]-chord[cno].theta[0],dtheta>0?0:-2*pi) ;
        if((q>=0&&q>=dtheta)||(q<0&&q<dtheta)) continue ; 
        if(nX==2) throw up(": %d:%d:%d",visit,eqno,cno) ; 
        X[nX++] = crossing(l.pt[j],cno) ; 
      }

    if(nX==0) // the line doesn't cross the polygon
    { if(l.satisfiedby(chord[0].pt[0])) continue ;
      else { free(chord,line) ; return 0 ; }
    }
    else if(nX!=2) throw up("* %d:%d:%d",nX,eqno,visit) ;

    // we have a chord which cuts the polygon at 2 points: revise the polygon
    cno = X[0].chordno ;
    if(cno==X[1].chordno&&l.satisfiedby(chord[cno].pt[0]))
    { for(i=nchord-1;i>=cno;i--) chord[i+2] = chord[i] ; 
      nchord += 2 ; 
      dtheta = chord[cno].dtheta ; 
      for(j=0;j<2;j++) phi[j] = 
        radians(l.theta[j]-chord[cno].theta[0],dtheta>0?0:-2*pi) ;
      if((dtheta>=0&&phi[0]>=phi[1])||(dtheta<0&&phi[0]<phi[1])) 
        l.swapends() ; 
      chord[cno].revisept(l.pt[0],l.theta[0],1) ; 
      chord[cno+2].revisept(l.pt[1],l.theta[1],0) ; 
      chord[cno+1] = chordty(&chord[cno],&chord[cno+2],eqno) ; 
    }
    else if(cno==X[1].chordno) nchord = chordsegment(chord,&l,eqno) ;
    else 
    { if(l.satisfiedby(chord[cno].pt[1]))
      { cno = X[1].chordno ; 
        for(i=cno;i<nchord;i++) chord[2*neq+i-cno] = chord[i] ;
        for(i=0;i<cno;i++) chord[2*neq+i+nchord-cno] = chord[i] ; 
        for(i=0;i<nchord;i++) chord[i] = chord[2*neq+i] ; 
        swap(X[0].pt,X[1].pt) ; 
        X[1].chordno = X[0].chordno + nchord-cno ; 
        cno = X[0].chordno = 0 ;  
      }
      ndel = X[1].chordno - cno - 2 ; // no of old chords to delete
      if(ndel>0) for(j=X[1].chordno;j<nchord;j++) chord[j-ndel] = chord[j] ; 
      else if(ndel<0) for(j=nchord-1;j-ndel>=X[1].chordno;j--)
        chord[j-ndel] = chord[j] ;
      nchord -= ndel ; 
      chord[cno+1] = chordty(X[0].pt,X[1].pt,-(eqno+1)) ; 
      chord[cno].setpt(chord[cno+1].pt[0],1,1) ; 
      chord[cno+2].setpt(chord[cno+1].pt[1],0,1) ; 
    }
  } // end loop over inequality conditions to add to polygon

  chord[nchord].pt[0] = chord[nchord-1].pt[1] ; // shoelace formula
  for(q=i=0;i<nchord;i++) q += chord[i].pt[0].x*chord[i+1].pt[0].y - 
                               chord[i].pt[0].y*chord[i+1].pt[0].x ;
  for(q=fabs(q),i=0;i<nchord;i++) if(chord[i].lineno<0) // circular segments
  { dtheta = fabs(chord[i].dtheta) ; q += dtheta - sin(dtheta) ; }
  free(chord,line) ;
  return q/2 ;
}


#include "memory.h"
#include <math.h>
#include "random.h"
double radians() ;
static double pi=radians()/2 ; 
double discint(double **le,int neq) ;

int main()
{ int i,j,testno,sampleno,n=100000,m ; 
  double **a=matrix(3,3),theta,h,p,sigma ; 
  for(testno=0;testno<10000;testno++) 
  { ranset(testno) ; 
    for(i=0;i<3;i++)
    { theta = 2*pi*ranf() ;
      h = 2.4 * (ranf()-0.5) ; 
      a[i][0] = cos(theta) ; 
      a[i][1] = sin(theta) ; 
      a[i][2] = cos(2*theta) ; 
    }
    for(m=sampleno=0;sampleno<n;sampleno++)
    { theta = 2 * pi * ranf() ;
      h = sqrt(ranf()) ; 
      for(i=0;i<3;i++) 
        if(a[i][0]*h*cos(theta)+a[i][1]*h*sin(theta)>=a[i][2]) break ; 
      if(i==3) m += 1 ; 
    }
    p = discint(a,3)/pi ;
    sigma = sqrt(n*p*(1-p)) ; 
    if(fabs(m-n*p)>3.5*sigma) printf("%3d: %5d %9.3f\n",testno,m,n*p) ; 
  }
}