/* rng/uni32.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough

 * 

 * 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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

/**

  This is a lagged Fibonacci generator which supposedly excellent

  statistical properties (I do not concur)

  I got it from the net and translated into C.

* ======================================================================

* NIST Guide to Available Math Software.

* Fullsource for module UNI from package CMLIB.

* Retrieved from CAMSUN on Tue Oct  8 14:04:10 1996.

* ======================================================================

C***BEGIN PROLOGUE  UNI

C***DATE WRITTEN   810915

C***REVISION DATE  830805

C***CATEGORY NO.  L6A21

C***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS

C***AUTHOR    BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS

C             KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS

C             MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV

C

C***PURPOSE  THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1

C             AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS

C             AT LEAST AS LARGE AS 32767.

C***DESCRIPTION

C

C       THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTER

C       [0,1).  IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS

C       INTEGERS AT LEAST AS LARGE AS 32767.

C

C

C   USE

C       FIRST TIME....

C                   Z = UNI(JD)

C                     HERE JD IS ANY  N O N - Z E R O  INTEGER.

C                     THIS CAUSES INITIALIZATION OF THE PROGRAM

C                     AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z.

C       SUBSEQUENT TIMES...

C                   Z = UNI(0)

C                     CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z.

C

C

C..................................................................

C   NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER

C         TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION.....

C

C   MACHINE DEPENDENCIES...

C      MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE

C              FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT.

C              THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED

C              IN LINE WITH REMARK A BELOW.

C

C   REMARKS...

C     A. THIS PROGRAM CAN BE USED IN TWO WAYS:

C        (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS,

C            SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR,

C        (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE

C            GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE

C            LARGEST POSSIBLE VALUE.

C     B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL

C          INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'.

C          IF MDIG=16 ONE SHOULD FIND THAT

   Editors Note: set the seed using 152 in order to get uni(305)

   -jt

C            THE FIRST EVALUATION

C              Z=UNI(305) GIVES Z=.027832881...

C            THE SECOND EVALUATION

C              Z=UNI(0) GIVES   Z=.56102176...

C            THE THIRD EVALUATION

C              Z=UNI(0) GIVES   Z=.41456343...

C            THE THOUSANDTH EVALUATION

C              Z=UNI(0) GIVES   Z=.19797357...

C

C***REFERENCES  MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM

C                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.

C***ROUTINES CALLED  I1MACH,XERROR

C***END PROLOGUE  UNI

  **/

#include <config.h>

#include <stdlib.h>

#include <gsl/gsl_rng.h>

static inline unsigned long int uni32_get (void *vstate);

static double uni32_get_double (void *vstate);

static void uni32_set (void *state, unsigned long int s);

static const unsigned long int MDIG = 32;       /* Machine digits in int */

static const unsigned long int m1 = 2147483647;         /* 2^(MDIG-1) - 1 */

static const unsigned long int m2 = 65536;      /* 2^(MDIG/2) */

typedef struct

  {

    int i, j;

    unsigned long m[17];

  }

uni32_state_t;

static inline unsigned long

uni32_get (void *vstate)

{

  uni32_state_t *state = (uni32_state_t *) vstate;

  const long int i = state->i;

  const long int j = state->j;

  /* important k not be unsigned */

  long int k = state->m[i] - state->m[j];

  if (k < 0)

    k += m1;

  state->m[j] = k;

  if (i == 0)

    {

      state->i = 16;

    }

  else

    {

      (state->i)--;

    }

  if (j == 0)

    {

      state->j = 16;

    }

  else

    {

      (state->j)--;

    }

  return k;

}

static double

uni32_get_double (void *vstate)

{

  return uni32_get (vstate) / 2147483647.0 ;

}

static void

uni32_set (void *vstate, unsigned long int s)

{

  long int seed, k0, k1, j0, j1;

  int i;

  uni32_state_t *state = (uni32_state_t *) vstate;

  /* For this routine, the seeding is very elaborate! */

  /* A flaw in this approach is that seeds 1,2 give exactly the

     same random number sequence!  */

  /*s = 2*s+1; *//* enforce seed be odd */

  seed = (s < m1 ? s : m1);     /* seed should be less than m1 */

  seed -= (seed % 2 == 0 ? 1 : 0);

  k0 = 9069 % m2;

  k1 = 9069 / m2;

  j0 = seed % m2;

  j1 = seed / m2;

  for (i = 0; i < 17; i++)

    {

      seed = j0 * k0;

      j1 = (seed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);

      j0 = seed % m2;

      state->m[i] = j0 + m2 * j1;

    }

  state->i = 4;

  state->j = 16;

  return;

}

static const gsl_rng_type uni32_type =

{"uni32",                       /* name */

 2147483646,                    /* RAND_MAX */

 0,                             /* RAND_MIN */

 sizeof (uni32_state_t),

 &uni32_set,

 &uni32_get,

 &uni32_get_double};

const gsl_rng_type *gsl_rng_uni32 = &uni32_type;