/* timegauss.c - measure the speed of random number generation
 *
 * Copyright (C) 2005  Jochen Voss.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <stdio.h>
#include <stdlib.h>
#include <signal.h>
#include <setjmp.h>
#include <sys/time.h>

#include <gsl/gsl_rng.h>


extern double gsl_ran_gaussian_ziggurat (gsl_rng *r, const double sigma);
extern double gsl_ran_gaussian2 (gsl_rng *r, const double sigma);
extern double gsl_ran_gaussian3 (gsl_rng *r, const double sigma);


static  sigjmp_buf  leave_loop;

static void
handle_alarm (int signum)
{
  siglongjmp (leave_loop, 1);
}

static unsigned long
count_calls (gsl_rng *rng, double (*f)(gsl_rng *, const double))
{
  volatile  unsigned long  n = 0;
  struct itimerval  pause;
  
  pause.it_interval.tv_usec = 0;
  pause.it_interval.tv_sec = 0;
  pause.it_value.tv_usec = 0;
  pause.it_value.tv_sec = 1;

  signal (SIGVTALRM, handle_alarm);
  if (sigsetjmp (leave_loop, 1) == 0) {
    if (setitimer (ITIMER_VIRTUAL, &pause, NULL) < 0)  abort ();
    while (1) {
      f (rng, 1);
      ++n;
    }
  }

  return  n;
}

static void
test_with_rng (const gsl_rng_type *T, FILE *res)
{
  gsl_rng *rng = gsl_rng_alloc (T);
  unsigned long  c1, c2, c3;	/* generated numbers / second */
  double  m1, m2, m3;		/* averages */
  double  q;
  int  i;

  printf ("testing with %s:\n", T->name);

  m1 = m2 = m3 = 0;
  for (i=0; i<=10; ++i) {
    c1 = count_calls (rng, gsl_ran_gaussian_ziggurat);
    c2 = count_calls (rng, gsl_ran_gaussian2);
    c3 = count_calls (rng, gsl_ran_gaussian3);
    printf ("  %lu %lu %lu numbers/second\n", c1, c2, c3);
    if (i>0) {
      m1 += c1;
      m2 += c2;
      m3 += c3;
    }
  }
  q = m1 / (m2>m3 ? m2 : m3);
  fprintf (res, "<tr><td>%s<td>%.3g<td>%.3g<td>%.3g<td>%.1f\n",
	   T->name, m1/10.0e6, m2/10.0e6, m3/10.0e6, q);
  fflush (NULL);
  gsl_rng_free (rng);
}

int
main ()
{
  const gsl_rng_type **t, **t0;
  FILE *res;

  res = fopen ("data.html", "w");
  gsl_rng_env_setup();


  t0 = gsl_rng_types_setup ();
  for (t = t0; *t != 0; t++) {
    test_with_rng (*t, res);
  }

  fclose (res);
  return  0;
}