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