/* ziggurat.c - construct the tables for "normal.c"
 *
 * For details see the article
 * George Marsaglia, Wai Wan Tsang
 * The Ziggurat Method for Generating Random Variables
 * Journal of Statistical Software, vol. 5 (2000), no. 8
 * http://www.jstatsoft.org/v05/i08/
 *
 * Copyright (C) 2005  Jochen Voss.
 *
 * $Id: ziggurat.c 6509 2005-07-07 18:31:10Z voss $
 */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

#include <stdio.h>
#include <math.h>


#define LEVELS 128


static  double  x[LEVELS];
static  double  v;


static double
f (double x)
/* the target density */
{
  return  exp(-x*x/2);
}

static double
f_inv (double y)
/* the inverse of the target density */
{
  return  sqrt(-2*log(y));
}

static double
try_r_value (double r)
{
  int  i;

  v = r*f(r) + exp(-0.5*r*r)/r;
  x[LEVELS-1] = r;
  for (i=LEVELS-1; i>1; --i) {
    x[i-1] = f_inv(v/x[i]+f(x[i]));
  }
  return  x[1]*(1-f(x[1]))/v;
}

int
main ()
{
  double  a, b, aa, bb, r;
  int  i;

  a=0;
  b=10;
  do {
    double  q;
    
    aa=a, bb=b;
    r=.5*(a+b);
    q=try_r_value(r);
    if (q>1) {
      b = r;
    } else {
      a = r;
    }
  } while (aa<r && r<bb);
  x[0] = 0;

  puts ("/* automatically generated by ziggurat.c */\n");
  
  puts ("/* number of levels in the ziggurat */");
  printf ("#define LEVELS %d\n\n", LEVELS);

  puts ("/* position of right-most step */");
  printf ("#define PARAM_R %.12g\n\n", r);
  
  puts ("/* level values */");
  printf ("static const double ytab[%d] = {\n", LEVELS);
  for (i=0; i<LEVELS; ++i) {
    printf ("  %.12g,\n", f(x[i]));
  }
  printf ("};\n\n");
  
  puts ("/* quick acceptance check */");
  printf ("static const unsigned long ktab[%d] = {\n", LEVELS);
  for (i=0; i<LEVELS-1; ++i) {
    printf ("  %lu,\n", (unsigned long)(pow(2,24)*x[i]/x[i+1]));
  }
  printf ("  %lu\n", (unsigned long)(pow(2,24)*r*f(r)/v));
  printf ("};\n\n");

  puts ("/* quick value conversion */");
  printf ("static const double wtab[%d] = {\n", LEVELS);
  for (i=0; i<LEVELS-1; ++i) {
    printf ("  %.12g,\n", x[i+1]/pow(2,24));
  }
  printf ("  %.12g\n", v/(pow(2,24)*f(r)));
  printf ("};\n");
  
  return 0;
}