Colin Champion
function | description |
functions provided : doc : source : test
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 doesnt 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 isnt 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.
Complex numbers: memory.h defines a complex data type (but not complex arithmetic, for which use unreal.h). This is so that it can provide allocation for the complex types. You could of course be a conformist and use <complex>.
Sorting: there are generic ascending and descending shell sorts, intended for data types of the users 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 youve 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 • settable • unset • double • print • free • cjcupalloc • cjcuprealloc • freename • swap • charvector • xivector • isortup • xysort • xisort • 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("%.1f") ; char *ptr = cjcbuf + cjcind ; cjcbuf[cjcind] = '(' ; if(cjcint) cjcind += 2 + sprintf(cjcbuf+cjcind+1,cjcbuf,(int)x) ; else cjcind += 2 + sprintf(cjcbuf+cjcind+1,cjcbuf,x) ; cjcbuf[cjcind-1] = ',' ; if(cjcint) cjcind += 2 + sprintf(cjcbuf+cjcind,cjcbuf,(int)y) ; else cjcind += 2 + sprintf(cjcbuf+cjcind,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("%.1f") ; char *ptr = cjcbuf + cjcind ; if(cjcint) cjcind += 1 + sprintf(cjcbuf+cjcind,cjcbuf,(int)x) ; else cjcind += 1 + sprintf(cjcbuf+cjcind,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) ; } } ; struct xi { double x ; int i ; xi() { x = i = 0 ; } xi(double a,int b) { x = a ; i = 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)) ; } \ #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(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 ; \ } \ } /* ----------------------------- integer sort ------------------------------- */ static void isortup(int *x,int n) { int i,j,inc,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 ; } } 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 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> #include "unreal.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) ; } }
functions provided : doc : source : test
unreal.h defines the usual set of arithmetic operations. You need to include both memory.h and <math.h> before unreal.h.
unreal.h defines the constant UNREAL_H as a sign of its presence.
• arg • norm • real • imag • conj • abs • exp • sqrt
#ifndef UNREAL_H #define UNREAL_H // addition inline complex operator+(complex z1,complex z2) { return complex(z1.re+z2.re,z1.im+z2.im) ; } inline complex operator +=(complex &z1,complex z2) { return z1 = complex(z1.re+z2.re,z1.im+z2.im) ; } inline complex operator+(complex z1,double z2) { return complex(z1.re+z2,z1.im) ; } inline complex operator+=(complex &z1,double z2) { return z1 = complex(z1.re+z2,z1.im) ; } inline complex operator+(double z1,complex z2) { return complex(z1+z2.re,z2.im) ; } // subtraction inline complex operator-(complex z1,complex z2) { return complex(z1.re-z2.re,z1.im-z2.im) ; } inline complex operator -=(complex &z1,complex z2) { return z1 = complex(z1.re-z2.re,z1.im-z2.im) ; } inline complex operator-(complex z1,double z2) { return complex(z1.re-z2,z1.im) ; } inline complex operator-=(complex &z1,double z2) { return z1 = complex(z1.re-z2,z1.im) ; } inline complex operator-(double z1,complex z2) { return complex(z1-z2.re,-z2.im) ; } // multiplication inline complex operator*(complex z1,complex z2) { return complex(z1.re*z2.re-z1.im*z2.im,z1.re*z2.im+z1.im*z2.re) ; } inline complex operator*=(complex &z1,complex z2) { return z1 = z1 * z2 ; } inline complex operator*(complex z1,double z2) { return complex(z1.re*z2,z1.im*z2) ; } inline complex operator*=(complex &z1,double z2) { return z1 = complex(z1.re*z2,z1.im*z2) ; } inline complex operator*(double z1,complex z2) { return complex(z1*z2.re,z1*z2.im) ; } // simple functions inline double arg(complex z) { return atan2(z.im,z.re) ; } inline double norm(complex z) { return z.re*z.re+z.im*z.im ; } inline double real(complex z) { return z.re ; } inline double imag(complex z) { return z.im ; } inline complex conj(complex z) { return complex(z.re,-z.im) ; } inline complex operator -(complex z) { return complex(-z.re,-z.im) ; } // abs inline double abs(complex z) { double q,qr=fabs(z.re),qi=fabs(z.im) ; if(qr>qi) { q = z.im/z.re ; return qr*sqrt(1+q*q) ; } else { if(qi==0) return 0 ; q = z.re/z.im ; return qi*sqrt(1+q*q) ; } } // powers inline complex exp(complex z) { double u=exp(z.re),v=z.im ; return complex(u*cos(v),u*sin(v)) ; } inline complex sqrt(complex z) { double r=sqrt(abs(z)),phi=arg(z)/2 ; return complex(r*cos(phi),r*sin(phi)) ; } // various forms of division inline complex operator/(complex z,double x) { return z * (1/x) ; } inline complex operator/=(complex &z,double x) { return z *= 1/x ; } inline complex operator/(complex z1,complex z2) { double q,qr=fabs(z2.re),qi=fabs(z2.im) ; if(qr>=qi) { q = z2.im / z2.re ; return complex(z1.re+z1.im*q,z1.im-z1.re*q) / (z2.re+q*z2.im) ; } else { q = z2.re / z2.im ; return complex(z1.re*q+z1.im,z1.im*q-z1.re) / (z2.re*q+z2.im) ; } } inline complex operator/(double x,complex z) { double q,qr=fabs(z.re),qi=fabs(z.im) ; if(qr>=qi) { q = z.im / z.re ; return complex(x,-x*q) / (z.re+q*z.im) ; } else { q = z.re / z.im ; return complex(x*q,-x) / (z.re*q+z.im) ; } } inline complex operator/=(complex &z1,complex z2) { return z1 = z1 / z2 ; } // test for equality/inequality with complex/double/int inline bool operator==(complex z1,complex z2) { return z1.re==z2.re && z1.im==z2.im ; } inline bool operator==(complex z,double y) { return z.re==y && z.im==0 ; } inline bool operator==(double y,complex z) { return z.re==y && z.im==0 ; } inline bool operator==(complex z,int y) { return z.re==y && z.im==0 ; } inline bool operator==(int y,complex z) { return z.re==y && z.im==0 ; } inline bool operator!=(complex z1,complex z2) { return z1.re!=z2.re || z1.im!=z2.im ; } inline bool operator!=(complex z,double y) { return z.re!=y || z.im!=0 ; } inline bool operator!=(double y,complex z) { return z.re!=y || z.im!=0 ; } inline bool operator!=(complex z,int y) { return z.re!=y || z.im!=0 ; } inline bool operator!=(int y,complex z) { return z.re!=y || z.im!=0 ; } #endif
The following program testunreal.c tests a few of the more complicated operators:
#include "memory.h" #include <math.h> #include "unreal.h" #include "random.h" #include <assert.h> double ranf() ; int main() { int t ; complex z,z1,z2 ; double q ; for(t=0;t<100;t++) { z = complex(ranf()-0.5,ranf()-0.5) ; q = abs(z) ; assert(norm(z)-q*q<1e-10) ; z1 = 1/z ; assert(abs(1-z*z1)<1e-10) ; z1 = complex(ranf()-0.5,ranf()-0.5) ; z1 *= z/z1 ; assert(abs(z1-z)<1e-10) ; z1 = z2 = complex(ranf()-0.5,ranf()-0.5) ; z1 /= z ; assert(abs(z1*z-z2)<1e-10) ; z1 = sqrt(z) ; assert(abs(z1*z1-z)<1e-10) ; } }
Compile and execute using
g++ -o testunreal testunreal.c ranf.c ; testunreal
Silent execution indicates that all is well; an assertion failure will be thrown in the case of error.
ij.h defines a structure for integer ordered pairs. memory.h needs to be included prior to it in your source code. There are two constructors
ij()
which returns a zero ij and
ij(i,j)
which returns an ij with the given components.
It defines two functions:
void swap(ij x,ij y) ; void ijsort(ij *x,int n) ;
the second of which sorts the array x in place into ascending order of its i component. (The algorithm is a shell sort with a relatively good choice of increment.)
Four vector allocators are also defined.
ij *ijvector(int n) ;
returns an empty array of n ijs.
ij *ijvector(int *x,int n) ;
returns an n-long array of ijs whose kth element has i component x[k] and j component k.
ij *ijvector(long long *x,int n) ;
does the same for long longs.
ij *ijvector(ij *x,int n) ;
reallocates x to length n and returns a pointer to it.
• ij • fprint • print • swap • ijsort • ijvector
#ifndef IJ_H #define IJ_H struct ij { long long i,j ; ij() { i = j = 0 ; } ij(long long a,long long b) { i = a ; j = b ; } int fprint(FILE *f) { return fprintf(f,"(%lld,%lld)",i,j) ; } int print() { return fprint(stdout) ; } } ; static inline void swap(ij &a,ij &b) { ij c=a ; a = b ; b = c ; } static void ijsort(ij *u,int n) { shellsortup(u,n,ij,i) ; } /* allocate vectors of ijs */ static ij *ijvector(int n) { return (ij *) cjcalloc(n,sizeof(ij)) ; } static ij *ijvector(int *x,int n) { ij *u=ijvector(n) ; for(int i=0;i<n;i++) u[i] = ij(x[i],i) ; return u ; } static ij *ijvector(long long *x,int n) { ij *u=ijvector(n) ; for(int i=0;i<n;i++) u[i] = ij(x[i],i) ; return u ; } static ij *ijvector(ij *a,int n) { return (ij *) cjcrealloc(a,n*sizeof(ij)) ; } #endif
The test program ijtest.c is very simple.
#include "memory.h" #include "ij.h" int main() { ij x ; int i ; for(i=0;i<10;i++) { x = ij(1ll<<(3*i),1ll<<5*i) ; printf("(%lld,%lld)\n",x.i,x.j) ; } }
Correct output:
(1,1) (8,32) (64,1024) (512,32768) (4096,1048576) (32768,33554432) (262144,1073741824) (2097152,34359738368) (16777216,1099511627776) (134217728,35184372088832)
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 ;
}
cholesky inverts a positive definite symmetric matrix, returning its determinant. logcholesky returns its log determinant. The prototypes are:
double cholesky(double **a,double **b,int n) ; double logcholesky(double **a,double **b,int n) ;
where a is the n×n matrix to invert and b is the n×n matrix to receive the result. b may overwrite a.
Cholesky decomposition works by first finding a lower-diagonal square root of a, by inverting the square root, and by then squaring it back up. The subroutines cholsqrt, cholinv and cholsqr perform the three components of the operation, and are available to be called separately.
If you just want to test for positive-definiteness, or to find the determinant, cholsqrt does all you need.
Note that cholesky looks only at those a[i][j] for which i≤j.
• cholesky • logcholesky • cholsqrt • logcholsqrt • cholinv • cholsqr
#include "memory.h" #include <math.h> double cholsqrt(double **,double **,int),logcholsqrt(double **,double **,int) ; void cholinv(double **,double **,int) ; void cholsqr(double **,double **,int) ; /* ------------------------------------------------------------------- cholesky inverts a positive-definite symmetric matrix. */ double cholesky(double **a,double **b,int n) { double qdet = cholsqrt(a,b,n) ; cholinv(b,b,n) ; cholsqr(b,b,n) ; return qdet ; } double logcholesky(double **a,double **b,int n) { double logqdet = logcholsqrt(a,b,n) ; cholinv(b,b,n) ; cholsqr(b,b,n) ; return logqdet ; } /* ------------------------------------------------------------------- cholsqrt looks at a[i][j] only for i<=j, and constructs b with b[i][j]=0 for i>j. it may be run in place if you wish, ie. with a=b. */ double cholsqrt(double **a,double **b,int n) { int i,j,k ; double q,qdet=1 ; for(j=0;j<n;j++) { for(q=a[j][j],i=0;i<j;i++) q -= b[i][j] * b[i][j] ; if(q<=0) return 0.0 ; qdet *= q ; b[j][j] = sqrt(q) ; for(k=j+1;k<n;k++) { for(q=i=0;i<j;i++) q += b[i][j] * b[i][k] ; b[j][k] = ( a[j][k] - q ) / b[j][j] ; } } for(j=0;j<n;j++) for(k=0;k<j;k++) b[j][k] = 0 ; return qdet ; } double logcholsqrt(double **a,double **b,int n) { int i,j,k ; double q,logqdet=0 ; for(j=0;j<n;j++) { for(q=a[j][j],i=0;i<j;i++) q -= b[i][j] * b[i][j] ; if(q<=0) throw up("cholesky doesn\'t think your input matrix is pd") ; logqdet += log(q) ; b[j][j] = sqrt(q) ; for(k=j+1;k<n;k++) { for(q=i=0;i<j;i++) q += b[i][j] * b[i][k] ; b[j][k] = ( a[j][k] - q ) / b[j][j] ; } } for(j=0;j<n;j++) for(k=0;k<j;k++) b[j][k] = 0 ; return logqdet ; } /* ------------------------------------------------------------------- cholinv looks at a[i][j] only for i<=j, and constructs b with b[i][j] zero for i<j. it may be run in place if you wish, ie. with a=b. */ void cholinv(double **a,double **b,int n) { int i,j,k ; double q ; for(j=0;j<n;j++) { b[j][j] = 1 / a[j][j] ; for(i=j-1;i>=0;i--) { for(q=0,k=i;k<j;k++) q += b[k][i] * a[k][j] ; b[j][i] = - q * b[j][j] ; } } for(j=0;j<n;j++) for(k=0;k<j;k++) b[k][j] = 0 ; } /* ------------------------------------------------------------------- cholsqr looks at a[i][j] only for i>=j, and constructs the complete matrix b. it may be run in place if you wish, ie. with a=b. */ void cholsqr(double **a,double **b,int n) { int i,j,k ; double q ; for(i=0;i<n;i++) for(j=i;j<n;j++) { for(q=0,k=j;k<n;k++) q += a[k][i] * a[k][j] ; b[i][j] = q ; } for(j=0;j<n;j++) for(k=0;k<j;k++) b[j][k] = b[k][j] ; }
void qr(double **a,int m,int n,double **Q,double **R) ;
where
See Wikipedia for the definition of QR decomposition. The subroutine uses the Householder method whose cost is quartic in the matrix dimension.
#include "memory.h" #include <math.h> void qr(double **a,int m,int n,double **Q,double **R) { int i,j,k,t ; double alpha,s,*x=vector(m),**aa,**z=matrix(m,m),**q=matrix(m,m) ; if(m<n) throw up("qr needs m>=n (m=%d, n=%d)",m,n) ; for(aa=a,k=0;k<n&&k<m-1;aa=R,k++) { for(i=0;i<m;i++) for(j=0;j<n;j++) z[i][j] = 0 ; for(i=0;i<k;i++) z[i][i] = 1 ; for(i=k;i<m;i++) for(j=k;j<n;j++) z[i][j] = aa[i][j] ; for(alpha=i=0;i<m;i++) { x[i] = z[i][k] ; alpha += x[i]*x[i] ; } if(a[k][k]>0) x[k] -= sqrt(alpha) ; else x[k] += sqrt(alpha) ; for(s=i=0;i<m;i++) s += x[i] * x[i] ; s = - 2/s ; for(i=0;i<m;i++) for(j=0;j<m;j++) Q[i][j] = s * x[i] * x[j] ; for(i=0;i<m;i++) Q[i][i] += 1 ; for(i=0;i<m;i++) for(j=0;j<n;j++) { for(s=t=0;t<m;t++) s += Q[i][t] * z[t][j] ; R[i][j] = s ; } if(k==0) for(i=0;i<m;i++) for(j=0;j<m;j++) q[i][j] = Q[i][j] ; else { for(i=0;i<m;i++) for(j=0;j<m;j++) { for(s=t=0;t<m;t++) s += q[i][t] * Q[t][j] ; z[i][j] = s ; } for(i=0;i<m;i++) for(j=0;j<m;j++) q[i][j] = z[i][j] ; } } for(i=0;i<m;i++) for(j=0;j<m;j++) Q[i][j] = q[i][j] ; for(i=0;i<m;i++) for(j=0;j<n;j++) { for(s=t=0;t<m;t++) s += Q[t][i] * a[t][j] ; R[i][j] = s ; } free(x) ; freematrix(z,q) ; }
To be stored as qrtest.c.
#include "memory.h"
void qr(double **a,int m,int n,double **Q,double **R) ;
int main(int argc,char **argv)
{ int i,j,m=5,n=3 ;
double in[5][3] =
{ { 12, -51, 4} ,
{ 6, 167, -68} ,
{ -4, 24, -41} ,
{ -1, 1, 0} ,
{ 2, 0, 3} ,
} ;
double **x=matrix(m,n),**Q=matrix(m,m),**R=matrix(m,n) ;
for(i=0;i<m;i++) for(j=0;j<n;j++) x[i][j] = in[i][j] ;
qr(x,m,n,Q,R) ;
for(i=0;i<m;i++)
{ for(j=0;j<m;j++) printf("%8.3f",Q[i][j]) ; printf(" |") ;
for(j=0;j<n;j++) printf("%8.3f",R[i][j]) ; printf("\n") ;
}
freematrix(x,Q,R) ;
return 0;
}
The correct output is:
0.846 -0.391 0.343 0.082 0.078 | 14.177 20.667 -13.402 0.423 0.904 -0.029 0.026 0.045 | 0.000 175.043 -70.080 -0.282 0.170 0.933 -0.047 -0.137 | -0.000 0.000 -35.202 -0.071 0.014 -0.001 0.980 -0.184 | -0.000 -0.000 -0.000 0.141 -0.017 -0.106 -0.171 -0.969 | 0.000 0.000 -0.000
Suppose that we fit a quadratic to the three points (1,a), (0,b), (1,c). Then
The inverse interpolate could of course take either of two values; the one returned is the one closer to the midpoint of the range over x (which is the one wanted for interpolation purposes).
The coefficients returned in co are respectively the constant, the coefficient of x, and the coefficient of x2.
The prototypes are included in algebra.h.
quadintcoefs: the calls are
void quadintcoefs(double x0,double x1,double x2,double *c) ; void quadintcoefs (double x0,double x1,double x2,double *c,double *lhs,double *rhs) ;
If a quadratic passes through (x0,y0), (x1,y1), and (x2,y2), then its integral from x0 to x1 will be
c0y0 + c1y1 + c2y2
The first call to quadintcoefs returns in c the coefficients corresponding to the given xi. The second does so unless the pointer c is 0.
The second call also optionally returns coefficients for the integrals from x0 to (x0+x1)/2 and from (x0+x1)/2 to x1.
There are companion routines in which the double arguments a, b, and c are replaced by xys A, B, and C specifying any three points through which the quadratic passes. And there are additional routines in which the quadratic is specified by its coefficients as returned by quadcoefs: in the case of invquadinterp, the return value is an xy (actually a pair of x-values) because there is no way of telling which ordinate the user prefers.
There is a further call to quadcoefs in which the quadratic is specified by its value y and gradient dy at ordinate x and by an additional point A.
There are also two linear interpolation routines which interpolate a linear function passing through the xy points A and B:
For cubics I provide cubcoefs and cubintcoefs, the latter integrating from x1 to x2.
• lininterp • invlininterp • quadinterp • quadcoefs • quadslope • invquadinterp • quadreach • quadmax
static double lininterp(double x,xy A,xy B) { return A.y + (x-A.x)*(B.y-A.y)/(B.x-A.x) ; } static double invlininterp(double y,xy A,xy B) { return A.x + (y-A.y)*(B.x-A.x)/(B.y-A.y) ; } /* ------------------------------ quadinterp -------------------------------- */ // if a parabola takes values a,b,c at -1,0,1, then quadinterp returns // its value at x static double quadinterp(double x,double a,double b,double c) { return ( x*(x-1)*a - 2*(x+1)*(x-1)*b + x*(x+1)*c ) / 2 ; } static double quadinterp(double x,xy A,xy B,xy C) { double y,lam,mindif,c[2],d[2] ; int ind; if(A.x==B.x||B.x==C.x||A.x==C.x) throw up("two of the ordinates " "for quadinterp are equal: %.3f, %.3f, %.3f",A.x,B.x,C.x) ; ind = 0 ; mindif = fabs(x-A.x) ; y = A.y ; if(fabs(x-B.x)<mindif) { ind = 1 ; mindif = fabs(x-B.x) ; y = B.y ; } if(fabs(x-C.x)<mindif) { ind = 2 ; mindif = fabs(x-C.x) ; y = C.y ; } lam = (B.y-A.y) / (A.x-B.x) ; d[0] = (B.x-x) * lam ; c[0] = (A.x-x) * lam ; lam = (C.y-B.y) / (B.x-C.x) ; d[1] = (C.x-x) * lam ; c[1] = (B.x-x) * lam ; if(ind==0) y += c[0] ; else { ind -= 1 ; y += d[ind] ; } lam = (c[1]-d[0]) / (A.x-C.x) ; if(ind==0) y += (A.x-x)*lam ; else y += (C.x-x)*lam ; return y ; } static double quadinterp(double x,double *co) { return co[0] + x*(co[1]+x*co[2]) ; } /* ------------------------------- quadcoefs -------------------------------- */ static void quadcoefs(xy A,xy B,xy C,double *co) { int i,j ; double p[] = { A.x*B.x*C.x , -(A.x*B.x+A.x*C.x+B.x*C.x) , A.x+B.x+C.x } ; double q,u ; if(A.x==B.x||B.x==C.x||A.x==C.x) throw up("two of the ordinates " "for quadcoefs are equal: %.3f, %.3f, %.3f",A.x,B.x,C.x) ; for(q=3,j=2;j>0;j--) q = q*A.x - j*p[j] ; q = A.y / q ; for(u=1,j=2;j>=0;j--) { co[j] = u*q ; u = u*A.x - p[j] ; } for(q=3,j=2;j>0;j--) q = q*B.x - j*p[j] ; q = B.y / q ; for(u=1,j=2;j>=0;j--) { co[j] += u*q ; u = u*B.x - p[j] ; } for(q=3,j=2;j>0;j--) q = q*C.x - j*p[j] ; q = C.y / q ; for(u=1,j=2;j>=0;j--) { co[j] += u*q ; u = u*C.x - p[j] ; } } // linear embedded in a quadratic static void quadcoefs(xy A,xy B,double *co) { if(A.x==B.x) throw up("the ordinates for quadcoefs are equal: %.3e",A.x) ; co[2] = 0 ; co[1] = (B.y-A.y) / (B.x-A.x) ; co[0] = A.y - co[1]*A.x ; } static void quadcoefs(double a,double b,double c,double *co) { return quadcoefs(xy(-1,a),xy(0,b),xy(1,c),co) ; } // value at one point and value and gradient at a second static void quadcoefs(double x,double y,double dy,xy A,double *co) { if(A.x==x) throw up("equal ordinates supplied to quadcoefs: %.3f, %.3f",A.x,x) ; A = xy(A.x-x,A.y-y) ; co[2] = (A.y-dy*A.x) / (A.x*A.x) ; co[1] = dy - 2*x*co[2] ; co[0] = y - dy*x + co[2]*x*x ; } // two points and gradient at a third ordinate - not well conditioned static void quadcoefs(xy A,xy B,double x,double dy,double *co) { if(A.x==B.x||A.x==x||B.x==x) throw up("two of the ordinates " "for quadcoefs are equal: %.3f, %.3f, %.3f",A.x,B.x,x) ; co[2] = ( (B.y-A.y)/(B.x-A.x) ) ; co[2] = ( (B.y-A.y)/(B.x-A.x) - dy ) / ( (A.x-x)+(B.x-x) ) ; co[1] = dy - 2*co[2]*x ; co[0] = A.y - co[1]*A.x - co[2]*A.x*A.x ; } /* ------------------------------- quadslope -------------------------------- */ // if a parabola takes values a,b,c at -1,0,1, then quadslope returns // its gradient at x static double quadslope(double x,double a,double b,double c) { return ( (2*x-1)*a - 4*x*b + (2*x+1)*c ) / 2 ; } static double quadslope(double x,double *co) { return co[1]+2*x*co[2] ; } static double quadslope(double x,xy A,xy B,xy C) { double co[3] ; quadcoefs(A,B,C,co) ; return quadslope(x,co) ; } /* -------------------------- invquadinterp --------------------------------- */ // if a parabola takes values a,b,c at -1,0,1, then invquadinterp // returns the argument x such the parabola's value at x is y // (and given that there are two such, it chooses the one closer to 0) static double invquadinterp(double y,double a,double b,double c) { double t=b-y,u=(a-c)/4,v=(a-2*b+c)/2,q ; if(v==0) return (2*y-(c+a))/(c-a) ; q = u*u - t*v ; if(q<0) throw up("parabola through (%.2f,%.2f,%.2f) doesn\'t reach %.2f",a,b,c,y) ; q = sqrt(q) ; if(u>0) q = u + q ; else q = u - q ; if(q*q<fabs(t*v)) return q/v ; else return t/q ; } static double invquadinterp(double y,xy A,xy B,xy C) { double co[3],q,t,u,v,xmin=A.x,xmax=A.x ; quadcoefs(A,B,C,co) ; t = co[0] - y ; u = -co[1]/2 ; v = co[2] ; if(v==0) return -t/co[1] ; q = u*u - t*v ; if(q<0) throw up("parabola doesn\'t reach %.2f",y) ; q = sqrt(q) ; if(u>0) q = u + q ; else q = u - q ; if(B.x<xmin) xmin = B.x ; else if(B.x>xmax) xmax = B.x ; if(C.x<xmin) xmin = C.x ; else if(C.x>xmax) xmax = C.x ; u = q * v * (xmin+xmax) / 2 ; if(fabs(q*q-u)<fabs(t*v-u)) return q/v ; else return t/q ; } static xy invquadinterp(double y,double *co) { double q,t,u,v ; if(co[2]==0) { if(co[1]==0) throw up("invquadinterp called on constant quadratic %.2e",co[0]) ; q = (y-co[0]) / co[1] ; return xy(q,q) ; } t = co[0] - y ; u = -co[1]/2 ; v = co[2] ; if(v==0) { v = -t/co[1] ; return xy(v,v) ; } q = u*u - t*v ; if(q<0) throw up("parabola doesn\'t reach %.2e",y) ; q = sqrt(q) ; if(u>0) q = u + q ; else q = u - q ; if(q/v<t/q) return xy(q/v,t/q) ; else return xy(t/q,q/v) ; } static double invquadinterp(double y,double *co,double flag) { xy Z = invquadinterp(y,co) ; return (flag>=0)?Z.y:Z.x ; } /* -------------------------------- quadreach ------------------------------- */ // if a parabola takes values a,b,c at -1,0,1, then quadreach // returns 1 if it ever reaches the value y, and 0 otherwise static int quadreach(double y,double *coef) { if(coef[2]==0) return (coef[1]!=0) ; else return coef[1]*coef[1] >= 4*(coef[0]-y)*coef[2] ; } static int quadreach(double y,xy A,xy B,xy C) { double coef[3] ; quadcoefs(A,B,C,coef) ; return quadreach(y,coef) ; } static int quadreach(double y,double a,double b,double c) { return quadreach(y,xy(-1,a),xy(0,b),xy(1,c)) ; } /* --------------------------------- quadmax -------------------------------- */ // if a parabola takes values a,b,c at -1,0,1, then quadmax returns // the (x,y) values at its turning point (which may be a min or a max) static xy quadmax(double a,double b,double c) { double u=(a-c)/4,v=(a+c-2*b)/2 ; return xy(u/v,b-u*u/v) ; } static xy quadmax(xy A,xy B,xy C) { double co[3],t,u,v,x ; quadcoefs(A,B,C,co) ; t = co[0] ; u = -co[1]/2 ; v = co[2] ; x = u/v ; return xy(x,t+x*(v*x-2*u)) ; } static xy quadmax(double *co) { double t,u,v,x ; if(co[2]==0) throw up("quadmax called for a linear quadratic") ; t = co[0] ; u = -co[1]/2 ; v = co[2] ; x = u/v ; return xy(x,t+x*(v*x-2*u)) ; }
The following program should be stored as testquadinterp.c.
#include "memory.h" #include <math.h> #include "quadinterp.h" #include "random.h" double a=fabs(gaussv()),b=gaussv(),c=gaussv() ; double f(double x) { return c+x*(b+x*a) ; } int main() { int i ; double x,y,co[3],eps=1e-6 ; xy A=xy(0.1,f),B=xy(0.4,f),C=xy(0.9,f),U,V ; quadcoefs(f(-1),f(0),f(1),co) ; printf("quadcoefs: %.1e %.1e %.1e\n", fabs(co[0]-c),fabs(co[1]-b),fabs(co[2]-a)) ; quadcoefs(B,C,A,co) ; printf("quadcoefs: %.1e %.1e %.1e\n", fabs(co[0]-c),fabs(co[1]-b),fabs(co[2]-a)) ; U = quadmax(f(-1),f(0),f(1)) ; V = quadmax(A,B,C) ; printf("quadmax : %.1e %.1e %.1e %.1e\n", fabs(U.x-V.x),fabs(U.y-V.y),fabs(U.y-f(U.x)),fabs(V.y-f(V.x))) ; for(i=-10;i<=10;i++) if(f(0.1*i)<U.y||f(0.1*i)<V.y) throw up("error") ; printf("quadinterp: ") ; for(i=-8;i<8;i+=2) printf("%.1e ",fabs(f(0.1*i)-quadinterp(0.1*i,f(-1),f(0),f(1)))) ; printf("\nquadinterp: ") ; for(i=0;i<8;i++) printf("%.1e ",fabs(f(0.1*i)-quadinterp(0.1*i,A,B,C))) ; printf("\n") ; printf("invquadinterp: ") ; for(i=0;i<8;i++) { y = U.y + 1e-8 + 0.1 * i ; printf("%.1e ",fabs(y-f(invquadinterp(y,f(-1),f(0),f(1))))) ; } printf("\ninvquadinterp: ") ; for(i=0;i<8;i++) { y = V.y + 1e-8 + 0.1 * i ; printf("%.1e ",fabs(y-f(invquadinterp(y,A,B,C)))) ; } printf("\n") ; printf("quadreach: ") ; for(i=-4;i<8;i++) { y = U.y + 1e-8 + 0.1*i ; printf("%d ",quadreach(y,f(-1),f(0),f(1))!=(i>=0)) ; } printf("\nquadreach: ") ; for(i=-4;i<8;i++) { y = V.y + 1e-8 + 0.1*i ; printf("%d ",quadreach(y,A,B,C)!=(i>=0)) ; } printf("\n") ; printf("quadslope: ") ; for(i=-8;i<8;i+=2) { x = 0.1*i ; printf("%.1e ",fabs(f(x+eps)-f(x)-eps*quadslope(x,f(-1),f(0),f(1)))) ; } printf("\nquadslope: ") ; for(i=-8;i<8;i+=2) { x = 0.1*i ; printf("%.1e ",fabs(f(x+eps)-f(x)-eps*quadslope(x,A,B,C))) ; } printf("\nquadcoefs: ") ; { quadcoefs(A.x,A.y,quadslope(A.x,A,B,C),B,co) ; printf("%.1e %.1e %.1e\n",fabs(A.y-quadinterp(A.x,co)), fabs(B.y-quadinterp(B.x,co)),fabs(C.y-quadinterp(C.x,co))) ; } }
Compile and execute as
g++ -o testquadinterp testquadinterp.c quadinterp.c ranf.c ; testquadinterp
and obtain the output
quadcoefs: 0.0e+00 1.1e-16 0.0e+00 quadcoefs: 5.6e-16 3.2e-15 3.1e-15 quadmax : 3.1e-15 2.4e-15 0.0e+00 2.3e-15 quadinterp: 0.0e+00 0.0e+00 2.2e-16 1.1e-16 0.0e+00 1.1e-16 4.4e-16 2.2e-16 quadinterp: 0.0e+00 0.0e+00 1.1e-16 1.1e-16 0.0e+00 2.2e-16 0.0e+00 0.0e+00 invquadinterp: 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 2.2e-16 2.2e-16 0.0e+00 invquadinterp: 2.3e-15 3.3e-16 1.1e-15 1.6e-15 1.8e-15 2.0e-15 2.0e-15 2.0e-15 quadreach: 0 0 0 0 0 0 0 0 0 0 0 0 quadreach: 0 0 0 0 0 0 0 0 0 0 0 0 quadslope: 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 quadslope: 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 2.7e-12 quadcoefs: 0.0e+00 0.0e+00 1.3e-14
doc : reference : source : test
The prototype is
double ***bezier(double **x,int nx) ;
where
#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 should be stored as beziertest.c. Compile and execute as
g++ -o -g beziertest beziertest.c bezier.c linsolve.c ; beziertest
and obtain the output
( 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)
#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 remaining test programs generate plots illustrative of Keynesian economics.
#include "memory.h" #include <math.h> double ***bezier(double **x,int nx) ; int main() { double **x=matrix(7,2),***b,q,h ; int i,j ; // green for(i=0;i<7;i++) x[i][0] = 60+30*i ; for(i=0;i<7;i++) x[i][1] = 150 - 150*exp(-(x[i][0]-50)/100) ; printf("y1=%d\n",(int)x[2][1]) ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<6;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int) b[i][j][0],b[i][j][1]) ; printf("\n ") ; } // blue for(i=0;i<7;i++) x[i][1] = 400 - 230*exp(-pow((x[i][0]-50)/100,0.8)) ; printf("y1=%d\n",(int)x[1][1]) ; h = 390 - x[2][1] ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<6;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int)b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; // red for(i=0;i<7;i++) x[i][1] = 390 - 120*exp((x[i][0]-50)/400) ; q = 390 - x[2][1] ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<5;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int)b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; for(i=0;i<7;i++) x[i][1] = 390 - 120*(h/q)*exp((x[i][0]-50)/400) ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<5;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int)b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; for(i=0;i<7;i++) x[i][1] = 390 - 45*exp((x[i][0]-50)/400) ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<6;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int)b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; }
#include "memory.h" #include <math.h> double ***bezier(double **x,int nx) ; int main() { double **x=matrix(11,2),***b,q,h ; int i,j ; // islm for(i=0;i<11;i++) x[i][0] = 40+20*i ; for(i=1;i<11;i++) x[i][1] = 180 - 150*exp(-pow((x[i][0]-50)/100,1.2)) ; q = x[6][1] ; printf("y1=%d\n",(int)x[6][1]) ; b = bezier(x+1,10) ; printf("M %5.1f %5.1f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<9;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%5.1f %5.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } for(i=0;i<11;i++) x[i][1] = 200 - 20*exp((x[i][0]-50)/90) ; h = x[6][1] ; for(i=0;i<11;i++) x[i][1] += q-h ; printf("\n") ; b = bezier(x,11) ; printf("M %5.1f %5.1f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<10;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%5.1f %5.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } }
#include "memory.h" #include <math.h> double ***bezier(double **x,int nx) ; int main() { double **x=matrix(11,2),***b,q,h ; int i,j ; // samuelson cross printf("\n") ; for(i=0;i<11;i++) x[i][0] = 40 + 20*i ; for(i=0;i<11;i++) x[i][1] = 20 + 120*exp(-(x[i][0]-50)/200) ; h = x[7][1] ; for(i=0;i<11;i++) x[i][1] += (70-h) ; b = bezier(x,11) ; printf("M %.1f %.2f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<10;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%2.1f %.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } }
#include "memory.h" #include <math.h> double ***bezier(double **x,int nx) ; int main() { double **x=matrix(50,2),***b,q,h ; int i,j ; // liquidity trap printf("\n") ; for(i=0;i<24;i++) x[i][0] = pow(4-0.125*i,-0.25) ; for(i=0;i<25;i++) x[24+i][0] = 1+0.125*i ; for(i=0;i<49;i++) x[i][1] = exp(-x[i][0]/3) + pow(x[i][0],-4) ; h = x[0][1] ; for(i=0;i<49;i++) { x[i][0] = 40 + 50*x[i][0] ; x[i][1] = 150 - 140*x[i][1]/h ; } h = x[7][1] ; printf("y1=%d\n",(int)h) ; b = bezier(x,49) ; printf("M %.2f %.2f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<48;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%.2f %.2f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } for(i=0;i<11;i++) x[i][0] = 40+20*i ; x[2][0] = 77 ; // blue for(i=0;i<11;i++) x[i][1] = 450 - 230*exp(-pow((x[i][0]-30)/200,0.8)) ; printf("y1=%d\n",(int)x[2][1]) ; h = 390 - x[2][1] ; b = bezier(x,11) ; printf("M %.1f %.1f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<10;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%2.1f %2.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; // red for(i=0;i<11;i++) x[i][1] = 390 - 140*exp((x[i][0]-50)/400) ; q = 390 - x[2][1] ; b = bezier(x,11) ; printf("M %.1f %.1f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<8;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%.1f %.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; for(i=0;i<11;i++) x[i][1] = 390 - (h/q)*(390-x[i][1]) ; b = bezier(x,11) ; printf("M %.1f %.1f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<8;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%.1f %.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; for(i=0;i<11;i++) x[i][1] = 390 - 0.9*(h/q)*(390-x[i][1]) ; b = bezier(x,11) ; printf("M %.1f %.1f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<10;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%.1f %.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; }
#include "memory.h" #include <math.h> double ***bezier(double **x,int nx) ; int main() { double **x=matrix(100,2),***b,q,h ; int i,j ; // blue for(i=0;i<10;i++) x[50+i][0] = 130 + i ; for(i=0;i<5;i++) x[60+i][0] = 140 + 2*i ; for(i=0;i<19;i++) x[65+i][0] = 150 + 5*i ; for(i=0;i<34;i++) x[50+i][1] = 210 - 2*pow((x[50+i][0]-120)/300,-1.2) ; q = x[50][1] ; printf("q=%.1f\n",q) ; x[46][0] = 105 ; x[46][1] = q - 80 ; x[47][0] = 115 ; x[47][1] = q - 69 ; x[48][0] = 123 ; x[48][1] = q - 52 ; x[49][0] = 127 ; x[49][1] = q - 30 ; b = bezier(&x[46],38) ; printf("M %6.2f %6.2f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<37;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%6.2f %6.2f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; // red for(i=0;i<7;i++) x[i][0] = 60+30*i ; x[2][0] = 130 ; for(i=0;i<7;i++) x[i][1] = 210 - 120*(0.75)*exp((x[i][0]-50)/400) ; q -= x[2][1] ; for(i=0;i<7;i++) x[i][1] += q ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<6;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int)b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; for(i=0;i<7;i++) x[i][1] = 210 - 120*(0.6)*exp((x[i][0]-50)/400) ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<6;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int)b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; for(i=0;i<7;i++) x[i][1] = 330 - 120*(0.75)*exp((x[i][0]-50)/400) ; b = bezier(x,7) ; printf("M %3d %5.1f ",(int) b[0][0][0],b[0][0][1]) ; for(i=0;i<6;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%3d %5.1f ",(int)b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; }
#include "memory.h" #include <math.h> double ***bezier(double **x,int nx) ; static double pi=3.141592653589793 ; int main() { double **x=matrix(18,2),***b,q,h,theta,gamma=2.8,alpha ; int i,j ; // lm for(i=0;i<18;i++) { x[i][0] = 40 + 10*i ; if(i==17) x[i][0] = 206 ; x[i][1] = 210 - 12 * ( 4 + 1/(1-(x[i][0]-40)/180) ) ; } b = bezier(x,18) ; printf("M %5.1f %5.1f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<17;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%5.1f %5.1f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; // asymptotic x[0][0] = 40 ; x[1][0] = 105 ; x[2][0] = 110 ; x[3][0] = 115 ; x[14][1] = 65 ; x[15][1] = 60 ; x[16][1] = 55 ; x[17][1] = 7 ; for(i=0;i<4;i++) { x[i][1] = 165 ; x[14+i][0] = 215 ; } for(i=0;i<10;i++) { theta = (1+i) * (pi/2) / 11 ; alpha = tan(theta) ; x[13-i][0] = 115 + pow(pow(100,gamma)/(1+pow(alpha,gamma)),1/gamma) ; x[13-i][1] = 65 + (x[13-i][0]-115)*alpha ; } for(i=0;i<18;i++) printf("%.2f %.2f\n",x[i][0],x[i][1]) ; b = bezier(x,18) ; printf("M %6.2f %6.2f ",b[0][0][0],b[0][0][1]) ; for(i=0;i<17;i++) { printf("C ") ; for(j=1;j<4;j++) printf("%6.2f %6.2f ",b[i][j][0],b[i][j][1]) ; printf("\n ") ; } printf("\n") ; }
doc : reference : source : test
double ranf(),gaussv(),gammav(double),betav(double,double) ;
int rani(int) ;
void ranset(int),ranperm(int*,int) ;
void direv(double*,int,double*),direv(double,int,double*) ;
void ranrot(double **q,int n) ;
struct multi { int n,*alias ; double *thresh ; } ;
multi multigen(double *p,int n) ;
int multiv(multi m) ;
void multifree(multi m),multichk(double *p,multi m) ;
void *mgvinit(double **,int),multigaussv(double *,void *) ;
void freemgvstate(void *) ;
functions provided : source : ranf test : gaussv 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. |
ranperm(a,n) | randomly permutes the supplied integer array a of length n (the version first provided instead generated a random permutation of 0…n1). |
gaussv() | returns a standard gaussian 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 • rani • ranperm • gaussv
#include <math.h> #define IA 16807 #define IM 2147483647 #define IQ 127773 #define IR 2836 #define NDIV (1+(IM-1)/32) static long val=0,tab[32],state=123467 ; static long step() { long k = state / IQ ; state = IA * (state-k*IQ) - IR*k ; if(state<0) state += IM ; return state ; } double ranf() { int i ; if(val==0) { for(i=63;i>=0;i--) tab[i&31] = step() ; val = tab[0] ; } i = val / NDIV ; val = tab[i] ; tab[i] = step() ; return val / (double) IM; } void ranset(int seed) { state = seed ; if(state>0) state *= 2 ; else state = 1 - 2*state ; if(state>=IM) state = 1 + state % (IM-1) ; val = 0 ; } int rani(int n) { return (int)(n*ranf()) ; } void ranperm(int *a,int n) { int i,j,u ; for(i=n-1;i>0;i--) { j = (int) ( (i+1)*ranf() ) ; u = a[i] ; a[i] = a[j] ; a[j] = u ; } } double sin(double),cos(double),log(double),sqrt(double) ; double gaussv() { static int toggle = 1 ; static double u1,u2,pi=3.141592653589793 ; if(0==(toggle=1-toggle)) { u1 = 0 ; while(u1==0) u1 = ranf() ; u1 = sqrt(-2*log(u1)) ; u2 = 2 * pi * ranf() ; return u1 * cos(u2) ; } else return u1 * sin(u2) ; }
To be stored as ranftest.c.
#include <stdio.h>
#include "random.h"
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"
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) ;
}
}
The results are given below.
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 # ranf and header files ranftest.c gaussvtest.c # test progs --- g++ -O -o ranftest ranftest.c ranf.c ranftest g++ -O -o gaussvtest gaussvtest.c ranf.c gaussvtest
functions provided : source : ranf test : gaussv test : download
void ranrot(double **q,int n) ;
fills n×n matrix q with a random rotation matrix.
#include "random.h"
#include "memory.h"
void qr(double **a,int m,int n,double **Q,double **R) ;
void ranrot(double **q,int n)
{ int i,j ;
double **Q=matrix(n,n),**R=matrix(n,n) ;
for(i=0;i<n;i++) for(j=0;j<n;j++) Q[i][j] = gaussv() ;
qr(Q,n,n,q,R) ;
for(j=0;j<n;j++) if(R[j][j]<0) for(i=0;i<n;i++) q[i][j] = -q[i][j] ;
freematrix(Q,R) ;
}