single-rnd.c

     
   1  //! @file single-rnd.c
   2  //! @author J. Marcel van der Veer
   3  //!
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
   8  //!
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! REAL pseudo-random number generator.
  25  
  26  #include "a68g.h"
  27  #include "a68g-genie.h"
  28  #include "a68g-prelude.h"
  29  #include "a68g-double.h"
  30  #include "a68g-numbers.h"
  31  
  32  // Next part is a "stand-alone" version of GNU Scientific Library (GSL)
  33  // random number generator "taus113", based on GSL file "rng/taus113.c".
  34  
  35  // rng/taus113.c
  36  // Copyright (C) 2002 Atakan Gurkan
  37  // Based on the file taus.c which has the notice
  38  // Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
  39  // 
  40  // This program is free software; you can redistribute it and/or modify
  41  // it under the terms of the GNU General Public License as published by
  42  // the Free Software Foundation; either version 3 of the License, or (at
  43  // your option) any later version.
  44  // 
  45  // This program is distributed in the hope that it will be useful, but
  46  // WITHOUT ANY WARRANTY; without even the implied warranty of
  47  // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  48  // General Public License for more details.
  49  // 
  50  // You should have received a copy of the GNU General Public License
  51  // along with this program; if not, write to the Free Software
  52  // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  53  //
  54  // This is a maximally equidistributed combined, collision free 
  55  // Tausworthe generator, with a period ~2^{113}. The sequence is,
  56  //
  57  //   x_n = (z1_n ^ z2_n ^ z3_n ^ z4_n)  
  58  //
  59  //   b = (((z1_n <<  6) ^ z1_n) >> 13)
  60  //   z1_{n+1} = (((z1_n & 4294967294) << 18) ^ b)
  61  //   b = (((z2_n <<  2) ^ z2_n) >> 27)
  62  //   z2_{n+1} = (((z2_n & 4294967288) <<  2) ^ b)
  63  //   b = (((z3_n << 13) ^ z3_n) >> 21)
  64  //   z3_{n+1} = (((z3_n & 4294967280) <<  7) ^ b)
  65  //   b = (((z4_n <<  3)  ^ z4_n) >> 12)
  66  //   z4_{n+1} = (((z4_n & 4294967168) << 13) ^ b)
  67  //
  68  // computed modulo 2^32. In the formulas above '^' means exclusive-or 
  69  // (C-notation), not exponentiation. 
  70  // The algorithm is for 32-bit integers, hence a bitmask is used to clear 
  71  // all but least significant 32 bits, after left shifts, to make the code 
  72  // work on architectures where integers are 64-bit.
  73  //
  74  // The generator is initialized with 
  75  // z{i+1} = (69069 * zi) MOD 2^32 where z0 is the seed provided
  76  // During initialization a check is done to make sure that the initial seeds 
  77  // have a required number of their most significant bits set.
  78  // After this, the state is passed through the RNG 10 times to ensure the
  79  // state satisfies a recurrence relation.
  80  //
  81  // References:
  82  //   P. L'Ecuyer, "Tables of Maximally-Equidistributed Combined LFSR Generators",
  83  //   Mathematics of Computation, 68, 225 (1999), 261--269.
  84  //     http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps
  85  //   P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe Generators", 
  86  //   Mathematics of Computation, 65, 213 (1996), 203--213.
  87  //     http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
  88  //   the online version of the latter contains corrections to the print version.
  89  
  90  #define LCG(n) ((69069UL * n) & 0xffffffffUL)
  91  #define MASK 0xffffffffUL
  92  
  93  unt taus113_get (void *vstate);
  94  double taus113_get_double (void *vstate);
  95  void taus113_set (void *state, unt long int s);
  96  
  97  typedef struct
  98  {
  99    unt long int z1, z2, z3, z4;
 100  }
 101  taus113_state_t;
 102  
 103  static taus113_state_t rng_state;
 104  
 105  unt taus113_get (void *vstate)
 106  {
 107    taus113_state_t *state = (taus113_state_t *) vstate;
 108    unt long b1, b2, b3, b4;
 109    b1 = ((((state->z1 << 6UL) & MASK) ^ state->z1) >> 13UL);
 110    state->z1 = ((((state->z1 & 4294967294UL) << 18UL) & MASK) ^ b1);
 111    b2 = ((((state->z2 << 2UL) & MASK) ^ state->z2) >> 27UL);
 112    state->z2 = ((((state->z2 & 4294967288UL) << 2UL) & MASK) ^ b2);
 113    b3 = ((((state->z3 << 13UL) & MASK) ^ state->z3) >> 21UL);
 114    state->z3 = ((((state->z3 & 4294967280UL) << 7UL) & MASK) ^ b3);
 115    b4 = ((((state->z4 << 3UL) & MASK) ^ state->z4) >> 12UL);
 116    state->z4 = ((((state->z4 & 4294967168UL) << 13UL) & MASK) ^ b4);
 117    return (state->z1 ^ state->z2 ^ state->z3 ^ state->z4);
 118  }
 119  
 120  double taus113_get_double (void *vstate)
 121  {
 122    return taus113_get (vstate) / 4294967296.0;
 123  }
 124  
 125  void taus113_set (void *vstate, unt long int s)
 126  {
 127    taus113_state_t *state = (taus113_state_t *) vstate;
 128    if (!s) {
 129      s = 1UL; // default seed is 1
 130    }
 131    state->z1 = LCG (s);
 132    if (state->z1 < 2UL) {
 133      state->z1 += 2UL;
 134    }
 135    state->z2 = LCG (state->z1);
 136    if (state->z2 < 8UL) {
 137      state->z2 += 8UL;
 138    }
 139    state->z3 = LCG (state->z2);
 140    if (state->z3 < 16UL) {
 141      state->z3 += 16UL;
 142    }
 143    state->z4 = LCG (state->z3);
 144    if (state->z4 < 128UL) {
 145      state->z4 += 128UL;
 146    }
 147  // Calling RNG ten times to satify recurrence condition
 148    taus113_get (state);
 149    taus113_get (state);
 150    taus113_get (state);
 151    taus113_get (state);
 152    taus113_get (state);
 153    taus113_get (state);
 154    taus113_get (state);
 155    taus113_get (state);
 156    taus113_get (state);
 157    taus113_get (state);
 158    return;
 159  }
 160  
 161  // Initialise rng.
 162  
 163  void init_rng (unt u)
 164  {
 165    taus113_set (&rng_state, u);
 166  }
 167  
 168  // A68G rng in R mathlib style.
 169  
 170  REAL_T a68_unif_rand (void)
 171  {
 172  // In [0, 1>
 173    return taus113_get_double (&rng_state);
 174  }
 175  
 176  static char *state_file = ".Random.seed";
 177  
 178  void GetRNGstate (void)
 179  {
 180    INT_T fd = open (state_file, A68_READ_ACCESS);
 181    if (fd != -1) {
 182      ASSERT (read (fd, &rng_state, sizeof (taus113_state_t)) != -1);
 183      close (fd);
 184    }
 185  }
 186  
 187  void PutRNGstate (void)
 188  {
 189    INT_T fd = open (state_file, A68_WRITE_ACCESS, A68_PROTECTION);
 190    if (fd != -1) {
 191      ASSERT (write (fd, &rng_state, sizeof (taus113_state_t)) != -1);
 192      close (fd);
 193    }
 194  }