Colin Champion, 19 Jan 2023

This page provides a short C function to perform cubature – i.e. quadrature over a space of arbitrary dimension.

doc : algorithm : source : tests : example : download : index

cubate is an adaptive cubature subroutine which I’ve written for my own interest, and which is available to anyone who cares to use it but is not recommended, still less guaranteed. cubature by Steven G. Johnson may be mentioned as an alternative which appears to be based on sounder knowledge of the topic. (I haven’t tried it.)

Four calling sequences are provided, all prototyped in quadrate.h. The integral is over the unit hypercube, i.e. the set of x such that 0≤x_i≤1.

double cubate(double (f)(double*),int dim,double qtol,int maxeval) ;

f is a function on a vector of dim arguments which returns a single value. This value is integrated over the unit cube to within the given tolerance and the integral returned. It crashes out on failure if maxeval function evaluations are exceeded. See quadrate for a description of the tolerance. f is assumed to be written to incorporate the value of dim.

double cubate(double (f)(double*,int),int dim,double qtol,int maxeval) ;

This call differs from the preceding one in that f has the dimension dim passed to it as a second parameter.

quadratereport cubate(void (f)(double*,double *),int dim,int n,double *res,
                      double qtol,int maxeval) ;

This call integrates the vector-valued function f which computes n values at every point, and whose prototype is void f(double *x,double *y) where x is the input vector and y the output. The n-dimensional result is returned in res The function is assumed to be written with knowledge of dim and n.

The returned structure is as for quadrate, but the nround field is unused whereas the leveval field now serves a purpose. Since the user now knows the estimated precision of the integral, cubate terminates normally rather than crashing out when the maximum number of function evaluations is reached.

leveval is set up as an array by cubate (nleveval is its length). leveval[i] is the number of function evaluations at level i of the recursion. These levels exclude the initial expansion of the cube into polytopes (during which no integrations are performed), and refer only to integration over simplices. You must free(r.leveval) after return from cubate.

quadratereport cubate(void (f)(double*,double *,int,int),int dim,int n,
                      double *res,double qtol,int maxeval) ;

This call differs from the preceding one in that f is told the dimensions of the spaces, having the prototype void f(double *x,double *y,int dim,int n).

The error estimate used by cubate is hopelessly pessimistic; this seems to me its main weakness. See below.

Let us start by stating Lipson’s formula [1]. Write μ_q for the average value of a given quadratic over a d-dimensional simplex; write μ_v for its average value at the vertices; and write μ_c for its value at the centroid. Then Lipson’s formula states that:

μq =	μv + (d+1)μc
	2+d

The integral of a quadratic over a simplex can then be computed as the product of μ_q with the volume of the simplex. (I don’t know how widely Lipson’s formula is known.)

If we additionally define μ_m as the average value of the quadratic at the midpoints of the edges, then a companion formula tells us that:

μq =	2dμm + (2–d)μv
	2+d

Notice that in two dimensions – i.e. integrating over triangles – μ_v drops out.

We can now summarise our algorithm as a sequence of steps:

• The first step is a reduction of the polytope defining the region of integration. Initially this is a hypercube, but it is reduced by successive subdivisions to polytopes in which the number of triangular faces increases and the number of square faces diminishes.
• Eventually the integration region has been subdivided into dim! simplices. These can be integrated over using Lipson’s formula; an estimate is made of the error; and if the error is too large, the simplex is subdivided by bisecting the longest edge.
• If we write I as our estimate of the integral over a simplex, then the error is estimated as |I–μ_v| multiplied by the volume.
• The simplex chosen for subdivision is the one whose estimated error is largest. Following Johnson I maintain a heap of simplices in order to efficiently find it.
• The process terminates when the total estimated error is less than a threshold, which is either the tolerance if absolute, or the product of the tolerance with sum of absolute integrals over simplices if the tolerance is relative.
• I keep all simplices in memory. To facilitate access to vertices I keep them (with their associated function values) in a separate list. This saves re-evaluating f at a given vertex. The structure for a vertex contains a list of connected edges.

I use the standard recursive subdivision of a hypercube into simplices. The hypercube is represented as the set of inequalities 0 ≤x₀, x₁, x₂, ... ≤ 1. This is then partitioned into polytopes given by the sets of inequalities

• 0 ≤x₀ ≤ x₁, x₂, ... ≤ 1;
• 0 ≤x₁ ≤ x₀, x₂, ... ≤ 1; ...

and these are split in turn into polytopes such as 0 ≤x₄ ≤ x₂ ≤ x₀, x₁, ... ≤ 1; and so forth until we obtain inequalities similar to 0 ≤x₄ ≤ x₂ ≤ x₀ ≤ x₃ ≤ x₁ ≤ 1 which define simplices.

