1// -*- C++ -*- 2// $Id: drand48.src,v 1.3 2010/06/16 17:24:53 garren Exp $ 3// --------------------------------------------------------------------------- 4// 5// This code is based on a code extracted from GNU C Library 2.1.3 with 6// a purpose to provide drand48() on Windows NT. 7// 8 9#include <stdlib.h> 10#include <string.h> 11#include <limits.h> 12 13#define __LITTLE_ENDIAN 1234 14#define __BIG_ENDIAN 4321 15 16#ifdef __BIG_ENDIAN__ 17#define __BYTE_ORDER __BIG_ENDIAN /* powerpc-apple is big-endian. */ 18#else 19#define __BYTE_ORDER __LITTLE_ENDIAN /* i386 is little-endian. */ 20#endif 21#define __FLOAT_WORD_ORDER __BYTE_ORDER 22 23#define IEEE754_DOUBLE_BIAS 0x3ff /* Added to exponent. */ 24 25//#ifdef WIN32 26//#include <wtypes.h> 27//typedef ULONGLONG u_int64_t; 28//#else 29typedef unsigned long long int u_int64_t; 30//#endif 31 32union ieee754_double { 33 double d; 34 35 /* This is the IEEE 754 double-precision format. */ 36 struct { 37#if __BYTE_ORDER == __BIG_ENDIAN 38 unsigned int negative:1; 39 unsigned int exponent:11; 40 /* Together these comprise the mantissa. */ 41 unsigned int mantissa0:20; 42 unsigned int mantissa1:32; 43#endif /* Big endian. */ 44#if __BYTE_ORDER == __LITTLE_ENDIAN 45#if __FLOAT_WORD_ORDER == BIG_ENDIAN 46 unsigned int mantissa0:20; 47 unsigned int exponent:11; 48 unsigned int negative:1; 49 unsigned int mantissa1:32; 50#else 51 /* Together these comprise the mantissa. */ 52 unsigned int mantissa1:32; 53 unsigned int mantissa0:20; 54 unsigned int exponent:11; 55 unsigned int negative:1; 56#endif 57#endif /* Little endian. */ 58 } ieee; 59}; 60 61/* Data structure for communication with thread safe versions. */ 62struct drand48_data { 63 unsigned short int x[3]; /* Current state. */ 64 unsigned short int a[3]; /* Factor in congruential formula. */ 65 unsigned short int c; /* Additive const. in congruential formula. */ 66 unsigned short int old_x[3]; /* Old state. */ 67 int init; /* Flag for initializing. */ 68}; 69 70extern "C" { 71 double drand48 (); 72 int erand48_r (unsigned short int[3], drand48_data *, double *result); 73 int drand48_iterate (unsigned short int[3], drand48_data *); 74 void srand48 (long int); 75 int srand48_r (long int, drand48_data *); 76 unsigned short int * seed48 (unsigned short int[3]); 77 int seed48_r (unsigned short int[3], drand48_data *); 78} 79 80/* Global state for non-reentrant functions. */ 81static drand48_data libc_drand48_data; 82 83// --------------------------------------------------------------------------- 84double drand48 () 85{ 86 double result; 87 (void) erand48_r (libc_drand48_data.x, &libc_drand48_data, &result); 88 return result; 89} 90 91// --------------------------------------------------------------------------- 92int erand48_r (unsigned short int xsubi[3], 93 drand48_data *buffer, 94 double *result) 95{ 96 union ieee754_double temp; 97 98 /* Compute next state. */ 99 if (drand48_iterate (xsubi, buffer) < 0) return -1; 100 101 /* Construct a positive double with the 48 random bits distributed over 102 its fractional part so the resulting FP number is [0.0,1.0). */ 103 104#if USHRT_MAX == 65535 105 temp.ieee.negative = 0; 106 temp.ieee.exponent = IEEE754_DOUBLE_BIAS; 107 temp.ieee.mantissa0 = (xsubi[2] << 4) | (xsubi[1] >> 12); 108 temp.ieee.mantissa1 = ((xsubi[1] & 0xfff) << 20) | (xsubi[0] << 4); 109#elif USHRT_MAX == 2147483647 110 temp.ieee.negative = 0; 111 temp.ieee.exponent = IEEE754_DOUBLE_BIAS; 112 temp.ieee.mantissa0 = (xsubi[1] << 4) | (xsubi[0] >> 28); 113 temp.ieee.mantissa1 = ((xsubi[0] & 0xfffffff) << 4); 114#else 115# error Unsupported size of short int 116#endif 117 118 /* Please note the lower 4 bits of mantissa1 are always 0. */ 119 *result = temp.d - 1.0; 120 121 return 0; 122} 123 124// --------------------------------------------------------------------------- 125int drand48_iterate (unsigned short int xsubi[3], drand48_data *buffer) 126{ 127 u_int64_t X, a, result; 128 129 /* Initialize buffer, if not yet done. */ 130 if (!buffer->init) { 131#if (USHRT_MAX == 0xffffU) 132 buffer->a[2] = 0x5; 133 buffer->a[1] = 0xdeec; 134 buffer->a[0] = 0xe66d; 135#else 136 buffer->a[2] = 0x5deecUL; 137 buffer->a[1] = 0xe66d0000UL; 138 buffer->a[0] = 0; 139#endif 140 buffer->c = 0xb; 141 buffer->init = 1; 142 } 143 144 /* Do the real work. We choose a data type which contains at least 145 48 bits. Because we compute the modulus it does not care how 146 many bits really are computed. */ 147 148 X = (u_int64_t)xsubi[2] << 32 | (u_int64_t)xsubi[1] << 16 | xsubi[0]; 149 a = ((u_int64_t)buffer->a[2] << 32 | (u_int64_t)buffer->a[1] << 16 150 | buffer->a[0]); 151 152 result = X * a + buffer->c; 153 154 xsubi[0] = result & 0xffff; 155 xsubi[1] = (result >> 16) & 0xffff; 156 xsubi[2] = (result >> 32) & 0xffff; 157 158 return 0; 159} 160 161// --------------------------------------------------------------------------- 162void srand48 (long int seedval) 163{ 164 (void) srand48_r (seedval, &libc_drand48_data); 165} 166 167// --------------------------------------------------------------------------- 168int srand48_r (long int seedval, drand48_data *buffer) 169{ 170 /* The standards say we only have 32 bits. */ 171 if (sizeof (long int) > 4) 172 seedval &= 0xffffffffl; 173 174#if USHRT_MAX == 0xffffU 175 buffer->x[2] = seedval >> 16; 176 buffer->x[1] = seedval & 0xffffl; 177 buffer->x[0] = 0x330e; 178 179 buffer->a[2] = 0x5; 180 buffer->a[1] = 0xdeec; 181 buffer->a[0] = 0xe66d; 182#else 183 buffer->x[2] = seedval; 184 buffer->x[1] = 0x330e0000UL; 185 buffer->x[0] = 0; 186 187 buffer->a[2] = 0x5deecUL; 188 buffer->a[1] = 0xe66d0000UL; 189 buffer->a[0] = 0; 190#endif 191 buffer->c = 0xb; 192 buffer->init = 1; 193 194 return 0; 195} 196 197// --------------------------------------------------------------------------- 198unsigned short int * seed48 (unsigned short int seed16v[3]) 199{ 200 (void) seed48_r (seed16v, &libc_drand48_data); 201 return libc_drand48_data.old_x; 202} 203 204// --------------------------------------------------------------------------- 205int seed48_r (unsigned short int seed16v[3], drand48_data *buffer) 206{ 207 /* Save old value at a private place to be used as return value. */ 208 memcpy (buffer->old_x, buffer->x, sizeof (buffer->x)); 209 210 /* Install new state. */ 211#if USHRT_MAX == 0xffffU 212 buffer->x[2] = seed16v[2]; 213 buffer->x[1] = seed16v[1]; 214 buffer->x[0] = seed16v[0]; 215 216 buffer->a[2] = 0x5; 217 buffer->a[1] = 0xdeec; 218 buffer->a[0] = 0xe66d; 219#else 220 buffer->x[2] = (seed16v[2] << 16) | seed16v[1]; 221 buffer->x[1] = seed16v[0] << 16; 222 buffer->x[0] = 0; 223 224 buffer->a[2] = 0x5deecUL; 225 buffer->a[1] = 0xe66d0000UL; 226 buffer->a[0] = 0; 227#endif 228 buffer->c = 0xb; 229 buffer->init = 1; 230 231 return 0; 232} 233