/* This file is part of the alrand library.
   alrand is Copyright by Vincent Penquerc'h <lyrian -at- kezako -dot- net>.

   The alrand library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The alrand library 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
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the alrand Library; 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 <time.h>
#include <math.h>
#include <sys/time.h>
#include <unistd.h>
#include "lsstdio.h"
#include "alrand.h"

/* some systems (eg win32) don't seem to bother defining this */
#ifndef M_PI
#define M_PI           3.14159265358979323846  /* pi */
#endif

static struct {
  char *name;
  alrand_generator_type type;
} types[]={
  {"Mersenne Twister",alrand_type_mersenne_twister},
  {"Linear congruential",alrand_type_linear_congruential},
};

static void test_mersenne_twister()
{
  int i;
  unsigned long seed=4357;
  alrand_state *mt=alrand_create(alrand_type_mersenne_twister,seed);

  fprintf(stderr,"\nThis test compares Mersenne Twister generated numbers against reference data\n");

  for (i=0;i<1000;i++) {
    printf("%10lu ", alrand_generate(mt));
    if (i%5==4) printf("\n");
  }

  alrand_destroy(mt);
}

static void approximate_pi()
{
  size_t n;
  unsigned long seed=time(NULL);

  fprintf(stderr,"\nThis test computes an approximate value of PI by checking the distribution\nof random points in a square and circle\n");
  for (n=0;n<sizeof(types)/sizeof(types[0]);++n) {
    alrand_state *mt=alrand_create(types[n].type,seed);
    int i,in=0;
    double pi;
    int num_points=4000000;
    for (i=0;i<num_points;i++) {
      double r1=(double)(alrand_generate(mt)&0xffff)/65536.0;
      double r2=(double)(alrand_generate(mt)&0xffff)/65536.0;
      double radius=r1*r1+r2*r2;
      if (radius<1.0) in++;
    }
    /* square area: 4 */
    /* circle area: pi */
    /* => pi/4 ~= in/num_points */
    pi=4.0*in/num_points;
    fprintf(stderr,"Using %s: PI==%f (%f variation)\n",types[n].name,pi,fabs(M_PI-pi));
    alrand_destroy(mt);
  }
}

static void test_persistence()
{
  size_t n,i;

  fprintf(stderr,"\nThis test checks that two states with the same parameters and seed yield\nthe same numbers after one of them is saved and restored\n");

  for (n=0;n<sizeof(types)/sizeof(types[0]);++n) {
    FILE *file;
    alrand_uint32_t seed=time(NULL);
    alrand_state *state1=alrand_create(types[n].type,seed);
    alrand_state *state2=alrand_create(types[n].type,seed);
    fprintf(stderr,"%s",types[n].name);
    for (i=0;i<100;++i) {
      int n1=alrand_generate(state1);
      int n2=alrand_generate(state2);
      if (n1!=n2) {
        fprintf(stderr,"Error: generated data differs\n");
        exit(1);
      }
    }
    file=tmpfile();
    if (!file) {
        fprintf(stderr,"Error: unable to create a file to save state\n");
        exit(1);
    }
    if (alrand_save_to_stdio_file(state2,file)) {
      fprintf(stderr,"Error: unable to save state\n");
      exit(1);
    }
    alrand_destroy(state2);
    fseek(file,0,SEEK_SET);
    state2=alrand_load_from_stdio_file(file);
    if (!state2) {
      fprintf(stderr,"Error: unable to load saved state\n");
      exit(1);
    }
    for (i=0;i<10000;++i) {
      int n1=alrand_generate(state1);
      int n2=alrand_generate(state2);
      if (n1!=n2) {
        fprintf(stderr,"Error: generated data differs\n");
        exit(1);
      }
    }
    fclose(file);
    fprintf(stderr," save/load test passed\n");
  }
}

static void test_speed()
{
#ifndef _WIN32
  size_t n,i;
  size_t num_points=40000000;
  struct timeval start,stop;
  fprintf(stderr,"\nThis test makes a very rough generation speed measurement\n");
  for (n=0;n<sizeof(types)/sizeof(types[0]);++n) {
    alrand_state *state=alrand_create(types[n].type,time(NULL));
    gettimeofday(&start,NULL);
    for (i=0;i<num_points;++i) {
      alrand_generate(state);
    }
    gettimeofday(&stop,NULL);
    fprintf(stderr,"%s generated %d numbers in %2.3f seconds\n",types[n].name,num_points,(float)(stop.tv_sec-start.tv_sec+(stop.tv_usec-start.tv_usec)/1000000.0f));
  }

  gettimeofday(&start,NULL);
  for (i=0;i<num_points;++i) {
    rand();
  }
  gettimeofday(&stop,NULL);
  fprintf(stderr,"libc rand generated %d numbers in %2.3f seconds\n",num_points,(float)(stop.tv_sec-start.tv_sec+(stop.tv_usec-start.tv_usec)/1000000.0f));
#else
  fprintf(stderr,"POSIXness is too tough for poor Win32, so you don't get timing information\n");
#endif
}

int main()
{
  fprintf(stderr,"\nTesting alrand...\n");
  approximate_pi();
  test_persistence();
  test_speed();
  test_mersenne_twister();
  return 0;
}
