196 lines
		
	
	
		
			5.7 KiB
		
	
	
	
		
			C++
		
	
	
	
	
	
			
		
		
	
	
			196 lines
		
	
	
		
			5.7 KiB
		
	
	
	
		
			C++
		
	
	
	
	
	
| 
 | |
| ///////////////////////////////////////////////////////////////////////////
 | |
| //
 | |
| // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
 | |
| // Digital Ltd. LLC
 | |
| //
 | |
| // All rights reserved.
 | |
| //
 | |
| // Redistribution and use in source and binary forms, with or without
 | |
| // modification, are permitted provided that the following conditions are
 | |
| // met:
 | |
| // *       Redistributions of source code must retain the above copyright
 | |
| // notice, this list of conditions and the following disclaimer.
 | |
| // *       Redistributions in binary form must reproduce the above
 | |
| // copyright notice, this list of conditions and the following disclaimer
 | |
| // in the documentation and/or other materials provided with the
 | |
| // distribution.
 | |
| // *       Neither the name of Industrial Light & Magic nor the names of
 | |
| // its contributors may be used to endorse or promote products derived
 | |
| // from this software without specific prior written permission.
 | |
| //
 | |
| // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 | |
| // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 | |
| // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 | |
| // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 | |
| // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 | |
| // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 | |
| // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 | |
| // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 | |
| // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 | |
| // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 | |
| // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 | |
| //
 | |
| ///////////////////////////////////////////////////////////////////////////
 | |
| 
 | |
| //-----------------------------------------------------------------------------
 | |
| //
 | |
| //	Routines that generate pseudo-random numbers compatible
 | |
| //	with the standard erand48(), nrand48(), etc. functions.
 | |
| //
 | |
| //-----------------------------------------------------------------------------
 | |
| 
 | |
| #include "ImathRandom.h"
 | |
| #include "ImathInt64.h"
 | |
| 
 | |