[1]. Andrew Lipson, personal communication, 2005 (“Simpson’s rule on simplices and numerical integration”).

The main weakness seems to me to lie in the estimate of the error in the integral over a simplex. This amounts to subdividing simplices until the error in a linear fit satisfies the tolerance, and then returning the quadratic fit.

I suspect it would be better (though more expensive) to look at the difference between Lipson’s formula and its companion.

At present I simply bisect a simplex’s longest edge, but giving preference to an edge which has already been bisected (because shared with another simplex). The following might be better: (i) if all edges span the entire hypercube, choose the longest; (ii) otherwise each edge must have an ‘external point’ – i.e. the function must have been evaluated at a point along the line external to the segment; choose the closest such point and perform a cubic fit using the vertices and the midpoint; and choose the edge for which the cubic fit differs most from a quadratic fit excluding the external point.

The hypercubic polytope is subdivided into dim! simplices before the integration commences. This entails a minimum amount of work which soon becomes infeasible. It would be possible to integrate over the polytopes at each level, and to curtail the expansion if the error is small enough.

To integrate over a polytope, one would evaluate the function at enough points to fit a quadratic. One would then use the subdivision into simplices to integrate the quadratic using Lipson’s formula. So the amount of work is factorial in the dimension, even if the number of function evaluations is quadratic. The points to choose would be centroids of constituent simplices, allowing the results to be reused. If there aren’t enough, there’s no point in attempting the short cut, and one should simply subdivide regardless.

Unlike quadrate, cubate makes no attempt to cope with spiky functions. It could benefit from some attention to round-off error. And there are also whatever faults arise from my lack of knowledge of the topic.

• factorial     • adj     • control     • freeslist     • eval     • addleveval     • vertex     • release     • addneighbour     • delneighbour     • addneighbour     • delneighbour     • addvertex     • release     • polytope     • release     • simplex     • release     • setcentroid     • integrate     • findedge     • increment     • extend     • fanout     • printneighbours     • deleter     • stub     • cubate    

#include "memory.h"
#include <math.h>
#include "quadrate.h"
#include "heap.h"

static int factorial(int m) 
{ int i,t ; for(t=i=1;i<=m;i++) t *= i ; return t ; }

/* ----------------------------------- adj ---------------------------------- */

struct vertex ; 
struct simplex ; 
struct polytope ; 

static simplex* debug = (simplex *) (0) ; 

struct adj
{ vertex *v,*midpt ; int count ; double d2 ; 
  adj() { v = midpt = 0 ; d2 = count = 0 ; }
  adj(vertex *vv,double q,int n) 
  { this[0] = adj() ; v = vv ; d2 = q ; count = n ; }
} ;
genvector(vertex*,vpvector) ;  // a vpvector is a vector of (vertex *)s
genvector(simplex*,spvector) ; // an spvector is a vector of (vertex *)s
genvector(void*,vvvector) ;    // a vvvector is a vector of (void *)s
genvector(adj,adjvector) ; 

/* --------------------------------- control -------------------------------- */

struct control 
{ int neval,maxeval,thinslices,mode,nv,ns,*leveval,nleveval ; 
  vertex *vlist ; 
  simplex **slist ; 
  polytope *plist ; 
  void *f ; 

  control() 
  { nleveval = maxeval = neval = thinslices = nv = ns = mode = 0 ; 
    f = 0 ; 
    vlist = 0 ; 
    plist = 0 ; 
    slist = 0 ; 
    leveval = 0 ; 
  }
  control(void *ff,int momo,int dim,int n,int mm)
  { this[0] = control() ; 
    f = ff ; 
    mode = momo ; 
    maxeval = mm ; 
    slist = spvector(factorial(dim)) ;
  }
  void release() ;
  void freeslist() { free(slist) ; slist = 0 ; ns = 0 ; }
  vertex *addvertex(int dd,int nn,double *xx) ;
  void eval(double *x,double *y,int dim,int n)
  { neval += 1 ; 
    if(mode<2)
    { if(mode==0) y[0] = ((double (*)(double *))f)(x) ; 
      else y[0] = ((double (*)(double *,int))f)(x,dim) ; 
    }
    else 
    { if(mode==2) ((void (*)(double *,double *))f)(x,y) ;
      else ((void (*)(double *,double *,int,int))f)(x,y,dim,n) ;
    }
  } 
  void addleveval(int lev)
  { int t=nleveval ; 
    if(lev>=nleveval)
    { nleveval = 2*nleveval + 20 ; 
      leveval = ivector(leveval,nleveval) ; 
      for(;t<nleveval;t++) leveval[t] = 0 ; 
    }
    leveval[lev] += 1 ; 
  }  
} ;
/* ---------------------------------- vertex -------------------------------- */

