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.
x = imatrix(m,n,l)	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.
x = imatrix(m,n,l)	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).
x = charvector(char *s)	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.

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 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     • 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

remarks : 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.

source : remarks : test

doc : source : test

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 ij’s.

ij *ijvector(int *x,int n) ;

returns an n-long array of ij’s whose k^th 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)

doc : source : test

doc : source

linsolve solves one or more sets of linear equations ax_i = b_i. The prototype is

double linsolve(double **a,double **b,double **x,int n,int nrhs) ;

where

• a is an n×n matrix of equation coefficients;
• b[i] is the i^th right hand side;
• x[i] is to receive the i^th solution;
• n is the dimension; and
• nrhs is the number of right-hand sides.

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 ; 
}

doc : source

doc : source

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] ; 
}

doc : source

doc : source : test

void qr(double **a,int m,int n,double **Q,double **R) ;

where

• a is the m by n matrix to decompose (m≥n);
• Q is an m by m matrix
• m to receive one component of the decomposition; and
• R is an m by n matrix
• m to receive the other.

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

doc : source : test

Suppose that we fit a quadratic to the three points (–1,a), (0,b), (1,c). Then

• quadinterp returns the value of the quadratic at x,
• quadslope returns the gradient at x,
• invquadinterp returns the argument for which the value is y,
• quadreach returns 1 if the value y is reachable, else 0,
• quadcoefs returns the quadratic coefficients in co,
• quadmax returns the (x,y) pair at the turning point (which may be a minimum or a maximum), and
• quadintcoefs returns coefficients for integrating the quadratic interpolate (see below).

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 x².

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 (x₀,y₀), (x₁,y₁), and (x₂,y₂), then its integral from x₀ to x₁ will be

c₀y₀ + c₁y₁ + c₂y₂

The first call to quadintcoefs returns in c the coefficients corresponding to the given x_i. The second does so unless the pointer c is 0.

The second call also optionally returns coefficients for the integrals from x₀ to (x₀+x₁)/2 and from (x₀+x₁)/2 to x₁.

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:

• lininterp returns the value of the function at x, and
• invquadinterp returns the argument for which the value is y.

For cubics I provide cubcoefs and cubintcoefs, the latter integrating from x₁ to x₂.

• 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 : source : test

doc : reference : source : test

The prototype is

double ***bezier(double **x,int nx) ; 

where

• x is an nx by 2 matrix of points in 2-dimensional space, and
• the return value is an (nx–1) by 4 by 2 matrix of bezier control points ready for use in SVG.

• Notes by Kirby A. Baker, UCLA.

#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…n–1).
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 x⁴ 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

doc : source

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) ; 
}

doc : source