ImathRandom.cpp 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194
  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. // Routines that generate pseudo-random numbers compatible
  37. // with the standard erand48(), nrand48(), etc. functions.
  38. //
  39. //-----------------------------------------------------------------------------
  40. #include "ImathRandom.h"
  41. #include "ImathInt64.h"
  42. IMATH_INTERNAL_NAMESPACE_SOURCE_ENTER
  43. namespace {
  44. //
  45. // Static state used by Imath::drand48(), Imath::lrand48() and Imath::srand48()
  46. //
  47. unsigned short staticState[3] = {0, 0, 0};
  48. void
  49. rand48Next (unsigned short state[3])
  50. {
  51. //
  52. // drand48() and friends are all based on a linear congruential
  53. // sequence,
  54. //
  55. // x[n+1] = (a * x[n] + c) % m,
  56. //
  57. // where a and c are as specified below, and m == (1 << 48)
  58. //
  59. static const Int64 a = Int64 (0x5deece66dLL);
  60. static const Int64 c = Int64 (0xbLL);
  61. //
  62. // Assemble the 48-bit value x[n] from the
  63. // three 16-bit values stored in state.
  64. //
  65. Int64 x = (Int64 (state[2]) << 32) |
  66. (Int64 (state[1]) << 16) |
  67. Int64 (state[0]);
  68. //
  69. // Compute x[n+1], except for the "modulo m" part.
  70. //
  71. x = a * x + c;
  72. //
  73. // Disassemble the 48 least significant bits of x[n+1] into
  74. // three 16-bit values. Discard the 16 most significant bits;
  75. // this takes care of the "modulo m" operation.
  76. //
  77. // We assume that sizeof (unsigned short) == 2.
  78. //
  79. state[2] = (unsigned short)(x >> 32);
  80. state[1] = (unsigned short)(x >> 16);
  81. state[0] = (unsigned short)(x);
  82. }
  83. } // namespace
  84. double
  85. erand48 (unsigned short state[3])
  86. {
  87. //
  88. // Generate double-precision floating-point values between 0.0 and 1.0:
  89. //
  90. // The exponent is set to 0x3ff, which indicates a value greater
  91. // than or equal to 1.0, and less than 2.0. The 48 most significant
  92. // bits of the significand (mantissa) are filled with pseudo-random
  93. // bits generated by rand48Next(). The remaining 4 bits are a copy
  94. // of the 4 most significant bits of the significand. This results
  95. // in bit patterns between 0x3ff0000000000000 and 0x3fffffffffffffff,
  96. // which correspond to uniformly distributed floating-point values
  97. // between 1.0 and 1.99999999999999978. Subtracting 1.0 from those
  98. // values produces numbers between 0.0 and 0.99999999999999978, that
  99. // is, between 0.0 and 1.0-DBL_EPSILON.
  100. //
  101. rand48Next (state);
  102. union {double d; Int64 i;} u;
  103. u.i = (Int64 (0x3ff) << 52) | // sign and exponent
  104. (Int64 (state[2]) << 36) | // significand
  105. (Int64 (state[1]) << 20) |
  106. (Int64 (state[0]) << 4) |
  107. (Int64 (state[2]) >> 12);
  108. return u.d - 1;
  109. }
  110. double
  111. drand48 ()
  112. {
  113. return IMATH_INTERNAL_NAMESPACE::erand48 (staticState);
  114. }
  115. long int
  116. nrand48 (unsigned short state[3])
  117. {
  118. //
  119. // Generate uniformly distributed integers between 0 and 0x7fffffff.
  120. //
  121. rand48Next (state);
  122. return ((long int) (state[2]) << 15) |
  123. ((long int) (state[1]) >> 1);
  124. }
  125. long int
  126. lrand48 ()
  127. {
  128. return IMATH_INTERNAL_NAMESPACE::nrand48 (staticState);
  129. }
  130. void
  131. srand48 (long int seed)
  132. {
  133. staticState[2] = (unsigned short)(seed >> 16);
  134. staticState[1] = (unsigned short)(seed);
  135. staticState[0] = 0x330e;
  136. }
  137. float
  138. Rand32::nextf ()
  139. {
  140. //
  141. // Generate single-precision floating-point values between 0.0 and 1.0:
  142. //
  143. // The exponent is set to 0x7f, which indicates a value greater than
  144. // or equal to 1.0, and less than 2.0. The 23 bits of the significand
  145. // (mantissa) are filled with pseudo-random bits generated by
  146. // Rand32::next(). This results in in bit patterns between 0x3f800000
  147. // and 0x3fffffff, which correspond to uniformly distributed floating-
  148. // point values between 1.0 and 1.99999988. Subtracting 1.0 from
  149. // those values produces numbers between 0.0 and 0.99999988, that is,
  150. // between 0.0 and 1.0-FLT_EPSILON.
  151. //
  152. next ();
  153. union {float f; unsigned int i;} u;
  154. u.i = 0x3f800000 | (_state & 0x7fffff);
  155. return u.f - 1;
  156. }
  157. IMATH_INTERNAL_NAMESPACE_SOURCE_EXIT