mingwrt/mingwex/math/random.c

   1 /*
   2  * random.c
   3  *
   4  * Implementation of a (mostly) POSIX.1-1990 conforming pseudo-random
   5  * number generating API, for use with MinGW applications.
   6  *
   7  * $Id$
   8  *
   9  * Written by Keith Marshall <keith@users.osdn.me>
  10  * Copyright (C) 2021, MinGW.org Project
  11  *
  12  *
  13  * This is free software.  Permission is granted to copy, modify and
  14  * redistribute this software, under the provisions of the GNU General
  15  * Public License, Version 3, (or, at your option, any later version),
  16  * as published by the Free Software Foundation; see the file COPYING
  17  * for licensing details.
  18  *
  19  * Note, in particular, that this software is provided "as is", in the
  20  * hope that it may prove useful, but WITHOUT WARRANTY OF ANY KIND; not
  21  * even an implied WARRANTY OF MERCHANTABILITY, nor of FITNESS FOR ANY
  22  * PARTICULAR PURPOSE.  Under no circumstances will the author, or the
  23  * MinGW Project, accept liability for any damages, however caused,
  24  * arising from the use of this software.
  25  *
  26  *
  27  * This pseudo-random number generation suite has been derived from a
  28  * discourse on the GLIBC implementation, by Peter Selinger:
  29  *
  30  *   https://www.mathstat.dal.ca/~selinger/random/
  31  *
  32  * It is believed that the output from the PRNG will closely mimic that
  33  * from the GLIBC implementation; however, neither the author of this
  34  * implementation, nor the MinGW.org Project, offer any assurance as to
  35  * the statistical quality of the generated number sequence.
  36  *
  37  *
  38  * Note: this implementation has been provided for optional compliance
  39  * with the POSIX.1-1990 Extended Systems Interface.
  40  */
  41 #define _XOPEN_SOURCE  500
  42
  43 #include <stdint.h>
  44 #include <stdlib.h>
  45 #include <errno.h>
  46
  47 struct state
  48 { /* Internal structure, used to map the PRNG state buffer from its
  49    * arbitrarily specified sequence of 8, 32, 64, 128, or 256 bytes,
  50    * to an array of int32_t entities; the first entry is used as a
  51    * control sturcture, comprising three 8-bit indices, while the
  52    * remaining entries are used to represent PRNG state.
  53    */
  54   struct
  55   { uint32_t    phase:8;        /* current state data cycle pointer  */
  56     uint32_t    shift:8;        /* associated offset data pointer    */
  57     uint32_t    limit:8;        /* total count of state data entries */
  58   };
  59   int32_t       data[];         /* array of state data entries       */
  60 };
  61
  62 /* Provide a local state buffer, of length 128 bytes, (equivalent
  63  * to 32 int32_t entries), for use when no alternative buffer has
  64  * been specified, prior to calling random().
  65  */
  66 static int32_t default_state[32] = { 0 };
  67 static struct state *state = (struct state *)(default_state);
  68
  69 /* The following inline local function is provided to facilitate
  70  * setting of "errno", on abnormal return from any other function.
  71  */
  72 static __inline__ __attribute__((__always_inline__))
  73 intptr_t errout( int code, intptr_t retval ){ errno = code; return retval; }
  74
  75 static char *assign_state( char *buf, char *prev )
  76 { /* A local helper function, to attach a new state buffer to the
  77    * PRNG; aborts the assignment, if the specified buffer is NULL,
  78    * but otherwise, performs no buffer validation.
  79    */
  80   if( buf == NULL ) return (char *)(errout( EINVAL, (intptr_t)(NULL) ));
  81   state = (struct state *)(buf);
  82   return prev;
  83 }
  84
  85 char *__mingw_setstate( char *buf )
  86 { /* Public API entry point, exhibiting POSIX.1-1990 semantics,
  87    * for accessing the preceding buffer assignment function; the
  88    * "buf" argument is assumed to point to an existing non-NULL
  89    * buffer, which had been previously initialized by a prior
  90    * initstate() call, (and possibly already used), but, other
  91    * that checking that it is not NULL, no validation of its
  92    * initialization state is performed.
  93    */
  94   return assign_state( buf, (char *)(state) );
  95 }
  96
  97 static size_t normalized_cycle( size_t len )
  98 { /* A local helper function, to establish the effective number
  99    * of entries in the PRNG state data array, as a dependency of
 100    * the total byte count specified for the state buffer; note
 101    * that the compiler may be able to compute the resultant at
 102    * compile time, when called with a constant "len" argument,
 103    * and so optimize away such calls.
 104    */
 105   size_t cycle_counter = (len >= 32) ? 256 : 8;
 106   while( cycle_counter > len ) cycle_counter >>= 1;
 107   return (cycle_counter >> 2) - 1;
 108 }
 109
 110 static long update_state( void )
 111 { /* A local helper function, to update the content of the active
 112    * state buffer with each successive call of random().
 113    */
 114   if( state->phase == state->limit ) state->phase = 0;
 115   else if( state->shift == state->limit ) state->shift = 0;
 116   return (state->data[state->shift++] += state->data[state->phase++]);
 117 }
 118
 119 static void initialize_state_data( unsigned int seed )
 120 { /* A local helper function, called after initialization of the
 121    * state buffer length count, to populate the state data array.
 122    *
 123    * In all cases, we record the "seed" value, as the first entry
 124    * in the state data array, (rejecting zero, for which we force
 125    * substitution of 1).
 126    */
 127   state->data[0] = (seed != 0) ? seed : 1;
 128
 129   /* For a minimum-length state buffer, of only 8 bytes, there is
 130    * no more state data to be recorded...
 131    */
 132   if( state->limit > 1 )
 133   { /* ...but for any of the larger permitted buffer sizes, we must
 134      * populate the buffer using a multiplicative sequence of 31-bit
 135      * integers; the multiplication is performed in a 64-bit space,
 136      * then reduced modulo the Mersenne prime of order 31; set up a
 137      * collection of constants, pertinent to the computation...
 138      */
 139     const unsigned long order = 31;
 140     const unsigned long divisor = (unsigned long)((1ULL << order) - 1);
 141     const unsigned long long multiplier = 16807ULL;
 142
 143     /* ...then apply to each data array entry, in turn...
 144      */
 145     for( seed = 1; seed < state->limit; seed++ )
 146     { /* ...computing the sequence of products, in 64-bit space...
 147        */
 148       unsigned long long tmp = multiplier * state->data[seed - 1];
 149
 150       /* ...then reduce for storage; note that, since the divisor
 151        * is a Mersenne prime, we can perform the modulo division by
 152        * use of the more computationally efficient mask, shift and
 153        * addition, followed by a conditional correction, in case
 154        * of overflow into the sign bit.
 155        */
 156       state->data[seed] = (int)(tmp = (tmp & divisor) + (tmp >> order));
 157       if( state->data[seed] < 0 ) state->data[seed] += divisor;
 158     }
 159     /* We must also initialize the phase, and shift indices, for
 160      * the PRNG feed-back loop...
 161      */
 162     state->phase = 0;
 163     switch( state->limit )
 164     { /* ...noting that the shift interval is chosen differently,
 165        * depending on the limiting length of the sequence.
 166        */
 167       case 63: case 15: state->shift = 1; break;
 168       case 31: case  7: state->shift = 3; break;
 169     }
 170     /* Finally, run the PRNG through 10 update cycles for each and
 171      * every entry in the state data array, but without generating
 172      * any actual pseudo-random number output.
 173      */
 174     seed = ((state->limit << 2) + state->limit) << 1;
 175     while( seed-- > 0 ) update_state();
 176   }
 177 }
 178
 179 char *__mingw_initstate( unsigned int seed, char *buf, size_t len )
 180 { /* Public API entry point, implementing the POSIX.1-1990 initstate()
 181    * function; aborts if "buf" is passed as a NULL pointer, or if "len"
 182    * is less than the minimum allowed 8 bytes; otherwise, assigns "buf"
 183    * as a new PRNG state buffer...
 184    */
 185   char *prev = assign_state( buf, (len >= 8) ? (char *)(state) : NULL );
 186   if( prev != NULL )
 187   {
 188     /* ...initializes it to the specified "len", and as if srandom()
 189      * has been called with "seed" as argument...
 190      */
 191     state->limit = normalized_cycle( len );
 192     initialize_state_data( seed );
 193   }
 194   /* ...then returns a pointer to the original state buffer.
 195    */
 196   return prev;
 197 }
 198
 199 void __mingw_srandom( unsigned int seed )
 200 { /* Public API entry point, implementing the POSIX-1.1990 srandom()
 201    * function; this sets the PRNG to a known initial state, which is
 202    * dependent on the "seed" value specified.
 203    */
 204   if( (state == (struct state *)(default_state)) && (state->limit == 0) )
 205     state->limit = normalized_cycle( 128 );
 206   initialize_state_data( seed );
 207 }
 208
 209 long __mingw_random( void )
 210 { /* Public API entry point, implementing the POSIX-1.1990 random()
 211    * function; this initially checks whether the state buffer has been
 212    * provided by the user, (in which case we assume that it was properly
 213    * initialized, by calling initstate()), or if the default buffer is
 214    * in use...
 215    */
 216   if( state == (struct state *)(default_state) )
 217   { /* ...and we anticipate that, if in its default application start-up
 218      * state, this buffer may require implicit initialization; when this
 219      * is the case...
 220      */
 221     if( state->limit == 0 )
 222     { /* ...as indicated by no buffer size limit having been assigned,
 223        * we initialize it now, as if by calling initstate() with buffer
 224        * size specified appropriately, as 128 bytes, and seed of 1.
 225        */
 226       state->limit = normalized_cycle( 128 );
 227       initialize_state_data( 1 );
 228     }
 229   }
 230   else if( state->limit == normalized_cycle( 8 ) )
 231   { /* The state buffer has been explicitly initialized already, using
 232      * the minimum allowed length of 8 bytes; in this case, we use a
 233      * simple linear congruential PRNG, with multiplier of 1103515245,
 234      * and increment of 12345, masked to ensure that the result fits
 235      * within the range of positive 32-bit signed integer values.
 236      */
 237     const unsigned long long increment = 12345ULL;
 238     const unsigned long long multiplier = 1103515245ULL;
 239     const unsigned int result_mask = ((1ULL << 31) - 1);
 240     unsigned long long tmp = multiplier * state->data[0] + increment;
 241     return (long)(state->data[0] = (int)(tmp & result_mask));
 242   }
 243   /* For any state buffer, longer than the 8-byte minimum, update the
 244    * buffer state, and extract the next pseudo-random number to return,
 245    * scaling to discard its least significant bit.
 246    */
 247   return (long)((unsigned long)(update_state()) >> 1);
 248 }
 249
 250 /* $RCSfile$: end of file */