diff options
author | Roberto Ierusalimschy <roberto@inf.puc-rio.br> | 2018-04-04 13:12:53 -0300 |
---|---|---|
committer | Roberto Ierusalimschy <roberto@inf.puc-rio.br> | 2018-04-04 13:12:53 -0300 |
commit | b44787652bd3a51dfc2ee92120abe05b28159d32 (patch) | |
tree | 277a6db09ab19f1bb3fab87d1358fcc285c21dc3 /lmathlib.c | |
parent | 03c6a05ec836c3a90a6b8d730120afdad39c092b (diff) | |
download | lua-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.c | 175 |
1 files changed, 88 insertions, 87 deletions
@@ -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 */ |
271 | typedef unsigned long long I; | 271 | typedef unsigned long long I; |
272 | 272 | ||
273 | static I xorshift128plus (I *state) { | 273 | |
274 | I x = state[0]; | 274 | /* rotate left 'x' by 'n' bits */ |
275 | I y = state[1]; | 275 | static 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; | 279 | static 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 | |||
291 | static lua_Number I2d (I x) { | 301 | static 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) */ |
332 | static I Ixorshl (I i, int n) { | 343 | static 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) */ | ||
337 | static 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 | |||
341 | static I Ixor (I i1, I i2) { | 348 | static 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 | */ | ||
363 | static 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 | */ |
356 | static I xorshift128plus (I *state) { | 373 | static 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 | |||
398 | static lua_Number I2d (I x) { | 415 | static 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 | |||
404 | static lua_Unsigned I2UInt (I x) { | 422 | static 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 | |||
415 | static I Int2I (lua_Unsigned n) { | 427 | static 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 | */ | ||
433 | static 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 | */ |
460 | static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n, | 452 | static 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 | ||