/* rng/uni32.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 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., 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;