src/contrib/random.c

   1 /*
   2  * Print a lookup table for Gaussian numbers with 4 entries on each
   3  * line, formatted for inclusion in this file. Size is 2^bits.
   4  */
   5
   6 void
   7 print_gaussian_table(int bits)
   8 {
   9   int n,nh,i,j;
  10   double invn,fac,x,invgauss,det,dx;
  11   real  *table;
  12
  13   n = 1 << bits;
  14   table = (real *)malloc(n*sizeof(real));
  15
  16   /* Fill a table of size n such that random draws from it
  17     * produce a Gaussian distribution.
  18     * We integrate the Gaussian distribution G approximating:
  19     *   integral(x->x+dx) G(y) dy
  20     * with:
  21     *   G(x) dx + G'(x) dx^2/2 = G(x) dx - G(x) x dx^2/2
  22     * Then we need to find dx such that the integral is 1/n.
  23     * The last step uses dx = 1/x as the approximation is not accurate enough.
  24     */
  25   invn = 1.0/n;
  26   fac = sqrt(2*M_PI);
  27   x = 0.5*fac*invn;
  28   nh = n/2;
  29   for(i=0; i<nh; i++) {
  30     if (i > 0) {
  31       if (i < nh-1) {
  32         invgauss = fac*exp(0.5*x*x);
  33         /* det is larger than 0 for all x, except for the last */
  34         det = 1 - 2*invn*x*invgauss;
  35         dx = (1 - sqrt(det))/x;
  36       } else {
  37         dx = 1/x;
  38       }
  39       x = x + dx;
  40     }
  41     table[nh-1-i] = -x;
  42     table[nh+i]   =  x;
  43   }
  44   printf("static const real *\ngaussian_table[%d] = {\n",n);
  45   for(i=0;i<n;i+=4) {
  46     printf("  ");
  47     for(j=0;j<4;j++) {
  48       printf("%14.7e",table[i+j]);
  49       if(i+j<(n-1))
  50         printf(",");
  51     }
  52     printf("\n");
  53   }
  54   printf("};\n");
  55   free(table);
  56 }
  57