From owner-freebsd-numerics@FreeBSD.ORG Tue Mar 17 17:07:39 2015 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id D33D3FA4; Tue, 17 Mar 2015 17:07:39 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id B47C3AE4; Tue, 17 Mar 2015 17:07:39 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id t2HH7bwk024868 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Tue, 17 Mar 2015 10:07:37 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id t2HH7bCq024867; Tue, 17 Mar 2015 10:07:37 -0700 (PDT) (envelope-from sgk) Date: Tue, 17 Mar 2015 10:07:37 -0700 From: Steve Kargl To: Pedro Giffuni Subject: Re: Random number generators Message-ID: <20150317170737.GA24682@troutmask.apl.washington.edu> References: <7CBD7758-9472-4A2E-8065-EC6E68EE8DAB@FreeBSD.org> <20150317060310.GA21975@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: User-Agent: Mutt/1.5.23 (2014-03-12) Cc: freebsd-numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 17 Mar 2015 17:07:39 -0000 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 #include #include #include #include /* * 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; }