struct vertex
{ int dim,n,na,alen,serialno ; 
  double *x,*y ; 
  adj *neighbour ;
  vertex *next ;
  control *c ; 
  vertex() 
  { dim = n = na = alen = 0 ; 
    x = y = 0 ; 
    neighbour = 0 ; 
    next = 0 ; 
    c = 0 ; 
    serialno = 0x7fff ;
  }
  void release() { free(x,neighbour) ; this[0] = vertex() ; }
  vertex(int dd,int nn,control *cc,double *xx) 
  { this[0] = vertex() ; 
    c = cc ; 
    serialno = c->nv ; 
    dim = dd ; 
    n = nn ;  
    x = vector(dim+n) ; 
    y = x + dim ; 
    neighbour = adjvector(alen=2*dim) ; 
    for(int t=0;t<dim;t++) x[t] = xx[t] ; 
    c->eval(x,y,dim,n) ; 
  }
  void addneighbour(vertex *v)
  { int t ; 
    double q ; 
    for(t=0;t<na;t++) if(neighbour[t].v==v) 
    { neighbour[t].count += 1 ; return ; }
    if(na==alen) { alen += 10+alen ; neighbour = adjvector(neighbour,alen) ; }
    for(q=t=0;t<dim;t++) q += (x[t]-v->x[t]) * (x[t]-v->x[t]) ;
    neighbour[na++] = adj(v,q,1) ; 
  }
  void delneighbour(vertex *v)
  { for(int t=0;t<na;t++) if(neighbour[t].v==v) 
    { neighbour[t].count -= 1 ; 
      if(neighbour[t].count==0)
      { if(t!=na-1) swap(neighbour[t],neighbour[na-1]) ; na -= 1 ; }
      return ; 
    }
    throw up("deleting a non-existent neighbour %d of %d",
             v->serialno,serialno) ; 
  }
} ; 
void addneighbour(vertex *u,vertex *v)
{ if(u>v) v->addneighbour(u) ; else u->addneighbour(v) ; }

void delneighbour(vertex *u,vertex *v)
{ if(u>v) v->delneighbour(u) ; else u->delneighbour(v) ; }

vertex *control::addvertex(int dd,int nn,double *xx)
{ vertex *v=new vertex(dd,nn,this,xx) ;
  nv += 1 ; 
  v->next = vlist ; 
  return vlist = v ; 
}
void control::release() 
{ vertex *v ; 
  while(vlist->next) 
  { v = vlist ; vlist = vlist->next ; v->release() ; delete v ; }
  vlist->release() ; 
  delete vlist ; 
  this[0] = control() ; 
}
/* --------------------------------- polytope ------------------------------- */

struct polytope
{ double vol,*hypervol,*err ; 
  int dim,n ; 
  int ndescendents,*fixed,nfixed ; 
  control *c ; 
  void **descendents ;
  polytope()
  { vol = 0 ; hypervol = err = 0 ; 
    dim = n = 0 ; 
    ndescendents = nfixed = 0 ; 
    fixed = 0 ; 
    c = 0 ; 
  }
  polytope(int dd,int nn,control *cc,int nf) 
  { this[0] = polytope() ; 
    dim = dd ; 
    n = nn ; 
    c = cc ; 
    if((nfixed=nf)>0) fixed = ivector(nfixed) ; 
  }
  void fanout(vertex **corners) ;
  void release()
  { if(nfixed<dim-2) for(int i=0;i<ndescendents;i++) 
      ((polytope *)descendents[i])->release() ;
    free(fixed,descendents,hypervol,err) ;
  }
} ;
/* --------------------------------- simplex -------------------------------- */

