rand(3): Replace implementation with one backed by random(3) algorithm

rand(3)'s standard C API is extremely limiting, but we can do better
than the historical 32-bit state Park-Miller LCG we've shipped since
2001: r73156.

The justification provided at the time for not using random(3) was that
rand_r(3) could not be made to use the same algorithm.  That is still
true.  However, the irrelevance of rand_r(3) is increasingly obvious.
Since that time, POSIX has marked the interface obsolescent.  rand_r(3)
never became part of the standard C library.  If not for API
compatibility reasons, I would just remove rand_r(3) entirely.

So, I do not believe it is a problem for rand_r(3) and rand(3) to
diverge.

The 12 ABI is maintained with compatibility definitions, but this
revision does subtly change the API of rand(3).  The sequences of
pseudorandom numbers produced in programs built against new versions of
libc will differ from programs built against prior versions of libc.

Reviewed by:	kevans, markm
MFC after:	no
Relnotes:	yes
Differential Revision:	https://reviews.freebsd.org/D23290

include/stdlib.h

@@ -73,7 +73,11 @@ typedef struct {
 #define	EXIT_FAILURE	1
 #define	EXIT_SUCCESS	0
 
-#define	RAND_MAX	0x7ffffffd
+/*
+ * I.e., INT_MAX; rand(3) returns a signed integer but must produce output in
+ * the range [0, RAND_MAX], so half of the possible output range is unused.
+ */
+#define	RAND_MAX	0x7fffffff
 
 __BEGIN_DECLS
 #ifdef _XLOCALE_H_

lib/libc/stdlib/Symbol.map

@@ -54,8 +54,6 @@ FBSD_1.0 {
 	radixsort;
 	sradixsort;
 	rand_r;
-	rand;
-	srand;
 	srandom;
 	srandomdev;
 	initstate;
@@ -125,6 +123,8 @@ FBSD_1.5 {
 
 FBSD_1.6 {
 	qsort_s;
+	rand;
+	srand;
 };
 
 FBSDprivate_1.0 {

lib/libc/stdlib/rand.3

@@ -32,7 +32,7 @@
 .\"     @(#)rand.3	8.1 (Berkeley) 6/4/93
 .\" $FreeBSD$
 .\"
-.Dd December 14, 2019
+.Dd February 1, 2020
 .Dt RAND 3
 .Os
 .Sh NAME
@@ -59,49 +59,52 @@ Applications which require unpredictable random numbers should use
 instead.
 .Ef
 .Pp
-These interfaces are obsoleted by
-.Xr random 3 .
-.Pp
 The
 .Fn rand
 function computes a sequence of pseudo-random integers in the range
 of 0 to
-.Dv RAND_MAX
-(as defined by the header file
-.In stdlib.h ) .
+.Dv RAND_MAX ,
+inclusive.
 .Pp
 The
 .Fn srand
-function sets its argument
+function seeds the algorithm with the
 .Fa seed
-as the seed for a new sequence of
-pseudo-random numbers to be returned by
-.Fn rand .
-These sequences are repeatable by calling
+parameter.
+Repeatable sequences of
+.Fn rand
+output may be obtained by calling
 .Fn srand
-with the same seed value.
+with the same
+.Fa seed .
+.Fn rand
+is implicitly initialized as if
+.Fn srand "1"
+had been invoked explicitly.
 .Pp
-If no
-.Fa seed
-value is provided, the functions are automatically
-seeded with a value of 1.
-.Pp
-The
+In
+.Fx 13 ,
+.Fn rand
+is implemented using the same 128-byte state LFSR generator algorithm as
+.Xr random 3 .
+However, the legacy
 .Fn rand_r
-function
-provides the same functionality as
-.Fn rand .
-A pointer to the context value
-.Fa ctx
-must be supplied by the caller.
-.Pp
-For better generator quality, use
-.Xr random 3
-or
-.Xr lrand48 3 .
+function is not (and can not be, because of its limited
+.Fa *ctx
+size).
+.Fn rand_r
+implements the historical, poor-quality Park-Miller 32-bit LCG and should not
+be used in new designs.
+.Sh IMPLEMENTATION NOTES
+Since
+.Fx 13 ,
+.Fn rand
+is implemented with the same generator as
+.Xr random 3 ,
+so the low-order bits should no longer be significantly worse than the
+high-order bits.
 .Sh SEE ALSO
 .Xr arc4random 3 ,
-.Xr lrand48 3 ,
 .Xr random 3 ,
 .Xr random 4
 .Sh STANDARDS
@@ -115,5 +118,32 @@ conform to
 .Pp
 The
 .Fn rand_r
-function is marked as obsolescent in POSIX and may be removed in a future
-revision of the standard.
+function is not part of
+.St -isoC
+and is marked obsolescent in
+.St -p1003.1-2008 .
+It may be removed in a future revision of POSIX.
+.Sh CAVEATS
+Prior to
+.Fx 13 ,
+.Fn rand
+used the historical Park-Miller generator with 32 bits of state and produced
+poor quality output, especially in the lower bits.
+.Fn rand
+in earlier versions of
+.Fx ,
+as well as other standards-conforming implementations, may continue to produce
+poor quality output.
+.Pp
+.Em These functions should not be used in portable applications that want a
+.Em high quality or high performance pseudorandom number generator .
+One possible replacement,
+.Xr random 3 ,
+is portable to Linux — but it is not especially fast, nor standardized.
+.Pp
+If broader portability or better performance is desired, any of the widely
+available and permissively licensed SFC64/32, JSF64/32, PCG64/32, or SplitMix64
+algorithm implementations may be embedded in your application.
+These algorithms have the benefit of requiring less space than
+.Xr random 3
+and being quite fast (in header inline implementations).

lib/libc/stdlib/rand.c

@@ -40,11 +40,60 @@ __FBSDID("$FreeBSD$");
 #include "namespace.h"
 #include <sys/param.h>
 #include <sys/sysctl.h>
+#include <assert.h>
 #include <stdbool.h>
 #include <stdlib.h>
 #include <syslog.h>
 #include "un-namespace.h"
 
+#include "random.h"
+
+/*
+ * Implement rand(3), the standard C PRNG API, using the non-standard but
+ * higher quality random(3) implementation and the same size 128-byte state
+ * LFSR as the random(3) default.
+ *
+ * It turns out there are portable applications that want a PRNG but are too
+ * lazy to use better-but-nonstandard interfaces like random(3), when
+ * available, and too lazy to import higher-quality and faster PRNGs into their
+ * codebase (such as any of SFC, JSF, 128-bit LCGs, PCG, or Splitmix64).
+ *
+ * Since we're stuck with rand(3) due to the C standard, we can at least have
+ * it produce a relatively good PRNG sequence using our existing random(3)
+ * LFSR.  The random(3) design is not particularly fast nor compact, but it has
+ * the advantage of being the one already in the tree.
+ */
+static struct __random_state *rand3_state;
+
+static void
+initialize_rand3(void)
+{
+	int error;
+
+	rand3_state = allocatestate(TYPE_3);
+	error = initstate_r(rand3_state, 1, rand3_state->rst_randtbl, BREAK_3);
+	assert(error == 0);
+}
+
+int
+rand(void)
+{
+	if (rand3_state == NULL)
+		initialize_rand3();
+	return ((int)random_r(rand3_state));
+}
+
+void
+srand(unsigned seed)
+{
+	if (rand3_state == NULL)
+		initialize_rand3();
+	srandom_r(rand3_state, seed);
+}
+
+/*
+ * FreeBSD 12 and prior compatibility implementation of rand(3).
+ */
 static int
 do_rand(unsigned long *ctx)
 {
@@ -71,7 +120,9 @@ do_rand(unsigned long *ctx)
 	return (x);
 }
 
-
+/*
+ * Can't fix this garbage; too little state.
+ */
 int
 rand_r(unsigned *ctx)
 {
@@ -84,21 +135,23 @@ rand_r(unsigned *ctx)
 	return (r);
 }
 
-
 static u_long next = 1;
 
+int __rand_fbsd12(void);
 int
-rand(void)
+__rand_fbsd12(void)
 {
 	return (do_rand(&next));
 }
+__sym_compat(rand, __rand_fbsd12, FBSD_1.0);
 
+void __srand_fbsd12(unsigned seed);
 void
-srand(unsigned seed)
+__srand_fbsd12(unsigned seed)
 {
 	next = seed;
 }
-
+__sym_compat(srand, __srand_fbsd12, FBSD_1.0);
 
 void __sranddev_fbsd12(void);
 void

lib/libc/stdlib/random.3

@@ -28,7 +28,7 @@
 .\"     @(#)random.3	8.1 (Berkeley) 6/4/93
 .\" $FreeBSD$
 .\"
-.Dd January 20, 2020
+.Dd February 1, 2020
 .Dt RANDOM 3
 .Os
 .Sh NAME
@@ -74,8 +74,7 @@ The period of this random number generator is very large, approximately
 .Pp
 If initialized with less than 32 bytes of state,
 .Fn random
-uses the same poor-quality Park-Miller LCG as
-.Xr rand 3 .
+uses the poor-quality 32-bit Park-Miller LCG.
 .Pp
 The
 .Fn random
@@ -85,9 +84,6 @@ functions are analagous to
 .Xr rand 3
 and
 .Xr srand 3 .
-The difference is that
-.Xr rand 3
-is a worse pseudo-random number generator.
 .Pp
 Like
 .Xr rand 3 ,

lib/libc/stdlib/random.c

@@ -104,48 +104,13 @@ __FBSDID("$FreeBSD$");
  * byte buffer it is about 5 percent faster.
  */
 
-/*
- * For each of the currently supported random number generators, we have a
- * break value on the amount of state information (you need at least this
- * many bytes of state info to support this random number generator), a degree
- * for the polynomial (actually a trinomial) that the R.N.G. is based on, and
- * the separation between the two lower order coefficients of the trinomial.
- */
-#define	TYPE_0		0		/* linear congruential */
-#define	BREAK_0		8
-#define	DEG_0		0
-#define	SEP_0		0
-
-#define	TYPE_1		1		/* x**7 + x**3 + 1 */
-#define	BREAK_1		32
-#define	DEG_1		7
-#define	SEP_1		3
-
-#define	TYPE_2		2		/* x**15 + x + 1 */
-#define	BREAK_2		64
-#define	DEG_2		15
-#define	SEP_2		1
-
-#define	TYPE_3		3		/* x**31 + x**3 + 1 */
-#define	BREAK_3		128
-#define	DEG_3		31
-#define	SEP_3		3
-
-#define	TYPE_4		4		/* x**63 + x + 1 */
-#define	BREAK_4		256
-#define	DEG_4		63
-#define	SEP_4		1
-
-/*
- * Array versions of the above information to make code run faster --
- * relies on fact that TYPE_i == i.
- */
-#define	MAX_TYPES	5		/* max number of types above */
-
 #define NSHUFF 50       /* to drop some "seed -> 1st value" linearity */
 
 static const int degrees[MAX_TYPES] =	{ DEG_0, DEG_1, DEG_2, DEG_3, DEG_4 };
-static const int seps [MAX_TYPES] =	{ SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 };
+static const int seps[MAX_TYPES] =	{ SEP_0, SEP_1, SEP_2, SEP_3, SEP_4 };
+static const int breaks[MAX_TYPES] = {
+	BREAK_0, BREAK_1, BREAK_2, BREAK_3, BREAK_4
+};
 
 /*
  * Initially, everything is set up as if from:
@@ -524,3 +489,19 @@ random(void)
 {
 	return (random_r(&implicit));
 }
+
+struct __random_state *
+allocatestate(unsigned type)
+{
+	size_t asize;
+
+	/* No point using this interface to get the Park-Miller LCG. */
+	if (type < TYPE_1)
+		abort();
+	/* Clamp to widest supported variant. */
+	if (type > (MAX_TYPES - 1))
+		type = (MAX_TYPES - 1);
+
+	asize = sizeof(struct __random_state) + (size_t)breaks[type];
+	return (malloc(asize));
+}

lib/libc/stdlib/random.h

@@ -27,6 +27,44 @@
 
 #pragma once
 
+/*
+ * For each of the currently supported random number generators, we have a
+ * break value on the amount of state information (you need at least this
+ * many bytes of state info to support this random number generator), a degree
+ * for the polynomial (actually a trinomial) that the R.N.G. is based on, and
+ * the separation between the two lower order coefficients of the trinomial.
+ */
+#define	TYPE_0		0		/* linear congruential */
+#define	BREAK_0		8
+#define	DEG_0		0
+#define	SEP_0		0
+
+#define	TYPE_1		1		/* x**7 + x**3 + 1 */
+#define	BREAK_1		32
+#define	DEG_1		7
+#define	SEP_1		3
+
+#define	TYPE_2		2		/* x**15 + x + 1 */
+#define	BREAK_2		64
+#define	DEG_2		15
+#define	SEP_2		1
+
+#define	TYPE_3		3		/* x**31 + x**3 + 1 */
+#define	BREAK_3		128
+#define	DEG_3		31
+#define	SEP_3		3
+
+#define	TYPE_4		4		/* x**63 + x + 1 */
+#define	BREAK_4		256
+#define	DEG_4		63
+#define	SEP_4		1
+
+/*
+ * Array versions of the above information to make code run faster --
+ * relies on fact that TYPE_i == i.
+ */
+#define	MAX_TYPES	5		/* max number of types above */
+
 /* A full instance of the random(3) generator. */
 struct __random_state {
 	uint32_t	*rst_fptr;
@@ -40,6 +78,7 @@ struct __random_state {
 	uint32_t	rst_randtbl[];
 };
 
+struct __random_state *allocatestate(unsigned type);
 int initstate_r(struct __random_state *, unsigned, uint32_t *, size_t);
 long random_r(struct __random_state *);
 void srandom_r(struct __random_state *, unsigned);