Date: Tue, 17 Mar 2015 10:07:37 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Pedro Giffuni <pfg@FreeBSD.org> Cc: freebsd-numerics@FreeBSD.org Subject: Re: Random number generators Message-ID: <20150317170737.GA24682@troutmask.apl.washington.edu> In-Reply-To: <F6137E2C-FDF2-46B3-BFC2-1975AFA40951@FreeBSD.org> References: <7CBD7758-9472-4A2E-8065-EC6E68EE8DAB@FreeBSD.org> <20150317060310.GA21975@troutmask.apl.washington.edu> <F6137E2C-FDF2-46B3-BFC2-1975AFA40951@FreeBSD.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Tue, Mar 17, 2015 at 07:10:46AM -0500, Pedro Giffuni wrote: > > > One big issue is saving internal state. IIRC, MT requires 623-bit > > of internal state. KISS requires 4 32-bit int. Thus, if > > you want to reseed the generator, KISS requires far less effort. > > > > Yes, the problem is the that libc requires a single 32 (or 31) bit seed. > Given that restriction, our existing generator is not bad. Enforcing > something better breaks the API and is not really practical to get > crypto-grade randomness for stuff like refreshing a slide in a > presentation anyways. > > The musl libc approach seemed reasonable but I haven?t looked at the > base random generator there (I?ve heard the glibc one is not good at all). > GM's kiss generator has a period of 2**123. The manpage for random(3) claims a period of about 2**35. The rand(3) manpage does not give a period, but I suspect it is well short of 2**123. The code that follows my sig uses. a Lehmer LCG generator to provide the 4 seeds needed for kiss. The Lehmer LCG takes a single 32 bit seed. If reseed() is not called prior to calling kiss(), then a default set of seeds is used. If the argument to reseed() is 0 or 1, then the default set of seeds is used. It should be straight forward to map reseed() to srand() and kiss() to rand(). Do note that range kiss() is (0,UINT_MAX], so one needs to subtract 1 to get [0,UINT_MAX-1) if 0 needs to be in the range. To use this code in preference to random(3) (and in violation of POSIX?), initstate() and setstate() would need to muck with the internal state of kiss(). This is certainly possible, but I don't do it below. -- Steve #include <sys/types.h> #include <sys/time.h> #include <sys/resource.h> #include <stdint.h> #include <stdio.h> /* * The KISS generator requires 4 seeds. Use a simple Lehmer linear * congruential generator (LCG) that requires only a single seed to * generate the 4 KISS seeds. */ static uint64_t lcg_seed = 0; /* Need to remember the state of LCG. */ static uint32_t lehmer_lcg(uint32_t a) { return ((uint64_t)a * 279470273UL) % 4294967291UL; } /* * If reseed() is not called or it is called with an argument of 0 or 1, * then use DEFAULT_SEED as the set of seeds for kiss(). */ #define DEFAULT_SEED {123456789, 362436069, 521288629, 916191069} static const uint32_t default_seed[4] = DEFAULT_SEED; static uint32_t seed[4] = DEFAULT_SEED; uint32_t reseed(uint32_t s) { int i; uint32_t kiss_seed; /* Copy seed given by user. */ lcg_seed = s; if (lcg_seed == 0 || lcg_seed == 1) { for (i = 0; i < 4; i++) seed[i] = default_seed[i]; } else { kiss_seed = lcg_seed; for (i = 0; i < 4; i++) { kiss_seed = lehmer_lcg(kiss_seed); seed[i] = kiss_seed; } } return (lcg_seed); } /* * kiss() returns an integer value in the range of (0, 2**32-1]. The * distribution of pseudorandom numbers should be uniform. The overall * perdiod for this generator exceeds 2**123. */ #define LEFT(k, n) ((k)^((k)<<(n))) #define RGHT(k, n) ((k)^((k)>>(n))) uint32_t kiss(void) { seed[0] = 69069 * seed[0] + 1327217885; seed[1] = LEFT(RGHT(LEFT(seed[1],13),17),5); seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16); seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16); return (seed[0] + seed[1] + (seed[2] << 16) + seed[3]); } #define NUM 5 int main(void) { int i; uint32_t n; struct rusage t0, t1; /* Default seeds. */ printf("Default seeds used.\n"); for (i = 0; i < NUM; i++) printf("%u\n", kiss()); n = reseed(42); printf("\nseed = %u\n", n); for (i = 0; i < NUM; i++) printf("%u\n", kiss()); /* Default seeds. */ n = reseed(0); printf("\nseed = %u\n", n); for (i = 0; i < NUM; i++) printf("%u\n", kiss()); n = reseed(12345); printf("\nseed = %u\n", n); for (i = 0; i < NUM; i++) printf("%u\n", kiss()); /* Default seeds. */ n = reseed(1); printf("\nseed = %u\n", n); for (i = 0; i < NUM; i++) printf("%u\n", kiss()); getrusage(RUSAGE_SELF, &t0); for (i = 0; i < 1000000; i++) n = kiss(); getrusage(RUSAGE_SELF, &t1); { double u; u = (t1.ru_utime.tv_sec - t0.ru_utime.tv_sec); u += (t1.ru_stime.tv_sec - t0.ru_stime.tv_sec); u += 1.e-6*(t1.ru_utime.tv_usec - t0.ru_utime.tv_usec); u += 1.e-6*(t1.ru_stime.tv_usec - t0.ru_stime.tv_usec); printf("\n%u\n", n); printf("%.2f ms for 1 million values\n", n, 1000 * u); } return 0; }
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20150317170737.GA24682>