struct simplex 
{ double vol,*hypervol,*err,maxerr,preverr[2] ; 
  vertex **v ; 
  int dim,n,lev ; 
  control *c ; 
  simplex() 
  { vol = lev = maxerr = preverr[1] = preverr[0] = 0 ; 
    err = hypervol = 0 ; 
    v = 0 ; 
    c = 0 ; 
  } 
  simplex(int dd,int nn,control *cc) 
  { this[0] = simplex() ; 
    dim = dd ; 
    v = vpvector(dim+2) ; // last element is centroid, not a vertex 
    n = nn ; 
    c = cc ; 
    err = vector(n) ; 
    hypervol = vector(n) ; 
  }
  simplex(simplex *s) 
  { int t ; 
    this[0] = s[0] ; 
    v = vpvector(dim+2) ; 
    for(t=0;t<dim+2;t++) v[t] = s->v[t] ; 
    err = vector(n) ; 
    for(t=0;t<n;t++) err[t] = s->err[t] ; 
    hypervol = vector(n) ; 
    for(t=0;t<n;t++) hypervol[t] = s->hypervol[t] ; 
  }
  void release() { free(err,hypervol,v) ; } 
  void setcentroid()
  { int i,t ; 
    double *x=vector(dim) ; 
    for(i=0;i<dim;i++) for(t=0;t<=dim;t++) x[i] += v[t]->x[i] / (dim+1) ; 
    v[dim+1] = c->addvertex(dim,n,x) ;
    free(x) ; 
  }
  void integrate(double *asum)
  { int t,i ; 
    double q,qm,w=(dim+1)*(dim+1),qq=vol/((dim+1)*(dim+2)) ; 
    preverr[0] = preverr[1] / 2 ;
    preverr[1] = maxerr / 2 ; 
    for(maxerr=i=0;i<n;i++) 
    { for(qm=t=0;t<=dim;t++) qm += v[t]->y[i] ; 
      hypervol[i] = (w*v[dim+1]->y[i]+qm) * qq ;
      err[i] = fabs(hypervol[i]-qm*qq*(dim+2)) ;
      if(asum[i]&&err[i]/asum[i]>maxerr) maxerr = err[i] / asum[i] ;
    }
    c->addleveval(lev) ; 
  }
  adj *findedge(int v0,int v1) 
  { vertex *vtx0=v[v0],*vtx1=v[v1] ; 
    int t ; 
    if(vtx0>vtx1) swap(vtx0,vtx1) ; 
    for(t=0;t<vtx0->na&&vtx0->neighbour[t].v!=vtx1;t++) ; 
    if(t==vtx0->na) throw up("missing neighbour") ; 
    return vtx0->neighbour + t ; 
  }
  void increment(double *sum,double *asum,double *aerr,int decr)
  { double q ; 
    for(int t=0;t<n;t++)
    { sum[t] += decr * ( q = hypervol[t] ) ; 
      asum[t] += decr * fabs(q) ; 
      aerr[t] += decr * err[t] ; 
    }
  }
  /* --------------------------------- extend ------------------------------- */

  simplex *extend(double *asum)
  { int t,v0,v1,m0,m1,mf,xf ; 
    double qmax,*x=vector(dim) ; 
    simplex *s ; 
    vertex *vtxm,*vtx0,*vtx1 ; 
    adj *edge,*me ;

    for(qmax=mf=v0=0;v0<dim;v0++) for(v1=v0+1;v1<dim+1;v1++)
    { edge = findedge(v0,v1) ; 
      xf = (edge->midpt!=0) ;
      if(xf>mf||(xf==mf&&edge->d2>qmax))
      { m0 = v0 ; m1 = v1 ; me = edge ; mf = xf ; qmax = edge->d2 ; }
    }
    if(qmax<1e-20) { c->thinslices += 1 ; maxerr = 0 ; return 0 ; }

    // extend the linked list of vertices
    vtx0 = v[v0=m0] ; 
    vtx1 = v[v1=m1] ;
    for(t=0;t<dim;t++) x[t] = ( vtx0->x[t] + vtx1->x[t] ) / 2 ; 
    if(me->midpt==0) me->midpt = c->addvertex(dim,n,x) ; 
    free(x) ; 
    vtxm = me->midpt ; 
    me->count += 1 ; 
    for(t=0;t<=dim;t++) addneighbour(v[t],vtxm) ; 
    delneighbour(vtx0,vtx1) ; 

    // create a new simplex
    vol /= 2 ; 
    lev += 1 ; 
    s = new simplex(this) ; 
    s->v[v1] = v[v0] = vtxm ;
    s->setcentroid() ; 
    setcentroid() ; 
    s->integrate(asum) ; 
    integrate(asum) ; 
    return s ; 
  }
} ;
/* ---------------------------------- fanout -------------------------------- */

void polytope::fanout(vertex **corners)
{ int i , j , t , tau , *ind=ivector(dim) , co ;
  polytope *p ; 
  simplex *s ; 

  for(t=0;t<dim;t++) ind[t] = t ; 
  for(t=0;t<nfixed;t++) ind[fixed[t]] = -1 ; 
  for(tau=t=0;t<dim;t++) if(ind[t]>=0) ind[tau++] = ind[t] ;
  descendents = vvvector(ndescendents=tau) ;
 
  if(nfixed<dim-2) for(t=0;t<tau;t++)
  { descendents[t] = (void *) ( p = new polytope(dim,n,c,nfixed+1) ) ; 
    for(i=0;i<nfixed;i++) p->fixed[i] = fixed[i] ; 
    p->fixed[nfixed] = ind[t] ; 
    p->fanout(corners) ; 
  }
  else 
  { if(dim>1) for(i=0;i<2;i++) ind[dim-1-i] = ind[1-i] ;
    for(i=0;i<dim-2;i++) ind[i] = fixed[i] ; 
    for(t=0;t<tau;t++) 
    { c->slist[c->ns++] = s = new simplex(dim,n,c)  ; 
      descendents[t] = (void *) s ; 
      for(i=0;i<=dim;i++) 
      { for(co=j=0;j<i;j++) co |= 1 << ind[j] ; s->v[i] = corners[co] ; }
      s->setcentroid() ; 
      for(i=0;i<dim;i++) for(j=i+1;j<=dim;j++) addneighbour(s->v[i],s->v[j]) ;
      if(dim>1) swap(ind[dim-1],ind[dim-2]) ;
    }
  }
  free(ind) ; 
}
/* ------------------------------ printneighbours --------------------------- */

