aboutsummaryrefslogtreecommitdiff
path: root/lmathlib.c
diff options
context:
space:
mode:
authorRoberto Ierusalimschy <roberto@inf.puc-rio.br>2018-04-04 13:12:53 -0300
committerRoberto Ierusalimschy <roberto@inf.puc-rio.br>2018-04-04 13:12:53 -0300
commitb44787652bd3a51dfc2ee92120abe05b28159d32 (patch)
tree277a6db09ab19f1bb3fab87d1358fcc285c21dc3 /lmathlib.c
parent03c6a05ec836c3a90a6b8d730120afdad39c092b (diff)
downloadlua-b44787652bd3a51dfc2ee92120abe05b28159d32.tar.gz
lua-b44787652bd3a51dfc2ee92120abe05b28159d32.tar.bz2
lua-b44787652bd3a51dfc2ee92120abe05b28159d32.zip
using 'xoroshiro128+' for PRNG
(plus a rotate at the final result to have better lower bits)
Diffstat (limited to 'lmathlib.c')
-rw-r--r--lmathlib.c175
1 files changed, 88 insertions, 87 deletions
diff --git a/lmathlib.c b/lmathlib.c
index 9432d405..3e5ceea4 100644
--- a/lmathlib.c
+++ b/lmathlib.c
@@ -1,5 +1,5 @@
1/* 1/*
2** $Id: lmathlib.c,v 1.127 2018/03/22 19:54:49 roberto Exp roberto $ 2** $Id: lmathlib.c,v 1.128 2018/03/26 19:48:46 roberto Exp $
3** Standard mathematical library 3** Standard mathematical library
4** See Copyright Notice in lua.h 4** See Copyright Notice in lua.h
5*/ 5*/
@@ -247,7 +247,7 @@ static int math_type (lua_State *L) {
247 247
248/* 248/*
249** {================================================================== 249** {==================================================================
250** Pseudo-Random Number Generator based on 'xorshift128+'. 250** Pseudo-Random Number Generator based on 'xoroshiro128+'.
251** =================================================================== 251** ===================================================================
252*/ 252*/
253 253
@@ -270,34 +270,45 @@ static int math_type (lua_State *L) {
270/* a 64-bit value */ 270/* a 64-bit value */
271typedef unsigned long long I; 271typedef unsigned long long I;
272 272
273static I xorshift128plus (I *state) { 273
274 I x = state[0]; 274/* rotate left 'x' by 'n' bits */
275 I y = state[1]; 275static I rotl (I x, int n) {
276 state[0] = y; 276 return (x << n) | (x >> (64 - n));
277 x ^= x << 23; 277}
278 state[1] = (x ^ (x >> 18)) ^ (y ^ (y >> 5)); 278
279 return state[1] + y; 279static I nextrand (I *state) {
280 I s0 = state[0];
281 I s1 = state[1];
282 I res = s0 + s1;
283 res = rotl(res, 41); /* extra step to change place of lower bits */
284 s1 = s1 ^ s0;
285 state[0] = rotl(s0, 55) ^ (s1 ^ (s1 << 14));
286 state[1] = rotl(s1, 36);
287 return res;
280} 288}
281 289
290
282/* must take care to not shift stuff by more than 63 slots */ 291/* must take care to not shift stuff by more than 63 slots */
283 292
284#define shiftI (64 - FIGS) /* leave FIGS bits */
285#define shiftF (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIG) */
286 293
287/* 294/*
288** Convert bits from a random integer into a float in the 295** Convert bits from a random integer into a float in the
289** interval [0,1). 296** interval [0,1).
290*/ 297*/
298#define maskFIG (~(~1LLU << (FIGS - 1))) /* use FIGS bits */
299#define shiftFIG (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIGS) */
300
291static lua_Number I2d (I x) { 301static lua_Number I2d (I x) {
292 return (lua_Number)(x >> shiftI) * shiftF; 302 return (lua_Number)(x & maskFIG) * shiftFIG;
293} 303}
294 304
295/* convert an 'I' to a lua_Unsigned (using higher bits) */ 305/* convert an 'I' to a lua_Unsigned */
296#define I2UInt(x) ((lua_Unsigned)((x) >> (64 - LUA_UNSIGNEDBITS))) 306#define I2UInt(x) ((lua_Unsigned)(x))
297 307
298/* convert a lua_Integer to an 'I' */ 308/* convert a lua_Unsigned to an 'I' */
299#define Int2I(x) ((I)(x)) 309#define Int2I(x) ((I)(x))
300 310
311
301#else /* no long long }{ */ 312#else /* no long long }{ */
302 313
303/* 314/*
@@ -330,14 +341,10 @@ static I packI (lu_int32 h, lu_int32 l) {
330 341
331/* i ^ (i << n) */ 342/* i ^ (i << n) */
332static I Ixorshl (I i, int n) { 343static I Ixorshl (I i, int n) {
344 lua_assert(n > 0 && n < 32);
333 return packI(i.h ^ ((i.h << n) | (i.l >> (32 - n))), i.l ^ (i.l << n)); 345 return packI(i.h ^ ((i.h << n) | (i.l >> (32 - n))), i.l ^ (i.l << n));
334} 346}
335 347
336/* i ^ (i >> n) */
337static I Ixorshr (I i, int n) {
338 return packI(i.h ^ (i.h >> n), i.l ^ ((i.l >> n) | (i.h << (32 - n))));
339}
340
341static I Ixor (I i1, I i2) { 348static I Ixor (I i1, I i2) {
342 return packI(i1.h ^ i2.h, i1.l ^ i2.l); 349 return packI(i1.h ^ i2.h, i1.l ^ i2.l);
343} 350}
@@ -349,18 +356,29 @@ static I Iadd (I i1, I i2) {
349 return result; 356 return result;
350} 357}
351 358
359/*
360** Rotate left. As all offsets here are larger than 32, do a rotate right
361** of 64 - offset.
362*/
363static I Irotli (I i, int n) {
364 n = 64 - n;
365 lua_assert(n > 0 && n < 32);
366 return packI((i.h >> n) | (i.l << (32 - n)),
367 (i.h << (32 - n)) | (i.l >> n));
368}
352 369
353/* 370/*
354** implementation of 'xorshift128+' algorithm on 'I' values 371** implementation of 'xoroshiro128+' algorithm on 'I' values
355*/ 372*/
356static I xorshift128plus (I *state) { 373static I nextrand (I *state) {
357 I x = state[0]; 374 I s0 = state[0];
358 I y = state[1]; 375 I s1 = state[1];
359 state[0] = y; 376 I res = Iadd(s0, s1);
360 x = Ixorshl(x, 23); /* x ^= x << 23; */ 377 res = Irotli(res, 41);
361 /* state[1] = (x ^ (x >> 18)) ^ (y ^ (y >> 5)); */ 378 s1 = Ixor(s1, s0);
362 state[1] = Ixor(Ixorshr(x, 18), Ixorshr(y, 5)); 379 state[0] = Ixor(Irotli(s0, 55), Ixorshl(s1, 14));
363 return Iadd(state[1], y); /* return state[1] + y; */ 380 state[1] = Irotli(s1, 36);
381 return res;
364} 382}
365 383
366 384
@@ -373,45 +391,39 @@ static I xorshift128plus (I *state) {
373 391
374#if FIGS <= 32 392#if FIGS <= 32
375 393
376#define maskLOW 0 /* do not need bits from lower half */ 394#define maskHI 0 /* do not need bits from higher half */
377#define maskHI (~(~(lu_int32)0 >> (FIGS - 1) >> 1)) /* use FIGS bits */ 395#define maskLOW (~(~UONE << (FIGS - 1))) /* use FIGS bits */
378#define shiftHI 1 /* no shift */ 396#define shiftFIG (l_mathop(0.5) / (UONE << (FIGS - 1))) /* 2^(-FIGS) */
379#define shiftF (1 / l_mathop(4294967296.0)) /* 2^(-32) */
380 397
381#else /* 32 < FIGS <= 64 */ 398#else /* 32 < FIGS <= 64 */
382 399
383/* must take care to not shift stuff by more than 31 slots */ 400/* must take care to not shift stuff by more than 31 slots */
384 401
385/* use FIGS - 32 bits from lower half */ 402/* use FIGS - 32 bits from higher half */
386#define maskLOW (~(~(lu_int32)0 >> (FIGS - 33) >> 1)) 403#define maskHI (~(~UONE << (FIGS - 33)))
387
388/* use all bits from higher half */
389#define maskHI (~(lu_int32)0)
390 404
391#define shiftHI l_mathop(4294967296.0) /* 2^32 */ 405/* use all bits from lower half */
406#define maskLOW (~(lu_int32)0)
392 407
393/* 2^(-64) */ 408/* 2^(-FIGS) == (1 / 2^33) / 2^(FIGS-33) */
394#define shiftF ((lua_Number)(1 / (shiftHI * shiftHI))) 409#define shiftFIG ((lua_Number)(1.0 / 8589934592.0) / (UONE << (FIGS - 33)))
395 410
396#endif 411#endif
397 412
413#define twoto32 l_mathop(4294967296.0) /* 2^32 */
414
398static lua_Number I2d (I x) { 415static lua_Number I2d (I x) {
399 lua_Number h = (lua_Number)(x.h & maskHI); 416 lua_Number h = (lua_Number)(x.h & maskHI);
400 lua_Number l = (lua_Number)(x.l & maskLOW); 417 lua_Number l = (lua_Number)(x.l & maskLOW);
401 return (h * shiftHI + l) * shiftF; 418 return (h * twoto32 + l) * shiftFIG;
402} 419}
403 420
421
404static lua_Unsigned I2UInt (I x) { 422static lua_Unsigned I2UInt (I x) {
405#if (LUA_MAXINTEGER >> 30) <= 1 423 return ((lua_Unsigned)x.h << 31 << 1) | (lua_Unsigned)x.l;
406/* at most 32 bits; use only high bits */
407 return ((lua_Unsigned)x.h);
408#else
409/* at least 33 bits */
410 return ((lua_Unsigned)x.h << (LUA_UNSIGNEDBITS - 32)) |
411 (lua_Unsigned)x.l >> (64 - LUA_UNSIGNEDBITS);
412#endif
413} 424}
414 425
426
415static I Int2I (lua_Unsigned n) { 427static I Int2I (lua_Unsigned n) {
416 return packI((lu_int32)(n >> 31 >> 1), (lu_int32)n & ~(lu_int32)0); 428 return packI((lu_int32)(n >> 31 >> 1), (lu_int32)n & ~(lu_int32)0);
417} 429}
@@ -428,46 +440,35 @@ typedef struct {
428 440
429 441
430/* 442/*
431** Return the higher bit set in 'x' (first bit is 1).
432*/
433static int higherbit (lua_Unsigned x) {
434 /* table of higher bits from 0 to 255 */
435 static const unsigned char hb[256] = {
436 0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
437 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
438 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
439 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
440 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
441 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
442 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
443 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8
444 };
445 int l = 0;
446 while (x >= 256) { l += 8; x >>= 8; }
447 return l + hb[x];
448}
449
450
451/*
452** Project the random integer 'ran' into the interval [0, n]. 443** Project the random integer 'ran' into the interval [0, n].
453** To get a uniform projection into [0,n], we first compute 'shf', the 444** Because 'ran' has 2^B possible values, the projection can only be
454** largest number that we can right-shift 'ran' and still get numbers 445** uniform when the size of the interval is a power of 2 (exact
455** as larger as 'n'. We then shift 'ran'; if the result is inside [0, n], 446** division). To get a uniform projection into [0, n], we first compute
456** we are done. Otherwise, we try with another 'ran' until we have a 447** 'lim', the smallest Mersenne number not smaller than 'n'. We then
457** result inside the interval. (We use right shifts to avoid the lowest 448** project 'ran' into the interval [0, lim]. If the result is inside
458** bits of 'ran', which has poorer distributions.) 449** [0, n], we are done. Otherwise, we try with another 'ran' until we
450** have a result inside the interval.
459*/ 451*/
460static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n, 452static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n,
461 RanState *state) { 453 RanState *state) {
462 if (n == 0) return 0; /* special case for the unit set */ 454 lua_Unsigned lim = n;
463 else { 455 if ((lim & (lim + 1)) > 0) { /* 'lim + 1' is not a power of 2? */
464 int shf = LUA_UNSIGNEDBITS - higherbit(n); 456 /* compute the smallest (2^b - 1) not smaller than 'n' */
465 lua_assert(~(lua_Unsigned)0 >> shf >= n && /* not larger */ 457 lim |= (lim >> 1);
466 ~(lua_Unsigned)0 >> shf >> 1 < n); /* largest */ 458 lim |= (lim >> 2);
467 while ((ran >>= shf) > n) 459 lim |= (lim >> 4);
468 ran = I2UInt(xorshift128plus(state->s)); 460 lim |= (lim >> 8);
469 return ran; 461 lim |= (lim >> 16);
462#if (LUA_MAXINTEGER >> 30 >> 2) > 0
463 lim |= (lim >> 32); /* integer type has more than 32 bits */
464#endif
470 } 465 }
466 lua_assert((lim & (lim + 1)) == 0 /* 'lim + 1' is a power of 2 */
467 && lim >= n /* not smaller than 'n' */
468 && (lim == 0 || (lim >> 1) < n)); /* it is the smallest one */
469 while ((ran &= lim) > n)
470 ran = I2UInt(nextrand(state->s));
471 return ran;
471} 472}
472 473
473 474
@@ -475,7 +476,7 @@ static int math_random (lua_State *L) {
475 lua_Integer low, up; 476 lua_Integer low, up;
476 lua_Unsigned p; 477 lua_Unsigned p;
477 RanState *state = (RanState *)lua_touserdata(L, lua_upvalueindex(1)); 478 RanState *state = (RanState *)lua_touserdata(L, lua_upvalueindex(1));
478 I rv = xorshift128plus(state->s); /* next pseudo-random value */ 479 I rv = nextrand(state->s); /* next pseudo-random value */
479 switch (lua_gettop(L)) { /* check number of arguments */ 480 switch (lua_gettop(L)) { /* check number of arguments */
480 case 0: { /* no arguments */ 481 case 0: { /* no arguments */
481 lua_pushnumber(L, I2d(rv)); /* float between 0 and 1 */ 482 lua_pushnumber(L, I2d(rv)); /* float between 0 and 1 */
@@ -511,7 +512,7 @@ static void setseed (I *state, lua_Unsigned n) {
511 state[0] = Int2I(n); 512 state[0] = Int2I(n);
512 state[1] = Int2I(0xff); /* avoid a zero state */ 513 state[1] = Int2I(0xff); /* avoid a zero state */
513 for (i = 0; i < 16; i++) 514 for (i = 0; i < 16; i++)
514 xorshift128plus(state); /* discard initial values to "spread" seed */ 515 nextrand(state); /* discard initial values to "spread" seed */
515} 516}
516 517
517 518