1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2002-2012, Industrial Light & Magic, a division of Lucas
4 // Digital Ltd. LLC
5 //
6 // All rights reserved.
7 //
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
10 // met:
11 // * Redistributions of source code must retain the above copyright
12 // notice, this list of conditions and the following disclaimer.
13 // * Redistributions in binary form must reproduce the above
14 // copyright notice, this list of conditions and the following disclaimer
15 // in the documentation and/or other materials provided with the
16 // distribution.
17 // * Neither the name of Industrial Light & Magic nor the names of
18 // its contributors may be used to endorse or promote products derived
19 // from this software without specific prior written permission.
20 //
21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 //
33 ///////////////////////////////////////////////////////////////////////////
34
35 //-----------------------------------------------------------------------------
36 //
37 // Routines that generate pseudo-random numbers compatible
38 // with the standard erand48(), nrand48(), etc. functions.
39 //
40 //-----------------------------------------------------------------------------
41
42 #include "ImathRandom.h"
43 #include "ImathInt64.h"
44
45 IMATH_INTERNAL_NAMESPACE_SOURCE_ENTER
46 namespace {
47
48 //
49 // Static state used by Imath::drand48(), Imath::lrand48() and Imath::srand48()
50 //
51
52 unsigned short staticState[3] = {0, 0, 0};
53
54
55 void
rand48Next(unsigned short state[3])56 rand48Next (unsigned short state[3])
57 {
58 //
59 // drand48() and friends are all based on a linear congruential
60 // sequence,
61 //
62 // x[n+1] = (a * x[n] + c) % m,
63 //
64 // where a and c are as specified below, and m == (1 << 48)
65 //
66
67 static const Int64 a = Int64 (0x5deece66dLL);
68 static const Int64 c = Int64 (0xbLL);
69
70 //
71 // Assemble the 48-bit value x[n] from the
72 // three 16-bit values stored in state.
73 //
74
75 Int64 x = (Int64 (state[2]) << 32) |
76 (Int64 (state[1]) << 16) |
77 Int64 (state[0]);
78
79 //
80 // Compute x[n+1], except for the "modulo m" part.
81 //
82
83 x = a * x + c;
84
85 //
86 // Disassemble the 48 least significant bits of x[n+1] into
87 // three 16-bit values. Discard the 16 most significant bits;
88 // this takes care of the "modulo m" operation.
89 //
90 // We assume that sizeof (unsigned short) == 2.
91 //
92
93 state[2] = (unsigned short)(x >> 32);
94 state[1] = (unsigned short)(x >> 16);
95 state[0] = (unsigned short)(x);
96 }
97
98 } // namespace
99
100
101 double
erand48(unsigned short state[3])102 erand48 (unsigned short state[3])
103 {
104 //
105 // Generate double-precision floating-point values between 0.0 and 1.0:
106 //
107 // The exponent is set to 0x3ff, which indicates a value greater
108 // than or equal to 1.0, and less than 2.0. The 48 most significant
109 // bits of the significand (mantissa) are filled with pseudo-random
110 // bits generated by rand48Next(). The remaining 4 bits are a copy
111 // of the 4 most significant bits of the significand. This results
112 // in bit patterns between 0x3ff0000000000000 and 0x3fffffffffffffff,
113 // which correspond to uniformly distributed floating-point values
114 // between 1.0 and 1.99999999999999978. Subtracting 1.0 from those
115 // values produces numbers between 0.0 and 0.99999999999999978, that
116 // is, between 0.0 and 1.0-DBL_EPSILON.
117 //
118
119 rand48Next (state);
120
121 union {double d; Int64 i;} u;
122
123 u.i = (Int64 (0x3ff) << 52) | // sign and exponent
124 (Int64 (state[2]) << 36) | // significand
125 (Int64 (state[1]) << 20) |
126 (Int64 (state[0]) << 4) |
127 (Int64 (state[2]) >> 12);
128
129 return u.d - 1;
130 }
131
132
133 double
drand48()134 drand48 ()
135 {
136 return IMATH_INTERNAL_NAMESPACE::erand48 (staticState);
137 }
138
139
140 long int
nrand48(unsigned short state[3])141 nrand48 (unsigned short state[3])
142 {
143 //
144 // Generate uniformly distributed integers between 0 and 0x7fffffff.
145 //
146
147 rand48Next (state);
148
149 return ((long int) (state[2]) << 15) |
150 ((long int) (state[1]) >> 1);
151 }
152
153
154 long int
lrand48()155 lrand48 ()
156 {
157 return IMATH_INTERNAL_NAMESPACE::nrand48 (staticState);
158 }
159
160
161 void
srand48(long int seed)162 srand48 (long int seed)
163 {
164 staticState[2] = (unsigned short)(seed >> 16);
165 staticState[1] = (unsigned short)(seed);
166 staticState[0] = 0x330e;
167 }
168
169
170 float
nextf()171 Rand32::nextf ()
172 {
173 //
174 // Generate single-precision floating-point values between 0.0 and 1.0:
175 //
176 // The exponent is set to 0x7f, which indicates a value greater than
177 // or equal to 1.0, and less than 2.0. The 23 bits of the significand
178 // (mantissa) are filled with pseudo-random bits generated by
179 // Rand32::next(). This results in in bit patterns between 0x3f800000
180 // and 0x3fffffff, which correspond to uniformly distributed floating-
181 // point values between 1.0 and 1.99999988. Subtracting 1.0 from
182 // those values produces numbers between 0.0 and 0.99999988, that is,
183 // between 0.0 and 1.0-FLT_EPSILON.
184 //
185
186 next ();
187
188 union {float f; unsigned int i;} u;
189
190 u.i = 0x3f800000 | (_state & 0x7fffff);
191 return u.f - 1;
192 }
193
194 IMATH_INTERNAL_NAMESPACE_SOURCE_EXIT
195