static void printneighbours(vertex *v,control *c,int lev)
{ int i,j,t,n=c->nv ; 
  vertex *vv,*vmax,**vtx=vpvector(n) ; 
  for(vv=v;vv;vv=vv->next) 
  { j = vv->serialno ;
    if(j<0||j>=n) throw up("illegal serial number %d(/%d)",j,n) ; 
    if(vtx[j]) throw up("serial number %d used twice",j) ; 
    vtx[j] = vv ;
  }
  printf("=== adjacency:\n ") ; 
  for(i=0;i<n;i++) if(i%5==0) printf("%4d ",i%10) ;
  printf("\n") ; 
  for(i=0;i<n;i++) 
  { if(i%5==0) printf("%3d ",i) ; else printf("    ") ;
    for(j=0;j<n;j++) 
    { if(j<=i) { if(j==i) printf("=") ; else printf(" ") ; continue ; }
      if(vtx[i]<vtx[j]) { vv = vtx[i] ; vmax = vtx[j] ; } 
      else { vv = vtx[j] ; vmax = vtx[i] ; } 
      for(t=0;t<vv->na&&vv->neighbour[t].v!=vmax;t++) ;
      if(t==vv->na) printf(".") ; else printf("%d",vv->neighbour[t].count) ; 
    }
    printf("\n") ; 
    if(i==lev) break ; 
  }
  free(vtx) ; 
}
/* -------------------------------------------------------------------------- */

static void deleter(void *v)
{ simplex *s=(simplex *) v ; s->release() ; delete s ; }

static quadratereport stub(void *f,int dim,int n,double *res,
                           double qtol,int maxeval,int mode)
{ int i,t,sno,*allzero=ivector(n) ;
  double q,*x,ftol=fabs(qtol),*err=vector(n),*sum=vector(n),*asum=vector(n) ;
  simplex *s , *ss , *sptr ; 
  vertex **corners,*v ; 
  quadratereport report ; 
  control c(f,mode,dim,n,maxeval) ; 
  polytope p(dim,n,&c,0) ; 
  heap h ;

  // set up a vertex at each corner of the hypercube
  corners = vpvector(1<<dim) ; 
  x = vector(dim) ; 
  for(i=0;i<(1<<dim);i++)
  { for(t=0;t<dim;t++) x[t] = 1 & (i>>t) ; corners[i] = c.addvertex(dim,n,x) ; }

  // partition the hypercube into simplices
  p.fanout(corners) ; 
  free(corners,x) ; 
  p.release() ; 
  // printneighbours(c.vlist,&c,-1) ; 

  for(q=1.0/factorial(dim),sno=0;sno<c.ns;sno++)
  { s = c.slist[sno] ; ; s->vol = q ; s->integrate(asum) ; }

  // find where the function evaluates to zero everywhere 
  for(i=0;i<n;i++) for(allzero[i]=1,sno=0;sno<c.ns;sno++)
  { for(s=c.slist[sno],t=0;t<=dim;t++) if(s->v[t]->y[i]) allzero[i] = 0 ; 
    if(s->hypervol[i]) allzero[i] = 0 ; 
  }

  // compute sum, asum and err
  for(sno=0;sno<c.ns;sno++) s->increment(sum,asum,err,1) ; 
  if(qtol<0) for(i=0;i<n;i++) asum[i] = 1 ; // i.e. if tolerance is absolute

  // compute maxerrs and put the vertices on a heap
  for(sno=0;sno<c.ns;sno++) 
  { for(s=c.slist[sno],q=i=0;i<n;i++) 
      if(!allzero[i]&&asum[i]&&s->err[i]/asum[i]>q) q = s->err[i] / asum[i] ;
    s->maxerr = q ; 
    h.push(s,q) ; 
  }
  c.freeslist() ; 

  // main loop
  while(c.neval<c.maxeval) 
  { for(t=0;t<n&&(allzero[t]||err[t]<ftol*asum[t]);t++) ;
    if(t==n) break ; 
    s = (simplex *) h.pop() ; 
    s->increment(sum,asum,err,-1) ;
    ss = s->extend(asum) ; 
    for(sptr=s;sptr;sptr=(sptr==s?ss:0))
    { h.push(sptr,sptr->maxerr) ; sptr->increment(sum,asum,err,1) ; }
  }

  // wrap up
  for(q=t=0;t<n;t++) 
  { res[t] = sum[t] ; if(!allzero[t]) q = max(q,err[t]/asum[t]) ; }
  report.condition = q/ftol ;
  report.neval = c.neval ; 
  report.thinslices = c.thinslices ;
  while(c.nleveval>0&&c.leveval[c.nleveval-1]==0) c.nleveval -= 1 ; 
  report.leveval = c.leveval ; 
  report.nleveval = c.nleveval ; 
  free(allzero,err,sum,asum) ; 
  h.recurse(deleter) ; 
  h.release() ; 
  if(mode<2&&report.neval>=maxeval) 
    throw up("too many function evals performed (%d)",report.neval) ; 
  return report ; 
}
/* -------------------------------------------------------------------------- */

