Re: Completeness status of musl

[Szabolcs Nagy <nsz@port70.net>]

[-- Attachment #1: Type: text/plain, Size: 3110 bytes --]

* Rich Felker <dalias@aerifal.cx> [2011-06-23 02:54:52 -0400]:

> On Tue, Jun 21, 2011 at 03:46:40AM +0200, Szabolcs Nagy wrote:
> > * Rich Felker <dalias@aerifal.cx> [2011-06-08 19:58:10 -0400]:
> > > 
> > > Importing external code, with major cleanups or rewrite:
> > > - PRNG (random)
> > 
> > i looked into this one
> > these functions are rather useless
> > imho the posix api definition is insane
> > http://pubs.opengroup.org/onlinepubs/9699919799/functions/random.html
> > 
> > i wrote something anyway
> > see comments in the code:
> > http://port70.net/~nsz/musl/prng
> 
> I like this. A couple questions/requests:
> 
> 1. Is this simple code really equivalent (aside from the missing
> non-default-size code) to the ugly BSD code that's 100x bigger and
> full of undefined behavior? :)
> 

the historical bsd implementation used another lcg
for both the <32byte statebuffer and seeding

"modern" implementations seems to use the
 x = (7^5 * x) mod (2^31 - 1)
park-miller generator for seeding
(it makes the generator less predictable,
with the lcg everything is mod 2^k so
the numbers 'linearly depend' on the seed)
(glibc freebsd openbsd.. all use the same code here)
(they use overflowing signed ints and signed >>1 for whatever reason)
(they discard the first 10*n rands which i forgot to do)

(freebsd seems to use the park-miller generator even for
random() in case of small state buffer, which is bad as it
will produce uniform random in [0,2^31-2] not in [0,2^31-1])

i used lcg seeding for simplicity
the lcg seeding is bad, but the generator
does not have good statistical properties
and the seeding method won't change that..

now i made all the compatibility changes (attached)
the following produces identical output with glibc vs my code:
#include <stdio.h>
#include <stdlib.h>
int main() {
        char s[256];
        int i,j;
        for (j = 8; j <= 256; j*=2) {
                initstate(1, s, j);
                for (i = 0; i < 10; i++)
                        printf("%d %d %ld\n", j, i, random());
        }
        return 0;
}

(if compatibility does not matter then a 64bit lcg
could do the same job without much fuss:
 static unsigned long long init[] = {1};
 unsigned long long *x = init;
 long random(void){
    *x = (6364136223846793005* *x + 1442695040888963407);
    return *x >> 33;
 }
etc)

> 2. If you'd like me to include this in musl, can you please license it
> (either LGPL v2.1+ or any LGPL-compatible/less-restrictive license
> such as BSD).
> 

ok
license is
"permission to use, copy, modify, and/or distribute this code
for any purpose with or without fee is hereby granted.
there is no warranty."

so feel free to alter the code
eg the parkmiller code might be better to do with uint64 multiplication
instead of the div/mod hack to avoid overflow:
 return (unsigned long long)16807*x % 0x7fffffff;

> 3. For future code contributions, please attach the code/patches to
> the email rather than just including a link. That way the code being
> referred to/discussed is visible even if the version on the website
> changes or the website it taken down.
> 

ok

[-- Attachment #2: random.c --]
[-- Type: text/x-csrc, Size: 1918 bytes --]

#include <stdlib.h>

/*
32bit unsigned is asssumed
traditional bsd random implementation with park-miller seeding
lagged fibonacci generator params: (7,3), (15,1), (31,3), (63,1)
does not have good statistical properties
*/

static unsigned init[] = {
0x00000000,0x991539b1,0x16a5bce3,0x6774a4cd,
0x3e01511e,0x4e508aaa,0x61048c05,0xf5500617,
0x846b7115,0x6a19892c,0x896a97af,0xdb48f936,
0x14898454,0x37ffd106,0xb58bff9c,0x59e17104,
0xcf918a49,0x09378c83,0x52c7a471,0x8d293ea9,
0x1f4fc301,0xc3db71be,0x39b44e1c,0xf8a44ef9,
0x4c8b80b1,0x19edc328,0x87bf4bdd,0xc9b240e5,
0xe9ee4b1b,0x4382aee7,0x535b6b41,0xf3bec5da};

static int n = 31;
static int i = 3;
static int j = 0;
static unsigned *x = init+1;

static unsigned lcg(unsigned x) {
	return (1103515245*x + 12345) & 0x7fffffff;
}

static unsigned parkmiller(unsigned x) {
	/* x = 7^5 x mod (2^31-1) */
	x = 16807*(x%127773) - 2836*(x/127773);
	if (x & 0x80000000)
		x += 0x7fffffff;
	return x;
}

static void *savestate() {
	x[-1] = ((unsigned)n<<16)|(i<<8)|j;
	return x-1;
}

static void loadstate(unsigned *state) {
	x = state+1;
	n = x[-1]>>16;
	i = (x[-1]>>8)&0xff;
	j = x[-1]&0xff;
}

void srandom(unsigned seed) {
	int k;

	i = n == 31 || n == 7 ? 3 : 1;
	j = 0;
	x[0] = seed ? seed : 1;
	for (k = 1; k < n; k++)
		x[k] = parkmiller(x[k-1]);
	for (k = 0; k < 10*n; k++)
		random();
}

char *initstate(unsigned seed, char *state, size_t size) {
	void *old = savestate();
	if (size < 8)
		return 0;
	else if (size < 32)
		n = 0;
	else if (size < 64)
		n = 7;
	else if (size < 128)
		n = 15;
	else if (size < 256)
		n = 31;
	else
		n = 63;
	x = (unsigned*)state + 1;
	srandom(seed);
	return old;
}

char *setstate(char *state) {
	void *old = savestate();
	loadstate((unsigned*)state);
	return old;
}

long random(void) {
	unsigned k;

	if (n == 0)
		return x[0] = lcg(x[0]);
	x[i] += x[j];
	k = x[i]>>1;
	if (++i == n)
		i = 0;
	if (++j == n)
		j = 0;
	return k;
}