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