double cubate(double (f)(double*),int dim,double qtol,int maxeval)
{ double res ;
  stub((void *)f,dim,1,&res,qtol,maxeval,0) ; 
  return res ; 
}
double cubate(double (f)(double*,int),int dim,double qtol,int maxeval)
{ double res ;
  stub((void *)f,dim,1,&res,qtol,maxeval,1) ; 
  return res ; 
}
quadratereport cubate(void (f)(double*,double *),int dim,int n,double *res,
            double qtol,int maxeval)
{ return stub((void *)f,dim,n,res,qtol,maxeval,2) ; }

quadratereport cubate(void (f)(double*,double *,int,int),int dim,int n,
                      double *res,double qtol,int maxeval)
{ return stub((void *)f,dim,n,res,qtol,maxeval,3) ; }

/* -------------------------------------------------------------------------- */

heap.h manages a binary max heap of unlimited size. It should be self-explanatory.

• heapitem     • push     • print     • recurse     • pop     • release     • heap     • size    

struct heapitem
{ void *data ; double key ; heapitem *l,*r ; int nl,nr ; 
  heapitem() {} ; 
  heapitem(void *d,double k)
  { data = d ; key = k ; l = r = 0 ; nl = nr = 0 ; }
} ;

struct heap
{ 
private : 
  heapitem *h ; 
  void push(void *d,double k,heapitem **loc)
  { void *dd ; 
    double kk ; 
    if(loc[0]==0) loc[0] = new heapitem(d,k) ; 
    else if(k>=loc[0]->key)
    { dd = loc[0]->data ;
      kk = loc[0]->key ;
      loc[0]->data = d ;
      loc[0]->key  = k ;
      if(loc[0]->nl>loc[0]->nr) { loc[0]->nr += 1 ; push(dd,kk,&loc[0]->r) ; }
      else { loc[0]->nl += 1 ; push(dd,kk,&loc[0]->l) ; }
    }
    else if(loc[0]->nl==0) { loc[0]->l = new heapitem(d,k) ; loc[0]->nl = 1 ; }
    else if(loc[0]->nr==0) { loc[0]->r = new heapitem(d,k) ; loc[0]->nr = 1 ; }
    else if(loc[0]->nl>loc[0]->nr) { loc[0]->nr += 1 ; push(d,k,&loc[0]->r) ; }
    else { loc[0]->nl += 1 ; push(d,k,&loc[0]->l) ; }
  }
  void print(heapitem *x)
  { if(x==0) return ;
    printf("%d",(int)x->key) ;
    if(x->l||x->r)
    { printf("{") ; if(x->l) print(x->l) ; else printf(".") ; printf("}") ; 
      printf("{") ; if(x->r) print(x->r) ; else printf(".") ; printf("}") ; 
    }
  }
  void recurse(heapitem *x,void(f)(void *))
  { if(x==0) return ;
    if(x->l) recurse(x->l,f) ;
    if(x->r) recurse(x->r,f) ;
    f(x->data) ; 
  }
  void pop(heapitem **loc)
  { if(loc==0) ; 
    else if(loc[0]->nl==0&&loc[0]->nr==0) { delete loc[0] ; loc[0] = 0 ; }
    else if(loc[0]->nl==0||(loc[0]->nr>0&&loc[0]->r->key>=loc[0]->l->key))
    { loc[0]->data = loc[0]->r->data ;
      loc[0]->key  = loc[0]->r->key ;
      loc[0]->nr -= 1 ; 
      pop(&loc[0]->r) ; 
    }
    else 
    { loc[0]->data = loc[0]->l->data ;
      loc[0]->key  = loc[0]->l->key ;
      loc[0]->nl -= 1 ; 
      pop(&loc[0]->l) ; 
    }
  }
  void release(heapitem **loc)
  { if(loc[0]->nl) release(&loc[0]->l) ; 
    if(loc[0]->nr) release(&loc[0]->r) ; 
    delete(loc[0]) ; 
    loc[0] = 0 ; 
  }
public : 
  heap() { h = 0 ; }
  void push(void *d,double k) { push(d,k,&h) ; }
  void print() { if(h) print(h) ; }
  void *pop() 
  { void *r ; 
    if(h==0) return 0 ; else { r = h[0].data ; pop(&h) ; return r ; }
  }
  void recurse(void(f)(void *)) { recurse(h,f) ; }
  int size() { if(h==0) return 0 ; else return 1+h->nl+h->nr ; }
  void release() { release(&h) ; }
} ;

