Re: The POSIX random-number generator

This is a multi-part message in MIME format.
--------------562AFB5433C889B886E770E9
Content-Type: text/plain; charset=utf-8
Content-Transfer-Encoding: 8bit

On 02/04/2021 22:36, Keith Marshall wrote:
> On 28/03/2021 15:43, Anton Shepelev wrote:
>> I cannot build the last version of the Netpbm package because it now
>> depends on the random() and srandom() functions from POSIX:>
>>              https://manned.org/random.3
>>
>> Will the maintainers consider providing them in MinGW to simplify the
>> building of code that depends on basic POSIX facilities?
> 
> Sorry, but my initial inclination is to say no ...

Notwithstanding my reluctance to support this, FWIW I did find

   https://www.mathstat.dal.ca/~selinger/random/

from which, with hints from the GNULIB random_r.c implementation[1],
I've derived the attached patch.

While I still will not attest to the statistical reliability of this, if
you would care to submit a formal feature request[2], and there is
sufficient interest expressed, (by way of up-votes), I may consider
incorporation into libmingwex.a, subject to a requirement for the caller
to define either _XOPEN_SOURCE >= 500, or _BSD_SOURCE, before including
<stdlib.h>, to expose the API declarations.

[1] https://git.savannah.gnu.org/cgit/gnulib.git/tree/lib/random_r.c
Although GNULIB modules are normally licensed under LGPL, (which is NOT
compatible with the licence requirements for libmingwrt.a), this file
also exhibits a BSD licence, (which IS compatible).  Thus, I believe
that I am justified in incorporating hints as to suitable constant
definitions, under the terms of the BSD licence.

[2] https://mingw.osdn.io/index.html?page=contact.html#feature-request

-- 
Regards,
Keith.

Public key available from keys.gnupg.net
Key fingerprint: C19E C018 1547 DE50 E1D4 8F53 C0AD 36C6 347E 5A3F

--------------562AFB5433C889B886E770E9
Content-Type: text/x-patch; charset=UTF-8;
 name="posix-prng-implementation.patch"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
 filename="posix-prng-implementation.patch"

# HG changeset patch
# Parent  dedaf4b5585b853b055f333d18a1d85fbacf6c5f
Implement a POSIX.1 compliant pseudo-random number generator.

* mingwex/math/random.c: New file; it implements...
(__mingw_random, __mingw_srandom, __mingw_initstate)
(__mingw_setstate): ...each of these API functions.
* include/stdlib.h [_XOPEN_SOURCE >=3D 500]: Declare them; map...
(random, srandom, initstate, setstate): ...these aliases.
[_BSD_SOURCE]: Likewise.

* Makefile.in (libmingwex.a): Require random.$OBJEXT

diff --git a/mingwrt/Makefile.in b/mingwrt/Makefile.in
--- a/mingwrt/Makefile.in
+++ b/mingwrt/Makefile.in
@@ -4,12 +4,12 @@
 # Makefile template for MinGW.org Runtime Library Package
=20
 PACKAGE_TARNAME :=3D @PACKAGE_TARNAME@
 PACKAGE_VERSION :=3D @PACKAGE_VERSION@
=20
-# Written by Keith Marshall <[email protected]>
-# Copyright (C) 2014-2020, MinGW.org Project
+# Written by Keith Marshall <[email protected]>
+# Copyright (C) 2014-2021, MinGW.org Project
 #
 #
 # Permission is hereby granted, free of charge, to any person obtaining =
a
 # copy of this software and associated documentation files (the "Softwar=
e"),
 # to deal in the Software without restriction, including without limitat=
ion
@@ -420,10 +420,15 @@ libmingwex.a: $(addsuffix .$(OBJEXT), co
   powf powl powi powif powil remainder remainderf remainderl remquo remq=
uof \
   remquol rint rintf rintl round roundf roundl scalbn scalbnf scalbnl si=
gnbit \
   signbitf signbitl sqrtf sqrtl tgamma tgammaf tgammal trunc truncf trun=
cl \
   x87cvt x87cvtf x87log x87log1p x87pow x87remquo)