| namespace Imath {
 | |
| namespace {
 | |
| 
 | |
| //
 | |
| // Static state used by Imath::drand48(), Imath::lrand48() and Imath::srand48()
 | |
| //
 | |
| 
 | |
| unsigned short staticState[3] = {0, 0, 0};
 | |
| 
 | |
| 
 | |
| void
 | |
| rand48Next (unsigned short state[3])
 | |
| {
 | |
|     //
 | |
|     // drand48() and friends are all based on a linear congruential
 | |
|     // sequence,
 | |
|     //
 | |
|     //   x[n+1] = (a * x[n] + c) % m,
 | |
|     //
 | |
|     // where a and c are as specified below, and m == (1 << 48)
 | |
|     //
 | |
| 
 | |
|     static const Int64 a = Int64 (0x5deece66dLL);
 | |
|     static const Int64 c = Int64 (0xbLL);
 | |
| 
 | |
|     //
 | |
|     // Assemble the 48-bit value x[n] from the
 | |
|     // three 16-bit values stored in state.
 | |
|     //
 | |
| 
 | |
|     Int64 x = (Int64 (state[2]) << 32) |
 | |
|           (Int64 (state[1]) << 16) |
 | |
|            Int64 (state[0]);
 | |
| 
 | |
|     //
 | |
|     // Compute x[n+1], except for the "modulo m" part.
 | |
|     //
 | |
| 
 | |
|     x = a * x + c;
 | |
| 
 | |
|     //
 | |
|     // Disassemble the 48 least significant bits of x[n+1] into
 | |
|     // three 16-bit values.  Discard the 16 most significant bits;
 | |
|     // this takes care of the "modulo m" operation.
 | |
|     //
 | |
|     // We assume that sizeof (unsigned short) == 2.
 | |
|     //
 | |
| 
 | |
|     state[2] = (unsigned short)(x >> 32);
 | |
|     state[1] = (unsigned short)(x >> 16);
 | |
|     state[0] = (unsigned short)(x);
 | |
| }
 | |
| 
 | |
| } // namespace
 | |
| 
 | |
| 
 | |
| double
 | |
| erand48 (unsigned short state[3])
 | |
| {
 | |
|     //
 | |
|     // Generate double-precision floating-point values between 0.0 and 1.0:
 | |
|     //
 | |
|     // The exponent is set to 0x3ff, which indicates a value greater
 | |
|     // than or equal to 1.0, and less than 2.0.  The 48 most significant
 | |
|     // bits of the significand (mantissa) are filled with pseudo-random
 | |
|     // bits generated by rand48Next().  The remaining 4 bits are a copy
 | |
|     // of the 4 most significant bits of the significand.  This results
 | |
|     // in bit patterns between 0x3ff0000000000000 and 0x3fffffffffffffff,
 | |
|     // which correspond to uniformly distributed floating-point values
 | |
|     // between 1.0 and 1.99999999999999978.  Subtracting 1.0 from those
 | |
|     // values produces numbers between 0.0 and 0.99999999999999978, that
 | |
|     // is, between 0.0 and 1.0-DBL_EPSILON.
 | |
|     //
 | |
| 
 | |
|     rand48Next (state);
 | |
| 
 | |
|     union {double d; Int64 i;} u;
 | |
| 
 | |
|     u.i = (Int64 (0x3ff)    << 52) |	// sign and exponent
 | |
|       (Int64 (state[2]) << 36) |	// significand
 | |
|       (Int64 (state[1]) << 20) |
 | |
|       (Int64 (state[0]) <<  4) |
 | |
|       (Int64 (state[2]) >> 12);
 | |
| 
 | |
|     return u.d - 1;
 | |
| }
 | |
| 
 | |
| 
 | |
| double
 | |
| drand48 ()
 | |
| {
 | |
|     return Imath::erand48 (staticState);
 | |
| }
 | |
| 
 | |
| 
 | |
| long int
 | |
| nrand48 (unsigned short state[3])
 | |
| {
 | |
|     //
 | |
|     // Generate uniformly distributed integers between 0 and 0x7fffffff.
 | |
|     //
 | |
| 
 | |
|     rand48Next (state);
 | |
| 
 | |
|     return ((long int) (state[2]) << 15) |
 | |
|        ((long int) (state[1]) >>  1);
 | |
| }
 | |
| 
 | |
| 
 | |
| long int
 | |
| lrand48 ()
 | |
| {
 | |
|     return Imath::nrand48 (staticState);
 | |
| }
 | |
| 
 | |
| 
 | |
| void
 | |
| srand48 (long int seed)
 | |
| {
 | |
|     staticState[2] = (unsigned short)(seed >> 16);
 | |
|     staticState[1] = (unsigned short)(seed);
 | |
|     staticState[0] = 0x330e;
 | |
| }
 | |
| 
 | |
| 
 | |
| float
 | |
| Rand32::nextf ()
 | |
| {
 | |
|     //
 | |
|     // Generate single-precision floating-point values between 0.0 and 1.0:
 | |
|     //
 | |
|     // The exponent is set to 0x7f, which indicates a value greater than
 | |
|     // or equal to 1.0, and less than 2.0.  The 23 bits of the significand
 | |
|     // (mantissa) are filled with pseudo-random bits generated by
 | |
|     // Rand32::next().  This results in in bit patterns between 0x3f800000
 | |
|     // and 0x3fffffff, which correspond to uniformly distributed floating-
 | |
|     // point values between 1.0 and 1.99999988.  Subtracting 1.0 from
 | |
|     // those values produces numbers between 0.0 and 0.99999988, that is,
 | |
|     // between 0.0 and 1.0-FLT_EPSILON.
 | |
|     //
 | |
| 
 | |
|     next ();
 | |
| 
 | |
|     union {float f; unsigned int i;} u;
 | |
| 
 | |
|     u.i = 0x3f800000 | (_state & 0x7fffff);
 | |
|     return u.f - 1;
 | |
| }
 | |
| 
 | |
| } // namespace Imath
 | 