#include "memory.h"
#include <math.h>
#include "quadrate.h"

static void f(double *x,double *y,int dim,int n)
{ int t ; for(y[0]=t=0;t<dim;t++) y[0] += x[t]*x[t]*x[t] ; }

int main()
{ double res ; 
  quadratereport rep ; 
  for(int t=1;t<11;t++) 
  { rep = cubate(f,t,1,&res,1,pow(10,t+1)) ; 
    printf("%d: integral = %.10f; neval=%-6d, cond=%.1f\n",
           t,res,rep.neval,rep.condition) ; 
  }
}

This tests the cubate software rather than validating the algorithm. The tolerance is set meaninglessly high, and the errors in the results are suspiciously low – perhaps something is wrong. I get the following output:

1: integral = 0.2500000000; neval=6     , cond=0.2
2: integral = 0.5000000000; neval=9     , cond=0.8
3: integral = 0.7500000000; neval=17    , cond=0.9
4: integral = 1.0000000000; neval=40    , cond=1.0
5: integral = 1.2500000000; neval=155   , cond=1.0
6: integral = 1.5000000000; neval=787   , cond=1.0
7: integral = 1.7500000000; neval=5171  , cond=1.0
8: integral = 2.0000000000; neval=40579 , cond=1.0
9: integral = 2.2500000000; neval=363392, cond=1.0
10: integral = 2.4999999999; neval=3629824, cond=1.0

#include "memory.h"
#include "random.h"
#include "heap.h"

void printint(void *x) { printf("%d ",((int *)x)[0]) ; }

int main()
{ int t,n,*d ;
  heap h ; 
  for(t=0;t<10;t++) 
  { d = new int ; 
    d[0] = n = rani(100) ; 
    printf("+%2d: ",n) ; 
    h.push(d,n) ; 
    h.print() ; 
    printf("\n") ; 
    if(t%3!=2) continue ;
    printf("-%2d: ",((int *)h.pop())[0]) ;
    h.print() ; 
    printf("\n") ; 
  }
  h.recurse(printint) ; 
  printf("\n") ; 
  h.release() ;
  h.print() ; 
}

I get the following output:

+ 1: 1
+ 4: 4{1}{.}
+74: 74{1}{4}
-74: 4{1}{.}
+15: 15{1}{4}
+69: 69{15{1}{.}}{4}
+53: 69{15{1}{.}}{53{4}{.}}
-69: 53{15{1}{.}}{4}
+ 4: 53{15{1}{.}}{4{4}{.}}
+ 6: 53{15{1}{6}}{4{4}{.}}
+92: 92{15{1}{6}}{53{4}{4}}
-92: 53{15{1}{6}}{4{4}{.}}
+82: 82{15{1}{6}}{53{4}{4}}
1 6 15 4 4 53 82 

The following program computes the average distance of a voter from the candidate she ranks 1st, 2nd, etc. in an election between m candidates. Voters and candidates are drawn from a standard bivariate Gaussian. We can normalise the average distances to use as weights (which might be termed Euclidean weights) for a Borda-like voting method.

If we run borda 10 (the argument being the number of candidates) then we get the output

dim=1, cond=1.0; neval=9403
by level: 2 4 8 12 12 22 42 74 138 242 426 728 1148 1524 1638 946 314 102
 0.2152  0.4045  0.5841  0.7616  0.9428  1.1340  1.3441  1.5884  1.9018  2.4074 | avge=1.1284 (should be 2/sqrt(pi)=1.1284)
9 8.22 7.49 6.76 6.01 5.23 4.37 3.36 2.08 0 
dim=2, cond=1.0; neval=10038
by level: 2 4 8 14 18 22 42 78 148 260 474 814 1304 1898 2078 656 78 4
 0.6530  0.9649  1.2055  1.4184  1.6214  1.8259  2.0434  2.2906  2.6022  3.0992 | avge=1.7725 (should be sqrt(pi)=1.7725)
9 7.85 6.97 6.18 5.44 4.68 3.88 2.97 1.83 0 

The output gives sets of weights for univariate then bivariate voter distributions. The first two lines tell us that 9403 function evaluations were performed for a univariate distribution, and show how they break down by recursive level. The overall average distance is computed as 1.1284 equal to its theoretical value of 2/√π.