=20
+# An experimental implementation of the POSIX.1-1990 random() API;
+# once again, this is delivered in libmingwex.a
+#
+libmingwex.a: random.$(OBJEXT)
+
 # Replacement I/O functions in libmingwex.a, providing better POSIX
 # compatibility than their Microsoft equivalents.
 #
 vpath %.c ${mingwrt_srcdir}/mingwex/stdio
 libmingwex.a: $(addsuffix .$(OBJEXT), btowc fprintf fseeki64 ftelli64 \
diff --git a/mingwrt/include/stdlib.h b/mingwrt/include/stdlib.h
--- a/mingwrt/include/stdlib.h
+++ b/mingwrt/include/stdlib.h
@@ -5,11 +5,12 @@
  * associated macros, and manifest constant definitions.
  *
  * $Id$
  *
  * Written by Colin Peters <[email protected]>
- * Copyright (C) 1997-2009, 2011, 2014-2016, 2018, 2020, MinGW.org Proje=
ct.
+ * Copyright (C) 1997-2009, 2011, 2014-2016, 2018, 2020, 2021,
+ *   MinGW.org Project.
  *
  *
  * Permission is hereby granted, free of charge, to any person obtaining=20
a
  * copy of this software and associated documentation files (the "Softwa=
re"),
  * to deal in the Software without restriction, including without limita=
tion
@@ -539,10 +540,37 @@ wchar_t *_wfullpath (wchar_t *, const wc
 _CRTIMP __cdecl __MINGW_NOTHROW  int rand_s (unsigned int *);
=20
 #endif	/* Win-Vista || MSVCR80.DLL || later */
 #endif	/* _CRT_RAND_S enabled */
=20
+#if _XOPEN_SOURCE >=3D 500 || defined _BSD_SOURCE
+/* The POSIX.1-1990 Extended Systems Interface specification defines ano=
ther
+ * alternative to rand().  While the statistical quality of this alterna=
tive
+ * cannot be assured, this implementation in libmingwex.a may be useful.=

+ */
+__cdecl __MINGW_NOTHROW  long __mingw_random (void);
+__cdecl __MINGW_NOTHROW  void __mingw_srandom (unsigned int);
+__cdecl __MINGW_NOTHROW  char *__mingw_initstate (unsigned int, char *, =
size_t);
+__cdecl __MINGW_NOTHROW  char *__mingw_setstate (char *);
+
+__CRT_ALIAS __JMPSTUB__(( FUNCTION =3D random ))
+__cdecl __MINGW_NOTHROW  long random (void){ return __mingw_random(); }
+
+__CRT_ALIAS __JMPSTUB__(( FUNCTION =3D srandom ))
+__cdecl __MINGW_NOTHROW  void srandom (unsigned int __seed)
+{ __mingw_srandom (__seed); }
+
+__CRT_ALIAS __JMPSTUB__(( FUNCTION =3D initstate ))
+__cdecl __MINGW_NOTHROW  char *initstate (unsigned int __seed, char *__b=
uf, size_t __len)
+{ return __mingw_initstate (__seed, __buf, __len); }
+
+__CRT_ALIAS __JMPSTUB__(( FUNCTION =3D setstate ))
+__cdecl __MINGW_NOTHROW  char *setstate (char *__buf)
+{ return __mingw_setstate (__buf); }
+
+#endif	/* _XOPEN_SOURCE >=3D 500 || _BSD_SOURCE */
+
 _CRTIMP __cdecl __MINGW_NOTHROW  void abort (void) __MINGW_ATTRIB_NORETU=
RN;
 _CRTIMP __cdecl __MINGW_NOTHROW  void exit (int) __MINGW_ATTRIB_NORETURN=
;
=20
 /* Note: this is in startup code, not imported directly from the runtime=20
DLL
  */
diff --git a/mingwrt/mingwex/math/random.c b/mingwrt/mingwex/math/random.=
c
new file mode 100644
--- /dev/null
+++ b/mingwrt/mingwex/math/random.c
@@ -0,0 +1,250 @@
+/*
+ * random.c
+ *
+ * Implementation of a (mostly) POSIX.1-1990 conforming pseudo-random
+ * number generating API, for use with MinGW applications.
+ *
+ * $Id$
+ *
+ * Written by Keith Marshall <[email protected]>
+ * Copyright (C) 2021, MinGW.org Project
+ *
+ *
+ * This is free software.  Permission is granted to copy, modify and
+ * redistribute this software, under the provisions of the GNU General
+ * Public License, Version 3, (or, at your option, any later version),
+ * as published by the Free Software Foundation; see the file COPYING
+ * for licensing details.
+ *
+ * Note, in particular, that this software is provided "as is", in the
+ * hope that it may prove useful, but WITHOUT WARRANTY OF ANY KIND; not
+ * even an implied WARRANTY OF MERCHANTABILITY, nor of FITNESS FOR ANY
+ * PARTICULAR PURPOSE.  Under no circumstances will the author, or the
+ * MinGW Project, accept liability for any damages, however caused,
+ * arising from the use of this software.
+ *
+ *
+ * This pseudo-random number generation suite has been derived from a
+ * discourse on the GLIBC implementation, by Peter Selinger:
+ *
+ *   https://www.mathstat.dal.ca/~selinger/random/
+ *
+ * It is believed that the output from the PRNG will closely mimic that
+ * from the GLIBC implementation; however, neither the author of this
+ * implementation, nor the MinGW.org Project, offer any assurance as to
+ * the statistical quality of the generated number sequence.
+ *
+ *
+ * Note: this implementation has been provided for optional compliance
+ * with the POSIX.1-1990 Extended Systems Interface.
+ */
+#define _XOPEN_SOURCE  500
+
+#include <stdint.h>
+#include <stdlib.h>
+#include <errno.h>
+
+struct state
+{ /* Internal structure, used to map the PRNG state buffer from its
+   * arbitrarily specified sequence of 8, 32, 64, 128, or 256 bytes,
+   * to an array of int32_t entities; the first entry is used as a
+   * control sturcture, comprising three 8-bit indices, while the
+   * remaining entries are used to represent PRNG state.
+   */
+  struct
+  { uint32_t	phase:8;	/* current state data cycle pointer  */
+    uint32_t	shift:8;	/* associated offset data pointer    */
+    uint32_t	limit:8;	/* total count of state data entries */
+  };
+  int32_t	data[]; 	/* array of state data entries	     */
+};
+
+/* Provide a local state buffer, of length 128 bytes, (equivalent
+ * to 32 int32_t entries), for use when no alternative buffer has
+ * been specified, prior to calling random().
+ */
+static int32_t default_state[32] =3D { 0 };
+static struct state *state =3D (struct state *)(default_state);
+
+/* The following inline local function is provided to facilitate
+ * setting of "errno", on abnormal return from any other function.
+ */
+static __inline__ __attribute__((__always_inline__))
+intptr_t errout( int code, intptr_t retval ){ errno =3D code; return ret=
val; }
+
+static char *assign_state( char *buf, char *prev )
+{ /* A local helper function, to attach a new state buffer to the
+   * PRNG; aborts the assignment, if the specified buffer is NULL,
+   * but otherwise, performs no buffer validation.
+   */
+  if( buf =3D=3D NULL ) return (char *)(errout( EINVAL, (intptr_t)(NULL)=20
));
+  state =3D (struct state *)(buf);
+  return prev;
+}
+
+char *__mingw_setstate( char *buf )
+{ /* Public API entry point, exhibiting POSIX.1-1990 semantics,
+   * for accessing the preceding buffer assignment function; the
+   * "buf" argument is assumed to point to an existing non-NULL
+   * buffer, which had been previously initialized by a prior
+   * initstate() call, (and possibly already used), but, other
+   * that checking that it is not NULL, no validation of its
+   * initialization state is performed.
+   */
+  return assign_state( buf, (char *)(state) );
+}
+
+static size_t normalized_cycle( size_t len )
+{ /* A local helper function, to establish the effective number
+   * of entries in the PRNG state data array, as a dependency of
+   * the total byte count specified for the state buffer; note
+   * that the compiler may be able to compute the resultant at
+   * compile time, when called with a constant "len" argument,
+   * and so optimize away such calls.
+   */
+  size_t cycle_counter =3D (len >=3D 32) ? 256 : 8;
+  while( cycle_counter > len ) cycle_counter >>=3D 1;
+  return (cycle_counter >> 2) - 1;
+}
+
+static long update_state( void )
+{ /* A local helper function, to update the content of the active
+   * state buffer with each successive call of random().
+   */
+  if( state->phase =3D=3D state->limit ) state->phase =3D 0;
+  else if( state->shift =3D=3D state->limit ) state->shift =3D 0;
+  return (state->data[state->shift++] +=3D state->data[state->phase++]);=

+}
+
+static void initialize_state_data( unsigned int seed )
+{ /* A local helper function, called after initialization of the
+   * state buffer length count, to populate the state data array.
+   *
+   * In all cases, we record the "seed" value, as the first entry
+   * in the state data array, (rejecting zero, for which we force
+   * substitution of 1).
+   */
+  state->data[0] =3D (seed !=3D 0) ? seed : 1;
+
+  /* For a minimum-length state buffer, of only 8 bytes, there is
+   * no more state data to be recorded...
+   */
+  if( state->limit > 1 )
+  { /* ...but for any of the larger permitted buffer sizes, we must
+     * populate the buffer using a multiplicative sequence of 31-bit
+     * integers; the multiplication is performed in a 64-bit space,
+     * then reduced modulo the Mersenne prime of order 31; set up a
+     * collection of constants, pertinent to the computation...
+     */
+    const unsigned long order =3D 31;
+    const unsigned long divisor =3D (unsigned long)((1ULL << order) - 1)=
;
+    const unsigned long long multiplier =3D 16807ULL;
+
+    /* ...then apply to each data array entry, in turn...
+     */
+    for( seed =3D 1; seed < state->limit; seed++ )
+    { /* ...computing the sequence of products, in 64-bit space...
+       */
+      unsigned long long tmp =3D multiplier * state->data[seed - 1];
+
+      /* ...then reduce for storage; note that, since the divisor
+       * is a Mersenne prime, we can perform the modulo division by
+       * use of the more computationally efficient mask, shift and
+       * addition, followed by a conditional correction, in case
+       * of overflow into the sign bit.
+       */
+      state->data[seed] =3D (int)(tmp =3D (tmp & divisor) + (tmp >> orde=
r));
+      if( state->data[seed] < 0 ) state->data[seed] +=3D divisor;
+    }
+    /* We must also initialize the phase, and shift indices, for
+     * the PRNG feed-back loop...
+     */
+    state->phase =3D 0;
+    switch( state->limit )
+    { /* ...noting that the shift interval is chosen differently,
+       * depending on the limiting length of the sequence.
+       */
+      case 63: case 15: state->shift =3D 1; break;
+      case 31: case  7: state->shift =3D 3; break;
+    }
+    /* Finally, run the PRNG through 10 update cycles for each and
+     * every entry in the state data array, but without generating
+     * any actual pseudo-random number output.
+     */
+    seed =3D ((state->limit << 2) + state->limit) << 1;
+    while( seed-- > 0 ) update_state();
+  }
+}
+
+char *__mingw_initstate( unsigned int seed, char *buf, size_t len )
+{ /* Public API entry point, implementing the POSIX.1-1990 initstate()
+   * function; aborts if "buf" is passed as a NULL pointer, or if "len"
+   * is less than the minimum allowed 8 bytes; otherwise, assigns "buf"
+   * as a new PRNG state buffer...
+   */
+  char *prev =3D assign_state( buf, (len >=3D 8) ? (char *)(state) : NUL=
L );
+  if( prev !=3D NULL )
+  {
+    /* ...initializes it to the specified "len", and as if srandom()
+     * has been called with "seed" as argument...
+     */
+    state->limit =3D normalized_cycle( len );
+    initialize_state_data( seed );
+  }
+  /* ...then returns a pointer to the original state buffer.
+   */
+  return prev;
+}
+
+void __mingw_srandom( unsigned int seed )
+{ /* Public API entry point, implementing the POSIX-1.1990 srandom()
+   * function; this sets the PRNG to a known initial state, which is
+   * dependent on the "seed" value specified.
+   */
+  if( (state =3D=3D (struct state *)(default_state)) && (state->limit =3D=
=3D 0) )
+    state->limit =3D normalized_cycle( 128 );
+  initialize_state_data( seed );
+}
+
+long __mingw_random( void )
+{ /* Public API entry point, implementing the POSIX-1.1990 random()
+   * function; this initially checks whether the state buffer has been
+   * provided by the user, (in which case we assume that it was properly=

+   * initialized, by calling initstate()), or if the default buffer is
+   * in use...
+   */
+  if( state =3D=3D (struct state *)(default_state) )
+  { /* ...and we anticipate that, if in its default application start-up=

+     * state, this buffer may require implicit initialization; when this=

+     * is the case...
+     */
+    if( state->limit =3D=3D 0 )
+    { /* ...as indicated by no buffer size limit having been assigned,
+       * we initialize it now, as if by calling initstate() with buffer
+       * size specified appropriately, as 128 bytes, and seed of 1.
+       */
+      state->limit =3D normalized_cycle( 128 );
+      initialize_state_data( 1 );
+    }
+  }
+  else if( state->limit =3D=3D normalized_cycle( 8 ) )
+  { /* The state buffer has been explicitly initialized already, using
+     * the minimum allowed length of 8 bytes; in this case, we use a
+     * simple linear congruential PRNG, with multiplier of 1103515245,
+     * and increment of 12345, masked to ensure that the result fits
+     * within the range of positive 32-bit signed integer values.
+     */
+    const unsigned long long increment =3D 12345ULL;
+    const unsigned long long multiplier =3D 1103515245ULL;
+    const unsigned int result_mask =3D ((1ULL << 31) - 1);
+    unsigned long long tmp =3D multiplier * state->data[0] + increment;
+    return (long)(state->data[0] =3D (int)(tmp & result_mask));
+  }
+  /* For any state buffer, longer than the 8-byte minimum, update the
+   * buffer state, and extract the next pseudo-random number to return,
+   * scaling to discard its least significant bit.
+   */
+  return (long)((unsigned long)(update_state()) >> 1);
+}
+
+/* $RCSfile$: end of file */

--------------562AFB5433C889B886E770E9
Content-Type: text/plain; charset="utf-8"
MIME-Version: 1.0
Content-Transfer-Encoding: base64
Content-Disposition: inline

X19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX18KTWluR1ctVXNl
cnMgbWFpbGluZyBsaXN0Ck1pbkdXLVVzZXJzQGxpc3RzLm9zZG4ubWUKClRoaXMgbGlzdCBvYnNl
cnZlcyB0aGUgUG9zdGluZyBFdGlxdWV0dGUsIGFzIGRlc2NyaWJlZCBhdCBodHRwczovL21pbmd3
Lm9zZG4uaW8vaW5kZXguaHRtbD9wYWdlPW1haWxpbmcuaHRtbCNsaXN0LWV0aXF1ZXR0ZS4KV2Ug
YXNrIHRoYXQgeW91IGJlIHBvbGl0ZSBhbmQgZG8gdGhlIHNhbWUuICBEaXNyZWdhcmQgZm9yIHRo
ZSBsaXN0IGV0aXF1ZXR0ZSBtYXkgY2F1c2UgeW91ciBhY2NvdW50IHRvIGJlIG1vZGVyYXRlZC4K
Cl9fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fX19fCllvdSBtYXkg
Y2hhbmdlIHlvdXIgTWluR1cgQWNjb3VudCBPcHRpb25zIG9yIHVuc3Vic2NyaWJlIGF0OgpodHRw
czovL2xpc3RzLm9zZG4ubWUvbWFpbG1hbi9saXN0aW5mby9taW5ndy11c2VycwpBbHNvOiBtYWls
dG86bWluZ3ctdXNlcnMtcmVxdWVzdEBsaXN0cy5vc2RuLm1lP3N1YmplY3Q9dW5zdWJzY3JpYmU=

--------------562AFB5433C889B886E770E9--