Colin Champion

 

functiondescription

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.

• cjcup     • cjcuplog     • cjcformat     • cjcprint2     • cjcprint1     • complex     • print     • xy     • print     • xi     • settable     • 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>

#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))

/* ---------------------------- 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) { set = 1 ; x = y ; }
    operator double() 
    { if(!set) throw up("accessing an unset settable") ; return x ; }
    settable &operator=(double 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() ; }
} ;
/* ------------------------------ 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 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 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)

doc : source : test

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

• 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

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

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

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