The next line says that a voter’s first choice candidate is on average at a distance of 0.2152 from her; her second choice is at 0.4045; and so forth. If we wanted to transform these weights for a Borda count in which the highest score was best, we could assign 9 points to a voter’s top candidate, 8.22 to her second, and so forth.

Then we get similar information for bivariate gaussians.

• gausscirc     • f     • main

#include "memory.h"
#include <math.h>
#include "quadrate.h"
double besseli(int n,double x),rtlnorm(double x),ltlnorm(double x) ;
double gaussdisc(double dd,double rr),gaussexdisc(double dd,double rr) ;
static double gausscirc(double d,double r) 
{ return r * exp(-(d*d+r*r)/2) * besseli(0,r*d) ; }

static int m,dim ;

// For a voter at distance d from the origin, and given a radius r from the 
// voter, multiply together
// i.   The probability of a voter V having that position;
// ii.  The probability of a candidate C then being at distance r from V;  
// iii. For each position t in V's ballot, the probability that V will place t  
//      candidates higher than C and m-t-1 lower; 
// iv.  And r and m themselves.

static void f(double *x,double *y)
{ double u=x[0],v=x[1],d,r,p,q,pq,bico,bicon,bicod,qk,qcirc,qd ; 
  static double q2pi = 1/sqrt(2*3.141592653589793) ;
  int t ; 
  if(v<=0||v>=1||u>=1) { for(t=0;t<m;t++) y[t] = 0 ; return ; }
  d = 10 * atanh(u) ; 
  r = 10 * atanh(v) ;
  if(dim==2)
  { qd = d * exp(-d*d/2) ;                                           // i
    qcirc = gausscirc(d,r) ;                                         // ii
    if(r>1+1/(2*d+1)) p = 1 - ( q = gaussexdisc(d,r) ) ;
    else q = 1 - ( p = gaussdisc(d,r) ) ;
  }
  else 
  { qd = 2 * q2pi * exp(-d*d/2) ;                                    // i
    qcirc = q2pi * ( exp(-(d-r)*(d-r)/2) + exp(-(d+r)*(d+r)/2) ) ;   // ii
    if(r>1+1/(2*d+1)) p = 1 - ( q = rtlnorm(d+r) + ltlnorm(d-r) ) ;
    else q = 1 - ( p = rtlnorm(d-r) - rtlnorm(d+r) ) ;
  }
  pq = p / q ;
  qk = (r*m) * qd * qcirc * 100 / ( (1-u*u) * (1-v*v) ) ;
  y[0] = qk * pow(q,m-1) ;
  for(t=1;t<m;t++) y[t] = y[t-1] * pq ; 
  for(bicon=m-1,bicod=1,bico=t=1;t<m;t++,bicon--,bicod++)
  { bico *= bicon/bicod ; y[t] *= bico ; }
}
int main(int argc,char **argv)
{ int t,i ;
  if(argc>=2) m = atoi(argv[1]) ; else m = 9 ; 
  double q,*y=vector(m),tol,pi=3.141592653589793 ;
  if(argc>=3) tol = pow(0.1,atof(argv[2])) ; else tol = 0.01 ; 

  for(dim=1;dim<=2;dim++)
  { quadratereport r = cubate(f,2,m,y,tol,100000) ;
    printf("dim=%d, cond=%.1f; neval=%d\nby level:",dim,r.condition,r.neval) ;
    for(i=0;i<r.nleveval;i++) printf(" %d",r.leveval[i]) ; 
    printf("\n") ; 

    for(q=t=0;t<m;t++) { q += y[t] ; printf("%7.4f ",y[t]) ; }
    if(dim==2) printf("| avge=%.4f (should be sqrt(pi)=%.4f)\n",q/m,sqrt(pi)) ; 
    else printf("| avge=%.4f (should be 2/sqrt(pi)=%.4f)\n",q/m,2/sqrt(pi)) ; 
    q = (m-1) / (y[m-1]-y[0]) ; 
    for(t=0;t<m;t++) y[t] = q * ( y[m-1] - y[t] ) ; 
    printf("%d ",m-1) ; 
    for(t=1;t<m-1;t++) printf("%.2f ",y[t]) ; 
    printf("0 \n") ; 
  }
}

cubate.c heap.h
cubatetest.c heaptest.c
condorcet: memory.h ranf.c random.h
gaussint: quadrate.h

Compile and run the tests using g++ as compiler:

g++ -o cubatetest -O -g cubatetest.c cubate.c ; cubatetest
g++ -o heaptest heaptest.c ranf.c ; heaptest

borda.c
gaussint: quadrate.c gaussint.c rtlnorm.c log1plus.c
spatial solver: besseli.c

Compile borda by typing

g++ -O -g -o borda borda.c cubate.c gaussint.c rtlnorm.c log1plus.c quadrate.c besseli.c