ucw/random-fast.c

   1 /*
   2  *      UCW Library -- Fast pseudo-random generator
   3  *
   4  *      (c) 2020--2022 Pavel Charvat <pchar@ucw.cz>
   5  *
   6  *      This software may be freely distributed and used according to the terms
   7  *      of the GNU Lesser General Public License.
   8  */
   9
  10 #undef LOCAL_DEBUG
  11
  12 #include <ucw/lib.h>
  13 #include <ucw/random.h>
  14 #include <ucw/threads.h>
  15 #include <ucw/unaligned.h>
  16
  17 #include <float.h>
  18 #include <limits.h>
  19 #include <math.h>
  20 #include <stdlib.h>
  21 #include <unistd.h>
  22 #include <sys/time.h>
  23
  24 // Select exactly one of these pseudo-random generators:
  25 //#define FASTRAND_USE_PCG
  26 //#define FASTRAND_USE_SPLITMIX
  27 #define FASTRAND_USE_XOROSHIRO
  28 //#define FASTRAND_USE_LIBC_SVID
  29
  30 /**
  31  * FASTRAND_USE_PCG
  32  *
  33  * References:
  34  * -- https://en.wikipedia.org/wiki/Permuted_congruential_generator
  35  * -- https://www.pcg-random.org/
  36  * -- "PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation",
  37  *    2014, Melissa E. O’Neill
  38  *
  39  * Notes:
  40  * -- Simple and fast generator with good statistical properties.
  41  * -- We use a constant increment and 64->32bit XSH RR output function.
  42  **/
  43
  44 #define FASTRAND_PCG_MUL 6364136223846793005LLU
  45 #define FASTRAND_PCG_INC 1442695040888963407LLU
  46
  47 /**
  48  * FASTRAND_USE_SPLITMIX
  49  *
  50  * References:
  51  * -- http://prng.di.unimi.it/splitmix64.c  (public, no copyright)
  52  * -- http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
  53  * -- "Fast Splittable Pseudorandom Number Generators",
  54  *    2014, Guy L. Steele Jr., Doug Lea, Christine H. Flood
  55  *
  56  * Notes:
  57  * -- Each 64 bit number is generated exactly once in a cycle,
  58  *    so we only use half of it -> 32 bit output.
  59  **/
  60
  61 /**
  62  * FASTRAND_USE_XOROSHIRO
  63  *
  64  * References:
  65  * -- http://prng.di.unimi.it/  (examples are public, no copyright)
  66  * -- "Scrambled Linear Pseudorandom Number Generators",
  67  *    2019, David Blackman, Sebastiano Vigna
  68  *
  69  * Notes:
  70  * -- A family of very fast generators. Good statistical properties,
  71  *    some variants just have problems with lowest bits.
  72  * -- We support more xo(ro)shiro variants, see the options below.
  73  **/
  74
  75 // Engine: xoroshiro 128, xoshiro 256
  76 #define FASTRAND_XOROSHIRO_SIZE 128
  77
  78 // Scrambler: 1 (+), 2 (++), 3 (**)
  79 #define FASTRAND_XOROSHIRO_SCRAMBLER 1
  80
  81 // If defined, we only use the highest 60 bits from each 64 bit number;
  82 // useful especially for the weakest '+' scrambler. Even if not
  83 // we prefer higher bits where trivial.
  84 #undef FASTRAND_XOROSHIRO_ONLY_HIGHER
  85
  86 /**
  87  * FASTRAND_USE_LIBC_SVID
  88  * -- mrand48_r()
  89  * -- libc SVID generator
  90  **/
  91
  92 struct fastrand {
  93 #ifdef FASTRAND_USE_PCG
  94   u64 pcg_state;                        // PCG generator's current state
  95 #endif
  96
  97 #ifdef FASTRAND_USE_SPLITMIX
  98   u64 splitmix_state;
  99 #endif
 100
 101 #ifdef FASTRAND_USE_XOROSHIRO
 102   u64 xoroshiro_state[FASTRAND_XOROSHIRO_SIZE / 64];
 103 #endif
 104
 105 #ifdef FASTRAND_USE_LIBC_SVID
 106   struct drand48_data svid_data;        // SVID generator's state structure
 107 #endif
 108
 109   bool inited;                          // Is seeded?
 110   fastrand_seed_t seed_value;           // Initial seed value; useful for reproducible debugging
 111 };
 112
 113 /*** Common helpers ***/
 114
 115 // GCC should optimize these to a single instruction
 116
 117 static inline u32 fastrand_ror32(u32 x, uint r)
 118 {
 119   return (x >> r) | (x << (-r & 31));
 120 }
 121
 122 static inline u64 fastrand_rol64(u64 x, uint r)
 123 {
 124   return (x << r) | (x >> (-r & 63));
 125 }
 126
 127 /*** PCG random generator ***/
 128
 129 #ifdef FASTRAND_USE_PCG
 130 #define FASTRAND_ONE_BITS 32
 131
 132 static inline u32 fastrand_one(struct fastrand *c)
 133 {
 134   u64 s = c->pcg_state;
 135   uint r = s >> 59;
 136   c->pcg_state = s * FASTRAND_PCG_MUL + FASTRAND_PCG_INC;
 137   u32 x = ((s >> 18) ^ s) >> 27;
 138   return fastrand_ror32(x, r);
 139 }
 140
 141 static void fastrand_init(struct fastrand *c)
 142 {
 143   c->pcg_state = c->seed_value + FASTRAND_PCG_INC;
 144   (void) fastrand_one(c);
 145   DBG("RANDOM: Initialized PCG random generator, seed=0x%jx", (uintmax_t)c->seed_value);
 146 }
 147 #endif
 148
 149 /*** splitmix64 generator ***/
 150
 151 // A helper; used also for seeding xoroshiro
 152 static inline u64 fastrand_splitmix64_next(u64 *state)
 153 {
 154   u64 z = (*state += 0x9e3779b97f4a7c15LLU);
 155   z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9LLU;
 156   z = (z ^ (z >> 27)) * 0x94d049bb133111ebLLU;
 157   return z ^ (z >> 31);
 158 }
 159
 160 #ifdef FASTRAND_USE_SPLITMIX
 161 #define FASTRAND_ONE_BITS 32
 162
 163 static void fastrand_init(struct fastrand *c)
 164 {
 165   c->splitmix_state = c->seed_value;
 166   DBG("RANDOM: Initialized splitmix64 random generator, seed=0x%jx", (uintmax_t)c->seed_value);
 167 }
 168
 169 static inline u32 fastrand_one(struct fastrand *c)
 170 {
 171   return fastrand_splitmix64_next(&c->splitmix_state);
 172 }
 173 #endif
 174
 175 /*** xoroshiro variants ***/
 176
 177 #ifdef FASTRAND_USE_XOROSHIRO
 178 #ifdef FASTRAND_XOROSHIRO_ONLY_HIGHER
 179 #  define FASTRAND_ONE_BITS 60
 180 #else
 181 #  define FASTRAND_ONE_BITS 64
 182 #  define FASTRAND_PREFER_HIGHER
 183 #endif
 184
 185 static void fastrand_init(struct fastrand *c)
 186 {
 187   u64 x = c->seed_value;
 188   for (uint i = 0; i < ARRAY_SIZE(c->xoroshiro_state); i++)
 189     c->xoroshiro_state[i] = fastrand_splitmix64_next(&x);
 190   if (unlikely(!c->xoroshiro_state[0]))
 191     c->xoroshiro_state[0] = 1;
 192   DBG("RANDOM: Initialized xoroshiro random generator, seed=0x%jx, size=%u, scrambler=%u",
 193     (uintmax_t)c->seed_value, FASTRAND_XOROSHIRO_SIZE, FASTRAND_XOROSHIRO_SCRAMBLER);
 194 }
 195
 196 static inline u64 fastrand_one(struct fastrand *c)
 197 {
 198   u64 result;
 199   u64 *s = c->xoroshiro_state;
 200
 201 #if FASTRAND_XOROSHIRO_SIZE == 128
 202   u64 s0 = s[0];
 203   u64 s1 = s[1];
 204 #if FASTRAND_XOROSHIRO_SCRAMBLER == 1
 205   result = s0 + s1;
 206   uint ra = 24, rb = 16, rc = 37;
 207 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 2
 208   result = fastrand_rol64(s0 + s1, 17) + s0;
 209   uint ra = 49, rb = 21, rc = 28;
 210 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 3
 211   result = fastrand_rol64(s0 * 5, 7) * 9;
 212   uint ra = 24, rb = 16, rc = 37;
 213 #else
 214 #  error "Unsupported xoroshiro parameters"
 215 #endif
 216   s[0] = fastrand_rol64(s0, ra) ^ s1 ^ (s1 << rb);
 217   s[1] = fastrand_rol64(s1, rc);
 218
 219 #elif FASTRAND_XOROSHIRO_SIZE == 256
 220 #if FASTRAND_XOROSHIRO_SCRAMBLER == 1
 221   result = s[0] + s[3];
 222 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 2
 223   result = fastrand_rol64(s[0] + s[3], 23) + s[0];
 224 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 3
 225   result = fastrand_rol64(s[1] * 5, 7) * 9;
 226 #else
 227 #  error "Unsupported xoroshiro parameters"
 228 #endif
 229   u64 t = s[1] << 17;
 230   s[2] ^= s[0];
 231   s[3] ^= s[1];
 232   s[1] ^= s[2];
 233   s[0] ^= s[3];
 234   s[2] ^= t;
 235   s[3] = fastrand_rol64(s[3], 45);
 236 #else
 237 #  error "Unsupported xoroshiro parameters"
 238 #endif
 239
 240 #ifdef FASTRAND_XOROSHIRO_ONLY_HIGHER
 241   result >>= 4;
 242 #endif
 243   return result;
 244 }
 245
 246 #endif
 247
 248 /*** Libc SVID generator ***/
 249
 250 #ifdef FASTRAND_USE_LIBC_SVID
 251 #define FASTRAND_ONE_BITS 32
 252
 253 static void fastrand_init(struct fastrand *c)
 254 {
 255   DBG("RANDOM: Initialized libc SVID random generator, seed=0x%jx", (uintmax_t)c->seed_value);
 256   if (srand48_r(c->seed_value, &c->svid_data) < 0)
 257     die("RANDOM: Failed to call srand48_r(): %m");
 258 }
 259
 260 #ifdef LOCAL_DEBUG
 261 static NONRET void fastrand_mrand48_r_failed(void)
 262 {
 263   die("RANDOM: Failed to call mrand48_r(): %m");
 264 }
 265 #endif
 266
 267 static inline u32 fastrand_one(struct fastrand *c)
 268 {
 269   long int r;
 270 #ifdef LOCAL_DEBUG
 271   if (unlikely(mrand48_r(&c->svid_data, &r) < 0))
 272     fastrand_mrand48_r_failed();
 273 #else
 274   (void) mrand48_r(&c->svid_data, &r);
 275 #endif
 276   return r;
 277 }
 278
 279 #endif
 280
 281 /*** Interface for creating/deleting/setting contexts ***/
 282
 283 struct fastrand *fastrand_new(void)
 284 {
 285   struct fastrand *c = xmalloc_zero(sizeof(*c));
 286   return c;
 287 }
 288
 289 void fastrand_delete(struct fastrand *c)
 290 {
 291   ASSERT(c);
 292   xfree(c);
 293 }
 294
 295 void fastrand_set_seed(struct fastrand *c, fastrand_seed_t seed)
 296 {
 297   c->seed_value = seed;
 298   c->inited = true;
 299   fastrand_init(c);
 300 }
 301
 302 fastrand_seed_t fastrand_gen_seed_value(void)
 303 {
 304   // Generate a seed value as a mixture of:
 305   // -- current time
 306   // -- process id to deal with frequent fork() or startups
 307   // -- counter to deal with multiple contexts, for example in threaded application
 308   u64 val;
 309   static u64 static_counter;
 310   ucwlib_lock();
 311   val = (static_counter += 0x1927a7639274162dLLU);
 312   ucwlib_unlock();
 313   val += getpid();
 314   struct timeval tv;
 315   if (gettimeofday(&tv, NULL) < 0)
 316     die("RANDOM: Failed to call gettimeofday(): %m");
 317   val += (u64)(tv.tv_sec * 1000000LL + tv.tv_usec) * 0x5a03173d83f78347LLU;
 318   return val;
 319 }
 320
 321 fastrand_seed_t fastrand_gen_seed(struct fastrand *c)
 322 {
 323   fastrand_seed_t seed = fastrand_gen_seed_value();
 324   fastrand_set_seed(c, seed);
 325   return seed;
 326 }
 327
 328 bool fastrand_has_seed(struct fastrand *c)
 329 {
 330   return c->inited;
 331 }
 332
 333 void fastrand_reset(struct fastrand *c)
 334 {
 335   c->inited = false;
 336 }
 337
 338 #if 1
 339 // XXX: Adds some overhead, but probably worth it
 340 #define FASTRAND_CHECK_SEED(c) do { if (unlikely(!(c)->inited)) fastrand_check_seed_failed(); } while (0)
 341 static NONRET void fastrand_check_seed_failed(void)
 342 {
 343   ASSERT(0);
 344 }
 345 #else
 346 #define FASTRAND_CHECK_SEED(c) do { } while (0)
 347 #endif
 348
 349 /*** Reading interface ***/
 350
 351 // We assume 32-64 bit output of fastrand_one().
 352 // For 32 bits fastrand_one() can return u32 type, otherwise it must be u64.
 353 COMPILE_ASSERT(fastrand_one_range, FASTRAND_ONE_BITS >= 32 && FASTRAND_ONE_BITS <= 64);
 354
 355 #if FASTRAND_ONE_BITS == 32
 356 #define FASTRAND_ONE_MAX UINT32_MAX
 357 #else
 358 #define FASTRAND_ONE_MAX (UINT64_MAX >> (64 - FASTRAND_ONE_BITS))
 359 #endif
 360
 361 static inline u32 fastrand_u32_helper(struct fastrand *c)
 362 {
 363 #if FASTRAND_ONE_BITS > 32 && defined(FASTRAND_PREFER_HIGHER)
 364   return fastrand_one(c) >> (FASTRAND_ONE_BITS - 32);
 365 #else
 366   return fastrand_one(c);
 367 #endif
 368 }
 369
 370 u32 fastrand_u32(struct fastrand *c)
 371 {
 372   FASTRAND_CHECK_SEED(c);
 373   return fastrand_u32_helper(c);
 374 }
 375
 376 u64 fastrand_u64(struct fastrand *c)
 377 {
 378   FASTRAND_CHECK_SEED(c);
 379 #if FASTRAND_ONE_BITS < 64
 380   return ((u64)fastrand_one(c) << FASTRAND_ONE_BITS) | fastrand_one(c);
 381 #else
 382   return fastrand_one(c);
 383 #endif
 384 }
 385
 386 u64 fastrand_bits_u64(struct fastrand *c, uint bits)
 387 {
 388 #ifdef LOCAL_DEBUG
 389   ASSERT(bits <= 64);
 390 #endif
 391   // Notice that we never rotate an integer by all its bits (forbidden by C standard).
 392   FASTRAND_CHECK_SEED(c);
 393 #if FASTRAND_ONE_BITS == 64
 394   if (!bits)
 395     return 0;
 396   return fastrand_one(c) >> (FASTRAND_ONE_BITS - bits);
 397 #else
 398   u64 r = fastrand_one(c);
 399   if (bits <= FASTRAND_ONE_BITS)
 400     return r >> (FASTRAND_ONE_BITS - bits);
 401   return (r << (bits - FASTRAND_ONE_BITS)) | (fastrand_one(c) >> (2 * FASTRAND_ONE_BITS - bits));
 402 #endif
 403 }
 404
 405 /*
 406  * fastrand_max() and variants
 407  *
 408  * Algorithm 1:
 409  *
 410  * We can return r % bound from uniformly generated i-bit number r,
 411  * but only if:
 412  *   r - r % bound + bound <= 2^i
 413  *   r - r % bound <= 2^i-1 - bound + 1    -- no risk of integer overflow
 414  *   r - r % bound <= (2^i-1 + 1) - bound  -- can temporarily overflow
 415  *
 416  * Then r fits to a complete interval in the i-bit range
 417  * and results should be uniformly distributed.
 418  *
 419  * If r is too high, we simply generate a different r and repeat
 420  * the whole process. The probability is < 50% for any bound <= 2^i
 421  * (with worst case bound == 2^i / 2 + 1), so the average number
 422  * of iterations is < 2.
 423  *
 424  *
 425  * Algorithm 2 (without any division in many cases, but needs big multiplication):
 426  * ("Fast Random Integer Generation in an Interval", 2019, Daniel Lemier)
 427  *
 428  * Let bound < 2^i. Then for j in [0, bound) each interval [j * 2^i + 2^i % bound, (j + 1) * 2^i)
 429  * contains the same number of multiples [0, 1 * bound, 2 * bound, ..., 2^i * bound),
 430  * exactly 2^i / bound (rounded down) of them.
 431  *
 432  * Therefore if we have a random i-bit number r and r * bound fits to one of these intervals,
 433  * we can immediately return (r * bound) / 2^i (index of that interval) and results should be
 434  * correctly distributed. If we are lucky and r fits to even smaller interval
 435  * [j * 2^i + bound, (j + 1) * 2^i], we can check that very quickly without any division.
 436  *
 437  * Probability of multiple iterations is again < 50%.
 438  * Probability of division-less check is (2^i - bound) / 2^i.
 439  *
 440  *
 441  * Algorithm 3 (always without divisions or multiplications, but needs more iterations):
 442  *
 443  * Let 2^(i-1) <= bound < 2^i. Then i-bit random number r has at least 50% probability
 444  * to be lower than bound, when we can simply return it.
 445  * Can be useful for very fast fastrand_one() like xoroshiro128+.
 446  */
 447
 448 #define FASTRAND_MAX_U32_ALGORITHM 2
 449 #define FASTRAND_MAX_U64_ALGORITHM 3
 450
 451 static inline u32 fastrand_max_u32_helper(struct fastrand *c, u32 bound)
 452 {
 453   // Variant for bound <= UINT32_MAX.
 454   // This is the most common case -> we optimize it and inline the code to both fastrand_max() and fastrand_max_u64().
 455
 456 #if FASTRAND_MAX_U32_ALGORITHM == 1
 457   while (1)
 458     {
 459       u32 r = fastrand_u32_helper(c);
 460       u32 x = r % bound;
 461       if (r - x <= (u32)-bound)
 462         return x;
 463 #if 1
 464       // Possible optimization for unlucky cases
 465       u32 y = r - x;
 466       do
 467         r = fastrand_u32_helper(c);
 468       while (r >= y);
 469       return r % bound;
 470 #endif
 471     }
 472
 473 #elif FASTRAND_MAX_U32_ALGORITHM == 2
 474   u32 r = fastrand_u32_helper(c);
 475   u64 m = (u64)r * bound;
 476   u32 l = (u32)m;
 477   if (l < bound)
 478     {
 479       u32 t = (u32)-bound;
 480 #if 1
 481       // This condition can decrease the probability of division to <= 50% even in worst case and also optimize bad cases.
 482       if (t < bound)
 483         {
 484           // Implies l < t % bound, so we don't need to check it before switching to completely different loop.
 485           do
 486             r = fastrand_u32_helper(c);
 487           while (r >= bound);
 488           return r;
 489         }
 490 #endif
 491       t = t % bound;
 492       while (l < t)
 493         {
 494           r = fastrand_u32_helper(c);
 495           m = (u64)r * bound;
 496           l = (u32)m;
 497         }
 498     }
 499   return m >> 32;
 500
 501 #elif FASTRAND_MAX_U32_ALGORITHM == 3
 502   // XXX: Probably could be optimized with __builtin_clz()
 503   u32 r, y = bound - 1;
 504   y |= y >> 1;
 505   y |= y >> 2;
 506   y |= y >> 4;
 507   y |= y >> 8;
 508   y |= y >> 16;
 509   do
 510     r = fastrand_u32_helper(c) & y;
 511   while (r >= bound);
 512   return r;
 513
 514 #else
 515 #  error Invalid FASTRAND_MAX_U32_ALGORITHM
 516 #endif
 517 }
 518
 519 static u64 fastrand_max_u64_helper(struct fastrand *c, u64 bound)
 520 {
 521   // Less common case -> compromise between speed and size
 522 #if FASTRAND_MAX_U64_ALGORITHM == 1
 523   while (1)
 524     {
 525       u64 r = fastrand_one(c);
 526       u64 m = FASTRAND_ONE_MAX;
 527 #if FASTRAND_ONE_BITS < 64
 528       if (bound > m + 1)
 529         {
 530           r = (r << FASTRAND_ONE_BITS) | fastrand_one(c);
 531           m = (m << FASTRAND_ONE_BITS) | FASTRAND_ONE_MAX;
 532         }
 533 #endif
 534       u64 x = r % bound;
 535       if (r - x <= m + 1 - bound)
 536         return x;
 537     }
 538
 539 #elif FASTRAND_MAX_U64_ALGORITHM == 3
 540   u64 r, y = bound - 1;
 541   y |= y >> 1;
 542   y |= y >> 2;
 543   y |= y >> 4;
 544   y |= y >> 8;
 545   y |= y >> 16;
 546   y |= y >> 32;
 547   do
 548     {
 549       r = fastrand_one(c);
 550 #if FASTRAND_ONE_BITS < 64
 551       if (y > FASTRAND_ONE_MAX)
 552         r = (r << FASTRAND_ONE_BITS) | fastrand_one(c);
 553 #endif
 554       r &= y;
 555     }
 556   while (r >= bound);
 557   return r;
 558 #else
 559 #  error Invalid FASTRAND_MAX_U64_ALGORITHM
 560 #endif
 561 }
 562
 563 uint fastrand_max(struct fastrand *c, uint bound)
 564 {
 565   FASTRAND_CHECK_SEED(c);
 566 #if UINT32_MAX < UINT_MAX
 567   if (bound <= UINT32_MAX)
 568     return fastrand_max_u32_helper(c, bound);
 569   else
 570     return fastrand_max_u64_helper(c, bound);
 571 #else
 572   return fastrand_max_u32_helper(c, bound);
 573 #endif
 574 }
 575
 576 u64 fastrand_max_u64(struct fastrand *c, u64 bound)
 577 {
 578   FASTRAND_CHECK_SEED(c);
 579   if (bound <= UINT32_MAX)
 580     return fastrand_max_u32_helper(c, bound);
 581   else
 582     return fastrand_max_u64_helper(c, bound);
 583 }
 584
 585 void fastrand_mem(struct fastrand *c, void *ptr, size_t size)
 586 {
 587   // XXX: Could be optimized to aligned assignments instead of put_u32/64, but it
 588   // would be quite complex to make it deterministic (for same seed and any @ptr)
 589   // and large sizes are not used often anyway.
 590   FASTRAND_CHECK_SEED(c);
 591   byte *p = ptr;
 592 #if FASTRAND_ONE_BITS < 64
 593   for (size_t n = size >> 2; n--; p += 4)
 594     put_u32(p, fastrand_u32_helper(c));
 595   if (size &= 3)
 596     {
 597       u32 x = fastrand_u32_helper(c);
 598 #else
 599   for (size_t n = size >> 3; n--; p += 8)
 600     put_u64(p, fastrand_one(c));
 601   if (size &= 7)
 602     {
 603       u64 x = fastrand_one(c);
 604 #endif
 605       while (size--)
 606         {
 607           *p++ = x;
 608           x >>= 8;
 609         }
 610     }
 611 }
 612
 613 double fastrand_double(struct fastrand *c)
 614 {
 615   // Generate 52 random bits
 616   FASTRAND_CHECK_SEED(c);
 617   u64 r = fastrand_one(c);
 618 #if FASTRAND_ONE_BITS >= 52
 619   r >>= FASTRAND_ONE_BITS - 52;
 620 #else
 621   // Do we want to generate doubles with full 52 precision?
 622 #if 0
 623   r = (r << (52 - FASTRAND_ONE_BITS)) | (fastrand_one(c) >> (2 * FASTRAND_ONE_BITS - 52));
 624 #else
 625   r <<= 52 - FASTRAND_ONE_BITS;
 626 #endif
 627 #endif
 628
 629   // With standard IEC 60559 / IEEE 754 format, we can directly store @r as a double
 630   // from range [1.0, 2.0) and subtract 1.0.
 631   // But the alternative with multiplication is often faster anyway.
 632 #if 0 && defined(__STDC_IEC_559__)
 633   union {
 634     u64 u;
 635     double d;
 636   } tmp;
 637   ASSERT(sizeof(double) == 8 && sizeof(tmp) == 8);
 638   tmp.u = (0x3ffLLU << 52) | r;
 639   return tmp.d - 1.0;
 640 #elif 1
 641   static const double mul = 1.0 / (1LLU << 52);
 642   return r * mul;
 643 #else
 644   return ldexp(r, -52);
 645 #endif
 646 }