--- rand48.h.orig 2009-11-07 14:51:38.000000000 -0800 +++ rand48.h 2009-11-07 14:51:39.000000000 -0800 @@ -19,8 +19,6 @@ #include <math.h> #include <stdlib.h> -void _dorand48(unsigned short[3]); - #define RAND48_SEED_0 (0x330e) #define RAND48_SEED_1 (0xabcd) #define RAND48_SEED_2 (0x1234) @@ -29,4 +27,55 @@ void _dorand48(unsigned short[3]); #define RAND48_MULT_2 (0x0005) #define RAND48_ADD (0x000b) +typedef unsigned long long uint48; + +extern uint48 _rand48_seed; +extern uint48 _rand48_mult; +extern uint48 _rand48_add; + +#define TOUINT48(x,y,z) \ + ((uint48)(x) + (((uint48)(y)) << 16) + (((uint48)(z)) << 32)) + +#define RAND48_SEED TOUINT48(RAND48_SEED_0, RAND48_SEED_1, RAND48_SEED_2) +#define RAND48_MULT TOUINT48(RAND48_MULT_0, RAND48_MULT_1, RAND48_MULT_2) + +#define LOADRAND48(l,x) \ + (l) = TOUINT48((x)[0], (x)[1], (x)[2]) + +#define STORERAND48(l,x) \ + (x)[0] = (unsigned short)(l); \ + (x)[1] = (unsigned short)((l) >> 16); \ + (x)[2] = (unsigned short)((l) >> 32) + +#define _DORAND48(l) \ + (l) = (l) * _rand48_mult + _rand48_add + +#define DORAND48(l,x) \ + LOADRAND48(l, x); \ + _DORAND48(l); \ + STORERAND48(l, x) + +#include "fpmath.h" + +/* + * Optimization for speed: avoid int-to-double conversion. Assume doubles + * are IEEE-754 and insert the bits directly. To normalize, the (1 << 52) + * is the hidden bit, which the first set bit is shifted to. + */ +#define ERAND48_BEGIN \ + union { \ + union IEEEd2bits ieee; \ + unsigned long long l; \ + } u; \ + int s + +#define ERAND48_END(x) \ + u.l = ((x) & 0xffffffffffffULL); \ + if (u.l == 0) \ + return 0.0; \ + u.l <<= 5; \ + for(s = 0; !(u.l & (1LL << 52)); s++, u.l <<= 1) {} \ + u.ieee.bits.exp = 1022 - s; \ + return u.ieee.d + #endif /* _RAND48_H_ */