1 /* 2 * Copyright (c) 1993 Martin Birgmeier 3 * All rights reserved. 4 * 5 * You may redistribute unmodified or modified versions of this source 6 * code provided that the above copyright notice and this and the 7 * following conditions are retained. 8 * 9 * This software is provided ``as is'', and comes with no warranties 10 * of any kind. I shall in no event be liable for anything that happens 11 * to anyone/anything when using this software. 12 */ 13 14 #ifndef _RAND48_H_ 15 #define _RAND48_H_ 16 17 #include <sys/types.h> 18 #include <math.h> 19 #include <stdlib.h> 20 21 #include "fpmath.h" 22 23 #define RAND48_SEED_0 (0x330e) 24 #define RAND48_SEED_1 (0xabcd) 25 #define RAND48_SEED_2 (0x1234) 26 #define RAND48_MULT_0 (0xe66d) 27 #define RAND48_MULT_1 (0xdeec) 28 #define RAND48_MULT_2 (0x0005) 29 #define RAND48_ADD (0x000b) 30 31 typedef uint64_t uint48; 32 33 extern uint48 _rand48_seed; 34 extern uint48 _rand48_mult; 35 extern uint48 _rand48_add; 36 37 #define TOUINT48(x, y, z) \ 38 ((uint48)(x) + (((uint48)(y)) << 16) + (((uint48)(z)) << 32)) 39 40 #define RAND48_SEED TOUINT48(RAND48_SEED_0, RAND48_SEED_1, RAND48_SEED_2) 41 #define RAND48_MULT TOUINT48(RAND48_MULT_0, RAND48_MULT_1, RAND48_MULT_2) 42 43 #define LOADRAND48(l, x) do { \ 44 (l) = TOUINT48((x)[0], (x)[1], (x)[2]); \ 45 } while (0) 46 47 #define STORERAND48(l, x) do { \ 48 (x)[0] = (unsigned short)(l); \ 49 (x)[1] = (unsigned short)((l) >> 16); \ 50 (x)[2] = (unsigned short)((l) >> 32); \ 51 } while (0) 52 53 #define _DORAND48(l) do { \ 54 (l) = (l) * _rand48_mult + _rand48_add; \ 55 } while (0) 56 57 #define DORAND48(l, x) do { \ 58 LOADRAND48(l, x); \ 59 _DORAND48(l); \ 60 STORERAND48(l, x); \ 61 } while (0) 62 63 #define ERAND48_BEGIN \ 64 union { \ 65 union IEEEd2bits ieee; \ 66 uint64_t u64; \ 67 } u; \ 68 int s 69 70 /* 71 * Optimization for speed: assume doubles are IEEE 754 and use bit fiddling 72 * rather than converting to double. Specifically, clamp the result to 48 bits 73 * and convert to a double in [0.0, 1.0) via division by 2^48. Normalize by 74 * shifting the most significant bit into the implicit one position and 75 * adjusting the exponent accordingly. The store to the exponent field 76 * overwrites the implicit one. 77 */ 78 #define ERAND48_END(x) do { \ 79 u.u64 = ((x) & 0xffffffffffffULL); \ 80 if (u.u64 == 0) \ 81 return (0.0); \ 82 u.u64 <<= 5; \ 83 for (s = 0; !(u.u64 & (1LL << 52)); s++, u.u64 <<= 1) \ 84 ; \ 85 u.ieee.bits.exp = 1022 - s; \ 86 return (u.ieee.d); \ 87 } while (0) 88 89 #endif /* _RAND48_H_ */ 90