Colin Champion
This page contains the utility functions I have written (or adapted) for use with voting software.
| function | description | |
| function | under | description |
functions provided : doc : source : test
memory.h is my standard header, used in nearly all my code.
| free(x,...) | frees up to 8 memory locations in one call. Saves typing. |
| x = vector(n) | allocates an n-long vector of doubles (which may be freed using free). |
| x = vector(x,n) | reallocates a vector x of doubles to length n. |
| x = matrix(m,n) | allocates an m×n two-dimensional matrix of doubles. |
| x = matrix(m,n,l) | allocates an m×n×l three-dimensional matrix of doubles. |
| freematrix(x) | frees a matrix thus generated. |
| freematrix(x,...) | frees up to 6 two-dimensional matrices in one call. Saves typing. |
| x = ivector(n) | allocates an n-long vector of ints (which may be freed using free). |
| x = ivector(x,n) | reallocates a vector x of ints to length n. |
| x = imatrix(m,n) | allocates an m×n two-dimensional matrix of ints. |
| allocates an m×n×l three-dimensional matrix of ints. | |
| freeimatrix(x) | frees a matrix thus generated. |
| x = cvector(n) | allocates an n-long vector of complexes (which may be freed using free). |
| x = cvector(x,n) | reallocates a vector x of complexes to length n. |
| x = cmatrix(m,n) | allocates an m×n two-dimensional matrix of complexes. |
| allocates an m×n×l three-dimensional matrix of complexes. | |
| freecmatrix(x) | frees a matrix thus generated. |
| x = charvector(n) | allocates an n-long vector of chars (which may be freed using free). |
| allocates a strlen(s)+1-long vector of chars, initialising it with s. | |
| x = charvector(x,n) | reallocates a vector x of chars to length n. |
| x = strvector(n) | allocates an n-long vector of char *s (which may be freed using free). |
| x = strvector(x,n) | reallocates a vector x of char *s to length n. |
| x = shortvector(n) | allocates an n-long vector of shorts (which may be freed using free). |
| x = shortvector(x,n) | reallocates a vector x of shorts to length n. |
| swap(x,y) | swaps x and y, which may be any mixture of ints and doubles, or else both be chars. |
| fopenread(char *s) | calls fopen(s,"r"), returning the file descriptor but erroring out informatively if the call fails. If s is the string ‘--’ then fopenread returns stdin. |
| fopenwrite(char *s) | the same for writing instead of reading, and with stdout as the default. |
| freadline(FILE *ifl) | returns a line read from ifl (like fgets). |
| readline() | returns a line read from stdin (like gets). |
memory.h provides a number of convenience functions to the user.
Error handling: it provides a notation for erroring out. In case of error you throw the uncaught exception ‘up’, whose arguments are precisely those of printf. The result will be to terminate with the error message specified by the arguments (but also identifying the location in your program at which the error occurred). Eg.
if(n<0) throw up("Your size %d is negative.",n) ;
There is no need for a terminating new-line: one is added automatically.
I have experimented with several alternative methods. There is an underlying problem which is that C doesn’t really have a semantics for errors, just for terminating with different values. But an uncaught exception has to be considered an error.
Memory allocation: the functions listed above allow you to allocate and deallocate vectors and matrices of various data types. Inevitably the list of options provided is a little arbitrary. In case of failure the program will error out (saving you from having to test against NULL every time). The location reported for the error will be in the function you called (eg. ivector) rather than being assigned to a location in your own program. The latter would be possible if I could use preprocessor substitutions, but these are incompatible with overloading and nearly all my functions are overloaded.
The memory allocation is performed by calling one of cjcalloc and cjcrealloc which are similar to alloc and realloc except for erroring out in case of failure. If you need to allocate memory in a way which isn’t handled by my vector and matrix functions you may wish to call them yourself. The prototypes are
void *cjcalloc(int nitems,int itemsize) ; void *cjcrealloc(void *ptr,int nbytes) ;
cjcrealloc differs from realloc in that a = cjcrealloc(a,0) returns NULL if a is supplied as NULL.
Matrices (and submatrices) are allocated as consecutive regions of memory, allowing them to be accessed by single or double indirection, or by triple indirection in the case of 3 dimensions.
Sorting: there are generic ascending and descending shell sorts, intended for data types of the user’s own defining (eg. ordered pairs) and an ascending integer sort
isortup(int *x,int n) ;
See ij.h for a use of the generic sort.
max/min: diadic operators are defined as preprocessor substitutions.
File opening: functions are provided which error out in case of failure.
Reading line by line: since C made a complete hash of this I provide readline and freadline. They return strings which need to be freed when you’ve finished with them. An empty line is returned as an empty string; a null return corresponds to end of file.
Standard header inclusion: memory.h includes all 4 of stdio.h, stdlib.h, string.h and stdarg.h, so logically it precedes the rest of your #include files.
• max • min • cjcup • cjcuplog • cjcformat • cjcprint2 • cjcprint1 • complex • print • xy • print • xi • ij • settable • unset • double • print • free • cjcupalloc • cjcuprealloc • swap • freename • charvector • xivector • isortup • realsort • realsortdown • xysort • xisort • ijsort • xisortdown • fupopenread • fupopenwrite • freadline • readline
#ifndef MEMORY_H #define MEMORY_H #include <stdio.h> #include <stdlib.h> #include <string.h> #include <stdarg.h> static double max(double a,double b) { if(a>b) return a ; else return b ; } static double min(double a,double b) { if(a<b) return a ; else return b ; } static double max(double a,int b) { if(a>b) return a ; else return b ; } static double min(double a,int b) { if(a<b) return a ; else return b ; } static double max(int a,double b) { if(a>b) return a ; else return b ; } static double min(int a,double b) { if(a<b) return a ; else return b ; } static int max(int a,int b) { if(a>b) return a ; else return b ; } static int min(int a,int b) { if(a<b) return a ; else return b ; } #define abnormal(x) (!!((x)-(x))) // x-x forced to boolean /* ---------------------------- define throwing up -------------------------- */ static int cjcupline=-1,cjcperror=0 ; static const char *cjcupfile="",*cjcupfunc="" ; static int cjcup(const char *m,...) { va_list vl ; fprintf(stderr,"*** Error at line %d of %s [function %s]:\n", cjcupline,cjcupfile,cjcupfunc) ; if(cjcperror) perror(0) ; va_start(vl,m) ; vfprintf(stderr,m,vl) ; va_end(vl) ; fprintf(stderr,"\n") ; fflush(0) ; return 0 ; } ; static void cjcuplog(const char *f,int l,const char *ff) { cjcupfile = f ; cjcupline = l ; cjcupfunc = ff ; } #define up (cjcuplog(__FILE__,__LINE__,__PRETTY_FUNCTION__),cjcup) /* -------------------------- simple ordered pairs -------------------------- */ static char cjcbuf[600]={0} ; static int cjcind = 40 , cjcint = 0 ; static void cjcformat(char *fmt) { int i ; if(strlen(fmt)>18) throw up("overlong cjcprint format %s") ; strcpy(cjcbuf,fmt) ; for(i=0;fmt[i]&&fmt[i]!='%';i++) ; for(;fmt[i]&&(fmt[i]=='l'||(fmt[i]<'a'||fmt[i]>'z'));i++) ; cjcint = (fmt[i]=='d') ; } static char *cjcprint2(double x,double y) { if(cjcbuf[0]==0) cjcformat((char *)"%.1f") ; char *ptr = cjcbuf + cjcind ; cjcbuf[cjcind] = '(' ; if(cjcint) cjcind += 2 + snprintf(cjcbuf+cjcind+1,100,cjcbuf,(int)x) ; else cjcind += 2 + snprintf(cjcbuf+cjcind+1,100,cjcbuf,x) ; cjcbuf[cjcind-1] = ',' ; if(cjcint) cjcind += 2 + snprintf(cjcbuf+cjcind,100,cjcbuf,(int)y) ; else cjcind += 2 + snprintf(cjcbuf+cjcind,100,cjcbuf,y) ; cjcbuf[cjcind-2] = ')' ; cjcbuf[cjcind-1] = 0 ; if(cjcind>600) throw up("data overflow on printing an ordered pair") ; if(cjcind>500) cjcind = 40 ; return ptr ; } static char *cjcprint2(char *fmt,double x,double y) { cjcformat(fmt) ; return cjcprint2(x,y) ; } static char *cjcprint1(double x) { if(cjcbuf[0]==0) cjcformat((char *)"%.1f") ; char *ptr = cjcbuf + cjcind ; if(cjcint) cjcind += 1 + snprintf(cjcbuf+cjcind,100,cjcbuf,(int)x) ; else cjcind += 1 + snprintf(cjcbuf+cjcind,100,cjcbuf,x) ; if(cjcind>600) throw up("data overflow on printing") ; if(cjcind>500) cjcind = 40 ; return ptr ; } static char *cjcprint1(char *fmt,double x) { cjcformat(fmt) ; return cjcprint1(x) ; } struct complex { double re,im ; complex() { re = im = 0 ; } complex(double x) { re = x ; im = 0 ; } complex(double x,double y) { re = x ; im = y ; } complex &operator=(double x) { re = x ; im = 0 ; return *this ; } char *print(char *fmt) { return cjcprint2(fmt,re,im) ; } char *print() { return cjcprint2(re,im) ; } } ; struct xy { double x,y ; xy() { x = y = 0 ; } xy(double xx) { x = xx ; y = 0 ; } xy(double xx,double yy) { x = xx ; y = yy ; } xy(double xx,double(*f)(double)) { x = xx ; y = f(x) ; } xy &operator=(double xx) { x = xx ; y = 0 ; return *this ; } char *print(char *fmt) { return cjcprint2(fmt,x,y) ; } char *print() { return cjcprint2(x,y) ; } } ; inline bool operator==(xy a,xy b) { return a.x==b.x&&a.y==b.y ; } struct xi { double x ; int i ; xi() { x = i = 0 ; } xi(double a,int b) { x = a ; i = b ; } } ; struct ij { int i,j ; ij() { j = i = 0 ; } ij(int a,int b) { i = a ; j = b ; } } ; /* --------------------------------- settables ------------------------------ */ struct settable { private: double x ; public: bool set ; settable() { set = 0 ; } settable(double y) { if(abnormal(y)) throw up("setting a settable to %.1e",y) ; set = 1 ; x = y ; } settable unset() { return this[0] = settable() ; } operator double() { if(!set) throw up("accessing an unset settable") ; return x ; } settable &operator=(double y) { if(abnormal(y)) throw up("assigning %.1e to a settable",y) ; set = 1 ; x = y ; return *this ; } char *print() { if(set) return cjcprint1(x) ; char *ptr=cjcbuf+cjcind ; strcpy(ptr,"undef") ; cjcind += 6 ; if(cjcind>500) cjcind = 40 ; return ptr ; } char *print(char *fmt) { cjcformat(fmt) ; return print() ; } } ; inline bool operator==(settable x,settable y) { return (x.set==0&&y.set==0)||((double)x)==(double(y)) ; } inline bool operator==(settable x,double y) { return x.set && y==(double)x ; } inline bool operator==(double x,settable y) { return (y==x) ; } inline bool operator!=(settable x,settable y) { return !(x==y) ; } inline settable operator+=(settable &x,double y) { if(!x.set) throw up("incrementing an unset settable") ; return x = settable(x+y) ; } inline settable operator-=(settable &x,double y) { if(!x.set) throw up("decrementing an unset settable") ; return x = settable(x-y) ; } inline settable operator*=(settable &x,double y) { if(y==0) return x = settable(0) ; if(!x.set) throw up("*= applied to an unset settable") ; return x = settable(x*y) ; } inline settable operator/=(settable &x,double y) { if(y==0) throw up("settable divided by 0") ; if(!x.set) throw up("/= applied to an unset settable") ; return x = settable(x/y) ; } inline settable max(settable a,settable b) { if(a.set&&(!b.set||a>b)) return a ; else return b ; } inline settable min(settable a,settable b) { if(a.set&&(!b.set||a<b)) return a ; else return b ; } inline settable max(settable a,double b) { if(a.set&&a>b) return a ; else return b ; } inline settable min(settable a,double b) { if(a.set&&a<b) return a ; else return b ; } inline settable max(double a,settable b) { if(b.set&&b>a) return b ; else return a ; } inline settable min(double a,settable b) { if(b.set&&b<a) return b ; else return a ; } inline settable max(settable a,int b) { if(a.set&&a>b) return a ; else return b ; } inline settable min(settable a,int b) { if(a.set&&a<b) return a ; else return b ; } inline settable max(int a,settable b) { if(b.set&&b>a) return b ; else return a ; } inline settable min(int a,settable b) { if(b.set&&b<a) return b ; else return a ; } /* ------------------------------ variadic free() --------------------------- */ static void free(void *a,void *b) { free(a) ; free(b) ; } static void free(void *a,void *b,void *c) { free(a) ; free(b) ; free(c) ; } static void free(void *a,void *b,void *c,void *d) { free(a) ; free(b) ; free(c) ; free(d) ; } static void free(void *a,void *b,void *c,void *d,void *e) { free(a) ; free(b) ; free(c) ; free(d) ; free(e) ; } static void free(void *a,void *b,void *c,void *d,void *e,void *f) { free(a) ; free(b) ; free(c) ; free(d) ; free(e) ; free(f) ; } static void free(void *a,void *b,void *c,void *d,void *e,void *f,void *g) { free(a) ; free(b) ; free(c) ; free(d) ; free(e) ; free(f) ; free(g) ; } static void free(void *a,void *b,void *c,void *d, void *e,void *f,void *g,void *h) { free(a) ; free(b) ; free(c) ; free(d) ; free(e) ; free(f) ; free(g) ; free(h) ; } /* ------------------------- robust allocs ---------------------------- */ static void *cjcupalloc(int a,int b) { if(a<0||b<0) throw cjcup("negative length %d requested from cjcalloc.",a*b) ; void *p=calloc(a,b) ; if(p==0) throw cjcup("cjcalloc unable to allocate %d bytes of memory.",b) ; memset(p,0,a*b) ; return p ; } static void *cjcuprealloc(void *a,int b) { if(b<0) throw cjcup("negative length %d requested from cjcrealloc.",b) ; if(a==0&&b==0) return 0 ; void *p=realloc(a,b) ; if(b>0&&p==0) throw cjcup("cjcrealloc unable to reallocate %x to %d bytes.",a,b) ; return p ; } #define cjcalloc (cjcuplog(__FILE__,__LINE__,__PRETTY_FUNCTION__),cjcupalloc) #define cjcrealloc (cjcuplog(__FILE__,__LINE__,__PRETTY_FUNCTION__),cjcuprealloc) /* ---------------------------- generic matrix ------------------------------ */ #define genvector(type,vecname) \ static type *vecname(int n) \ { return (type*) cjcalloc(n,sizeof(type)) ; } \ static type *vecname(type *x,int n) \ { return (type *) cjcrealloc(x,n*sizeof(type)) ; } \ static void swap(type &a,type &b) { type c=b ; b = a ; a = c ; } #define genmatrix(type,vecname,name,freename) \ static type *vecname(int n) \ { return (type*) cjcalloc(n,sizeof(type)) ; } \ static type *vecname(type *x,int n) \ { return (type *) cjcrealloc(x,n*sizeof(type)) ; } \ static type **name(int m,int n) \ { type **a = (type **) cjcalloc(m,sizeof(type *)) ; \ a[0] = vecname(m*n) ; \ for(int i=1;i<m;i++) a[i] = a[i-1] + n ; \ return a ; \ } \ static type ***name(int m,int n,int l) \ { int i ; \ type ***a = (type ***) cjcalloc(m,sizeof(type **)) ; \ a[0] = (type **) cjcalloc(m*n,sizeof(type *)) ; \ for(i=1;i<m;i++) a[i] = a[i-1] + n ; \ a[0][0] = vecname(m*n*l) ; \ for(i=1;i<m*n;i++) a[0][i] = a[0][i-1] + l ; \ return a ; \ } \ static void freename(type **a) { if(a) free(a[0],a) ; } \ static void freename(type **a,type **b) \ { freename(a) ; freename(b) ; } \ static void freename(type **a,type **b,type **c) \ { freename(a) ; freename(b) ; freename(c) ; } \ static void freename(type **a,type **b,type **c,type **d) \ { freename(a) ; freename(b) ; freename(c) ; freename(d) ; } \ static void freename(type **a,type **b,type **c,type **d, \ type **e) \ { freename(a) ; freename(b) ; freename(c) ; freename(d) ; \ freename(e) ; } \ static void freename(type **a,type **b,type **c,type **d, \ type **e,type **f) \ { freename(a) ; freename(b) ; freename(c) ; freename(d) ; \ freename(e) ; freename(f) ; } \ static void freename(type ***a) { if(a) free(a[0][0],a[0],a) ; } \ static void swap(type &a,type &b) { type c=b ; b = a ; a = c ; } /* --------------------------- matrix instances ----------------------------- */ genmatrix(double,vector,matrix,freematrix) ; genmatrix(int,ivector,imatrix,freeimatrix) ; genmatrix(xi,xivector,ximatrix,freeximatrix) ; genmatrix(ij,ijvector,ijmatrix,freeijmatrix) ; genmatrix(xy,xyvector,xymatrix,freexymatrix) ; genmatrix(char,charvector,charmatrix,freecharmatrix) ; genmatrix(char*,strvector,strmatrix,freestrmatrix) ; genmatrix(short,shortvector,shortmatrix,freeshortmatrix) ; genmatrix(complex,cvector,cmatrix,freecmatrix) ; genmatrix(settable,setvector,setmatrix,freesetmatrix) ; // a couple of specials static char *charvector(char *c) { char *ret=charvector(1+strlen(c)) ; strcpy(ret,c) ; return ret ; } static xi *xivector(double *x,int n) { xi *y=xivector(n) ; for(int i=0;i<n;i++) y[i] = xi(x[i],i) ; return y ; } /* ----------------------------- generic sorts ------------------------------ */ #define shellsortup(x,n,type,field) \ { int i,j,inc ; type y ; \ for(inc=1;1+3*inc<(n);inc=1+3*inc) ; \ for(;inc>0;inc/=3) for(i=inc;i<(n);i++) \ { for(y=x[i],j=i;j>=inc;j-=inc) \ { if(y.field<x[j-inc].field) x[j] = x[j-inc] ; else break ; } \ x[j] = y ; \ } \ } #define shellsortdown(x,n,type,field) \ { int i,j,inc ; type y ; \ for(inc=1;1+3*inc<(n);inc=1+3*inc) ; \ for(;inc>0;inc/=3) for(i=inc;i<(n);i++) \ { for(y=x[i],j=i;j>=inc;j-=inc) \ { if(y.field>x[j-inc].field) x[j] = x[j-inc] ; else break ; } \ x[j] = y ; \ } \ } #define shsortup(x,n,type) \ { int i,j,inc ; type y ; \ for(inc=1;1+3*inc<(n);inc=1+3*inc) ; \ for(;inc>0;inc/=3) for(i=inc;i<(n);i++) \ { for(y=x[i],j=i;j>=inc;j-=inc) \ { if(y<x[j-inc]) x[j] = x[j-inc] ; else break ; } \ x[j] = y ; \ } \ } #define shsortdown(x,n,type) \ { int i,j,inc ; type y ; \ for(inc=1;1+3*inc<(n);inc=1+3*inc) ; \ for(;inc>0;inc/=3) for(i=inc;i<(n);i++) \ { for(y=x[i],j=i;j>=inc;j-=inc) \ { if(y>x[j-inc]) x[j] = x[j-inc] ; else break ; } \ x[j] = y ; \ } \ } /* ----------------------------- sundry sorts ------------------------------- */ static void isortup(int *u,int n) { shsortup(u,n,int) ; } static void realsort(double *u,int n) { shsortup(u,n,double) ; } static void realsortdown(double *u,int n) { shsortdown(u,n,double) ; } static void xysort(xy *u,int n) { shellsortup(u,n,xy,x) ; } static void xisort(xi *u,int n) { shellsortup(u,n,xi,x) ; } static void ijsort(ij *u,int n) { shellsortup(u,n,ij,i) ; } static void xisortdown(xi *u,int n) { shellsortdown(u,n,xi,x) ; } /* -------------------------- define robust fopens -------------------------- */ static FILE *fupopenread(char *name) { FILE *f ; if(name[0]=='-'&&name[1]=='-'&&name[2]==0) f = stdin ; else f = fopen(name,"r") ; if(f==0) { cjcperror = 1 ; throw cjcup("Your input file %s could not be found.",name) ; } return f ; } static FILE *fupopenwrite(char *name) { FILE *f ; if(name[0]=='-'&&name[1]=='-'&&name[2]==0) f = stdout ; else f = fopen(name,"w") ; if(f==0) { cjcperror = 1 ; throw cjcup("Unable to write to your file %s.",name) ; } return f ; } #define fopenread (cjcuplog(__FILE__,__LINE__,__PRETTY_FUNCTION__),fupopenread) #define fopenwrite (cjcuplog(__FILE__,__LINE__,__PRETTY_FUNCTION__),fupopenwrite) static char *freadline(FILE *ifl) { char *s=0 ; int slen,ns,c,i ; for(slen=ns=0;;) { c = fgetc(ifl) ; if(c==EOF||c=='\n') { if(slen>ns+1||(ns==0&&c=='\n')) s = charvector(s,ns+1) ; return s ; } if(ns>=slen-1) { slen += 20 + slen/2 ; s = charvector(s,slen) ; for(i=ns;i<slen;i++) s[i] = 0 ; } s[ns++] = (char) c ; } } static char *readline() { return freadline(stdin) ; } #endif
The rather rudimentary test program memorytest.c verifies some simple operations.
#include "memory.h"
#include <math.h>
int main()
{ int i,j,k ;
{
double **a,***b ;
a = matrix(3,4) ;
for(i=0;i<3;i++) for(j=0;j<4;j++) a[i][j] = 10*i + j ;
for(i=0;i<12;i++) printf("%d ",(int) a[0][i]) ;
printf("\n") ;
freematrix(a) ;
b = matrix(2,3,4) ;
for(i=0;i<2;i++) for(j=0;j<3;j++) for(k=0;k<4;k++)
b[i][j][k] = 100*i + 10*j + k ;
for(i=0;i<24;i++) printf("%d ",(int) b[0][0][i]) ;
printf("\n") ;
freematrix(b) ;
}
{
int **a,***b ;
a = imatrix(3,4) ;
for(i=0;i<3;i++) for(j=0;j<4;j++) a[i][j] = 10*i + j ;
for(i=0;i<12;i++) printf("%d ",(int) a[0][i]) ;
printf("\n") ;
freeimatrix(a) ;
b = imatrix(2,3,4) ;
for(i=0;i<2;i++) for(j=0;j<3;j++) for(k=0;k<4;k++)
b[i][j][k] = 100*i + 10*j + k ;
for(i=0;i<24;i++) printf("%d ",(int) b[0][0][i]) ;
printf("\n") ;
freeimatrix(b) ;
}
{
complex **a,***b ;
a = cmatrix(3,4) ;
for(i=0;i<3;i++) for(j=0;j<4;j++) a[i][j] = 10*i + j ;
for(i=0;i<12;i++) printf("%d ",(int) a[0][i].re) ;
printf("\n") ;
freecmatrix(a) ;
b = cmatrix(2,3,4) ;
for(i=0;i<2;i++) for(j=0;j<3;j++) for(k=0;k<4;k++)
b[i][j][k] = 100*i + 10*j + k ;
for(i=0;i<24;i++) printf("%d ",(int) b[0][0][i].re) ;
printf("\n") ;
freecmatrix(b) ;
}
}
memory.h memorytest.c
g++ -g -O -o memorytest memorytest.c ; memorytest
The correct output is
0 1 2 3 10 11 12 13 20 21 22 23 0 1 2 3 10 11 12 13 20 21 22 23 100 101 102 103 110 111 112 113 120 121 122 123 0 1 2 3 10 11 12 13 20 21 22 23 0 1 2 3 10 11 12 13 20 21 22 23 100 101 102 103 110 111 112 113 120 121 122 123 0 1 2 3 10 11 12 13 20 21 22 23 0 1 2 3 10 11 12 13 20 21 22 23 100 101 102 103 110 111 112 113 120 121 122 123
functions provided : doc : source : test
linsolve solves one or more sets of linear equations axi = bi. The prototype is
double linsolve(double **a,double **b,double **x,int n,int nrhs) ;
where
a is left unchanged, and the solution x may overwrite b, otherwise b too is unchanged.
The return value is the arithmetic precision of the result. Zero implies rubbish; very small (eg. 1e-12) implies that roundoff error will have contaminated the result; close to 1 implies that the result is accurate.
#include <math.h>
#include "memory.h"
double linsolve(double **adash,double **bdash,double **x,int n,int m)
{ int i,j,k,imax,idash,jdash ;
double **a=matrix(n,n),**b=matrix(m,n),*cmax=vector(n),q,qerr ;
for(i=0;i<n;i++)
{ for(k=0;k<m;k++) b[k][i] = bdash[k][i] ;
for(j=0;j<n;j++)
{ q = a[i][j] = adash[i][j] ; cmax[j] = max(cmax[j],fabs(q)) ; }
}
for(j=0;j<n;j++)
{ for(q=0,imax=i=j;i<n;i++) if(fabs(a[i][j])>q)
{ imax = i ; q = fabs(a[i][j]) ; }
if(imax>j)
{ for(k=0;k<m;k++) swap(b[k][imax],b[k][j]) ;
for(jdash=j;jdash<n;jdash++) swap(a[imax][jdash],a[j][jdash]) ;
}
if(j==0||fabs(a[j][j])<cmax[j]*qerr) qerr = fabs(a[j][j]) / cmax[j] ;
for(idash=j+1;idash<n;idash++)
{ q = a[idash][j] / a[j][j] ;
for(k=0;k<m;k++) b[k][idash] -= q * b[k][j] ;
for(jdash=j;jdash<n;jdash++) a[idash][jdash] -= q * a[j][jdash] ;
}
}
if(qerr) for(j=n-1;j>=0;j--) for(k=0;k<m;k++)
{ x[k][j] = b[k][j] / a[j][j] ;
for(i=j-1;i>=0;i--) b[k][i] -= a[i][j] * x[k][j] ;
}
freematrix(a,b) ; free(cmax) ;
return qerr ;
}
lu performs LU decomposition on an n×n matrix. It corresponds to ludcmp in NR. I only use it in order to find a matrix’s determinant. The prototypes are
int lu(double **x,int n,int *ind) ; int lu(double **x,double **y,int n,int *ind) ; double det(double **x,int n) ;
where x is the input matrix; y is the output matrix; ind is an output array of row indices; and the return value of lu is the parity of the count of row interchanges. The call with separate matrices leaves x unchanged. Pass a null pointer in place of ind if you don’t want it. det() returns the determinant of a matrix (which is left unchanged).
Dependency: memory.h.
#include "memory.h" #include <math.h> int lu(double **x,int n,int *ind) { int i,imax,j,k,parity ; double q,qmax,*v=vector(n) ; for(parity=1,i=0;parity&&i<n;i++) { for(qmax=j=0;j<n;j++) qmax = max(qmax,fabs(x[i][j])) ; if(qmax==0) parity = 0 ; else v[i] = 1/qmax ; } for(j=0;parity&&j<n;j++) { for(i=0;i<j;i++) { for(q=x[i][j],k=0;k<i;k++) q -= x[i][k]*x[k][j] ; x[i][j] = q ; } for(qmax=imax=0,i=j;i<n;i++) { for(q=x[i][j],k=0;k<j;k++) q -= x[i][k]*x[k][j] ; x[i][j] = q ; if(v[i]*fabs(q)>=qmax) { qmax = v[i]*fabs(q) ; imax = i ; } } if(imax!=j) { for(k=0;k<n;k++) swap(x[imax][k],x[j][k]) ; v[imax] = v[j] ; parity = -parity ; } if(ind) ind[j] = imax ; if(x[j][j]==0) parity = 0 ; else for(i=j+1;i<n;i++) x[i][j] /= x[j][j] ; } free(v) ; return parity ; } int lu(double **x,double **y,int n,int *ind) { if(x!=y) for(int i=0;i<n;i++) for(int j=0;j<n;j++) y[i][j] = x[i][j] ; return lu(y,n,ind) ; } double det(double **x,int n) { double **y=matrix(n,n),q=lu(x,y,n,0) ; for(int i=0;i<n;i++) q *= y[i][i] ; freematrix(y) ; return q ; }
#include "memory.h"
double det(double **x,int n) ;
int main()
{ double a[3][3] = { { -2,-1,2} , {2,1,4} , {-3,3,-1} },**x=matrix(3,3) ;
int i,j,n=3 ;
for(i=0;i<n;i++) for(j=0;j<n;j++) x[i][j] = a[i][j] ;
printf("det=%.1f\n",det(x,n)) ;
}
memory.h lu.c lutest.c
Compile and execute as:
g++ -o lutest -g -O lu.c lutest.c ; lutest
The example comes from the Wikipedia article. The correct output is
det=54.0
gramschmidt orthonormalises an n×n matrix, i.e. replaces it by a matrix whose rows comprise an orthonormal basis. It works from the end, so that if the last few rows are already orthonormal, they will be left unchanged. The input rows must span the space. (You will usually fill the earlier ones with unit vectors.) By construction, the resulting matrix has the property that its transpose is its inverse. The prototype is
void gramschmidt(double **x,int n,) ;
where x is the matrix.
Dependency: memory.h.
#include "memory.h"
#include <math.h>
void gramschmidt(double **x,int n)
{ int i,j,ii ;
double q ;
for(i=n-1;i>=0;i--)
{ for(q=j=0;j<n;j++) q += x[i][j] * x[i][j] ;
if(q==0) throw up("matrix supplied to gramschmidt is singular") ;
for(q=1/sqrt(q),j=0;j<n;j++) x[i][j] *= q ;
for(ii=i-1;ii>=0;ii--)
{ for(q=j=0;j<n;j++) q += x[ii][j] * x[i][j] ;
for(j=0;j<n;j++) x[ii][j] -= q*x[i][j] ;
}
}
}
#include "memory.h"
void gramschmidt(double **x,int n) ;
int main()
{ int i,ii,j,n=5 ;
double q,**x=matrix(n,n) ;
for(i=0;i<n-1;i++) x[i][i+1] = i+1 ;
for(i=0;i<n;i++) x[n-1][i] = 1 ;
gramschmidt(x,n) ;
for(i=0;i<n;i++) { for(j=0;j<n;j++) printf("%6.3f",x[i][j]) ; printf("\n") ; }
for(i=0;i<n-1;i++) for(ii=i+1;ii<n;ii++)
{ for(q=j=0;j<n;j++) q += x[i][j] * x[ii][j] ; printf("%6.3f",q-(i==ii)) ; }
printf("\n") ;
}
memory.h gramschmidt.c gramschmidttest.c
Compile and execute as:
g++ -o gramschmidttest -g -O gramschmidt.c gramschmidttest.c ; gramschmidttest
The correct output is
-0.707 0.707-0.000 0.000 0.000 -0.408-0.408 0.816 0.000 0.000 -0.289-0.289-0.289 0.866 0.000 -0.224-0.224-0.224-0.224 0.894 0.447 0.447 0.447 0.447 0.447 -0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
in which the last line is a list of errors in row pairs.
doc : reference : source : test : download
The prototype is
double ***bezier(double **x,int nx) ;
where
Dependencies: memory.h , linsolve.c
#include "memory.h" double linsolve(double **a,double **b,double **x,int n,int nrhs) ; double ***bezier(double **x,int nx) { int t,k,n=2*(nx-2) ; if(nx<2) throw up("%d (<2) points supplied for bezier interpolation",nx) ; double **a=matrix(n,n),*b=vector(n),*y=vector(n+4),*y2=y+2 ; double ***bez=matrix(nx-1,4,2) ; // special handling of end points for(k=0;k<2;k++) { b[k] = 6*x[1][k] - x[0][k] ; b[n-2+k] = 6*x[nx-2][k] - x[nx-1][k] ; y[k] = x[0][k] ; y[2*(nx-1)+k] = x[nx-1][k] ; bez[0][0][k] = x[0][k] ; bez[nx-2][3][k] = x[nx-1][k] ; } // set up and solve equations for(t=0;t<n;t+=2) { if(t>0) a[t][t-2] = a[t+1][t-1] = 1 ; a[t][t] = a[t+1][t+1] = 4 ; if(t<n-2) a[t][t+2] = a[t+1][t+3] = 1 ; if(t>0&&t<n-2) { b[t] = 6*x[t/2+1][0] ; b[t+1] = 6*x[t/2+1][1] ; } } linsolve(a,&b,&y2,n,1) ; // y now holds the B-spline control points - compute the Bezier control points for(k=0;k<2;k++) for(t=0;t<nx-1;t++) { if(t>0) bez[t][0][k] = (y[2*(t-1)+k]+y[2*(t+1)+k])/6 + 2*y[2*t+k]/3 ; bez[t][1][k] = (2*y[2*t+k]+y[2*(t+1)+k]) / 3 ; bez[t][2][k] = (y[2*t+k]+2*y[2*(t+1)+k]) / 3 ; if(t<nx-2) bez[t][3][k] = (y[2*t+k]+y[2*(t+2)+k])/6 + 2*y[2*(t+1)+k]/3 ; } freematrix(a) ; free(b,y) ; return bez ; }
The first program is a simple workout.
#include "memory.h"
#include <math.h>
double ***bezier(double **x,int nx) ;
int main()
{ double **x=matrix(5,2),***b,q,h ;
int i,j ;
x[0][0] = 1 ; x[0][1] = -1 ;
x[1][0] = -1 ; x[1][1] = 2 ;
x[2][0] = 1 ; x[2][1] = 4 ;
x[3][0] = 4 ; x[3][1] = 3 ;
x[4][0] = 7 ; x[4][1] = 5 ;
b = bezier(x,5) ;
for(i=0;i<4;i++)
{ for(j=0;j<4;j++) printf("(%2d,%2d) ",(int) b[i][j][0],(int) b[i][j][1]) ;
printf("\n") ;
}
}
The second test program generates a plot illustrative of the spatial Borda count.
#include "memory.h"
#include <math.h>
double ***bezier(double **x,int nx) ;
static double w1[10]={9,8.22,7.49,6.76,6.01,5.23,4.37,3.36,2.08,0} ;
static double w2[10]={9,7.85,6.97,6.18,5.44,4.68,3.88,2.97,1.83,0} ;
static double w[10]={9,7.88,6.68,5.46,4.28,3.20,2.23,1.40,0.67,0} ;
static double ww[10]={9,7.65,6.38,5.27,4.28,3.39,2.58,1.81,0.96,0} ;
double xmap(double x) { return 20+600*x/9 ; }
double ymap(double y) { return 620-600*y/9 ; }
void drawX(double **X,int n,char *col)
{ int i,j ;
double **Y=matrix(n,2) ;
for(i=0;i<n;i++)
{ Y[i][0] = xmap(X[i][0]) ; Y[i][1] = ymap(X[i][1]) ; }
double ***b=bezier(Y,n) ;
printf("<path d=\"M %8.3f %8.3f ",b[0][0][0],b[0][0][1]) ;
for(i=0;i<n-1;i++)
{ printf("C ") ;
for(j=1;j<4;j++) printf("%11.6f %11.6f ",b[i][j][0],b[i][j][1]) ;
if(i==n-2) printf("\"") ;
printf("\n ") ;
}
printf("fill-opacity=\"0\" ") ;
printf(" style=\"stroke:%s;stroke-width:1.5\"/>\n",col) ;
freematrix(Y) ; freematrix(b) ;
}
int main()
{ double **x=matrix(20,2),***b,q,h,*p ;
int i,j,lino,n ;
char *col ;
printf("<?xml version=\"1.0\" encoding=\"utf-8\"?>\n"
"<svg xmlns=\"http://www.w3.org/2000/svg\"\n"
" xmlns:xlink=\"http://www.w3.org/1999/xlink\"\n"
" width=\"640\" height=\"640\">\n"
"<title>Spatial Borda weights</title>\n"
"<defs><marker id=\"arrow\" markerWidth=\"10\" markerHeight=\"10\" "
" refX=\"9\" refY=\"3\" orient=\"auto\" markerUnits=\"strokeWidth\">\n"
"<path d=\"M0,0 L0,6 L9,3 z\" fill=\"#000\"/>\n</marker>\n</defs>\n") ;
printf("<line x1=\"20\" y1=\"620\" x2=\"20\" y2=\"20\" style="
"\"stroke:black;stroke-width:1\" marker-end=\"url(#arrow)\"/>\n") ;
printf("<line x1=\"20\" y1=\"620\" x2=\"620\" y2=\"620\" style="
"\"stroke:black;stroke-width:1\" marker-end=\"url(#arrow)\"/>\n") ;
for(i=1;i<10;i++)
printf("<line x1=\"%.1f\" y1=\"616\" x2=\"%.1f\" y2=\"620\" style="
"\"stroke:black;stroke-width:1\"/>\n",xmap(i),xmap(i)) ;
for(lino=0;lino<6;lino++)
{ n = 10 ;
for(p=0,i=0;i<10;i++) x[i][0] = i ;
if(lino==0) for(col="red",i=0;i<10;i++) x[i][1] = i ;
else if(lino==1)
{ for(i=0;i<19;i++) { x[i][0] = 0.5*i ; x[i][1] = 10/(10.0-0.5*i)-1 ; }
col = "orange" ;
n = 19 ;
}
else if(lino==2) { p = w1 ; col = "mediumblue" ; }
else if(lino==3) { p = w2 ; col = "dodgerblue" ; }
else if(lino==4) { p = w ; col = "darkseagreen" ; }
else if(lino==5) { p = ww ; col = "seagreen" ; }
if(p) for(i=0;i<10;i++) x[i][1] = p[9-i] ;
drawX(x,n,col) ;
for(i=0;i<10;i++)
{ if(lino==1) j = 2*i ; else j = i ;
printf("<circle cx=\"%.1f\" cy=\"%.1f\" r=\"2.2\" fill=\"%s\" "
"style=\"stroke-width:0;fill-opacity:1\"/>\n",
xmap(x[j][0]),ymap(x[j][1]),col) ;
}
if(lino==0) printf("<text x=\"220\" y=\"412\" font-size=\"17\" "
"transform=\"rotate(-45 220 412)\" fill=\"%s\">borda</text>\n",col) ;
else if(lino==1) printf("<text x=\"472\" y=\"470\" font-size=\"17\" "
"transform=\"rotate(-50 472 470)\" fill=\"%s\">nauru</text>\n",col) ;
else if(lino==2) printf("<text x=\"240\" y=\"331\" font-size=\"17\" "
"transform=\"rotate(-45 240 331)\" fill=\"%s\">ebc<tspan dy=\"3\" "
"font-size=\"13\">1</tspan></text>\n",col) ;
else if(lino==3) printf("<text x=\"240\" y=\"366\" font-size=\"17\" "
"transform=\"rotate(-45 240 366)\" fill=\"%s\">ebc<tspan dy=\"3\" "
"font-size=\"13\">2</tspan></text>\n",col) ;
else if(lino==4) printf("<text x=\"258\" y=\"457\" font-size=\"17\" "
"transform=\"rotate(-45 258 457)\" fill=\"%s\">sbc<tspan dy=\"3\" "
"font-size=\"13\">1</tspan></text>\n",col) ;
else if(lino==5) printf("<text x=\"438\" y=\"272\" font-size=\"17\" "
"transform=\"rotate(-45 438 272)\" fill=\"%s\">sbc<tspan dy=\"3\" "
"font-size=\"13\">2</tspan></text>\n",col) ;
}
printf("</svg>\n") ;
}
memory.h linsolve.c bezier.c beziertest.c sbcplot.c
The following call lines compile and execute the tests:
g++ -g -O -o beziertest beziertest.c bezier.c linsolve.c ; beziertest g++ -g -O -o sbcplot sbcplot.c linsolve.c bezier.c ; sbcplot
The correct output from beziertest is
( 1,-1) ( 0, 0) (-1, 1) (-1, 2) (-1, 2) (-1, 3) ( 0, 4) ( 1, 4) ( 1, 4) ( 2, 4) ( 3, 3) ( 4, 3) ( 4, 3) ( 5, 3) ( 6, 4) ( 7, 5)
doc : reference : source : test : download
doc : algorithm : source : quadrate.h : test : download
To use quadrate you should
#include "quadrate.h"
which defines a data structure and the function prototypes.
The main prototypes are as follows:
quadratereport quadraten(void (*f)(double,double *),int n,double *res,
double **pt,int npt,double qtol,int maxeval) ;
double quadrate(double (*f)(double),xy X,xy Y,xy Z,double qtol) ;
The second of these (quadrate) can be used for a function which computes (and returns) a single value.
The return value is the integral as evaluated.
The first call (quadraten) is for functions which compute some number n of values, all of which are integrated independently and the results returned in res. It also provides more control for advanced uses. You may supply a grid comprising any number n of initial ordinands (n≥2) with their associated function values.
The grid must be supplied in increasing order of ordinate (duplicates are ignored), so pt[0] and pt[npt-1] are the limits of the integration.
A third prototype is provided for convenience.
quadratereport quadraten(void (*f)(double *,double *),int n,double *res,
double **pt,int npt,double qtol,int maxeval) ;
in which f is passed a pointer to its argument rather than the argument itself. This is for compatibility with cubate.
quadrate may apply a relative or an absolute tolerance. If the value you supply is positive, eg. 10–6, then quadrate will try to find an answer to within one part in a million. If you supply a negative tolerance, eg. –10–6, then quadrate will try to find an answer which is correct to within 0.000001.
If it is unable to obtain the requested tolerance within the specified number of function evaluations, quadrate terminates in in one of two ways. If the first call (quadraten) is used, then it returns normally putting in res whatever results it has obtained so far. The user must consult the returned report to find out whether the result can be trusted. If the second call (quadrate) is used, then it simply crashes out if it is unable to satisfy the requested tolerance. At least you know that any value returned can be trusted.
The return value from the quadraten call is a report giving various statistics about the integration.
struct quadratereport
{ int neval,nround,thinslices,*leveval,nleveval,nzero ; double condition ; } ;
The condition number will be <1 (and ≥0) if the integral converged, and >1 otherwise.
To be specific, quadrate attempts to reduce the estimated error to below qtol (relatively or absolutely), and the condition number is the ratio of the estimated error to qtol on completion. Thus something slightly greater than 1 doesn’t matter, but when convergence is unattainable the condition number is likely to be huge (and the returned integral wholly unreliable).
Dependencies: memory.h , quadrate.h
I don’t know any good algorithms for numerical integration in one dimension (those in NR are hopeless). Probability functions are intrinsically difficult because in limiting cases they will be spiky – any non-adaptive algorithm will waste its effort in areas where the integrand is effectively zero. (Monte Carlo methods are adaptive by construction, but inefficient in spaces of low dimension.)
quadrate employs a simple-minded method of my own which is specifically designed to cope with for probability functions. If you know the integrand to be spiky you can supply an internal point, or even more than one, to guide the integration towards regions in which the integrand doesn’t vanish. You can find a suitable location by hillclimbing in the log domain. This saves quadrate from evaluating the function as zero everywhere it looks, and abandoning the task without noticing a spike in the interstices of the grid.
In my experience, it’s almost always worth preceding a numerical integration by a hillclimb when dealing with probability functions – the risk of missing a spiky peak is otherwise too high. The calling sequence for quadraten allows me to save the function evaluations from the hillclimb and input them to the integration, though the amount of work saved is negligible (say 2%). Hillclimbs are much cheaper than integrations.
I describe the case for a single function; the generalisation to several is straightforward.
A small initial grid is set up and the integration then takes place in a series of rounds. In each round a threshold is computed, and intervals are split by bisection (recursively) until no error exceeds the threshold. The integral over an interval is estimated by cubic interpolation, using the values at the end points and at the ends of the adjoining intervals. The integral of a cubic is actually computed as a left and right half-integral, separated by the mid point of the interval. Then, when an interval is split, the error in the left-hand half is estimated as the absolute difference between the new integral and the left-hand half integral of the interval before splitting; and likewise for the right hand side.
To ensure that the cubic approximation is reasonably good, an interval is not made much smaller than its neighbours: to prevent this from happening, the neighbours may be split at the same time.
The threshold is actually whichever is tighter of two thresholds. One is the maximum permitted total error, since no interval may have an error larger than this. The second, which becomes effective in the second round, is computed in such a way that if the error is made to vanish in all intervals whose current error is greater than this value, then the total error will satisfy the threshold. And since both thresholds are a little conservative, they are multiplied by ¾ for good measure.
The recursive splitting is terminated if the error for an interval is less than the total permitted error, and the last split didn’t appreciably reduce the error (even if the second threshold isn’t satisfied). This is to save quadrate from thrashing in the vicinity of a knot if the termination criterion can be satisfied by improving accuracy elsewhere.
A relative tolerance is computed proportional to the estimated integral of the absolute value of the integrand; thus one can meaningfully integrate sin x from 0 to 2π with a relative tolerance.
• radixsort • point • ~point • span • midpt • ship • err • quadratic • cubic • integrate • eval • split • splinter • reduce • stub • quadraten • quadrate
#include "memory.h" #include <math.h> #include "quadrate.h" static void radixsort(unsigned short *a,int n) // sort ascending { unsigned short *b=(unsigned short *) cjcalloc(n,sizeof(unsigned short)),k ; int i,bytepos,count[256],offs[256] ; for(i=0;i<256;i++) count[i] = 0 ; for(i=0;i<n;i++) count[a[i]&0xff] += 1 ; for(offs[0]=0,i=1;i<256;i++) offs[i] = offs[i-1] + count[i-1] ; for(i=0;i<n;i++) { k = a[i]&0xff ; b[offs[k]++] = a[i] ; } for(i=0;i<256;i++) count[i] = 0 ; for(i=0;i<n;i++) count[(b[i]&0xffff)>>8] += 1 ; for(offs[0]=0,i=1;i<256;i++) offs[i] = offs[i-1] + count[i-1] ; for(i=0;i<n;i++) { k = (b[i]&0xffff)>>8 ; a[offs[k]++] = b[i] ; } free(b) ; } /* -------------------------------------------------------------------------- */ struct point /* the two half-areas lie between this point and the next */ { double x,*y,*lhs,*rhs,*prevarea,w[4] ; point *prev,*next ; unsigned char integrated ; point(double xval,point *pval,point *nval,int n) { x = xval ; prev = pval ; next = nval ; integrated = 0 ; if(n==1) { y = w ; for(int i=0;i<4;i++) y[i] = 0 ; } else y = vector(4*n) ; lhs = y+n ; rhs = y+2*n ; prevarea = y+3*n ; } ~point() { if(y!=w) free(y) ; } double span() { return next->x - x ; } double midpt() { return (x+next->x)/2 ; } double ship(double *yy,int n) { for(int i=0;i<n;i++) yy[i] = y[i] ; return x ; } double err(int i) { return fabs(lhs[i]+rhs[i]-prevarea[i]) ; } void quadratic(double *xx,double **yy,int n,int offs) { double qa,qb,qa2,qa1,qy0,qy1,qxp[3],qx[3],x[3] ; int i,j ; x[0] = xx[offs] ; for(j=1;j<3;j++) x[j] = xx[offs+j] - x[0] ; if(offs==0) { qx[0] = x[1] ; qx[2] = x[2] ; } else { qx[0] = 0 ; qx[2] = x[1] ; } qx[1] = (qx[0]+qx[2]) / 2 ; qa2 = 1 / (x[2]*(x[2]-x[1])) ; qa1 = -1 / (x[1]*(x[2]-x[1])) ; for(i=0;i<n;i++) { qy0 = yy[offs][i] ; qy1 = yy[offs+1][i] - qy0 ; qa = qa2*(yy[offs+2][i]-qy0) + qa1*qy1 ; qb = qy1/x[1] - qa*x[1] ; for(j=0;j<3;j++) qxp[j] = qx[j] ; lhs[i] = rhs[i] = qy0 * (qxp[1]-qxp[0]) ; for(j=0;j<3;j++) qxp[j] *= qx[j] ; lhs[i] += qb*(qxp[1]-qxp[0])/2 ; rhs[i] += qb*(qxp[2]-qxp[1])/2 ; for(j=0;j<3;j++) qxp[j] *= qx[j] ; lhs[i] += qa*(qxp[1]-qxp[0])/3 ; rhs[i] += qa*(qxp[2]-qxp[1])/3 ; } } void cubic(double *xx,double **yy,int n) { double q12,q13,q23,x[4],qy[4],qx[3],qxp[3],qa,qb,qc,qa3,qa2,qa1,qb2,qb1 ; int i,j ; x[0] = xx[0] ; for(j=1;j<4;j++) x[j] = xx[j] - xx[0] ; qx[0] = x[1] ; qx[1] = (x[1]+x[2]) / 2 ; qx[2] = x[2] ; q12 = 1 / (x[2]-x[1]) ; q13 = 1 / (x[3]-x[1]) ; q23 = 1 / (x[3]-x[2]) ; qa3 = q13*q23/x[3] ; qa2 = -q12*q23/x[2] ; qa1 = q12*q13/x[1] ; qb2 = q12/x[2] ; qb1 = -q12/x[1] ; for(i=0;i<n;i++) { qy[0] = yy[0][i] ; for(j=1;j<4;j++) qy[j] = yy[j][i] - qy[0] ; for(j=0;j<3;j++) qxp[j] = qx[j] ; qa = qa3*qy[3] + qa2*qy[2] + qa1*qy[1] ; qb = qy[2]*qb2 + qy[1]*qb1 - qa*(x[1]+x[2]) ; qc = qy[2]/x[2] - x[2]*x[2]*qa - x[2]*qb ; lhs[i] = rhs[i] = qy[0] * (qxp[1]-qxp[0]) ; for(j=0;j<3;j++) qxp[j] *= qx[j] ; lhs[i] += qc*(qxp[1]-qxp[0])/2 ; rhs[i] += qc*(qxp[2]-qxp[1])/2 ; for(j=0;j<3;j++) qxp[j] *= qx[j] ; lhs[i] += qb*(qxp[1]-qxp[0])/3 ; rhs[i] += qb*(qxp[2]-qxp[1])/3 ; for(j=0;j<3;j++) qxp[j] *= qx[j] ; lhs[i] += qa*(qxp[1]-qxp[0])/4 ; rhs[i] += qa*(qxp[2]-qxp[1])/4 ; } } void integrate(double **yy,int n) // yy is work space { double xx[4] ; xx[1] = ship(yy[1],n) ; xx[2] = next->ship(yy[2],n) ; if(prev==0) { xx[3] = next->next->ship(yy[3],n) ; quadratic(xx,yy,n,1) ; } else { xx[0] = prev->ship(yy[0],n) ; if(next->next==0) quadratic(xx,yy,n,0) ; else { xx[3] = next->next->ship(yy[3],n) ; cubic(xx,yy,n) ; } } integrated = 1 ; } void eval(void *f,int mode,int &neval) { if(mode==2) ((void (*)(double *,double *))f)(&x,y) ; else if(mode) ((void (*)(double,double *))f)(x,y) ; else y[0] = ((double (*)(double))f)(x) ; neval += 1 ; } point *split(void *f,int n,int mode,int &neval) { point *p=next ; next = p->prev = new point(midpt(),this,next,n) ; next->eval(f,mode,neval) ; return p ; } double splinter(void *f,int n,int mode,double **yy,int &neval,double *asum) { int hardest,i ; double preverr,maxerr,thiserr,xmid=midpt() ; if(xmid==x||xmid==next->x) return -1 ; for(preverr=maxerr=hardest=i=0;i<n;i++) if(asum[i]>0&&err(i)/asum[i]>preverr) { maxerr = err(i) ; preverr = maxerr/asum[i] ; hardest = i ; } if(prev) prev->integrated = 0 ; point *p = split(f,n,mode,neval) ; p->integrated = 0 ; for(i=0;i<n;i++) { prevarea[i] = lhs[i] ; next->prevarea[i] = rhs[i] ; } integrate(yy,n) ; next->integrate(yy,n) ; thiserr = err(hardest) + next->err(hardest) ; if(maxerr==0) return 1 ; else if(thiserr>asum[hardest]) return 0 ; // always recurse if error... else return thiserr / maxerr ; // >max permitted total error } point *reduce(void *f,int n,int mode,double *thr,double *asum,double **yy, quadratereport &report,int maxeval) { int i ; double q ; point *p=next ; for(i=0;i<n&&(thr[i]<0||(err(i)<thr[i]&&err(i)<0.75*asum[i]));i++) ; if(i==n||report.neval>=maxeval) return p ; if(prev) while(prev->span()>span()) if(0>prev->splinter(f,n,mode,yy,report.neval,asum)) { report.thinslices += 1 ; break ; } if(next->next) while(next->span()>span()) if(0>next->splinter(f,n,mode,yy,report.neval,asum)) { report.thinslices += 1 ; break ; } if(0>(q=splinter(f,n,mode,yy,report.neval,asum))) { report.thinslices += 1 ; return p ; } if(q<0.5) // don't recurse unless beneficial { next->reduce(f,n,mode,thr,asum,yy,report,maxeval) ; reduce(f,n,mode,thr,asum,yy,report,maxeval) ; } return p ; } } ; /* -------------------------------------------------------------------------- */ static quadratereport stub(void *f,int n,int mode,double *res,double qtol, double **pt,int npt,int maxeval) { int i,k,flip=0,ngrid,undone,*allzero=ivector(n),ind ; double q,qq,maxerr,*sum=vector(n),*asum=vector(n),ftol=fabs(qtol) ; double **yy=matrix(4,n),*thr=vector(n),*thr2=vector(n),*err=vector(n) ; point *first,*p ; unsigned short *lerr ; quadratereport report ; // set up a grid for(p=0,ngrid=i=0;i<npt;i++) { if(i>0&&pt[i][0]==pt[i-1][0]) continue ; if(i>0&&pt[i][0]<pt[i-1][0]) throw up("quadrate error: points %d/%d out of order (%.6f,%.6f)", i-1,i,pt[i-1][0],pt[i][0]) ; p = new point(pt[i][0],p,0,n) ; if(p->prev) p->prev->next = p ; else first = p ; for(k=0;k<n;k++) p->y[k] = pt[i][1+k] ; ngrid += 1 ; } while(ngrid<5) for(p=first;p->next;ngrid++) p = p->split(f,n,mode,report.neval) ; for(p=first;p->next;p=p->next) p->integrate(yy,n) ; // find where the function evaluates to zero everywhere for(i=0;i<n;i++) { for(p=first;p&&p->y[i]==0;p=p->next) ; allzero[i] = !p ; } for(report.nround=0;;report.nround++) { // compute sum and asum for(i=0;i<n;i++) if(!allzero[i]) for(err[i]=sum[i]=asum[i]=0,p=first;p->next;p=p->next) { q = p->lhs[i] + p->rhs[i] ; sum[i] += q ; asum[i] += fabs(q) ; err[i] += fabs(q-p->prevarea[i]) ; } if(qtol<0) for(i=0;i<n;i++) asum[i] = 1 ; for(maxerr=k=i=0;i<n;i++) if(!allzero[i]) { if(qtol>0&&asum[i]<=0) { allzero[i] = 1 ; continue ; } err[i] /= asum[i] ; if(err[i]>maxerr) maxerr = err[i] ; if(err[i]>ftol&&report.neval<maxeval) k = 1 ; } if(k==0) break ; for(maxerr=ngrid=0,p=first;p->next;p=p->next,ngrid++) for(i=0;i<n;i++) if(!allzero[i]) maxerr = max(p->err(i)/asum[i],maxerr) ; // now list the errors lerr = (unsigned short *) cjcalloc(ngrid,sizeof(unsigned short)) ; for(ind=0,p=first;p->next;p=p->next,ind++) { for(q=i=0;i<n;i++) if(!allzero[i]) q = max(p->err(i)/asum[i],q) ; lerr[ind] = (unsigned short) (0.5+q*(64000/maxerr)) ; } radixsort(lerr,ngrid) ; for(q=i=0;q<ftol&&i<ngrid;i++) q += lerr[i] * (maxerr/64000) ; if(i==0) q = lerr[0] ; else if(i==ngrid) q = lerr[i-1] ; else q = ( lerr[i-1] + lerr[i] ) / 2 ; q *= maxerr/64000 ; for(i=0;i<n;i++) if(allzero[i]) thr[i] = thr2[i] = -1 ; else { thr[i] = 0.75 * q * asum[i] ; thr2[i] = ftol * asum[i] ; } free(lerr) ; for(p=first;p->next&&report.neval<maxeval;) p = p->reduce(f,n,mode,thr,thr2,yy,report,maxeval) ; for(p=first;p->next;p=p->next) if(p->integrated==0) p->integrate(yy,n) ; } while(first->next) { first = first->next ; delete first->prev ; } delete first ; for(i=0;i<n;i++) res[i] = sum[i] ; free(err,sum,asum,allzero,thr,thr2) ; freematrix(yy) ; if(mode==0&&report.neval>=maxeval) throw up("too many function evals performed (%d)",report.neval) ; report.condition = maxerr/ftol ; return report ; } /* -------------------------------------------------------------------------- */ quadratereport quadraten(void (*f)(double,double *),int n,double *res, double **pt,int npt,double qtol,int maxeval) { return stub((void *)f,n,1,res,qtol,pt,npt,maxeval) ; } quadratereport quadraten(void (*f)(double *,double *),int n,double *res, double **pt,int npt,double qtol,int maxeval) { return stub((void *)f,n,2,res,qtol,pt,npt,maxeval) ; } double quadrate(double (*f)(double),xy X,xy Y,xy Z,double qtol) { double res,**pt=matrix(3,2) ; pt[0][0] = X.x ; pt[0][1] = X.y ; pt[1][0] = Y.x ; pt[1][1] = Y.y ; pt[2][0] = Z.x ; pt[2][1] = Z.y ; stub((void *)f,1,0,&res,qtol,pt,3,20000) ; freematrix(pt) ; return res ; }
struct quadratereport { int neval,nround,thinslices,*leveval,nleveval,nzero ; double condition ; quadratereport() { neval = nround = thinslices = nleveval = nzero = 0 ; leveval = 0 ; condition = -1 ; } void release() { free(leveval) ; this[0] = quadratereport() ; } } ; // quadrate prototypes quadratereport quadraten(void (*f)(double,double *),int n,double *res, double **pt,int npt,double qtol,int maxeval) ; quadratereport quadraten(void (*f)(double *,double *),int n,double *res, double **pt,int npt,double qtol,int maxeval) ; double quadrate(double (*f)(double),xy X,xy Y,xy Z,double qtol) ; // cubate prototypes double cubate(double (f)(double*),int dim,double qtol,int maxeval) ; double cubate(double (f)(double*,int),int dim,double qtol,int maxeval) ; quadratereport cubate(void (f)(double*,double *),int dim,int n, double *res,double qtol,int maxeval) ; quadratereport cubate(void (f)(double*,double *,int,int),int dim,int n, double *res,double qtol,int maxeval) ; // cubate1 prototypes double cubate1(double (f)(double*),int dim,double qtol,int maxeval) ; double cubate1(double (f)(double*,int),int dim,double qtol,int maxeval) ; quadratereport cubate1(void (f)(double*,double *),int dim,int n, double *res,double qtol,int maxeval) ; quadratereport cubate1(void (f)(double*,double *,int,int),int dim,int n, double *res,double qtol,int maxeval) ; // the special call for nested tolerances quadratereport cubate1(void (f)(double*,double *,double *),int dim,int n, double *res,double qtol,int maxeval) ;
There are two test programs – quadratetest1 and quadratetestn. The first integrates a single function using both the shorter and longer call, and the second integrates 3 functions to find the first two moments of a beta distribution. I get the following outputs respectively:
tol = 1.0e-02 ; 15 fn evals (2 rounds); relative error = 1.5e-03; cond=0.6 tol = 1.0e-03 ; 26 fn evals (3 rounds); relative error = 5.4e-05; cond=1.0 tol = 1.0e-04 ; 50 fn evals (3 rounds); relative error = 1.2e-06; cond=0.7 tol = 1.0e-05 ; 82 fn evals (3 rounds); relative error = 4.4e-08; cond=0.9 tol = 1.0e-06 ; 139 fn evals (4 rounds); relative error = 2.7e-08; cond=0.8 tol = 1.0e-07 ; 246 fn evals (4 rounds); relative error = 3.9e-09; cond=0.8 tol = 1.0e-08 ; 427 fn evals (3 rounds); relative error = 2.2e-10; cond=0.8 tol = 1.0e-09 ; 748 fn evals (3 rounds); relative error = 4.2e-11; cond=0.9 tol = 1.0e-10 ; 1312 fn evals (3 rounds); relative error = 5.1e-12; cond=0.9 tol = 1.0e-02 ; 13 fn evals (1 rounds); absolute error = 3.4e-03; cond=1.0 tol = 1.0e-03 ; 24 fn evals (3 rounds); absolute error = 3.5e-05; cond=0.7 tol = 1.0e-04 ; 44 fn evals (4 rounds); absolute error = 3.9e-06; cond=0.8 tol = 1.0e-05 ; 78 fn evals (3 rounds); absolute error = 4.0e-08; cond=0.9 tol = 1.0e-06 ; 132 fn evals (4 rounds); absolute error = 2.8e-08; cond=0.8 tol = 1.0e-07 ; 221 fn evals (4 rounds); absolute error = 3.9e-09; cond=1.0 tol = 1.0e-08 ; 401 fn evals (3 rounds); absolute error = 3.0e-10; cond=0.8 tol = 1.0e-09 ; 700 fn evals (3 rounds); absolute error = 4.6e-11; cond=0.9 tol = 1.0e-10 ; 1219 fn evals (3 rounds); absolute error = 4.6e-12; cond=0.9
tol = 1.0e-02 ; 25 fn evals; moments ~= 0.450017518609 0.022492992249 tol = 1.0e-03 ; 49 fn evals; moments ~= 0.450001373728 0.022499450507 tol = 1.0e-04 ; 86 fn evals; moments ~= 0.450000165812 0.022499961119 tol = 1.0e-05 ; 151 fn evals; moments ~= 0.450000011228 0.022500007324 tol = 1.0e-06 ; 271 fn evals; moments ~= 0.450000001100 0.022499999792 tol = 1.0e-07 ; 386 fn evals; moments ~= 0.449999999824 0.022500000091 tol = 1.0e-08 ; 751 fn evals; moments ~= 0.450000000005 0.022500000015 tol = 1.0e-09 ; 1463 fn evals; moments ~= 0.449999999999 0.022499999999 tol = 1.0e-10 ; 2817 fn evals; moments ~= 0.450000000000 0.022500000000
#include "memory.h" #include <math.h> #include "quadrate.h" static double f(double x) { double u=1-x*x ; if(u>0) return sqrt(u) ; else return 0 ; } static void g(double x,double *yy) { yy[0] = f(x) ; } static void h(double *x,double *yy) { yy[0] = f(x[0]) ; } int main() { double tol,u,v,qerr,pi=3.1415926535897932384626433832795029,**pt=matrix(3,2) ; xy X=xy(0,f),Y=xy(0.5,f),Z=xy(1,f) ; quadratereport q ; pt[0][0] = X.x ; pt[0][1] = X.y ; pt[1][0] = Y.x ; pt[1][1] = Y.y ; pt[2][0] = Z.x ; pt[2][1] = Z.y ; for(int sign=1;sign>-2;sign-=2) for(tol=1e-2;tol>5e-11;tol/=10) { q = quadraten(g,1,&v,pt,3,sign*tol,20000) ; quadraten(h,1,&u,pt,3,sign*tol,20000) ; if(u!=v) throw up("discrepancy between quadraten prototypes") ; u = quadrate(f,X,Y,Z,sign*tol) ; if(sign>0) qerr = fabs(4*v/pi-1) ; else qerr = fabs(v-pi/4) ; printf("tol = %.1e ;%5d fn evals (%d rounds); %s error = " "%.1e; cond=%.1f\n",tol,q.neval+3,q.nround, sign>0?"relative":"absolute",qerr,q.condition) ; if(u!=v) printf("^^^ unequal results (%.12f,%.12f) ^^^\n",u,v) ; } }
#include "memory.h" #include <math.h> #include "quadrate.h" void f(double x,double *y) { y[2] = x * ( y[1] = x * ( y[0] = pow(x,3.5)*pow(1-x,4.5) ) ) ; } int main() { double **pt=matrix(3,4),tol,v[3] ; quadratereport report ; for(int i=0;i<3;i++) { pt[i][0] = 0.5*i ; f(pt[i][0],&pt[i][1]) ; } for(tol=1e-2;tol>5e-11;tol/=10) { report = quadraten(f,3,v,pt,3,tol,20000) ; v[1] /= v[0] ; v[2] = v[2]/v[0] - v[1]*v[1] ; printf("tol = %.1e ;%5d fn evals; moments ~= %.12f %.12f\n", tol,report.neval,v[1],v[2]) ; } }
memory.h quadrate.c quadrate.h quadratetest1.c quadratetestn.c
The following call lines compile and execute the two tests, using g++ as compiler:
g++ -o quadratetest1 -g -O quadratetest1.c quadrate.c ; quadratetest1 g++ -o quadratetestn -g -O quadratetestn.c quadrate.c ; quadratetestn
doc : algorithm : source : quadrate.h : test : download
doc : algorithm : weaknesses : source : test : download
cubate is an experimental cubature subroutine which comes with no guarantees. 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≤xi≤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 (note: dim > 1). 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. Since the user who invokes cubate in this way receives an estimate of the precision of the integral, the subroutine terminates normally rather than crashing out when the maximum number of function evaluations is reached.
The returned structure is as for quadrate, but the nround field is unused whereas the leveval and nzero fields now serve a purpose.
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 new fields of the return structure are as follows:
leveval is set up as an array by cubate (nleveval is its length). leveval[i] is the number of integrations performed 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 release the memory allocated to leveval after return from cubate; the best way to do this is through the release() method of the quadratereport structure.
nzero is returned as the number of simplices for which, for at least one integrand, the function has been evaluated as zero at all points computed. This should not happen unless your integrand is spiky; but when it is spiky, quadrature becomes a difficult task which I am not sure how to address, and which might require a lot of work. So for now I’m contenting myself with providing a diagnostic return which allows the problem to be identified. You can always fall back on Monte Carlo integration if there is only one spike.
Dependencies: memory.h , quadrate.h , infheap.h
The algorithm has several components:
I use the standard recursive subdivision of a hypercube into simplices. The hypercube is represented as the set of inequalities 0 ≤x0, x1, x2, ... ≤ 1. This is then partitioned into polytopes given by the sets of inequalities
and these are split in turn into polytopes such as 0 ≤x4 ≤ x2 ≤ x0, x1, ... ≤ 1; and so forth until we obtain inequalities similar to 0 ≤x4 ≤ x2 ≤ x0 ≤ x3 ≤ x1 ≤ 1 which define simplices.
I don’t know how widely Lipson’s formula is known [1]. Start by writing μ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 the formula states that:
| μq = | μv + (d+1)μc |
| d+2 |
The integral of a quadratic over a simplex can then be computed as the product of μq with the volume of the simplex.
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 – (d–2)μv |
| d+2 |
Notice that in two dimensions – i.e. integrating over triangles – μv drops out.
I previously used a cheaper criterion, namely the difference between the integral as computed by Lipson’s formula and the integral as implied by a linear fit, i.e. vol×μv. This amounts to subdividing simplices until the error in a linear fit satisfies the tolerance, and then returning the quadratic fit.
The new criterion is more work (because function evaluations are needed at all midpoints; there is a lot of reuse but still a lot of work whose only contribution is to deciding when to terminate). No real improvement is seen from the test program, but that’s not what it’s for.
We learn more when we consider the performance of the termination criterion in euclideanborda. With the old criterion the numbers of integrations by level were reported as follows:
dim=1, cond=1.0; neval=9403 dim=2, cond=1.0; neval=10038
With the improved criterion I get:
dim=1, cond=1.0; neval=1026 dim=2, cond=1.0; neval=1121
The accuracies are not as good, but they are within the requested tolerance. The saving in time is nearly a factor of 10.
[1]. Andrew Lipson, personal communication, 2005 (“Simpson’s rule on simplices and numerical integration”).
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. However (½,½,½,... ) is guaranteed to be chosen early on as an argument to the integrand, so if you transform your space to have its peak here, it will not be missed by the algorithm.
Some attention to round-off error might increase the accuracy; and no doubt there are casual faults arising from my own ignorance.
• factorial • adj • control • freeslist • eval • addleveval • vertex • release • addneighbour • delneighbour • addneighbour • delneighbour • addvertex • release • polytope • release • getcentroid • setmidpts • setvol • integrate • findedge • increment • extend • fanout • deleter • stub • cubate • cubate1
#include "memory.h" #include <math.h> #include "quadrate.h" #include "infheap.h" int lu(double **x,double **y,int n,int *ind) ; static double factorial(int m) { int i ; double q ; for(q=i=1;i<=m;i++) q *= i ; return q ; } /* ----------------------------------- adj ---------------------------------- */ struct vertex ; struct polytope ; 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(polytope*,ppvector) ; // a ppvector is a vector of (polytope *)s genvector(adj,adjvector) ; /* --------------------------------- control -------------------------------- */ struct control { int neval,maxeval,thinslices,mode,nv,np,*leveval,nleveval,nzero ; vertex *vlist ; polytope **plist ; void *f ; control() { nleveval = maxeval = neval = thinslices = nv = np = mode = nzero = 0 ; f = 0 ; vlist = 0 ; plist = 0 ; leveval = 0 ; } control(void *ff,int momo,int dim,int n,int mm) { this[0] = control() ; f = ff ; mode = momo ; maxeval = mm ; plist = ppvector((int) factorial(dim)) ; } void release() ; void freeslist() { free(plist) ; plist = 0 ; np = 0 ; } vertex *addvertex(int dd,int nn,double *xx,double *yguess) ; void eval(double *x,double *y,int dim,int n,double *yguess) { neval += 1 ; if((mode&3)<2) { if((mode&3)==0) y[0] = ((double (*)(double *))f)(x) ; else y[0] = ((double (*)(double *,int))f)(x,dim) ; } else { if(mode==14) ((void (*)(double *,double *,double *))f)(x,y,yguess) ; else if((mode&3)==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*lev + 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,double *yguess) { 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,yguess) ; } 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,double *yguess) { vertex *v=new vertex(dd,nn,this,xx,yguess) ; nv += 1 ; v->next = vlist ; return vlist = v ; } void control::release() { vertex *v ; polytope *p ; while(vlist->next) { v = vlist ; vlist = vlist->next ; v->release() ; delete v ; } vlist->release() ; delete vlist ; this[0] = control() ; } /* --------------------------------- polytope ------------------------------- */ struct polytope // a polytope of level dim-1 or more is a simplex { double vol,*hypervol,*err,maxerr ; int dim,n,lev,ndescendents,*fixed ; vertex **v ; control *c ; polytope **descendents ; polytope() { vol = maxerr = ndescendents = lev = dim = n = 0 ; err = hypervol = 0 ; v = 0 ; c = 0 ; fixed = 0 ; descendents = 0 ; } polytope(int dd,int nn,control *cc,int l) { this[0] = polytope() ; dim = dd ; n = nn ; c = cc ; lev = l ; if(lev<=0) return ; if(lev<dim-1) { fixed = ivector(lev) ; return ; } err = vector(n) ; hypervol = vector(n) ; v = vpvector(dim+1) ; } polytope(polytope *p) { int t ; this[0] = p[0] ; if(lev<dim-1) throw up("can\'t clone a polytope at level %d/%d",lev,dim) ; v = vpvector(dim+1) ; for(t=0;t<dim+1;t++) v[t] = p->v[t] ; err = vector(n) ; for(t=0;t<n;t++) err[t] = p->err[t] ; hypervol = vector(n) ; for(t=0;t<n;t++) hypervol[t] = p->hypervol[t] ; } void release() { if(lev<dim-2) for(int i=0;i<ndescendents;i++) { descendents[i]->release() ; delete descendents[i] ; } free(fixed,descendents,hypervol,err,v) ; } void getcentroid(double *y,double *yguess) { 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) ; c->eval(x,y,dim,n,yguess) ; free(x) ; } void setmidpts(double *y,double *yguess) { int v0,v1,t ; adj *edge ; double *x=vector(dim) ; for(t=0;t<n;t++) y[t] = 0 ; for(v0=0;v0<dim;v0++) for(v1=v0+1;v1<dim+1;v1++) { edge = findedge(v0,v1) ; if(!edge->midpt) { for(t=0;t<dim;t++) x[t] = ( v[v0]->x[t] + v[v1]->x[t] ) / 2 ; edge->midpt = c->addvertex(dim,n,x,yguess) ; } for(t=0;t<n;t++) y[t] += edge->midpt->y[t] ; } free(x) ; for(t=0;t<n;t++) y[t] /= ( dim*(dim+1) ) / 2 ; } double setvol() { int i,j ; double q , **x = matrix(dim,dim) ; for(i=0;i<dim;i++) for(j=0;j<dim;j++) x[i][j] = v[i]->x[j] - v[dim]->x[j] ; lu(x,x,dim,0) ; // often the volume is just halved... but when? for(q=1,i=0;i<dim;i++) q *= x[i][i] ; // determinant freematrix(x) ; return vol = fabs(q) / factorial(dim) ; } void integrate(double *asum,double *yguess) { int t,i ; double q,qv,ql,qq=vol/(dim+2),*qm=vector(n),*qc=vector(n) ; vertex *vtx0,*vtxm ; if((c->mode&4)==0) setmidpts(qm,0) ; // not cubate1 getcentroid(qc,lev>dim-1?yguess:0) ; for(maxerr=i=0;i<n;i++) { for(qv=t=0;t<=dim;t++) qv += v[t]->y[i] / (dim+1) ; ql = ( (dim+1)*qc[i] + qv ) / (dim+2) ; if((c->mode&4)==0) { hypervol[i] = vol * ql ; err[i] = vol * fabs( ql - (2*dim*qm[i]-(dim-2)*qv)/(dim+2) ) ; } else // cubate1 { hypervol[i] = vol * ((dim+1)*qv+qc[i]) / (dim+2) ; err[i] = vol * fabs(ql-qv) ; } if(asum[i]&&err[i]/asum[i]>maxerr) maxerr = err[i] / asum[i] ; } c->addleveval(lev) ; free(qm,qc) ; } 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 ------------------------------- */ polytope *extend(double *asum) { int t,v0,v1,m0,m1 ; double qmax,*yguess=0,*x ; polytope *p ; vertex *vtxm,*vtx0,*vtx1 ; adj *edge,*me ; for(qmax=v0=0;v0<dim;v0++) for(v1=v0+1;v1<dim+1;v1++) { edge = findedge(v0,v1) ; if(edge->d2>qmax) { m0 = v0 ; m1 = v1 ; me = edge ; qmax = edge->d2 ; } } if(qmax<1e-20) { c->thinslices += 1 ; maxerr = 0 ; return 0 ; } if(c->mode&8) for(yguess=vector(n),t=0;t<n;t++) yguess[t] = hypervol[t] / vol ; // extend the linked list of vertices vtx0 = v[v0=m0] ; vtx1 = v[v1=m1] ; if(me->midpt==0) { if((c->mode&4)==0) throw up("no midpt") ; // not cubate1 x = vector(dim) ; for(t=0;t<dim;t++) x[t] = ( v[v0]->x[t] + v[v1]->x[t] ) / 2 ; me->midpt = c->addvertex(dim,n,x,yguess) ; free(x) ; } vtxm = me->midpt ; me->count += 1 ; for(t=0;t<=dim;t++) addneighbour(v[t],vtxm) ; delneighbour(vtx0,vtx1) ; // create a new polytope lev += 1 ; p = new polytope(this) ; p->v[v1] = v[v0] = vtxm ; vol -= p->setvol() ; p->integrate(asum,yguess) ; integrate(asum,yguess) ; if(c->mode&8) free(yguess) ; return p ; } /* --------------------------------- fanout ------------------------------- */ void fanout(vertex **corners) { int i,j,t,tau,co , *ind=ivector(dim) ; polytope *p ; for(t=0;t<dim;t++) ind[t] = t ; for(t=0;t<lev;t++) ind[fixed[t]] = -1 ; for(tau=t=0;t<dim;t++) if(ind[t]>=0) ind[tau++] = ind[t] ; descendents = ppvector(ndescendents=tau) ; for(t=0;t<tau;t++) descendents[t] = new polytope(dim,n,c,lev+1) ; if(lev<dim-2) for(t=0;t<tau;t++) { p = descendents[t] ; for(i=0;i<lev;i++) p->fixed[i] = fixed[i] ; p->fixed[lev] = 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->plist[c->np++] = p = descendents[t] ; for(i=0;i<=dim;i++) { for(co=j=0;j<i;j++) co |= 1 << ind[j] ; p->v[i] = corners[co] ; } p->vol = 1 / factorial(dim) ; for(i=0;i<dim;i++) for(j=i+1;j<=dim;j++) addneighbour(p->v[i],p->v[j]) ; if(dim>1) swap(ind[dim-1],ind[dim-2]) ; } } free(ind) ; } double triang(double x[3][2],int lev,int cpt) ; } ; static void deleter(void *v) { polytope *p=(polytope *) v ; p->release() ; delete p ; } // #include "cubatedebug.c" /* -------------------------------------------------------------------------- */ static quadratereport stub(void *f,int dim,int n,double *res, double qtol,int maxeval,int mode) { if(dim<2) throw up("dimension %d(<2) not permitted for cubate",dim) ; int i,t,sno,*allzero=ivector(n) ; double q,*x,ftol=fabs(qtol),*err=vector(n),*sum=vector(n),*asum=vector(n) ; vertex **corners,*v ; quadratereport report ; control c(f,mode,dim,n,maxeval) ; polytope *s,*ss,*sptr , p(dim,n,&c,0) ; infheap 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,0) ; } // partition the hypercube into simplices p.fanout(corners) ; free(corners,x) ; p.release() ; // printneighbours(c.vlist,&c,-1) ; for(sno=0;sno<c.np;sno++) // integrate the simplices { s = c.plist[sno] ; s->integrate(asum,0) ; s->increment(sum,asum,err,1) ; } if(qtol<0) for(i=0;i<n;i++) asum[i] = 1 ; // i.e. if tolerance is absolute // find where the function evaluates to zero everywhere for(i=0;i<n;i++) for(allzero[i]=1,sno=0;sno<c.np;sno++) { for(s=c.plist[sno],t=0;t<=dim;t++) if(s->v[t]->y[i]) allzero[i] = 0 ; if(s->hypervol[i]) allzero[i] = 0 ; } // compute maxerrs and put the vertices on a heap for(sno=0;sno<c.np;sno++) { for(s=c.plist[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 = (polytope *) 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.nzero = c.nzero ; 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() ; c.release() ; if((mode&3)<2&&report.neval>=maxeval) throw up("too many function evals " "performed (%d; tol=%.1e)",report.neval,qtol) ; 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) ; } double cubate1(double (f)(double*),int dim,double qtol,int maxeval) { double res ; stub((void *)f,dim,1,&res,qtol,maxeval,4) ; return res ; } double cubate1(double (f)(double*,int),int dim,double qtol,int maxeval) { double res ; stub((void *)f,dim,1,&res,qtol,maxeval,5) ; return res ; } quadratereport cubate1(void (f)(double*,double *),int dim,int n, double *res,double qtol,int maxeval) { return stub((void *)f,dim,n,res,qtol,maxeval,6) ; } quadratereport cubate1(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,7) ; } quadratereport cubate1(void (f)(double*,double *,double *),int dim,int n, double *res,double qtol,int maxeval) { return stub((void *)f,dim,n,res,qtol,maxeval,14) ; } /* -------------------------------------------------------------------------- */
This is code I wrote while debugging and don’t want to throw away.
• triang • printneighbours • recorder
// this integrates component cpt of the objective function over a triangle // by recursive division into pow(4,lev) subtriangles and averaging the // function values at the centroids of the subtriangles. the result needs to be // multiplied by the volume. double simplex::triang(double x[3][2],int lev,int cpt) { double m[3][2],*y,q,z[2] ; int i,j ; if(lev==0) { for(z[0]=z[1]=i=0;i<3;i++) { z[0] += x[i][0]/3 ; z[1] += x[i][1]/3 ; } y = vector(n) ; c->eval(z,y,2,n) ; q = y[cpt] ; free(y) ; return q ; } for(i=0;i<3;i++) for(m[i][0]=m[i][1]=j=0;j<3;j++) if(j!=i) { m[i][0] += x[j][0]/2 ; m[i][1] += x[j][1]/2 ; } for(q=triang(m,lev-1,cpt),i=0;i<3;i++) { for(j=0;j<2;j++) { z[j] = m[i][j] ; m[i][j] = x[i][j] ; } q += triang(m,lev-1,cpt) ; for(j=0;j<2;j++) m[i][j] = z[j] ; } return q/4 ; } /* ------------------------------ printneighbours --------------------------- */ // this was written to help debug the code supporting the adj structures 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) ; } /* ------------------------------- recorder --------------------------------- */ // you can recurse over recorder to validate the result in 2 dimensions static void recorder(void *v) { simplex *s=(simplex *) v ; int i,t ; double x[3][2] , q , qer = s->err[0]/s->vol , qav , qr ; for(i=0;i<3;i++) for(t=0;t<2;t++) x[i][t] = s->v[i]->x[t] ; q = s->triang(x,3,0) ; qav = s->hypervol[0]/s->vol ; qr = fabs(qav-q)/qer ; for(t=0;t<s->dim+2;t++) printf("(%6.3f,%6.3f:%6.3f) ", s->v[t]->x[0],s->v[t]->x[1],s->v[t]->y[0]) ; printf(": area=%6.3f, avge=%6.3f/%6.3f, diff=%7.4f, err=%7.4f, ratio=%.1f\n", s->vol,qav,q,fabs(qav-q),qer,qr) ; }
The following program 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. Download as cubatetest.c.
euclideanborda is a program which uses cubate, and which may be used to verify its correctness.
#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=2;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) ; } }
I get the following output:
2: integral = 0.5000000000; neval=11 , cond=0.0 3: integral = 0.7500000000; neval=33 , cond=0.0 4: integral = 1.0000000000; neval=105 , cond=0.0 5: integral = 1.2500000000; neval=363 , cond=0.0 6: integral = 1.5000000000; neval=1449 , cond=0.0 7: integral = 1.7500000000; neval=7227 , cond=0.0 8: integral = 2.0000000000; neval=46881 , cond=0.0 9: integral = 2.2500000000; neval=382563, cond=0.0 10: integral = 2.4999999999; neval=3687849, cond=0.0
Full download:
memory.h cubate.c infheap.h quadrate.h cubatetest.c lu.c
Compile and run the test as follows:
g++ -o cubatetest -O -g cubatetest.c cubate.c lu.c ; cubatetest
doc : algorithm : weaknesses : source : test : download
double log1plus(double x),log1minus(double x) ;
log1plus(x) is correct over a
wider range of values than the call log(1+x), which for
instance yields 0 when
#include <math.h>
static double a1=0.49999999193204,a2=0.33333332236384 ,
a3=0.25008221234609,a4=0.20008290182233 ;
double log1plus(double x)
{ if(fabs(x)>0.02) return log(1+x) ;
else return x * ( 1 - x*(a1-x*(a2-x*(a3-x*a4))) ) ;
}
double log1minus(double x)
{ if(fabs(x)>0.02) return log(1-x) ;
else return -x * ( 1 + x*(a1+x*(a2+x*(a3+x*a4))) ) ;
}
doc : source : test : download
double rtlnorm(double x),logrtlnorm(double) ;
rtlnorm(x) returns the area to the right of
x under the pdf of a standard gaussian distribution;
ie. it returns the probability that a standard normal variate will have
Hence
rtlnorm(x) = erfc(x/sqrt(2))/2
logrtlnorm(x) returns log(rtlnorm(x)).
Dependencies: memory.h , log1plus.c
• rtlnorm • ltlnorm • logrtlnorm • logltlnorm • invrtlnorm • invltlnorm
#include "memory.h"
#include <math.h>
static double root = 1.0/sqrt(2.0) ;
static double a1=0.49999999193204,a2=0.33333332236384 ,
a3=0.25008221234609,a4=0.20008290182233 ;
double log1plus(double x),log1minus(double x) ;
double rtlnorm(double x)
{ double t,z=fabs(root*x),ans ;
t = 1 / ( 1 + z/2 ) ;
ans = z*z+1.26551223-t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277)))))))) ;
ans = t * exp(-ans) / 2 ;
if(x>0) return ans ; else return 1-ans ;
}
double ltlnorm(double x)
{ double t,z=fabs(root*x),ans ;
t = 1 / ( 1 + z/2 ) ;
ans = z*z+1.26551223-t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277)))))))) ;
ans = t * exp(-ans) / 2 ;
if(x>0) return 1-ans ; else return ans ;
}
double logrtlnorm(double x)
{ double t,z=fabs(root*x),ans ;
t = 1 / ( 1 + z/2 ) ;
ans = z*z+1.26551223-t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277)))))))) ;
if(x>0) return log(t/2)-ans ; else return log1minus(t*exp(-ans)/2) ;
}
double logltlnorm(double x)
{ double t,z=fabs(root*x),ans ;
t = 1 / ( 1 + z/2 ) ;
ans = z*z+1.26551223-t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277)))))))) ;
if(x>0) return log1minus(t*exp(-ans)/2) ; else return log(t/2)-ans ;
}
/*
* Upper tail quantile for standard normal distribution function.
*
* This function returns an approximation of the inverse cumulative
* standard normal distribution function. I.e., given P, it returns
* an approximation to the X satisfying P = Pr{Z >= X} where Z is a
* random variable from the standard normal distribution.
*
* The algorithm uses a minimax approximation by rational functions
* and the result has a relative error whose absolute value is less
* than 1.15e-9.
*
* Author: Peter John Acklam
* Time-stamp: 2002-06-09 18:45:44 +0200
* E-mail: jacklam@math.uio.no
* WWW URL: http://www.math.uio.no/~jacklam
*
* C implementation adapted from Peter's Perl version
*/
/* Coefficients in rational approximations. */
static const double a[] =
{
-3.969683028665376e+01,
2.209460984245205e+02,
-2.759285104469687e+02,
1.383577518672690e+02,
-3.066479806614716e+01,
2.506628277459239e+00
};
static const double b[] =
{
-5.447609879822406e+01,
1.615858368580409e+02,
-1.556989798598866e+02,
6.680131188771972e+01,
-1.328068155288572e+01
};
static const double c[] =
{
-7.784894002430293e-03,
-3.223964580411365e-01,
-2.400758277161838e+00,
-2.549732539343734e+00,
4.374664141464968e+00,
2.938163982698783e+00
};
static const double d[] =
{
7.784695709041462e-03,
3.224671290700398e-01,
2.445134137142996e+00,
3.754408661907416e+00
};
#define LOW 0.02425
#define HIGH (1-LOW)
double invrtlnorm(double p)
{ double q, r;
if(p<=0||p>=1) throw up("illegal argument %.2f to invrtlnorm",p) ;
if(p<LOW)
{ /* Rational approximation for lower region */
q = sqrt(-2*log(p));
return - (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) ;
}
if(p>HIGH)
{ /* Rational approximation for upper region */
q = sqrt(-2*log(1-p));
return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) ;
}
/* Rational approximation for central region */
q = p - 0.5 ;
r = q * q ;
return - (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q /
(((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1) ;
}
double invltlnorm(double p) { return -invrtlnorm(p) ; }
Store the following as rtlnormtest.c.
#include <math.h>
#include <stdio.h>
double rtlnorm(double x) ;
int main()
{ for(int i=0;i<21;i++)
{ printf("%.4f ",rtlnorm(0.2*i)) ; if(i%10==9) printf("\n") ; }
printf("\n") ;
}
Execution yields
0.5000 0.4207 0.3446 0.2743 0.2119 0.1587 0.1151 0.0808 0.0548 0.0359 0.0228 0.0139 0.0082 0.0047 0.0026 0.0013 0.0007 0.0003 0.0002 0.0001 0.0000
memory.h rtlnorm.c log1plus.c rtlnormtest.c
Compile and run by;
g++ -o rtlnormtest rtlnormtest.c rtlnorm.c log1plus.c ; rtlnormtest
doc : source : test : download
doc : source : test : download
double besseli(int n,double x) ;
computes the modified Bessel function In(x).
Dependency: memory.h
• besseli0 • besseli1 • besseli
#include <math.h> #include "memory.h" static double besseli0(double x) { double y ; x = fabs(x) ; if(x<3.75) { x *= x/(3.75*3.75) ; return 1 + x*( 3.5156229+x*(3.0899424+x*(1.2067492 + x*(0.2659732+x*(0.0360768+x*0.0045813))))); } y = exp(x) / sqrt(x) ; x = 3.75 / x ; return y * ( 0.39894228+x*(0.01328592+x*(0.00225319 + x*( -0.00157565+x*(0.00916281+x*(-0.02057706 + x*(0.02635537+x*(-0.01647633+x*0.00392377)))))))); } static double besseli1(double x) { double y,ax=fabs(x) ; if(ax<3.75) { y = x*x / (3.75*3.75) ; return x*( 0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934 + y*(0.02658733+y*(0.00301532+y*0.00032411)))))) ; } y = exp(ax) / sqrt(ax) ; if(x<0) y = -y ; x = 3.75 / ax ; return y * ( 0.39894228+x*(-0.03988024+x*(-0.00362018 + x*( 0.00163801+x*(-0.01031555+x*(0.02282967 + x*(-0.02895312+x*(0.01787654-x*0.00420059)))))))) ; } double besseli(int n,double x) { int t; double x1,x0,x2,qx,res,y; if(n<0) throw up("n=%d not allowed for besseli",n) ; if(n==0) return besseli0(x) ; else if(n==1) return besseli1(x) ; if(x==0) return 0 ; qx = 2 / fabs(x) ; for(res=x2=0,x1=1,t=2*(n+(int) sqrt(40*n));t>0;t--) { x0 = x2 + t*qx*x1 ; x2 = x1 ; x1 = x0; if(fabs(x1)>1e10) { res /= 1e10 ; x1 /= 1e10 ; x2 /= 1e10 ; } if(t==n) res = x2 ; } y = res * besseli0(x) / x1 ; if(x<0&&(n&1)) return -y ; else return y ; }
The following program testbesseli.c generates alternate entries for Table 9.8 of Abramowitz & Stegun.
#include <math.h>
#include "memory.h"
double besseli(int n,double x) ;
int main()
{ for(double x=0;x<5.00001;x+=0.2)
printf("%.1f %.10f %.10f %.10f\n",
x,exp(-x)*besseli(0,x),exp(-x)*besseli(1,x),
x?besseli(2,x)/(x*x):0.125) ;
}
Correct output:
0.0 1.0000000000 0.0000000000 0.1250000000 0.2 0.8269385470 0.0822831234 0.1254171871 0.4 0.6974021576 0.1367632235 0.1266750199 0.6 0.5993271850 0.1721644178 0.1287924377 0.8 0.5241489250 0.1944986909 0.1318014275 1.0 0.4657595967 0.2079104130 0.1357476666 1.2 0.4197820753 0.2152568578 0.1406914443 1.4 0.3830625169 0.2185075919 0.1467088842 1.6 0.3533150009 0.2190194903 0.1538934958 1.8 0.3288719519 0.2177262794 0.1623580911 2.0 0.3085083232 0.2152692896 0.1722371123 2.2 0.2913173335 0.2120877326 0.1836894254 2.4 0.2766223247 0.2084810884 0.1969016471 2.6 0.2639139986 0.2046522541 0.2120920863 2.8 0.2528055356 0.2007374113 0.2295153956 3.0 0.2430003531 0.1968267133 0.2494680479 3.2 0.2342688282 0.1929786227 0.2722947716 3.4 0.2264313998 0.1892298511 0.2983960992 3.6 0.2193462282 0.1856022491 0.3282372134 3.8 0.2129001300 0.1821075834 0.3623583114 4.0 0.2070019252 0.1787508396 0.4013868437 4.2 0.2015773894 0.1755325255 0.4460516746 4.4 0.1965655607 0.1724502351 0.4971998735 4.6 0.1919159144 0.1694997330 0.5558163299 4.8 0.1875862034 0.1666757060 0.6230467382 5.0 0.1835408131 0.1639722649 0.7002246005
memory.h besseli.c testbesseli.c
Compile and run by;
g++ -o testbesseli testbesseli.c besseli.c ; testbesseli
doc : source : test : download
double ranf(),gaussv(),gammav(double),betav(double,double) ; int rani(int) ; void ranset(int) ;
functions provided : source : ranf test : gaussv test : gammav test : download
| ranf() | returns a random uniform deviate x in the range 0≤x<1. |
| ranset(seed) | initialises the generator to an integer seed of your choice. |
| rani(n) | returns a random integer i in the range 0≤i<n. |
| gaussv() | returns a standard gaussian deviate. |
| gammav(α) | returns a Γ(α,1) deviate. |
| betav(r,s) | returns a Beta(r,s) deviate. |
random.h is a useful header file if you are using these routines or others which depend on them.
The random number generator is the one recommended by NR.
• step • ranf • ranset • gaussv • rani • gammav • betav
#include <math.h> #define IA 16807 #define IM 2147483647 #define IQ 127773 #define IR 2836 #define NDIV (1+(IM-1)/32) static long tabx=0,tab[32],state=123467 ; static int toggle=1 ; static long step() { long k = state / IQ ; state = IA * (state-k*IQ) - IR*k ; return state<0?state+IM:state ; } double ranf() { int i ; if(tabx==0) { for(i=63;i>=0;i--) tab[i&31] = step() ; tabx = tab[0] ; } i = tabx / NDIV ; tabx = tab[i] ; tab[i] = step() ; return tabx / (double) IM; } void ranset(int seed) { if(seed>0) state = 2*seed ; else state = 1 - 2*seed ; if(state>=IM) state = 1 + state%(IM-1) ; tabx = 0 ; toggle = 1 ; } double gaussv() { static double u1,u2,pi=3.141592653589793 ; if((toggle=1-toggle)) return u1 * sin(u2) ; for(u1=0;u1==0;u1=ranf()) ; u1 = sqrt(-2*log(u1)) ; u2 = 2*pi*ranf() ; return u1 * cos(u2) ; } int rani(int n) { return (int)(n*ranf()) ; } double gammav(double alpha) // marsaglia and tsang, 2000 { double d=alpha-1.0/3,c=1/sqrt(9*d),x,v,u ; while(1) { for(v=0;v<=0;) { x = gaussv() ; v = 1 + c*x ; } v *= v*v ; x *= x ; u = ranf() ; if(u<1-0.0331*x*x||log(u)<(x/2+d*(1-v+log(v)))) return d*v ; } return 0 ; } double betav(double r,double s) { double x=gammav(r),y=gammav(s) ; return x/(x+y) ; }
To be stored as ranftest.c.
#include <stdio.h>
#include "random.h"
int main()
{ int i,ind,histo[100] ;
for(i=0;i<100;i++) histo[i] = 0 ;
for(i=0;i<1000000;i++)
{ ind = (int)(100*ranf()) ; if(ind>99) ind = 99 ; histo[ind] += 1 ; }
for(i=0;i<100;i++)
{ printf("%6d",histo[i]) ; if(9==i%10) printf("\n") ; }
}
9666 10068 9914 10092 10091 10066 10086 10085 9909 9965 10155 10086 10053 9845 10032 10144 10141 10001 10140 10045 10102 10041 10004 10134 10126 10090 10078 9948 9949 9898 10104 10000 10108 9919 10002 9907 10047 10161 10112 10105 9910 10030 10053 9847 10063 10006 9849 10179 10022 9992 10086 10027 9817 9887 9936 9907 9922 9934 10058 9974 9976 10016 9960 9885 9860 9732 9804 9871 9989 9936 10084 9873 10058 9905 9998 9960 10104 10224 9868 10082 10155 10078 10022 10034 9845 10083 10027 9969 9960 9996 10154 9864 9841 9964 9957 10008 9950 9978 10117 9895
The following test program gaussvtest.c computes the mean value of x4 for a standard gaussian distribution over samples of increasing size: the true value is 3.
#include <stdio.h>
#include "random.h"
int main()
{ int i,lim ;
double q,u ;
for(u=i=0,lim=10;lim<10000001;lim*=10)
{ for(;i<lim;i++) { q = gaussv() ; q *= q ; u += q*q ; }
printf("%8d obs, E(x*x*x*x) = %.6f\n",lim,u/lim) ;
}
}
Results:
10 obs, E(x*x*x*x) = 8.396751
100 obs, E(x*x*x*x) = 3.862880
1000 obs, E(x*x*x*x) = 2.679729
10000 obs, E(x*x*x*x) = 3.055467
100000 obs, E(x*x*x*x) = 2.984881
1000000 obs, E(x*x*x*x) = 2.996360
10000000 obs, E(x*x*x*x) = 2.997730
Simple test for gammav and betav.
#include "memory.h"
#include <math.h>
#include "random.h"
int main()
{ int i,n=1000 ;
double x,sx,sxx,mu ;
for(sxx=sx=i=0;i<n;i++)
{ x = gammav(1.5) ; sx += x ; sxx += x*x ; }
mu = sx / n ;
printf("%.3f %.3f\n",mu,sxx/n-mu*mu) ;
for(sxx=sx=i=0;i<n;i++)
{ x = betav(3,1) ; sx += x ; sxx += x*x ; }
mu = sx / n ;
printf("%.3f %.3f\n",mu,sxx/n-mu*mu) ;
}
Results:
10 obs, E(x*x*x*x) = 8.396751
100 obs, E(x*x*x*x) = 3.862880
1000 obs, E(x*x*x*x) = 2.679729
10000 obs, E(x*x*x*x) = 3.055467
100000 obs, E(x*x*x*x) = 2.984881
1000000 obs, E(x*x*x*x) = 2.996360
10000000 obs, E(x*x*x*x) = 2.997730
ranf.c random.h memory.h ranftest.c gaussvtest.c
Compile and run:
g++ -O -g -o ranftest ranftest.c ranf.c ; ranftest g++ -O -g -o gaussvtest gaussvtest.c ranf.c ; gaussvtest
functions provided : source : ranf test : gaussv test : gammav test : download
doc : radixsort.h : source : test : download
uiradixsort sorts an array of unsigned ints; ui64radixsort sorts an array of unsigned long long ints. The options available are listed in the header which immediately follows. Where two arrays are provided as arguments, the first is sorted into the second with the original unaltered; where there is a single argument, the sort is in place. radixsort falls back on a shell sort for short arrays. It has not been exhaustively tested. It would be easy to extend to other integer data types.
void uiradixsort(unsigned int *x,unsigned int *y,int nn) ; void uiradixsort(unsigned int *x,int n) ; void uiradixsortdown(unsigned int *x,unsigned int *y,int n) ; void uiradixsortdown(unsigned int *x,int n) ; void ui64radixsort(unsigned long long *x,unsigned long long *y,int nn) ; void ui64radixsort(unsigned long long *x,int n) ; void ui64radixsortdown(unsigned long long *x,unsigned long long *y,int n) ; void ui64radixsortdown(unsigned long long *x,int n) ;
• uisortup • uisortdown • uiradixsort • uiradixsortdown • ui64sortup • ui64sortdown • ui64radixsort • ui64radixsortdown
#include "memory.h" genvector(unsigned int,uivector) ; static void uisortup(unsigned int *u,int n) { shsortup(u,n,unsigned int) ; } static void uisortdown(unsigned int *u,int n) { shsortdown(u,n,unsigned int) ; } void uiradixsort(unsigned int *x,unsigned int *y,int nn) { int i,j , n = abs(nn) ; if(n<150) { if(y==x) { if(nn>=0) uisortup(x,n) ; else uisortdown(x,n) ; return ; } for(i=0;i<n;i++) y[i] = x[i] ; if(nn>=0) uisortup(y,n) ; else uisortdown(y,n) ; return ; } unsigned int *a=x,*b,*aa,*c=uivector(n),k,*u,*v ; int size=sizeof(unsigned int),**count=imatrix(size,256),offs[256],pos ; // we sort from a to b, then between b and aa if(x==y) { b = c ; aa = a ; } else { b = y ; aa = c ; } // first pass through the data makes counts for(i=0;i<n;i++) for(k=x[i],j=0;j<size;j++) count[j][255&(k>>(8*j))] += 1 ; // now an additional pass per byte for(pos=0;pos<size;pos++) { if(pos==0) { u = a ; v = b ; } for(j=i=0;i<256&&j<2;i++) if(count[pos][i]) j += 1 ; if(j<2) continue ; // all values equal in this position if(nn>=0) for(offs[0]=0,i=1;i<256;i++) offs[i] = offs[i-1]+count[pos][i-1] ; else for(offs[255]=0,i=254;i>=0;i--) offs[i] = offs[i+1] + count[pos][i+1] ; for(i=0;i<n;i++) { k = (u[i]>>(8*pos))&255 ; v[offs[k]++] = u[i] ; } if(u==b) { v = b ; u = aa ; } else { u = b ; v = aa ; } } if(u!=y) for(i=0;i<n;i++) y[i] = u[i] ; free(c) ; freeimatrix(count) ; } void uiradixsort(unsigned int *x,int n) { uiradixsort(x,x,n) ; } void uiradixsortdown(unsigned int *x,unsigned int *y,int n) { uiradixsort(x,y,-n) ; } void uiradixsortdown(unsigned int *x,int n) { uiradixsort(x,x,-n) ; } genvector(unsigned long long int,ui64vector) ; static void ui64sortup(unsigned long long *u,int n) { shsortup(u,n,unsigned long long) ; } static void ui64sortdown(unsigned long long *u,int n) { shsortdown(u,n,unsigned long long) ; } void ui64radixsort(unsigned long long *x,unsigned long long *y,int nn) { int i,j , n = abs(nn) ; if(n<150) { if(y==x) { if(nn>=0) ui64sortup(x,n) ; else ui64sortdown(x,n) ; return ; } for(i=0;i<n;i++) y[i] = x[i] ; if(nn>=0) ui64sortup(y,n) ; else ui64sortdown(y,n) ; return ; } unsigned long long *a=x,*b,*aa,*c=ui64vector(n),k,*u,*v ; int size=sizeof(unsigned long long),**count=imatrix(size,256),offs[256],pos ; // we sort from a to b, then between b and aa if(x==y) { b = c ; aa = a ; } else { b = y ; aa = c ; } // first pass through the data makes counts for(i=0;i<n;i++) for(k=x[i],j=0;j<size;j++) count[j][255&(k>>(8*j))] += 1 ; // now an additional pass per byte for(pos=0;pos<size;pos++) { if(pos==0) { u = a ; v = b ; } for(j=i=0;i<256&&j<2;i++) if(count[pos][i]) j += 1 ; if(j<2) continue ; // all values equal in this position if(nn>=0) for(offs[0]=0,i=1;i<256;i++) offs[i] = offs[i-1]+count[pos][i-1] ; else for(offs[255]=0,i=254;i>=0;i--) offs[i] = offs[i+1] + count[pos][i+1] ; for(i=0;i<n;i++) { k = (u[i]>>(8*pos))&255 ; v[offs[k]++] = u[i] ; } if(u==b) { v = b ; u = aa ; } else { u = b ; v = aa ; } } if(u!=y) for(i=0;i<n;i++) y[i] = u[i] ; free(c) ; freeimatrix(count) ; } void ui64radixsort(unsigned long long *x,int n) { ui64radixsort(x,x,n) ; } void ui64radixsortdown(unsigned long long *x,unsigned long long *y,int n) { ui64radixsort(x,y,-n) ; } void ui64radixsortdown(unsigned long long *x,int n) { ui64radixsort(x,x,-n) ; }
#include "memory.h"
#include "random.h"
#include "radixsort.h"
genvector(unsigned int,uivector) ;
int main()
{ int n,i ;
unsigned int *a ;
for(n=10;n<=10000;n*=10)
{ a = uivector(n) ;
for(i=0;i<n;i++) a[i] = 10203 * rani(10203) ;
uiradixsort(a,n) ;
for(i=0;i<n-1;i++) if(a[i]>a[i+1]) throw up("error: %d|%d",a[i],a[i+1]) ;
free(a) ;
}
}
There will be no output if all is well.
ranf.c random.h memory.h radixsort.h radixsort.c radixsorttest.c
Compile and run:
g++ -O -g -o radixsorttest radixsorttest.c radixsort.c ranf.c ; radixsorttest
doc : radixsort.h : source : test : download
doc : source : test : download
infheap.h manages a binary max heap of unlimited size. It should be self-explanatory. It is used by cubate.
• infheapitem • push • print • recurse • pop • release • infheap • size
struct infheapitem { void *data ; double key ; infheapitem *l,*r ; int nl,nr ; infheapitem() {} ; infheapitem(void *d,double k) { data = d ; key = k ; l = r = 0 ; nl = nr = 0 ; } } ; struct infheap { private : infheapitem *h ; void push(void *d,double k,infheapitem **loc) { void *dd ; double kk ; if(loc[0]==0) loc[0] = new infheapitem(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 infheapitem(d,k) ; loc[0]->nl = 1 ; } else if(loc[0]->nr==0) { loc[0]->r = new infheapitem(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(infheapitem *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(infheapitem *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(infheapitem **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(infheapitem **loc) { if(loc[0]->nl) release(&loc[0]->l) ; if(loc[0]->nr) release(&loc[0]->r) ; delete(loc[0]) ; loc[0] = 0 ; } public : infheap() { 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 "random.h" #include "infheap.h" void printint(void *x) { printf("%d ",((int *)x)[0]) ; } int main() { int t,n,*d ; infheap 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
ranf.c random.h memory.h infheap.h infheaptest.c
Compile and run:
g++ -O -g -o infheaptest infheaptest.c ranf.c ; infheaptest