ImathRandom.h 9.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401
  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. #ifndef INCLUDED_IMATHRANDOM_H
  35. #define INCLUDED_IMATHRANDOM_H
  36. //-----------------------------------------------------------------------------
  37. //
  38. // Generators for uniformly distributed pseudo-random numbers and
  39. // functions that use those generators to generate numbers with
  40. // non-uniform distributions:
  41. //
  42. // class Rand32
  43. // class Rand48
  44. // solidSphereRand()
  45. // hollowSphereRand()
  46. // gaussRand()
  47. // gaussSphereRand()
  48. //
  49. // Note: class Rand48() calls erand48() and nrand48(), which are not
  50. // available on all operating systems. For compatibility we include
  51. // our own versions of erand48() and nrand48(). Our functions have
  52. // been reverse-engineered from the corresponding Unix/Linux man page.
  53. //
  54. //-----------------------------------------------------------------------------
  55. #include "ImathNamespace.h"
  56. #include "ImathExport.h"
  57. #include <stdlib.h>
  58. #include <math.h>
  59. IMATH_INTERNAL_NAMESPACE_HEADER_ENTER
  60. //-----------------------------------------------
  61. // Fast random-number generator that generates
  62. // a uniformly distributed sequence with a period
  63. // length of 2^32.
  64. //-----------------------------------------------
  65. class IMATH_EXPORT Rand32
  66. {
  67. public:
  68. //------------
  69. // Constructor
  70. //------------
  71. Rand32 (unsigned long int seed = 0);
  72. //--------------------------------
  73. // Re-initialize with a given seed
  74. //--------------------------------
  75. void init (unsigned long int seed);
  76. //----------------------------------------------------------
  77. // Get the next value in the sequence (range: [false, true])
  78. //----------------------------------------------------------
  79. bool nextb ();
  80. //---------------------------------------------------------------
  81. // Get the next value in the sequence (range: [0 ... 0xffffffff])
  82. //---------------------------------------------------------------
  83. unsigned long int nexti ();
  84. //------------------------------------------------------
  85. // Get the next value in the sequence (range: [0 ... 1[)
  86. //------------------------------------------------------
  87. float nextf ();
  88. //-------------------------------------------------------------------
  89. // Get the next value in the sequence (range [rangeMin ... rangeMax[)
  90. //-------------------------------------------------------------------
  91. float nextf (float rangeMin, float rangeMax);
  92. private:
  93. void next ();
  94. unsigned long int _state;
  95. };
  96. //--------------------------------------------------------
  97. // Random-number generator based on the C Standard Library
  98. // functions erand48(), nrand48() & company; generates a
  99. // uniformly distributed sequence.
  100. //--------------------------------------------------------
  101. class Rand48
  102. {
  103. public:
  104. //------------
  105. // Constructor
  106. //------------
  107. Rand48 (unsigned long int seed = 0);
  108. //--------------------------------
  109. // Re-initialize with a given seed
  110. //--------------------------------
  111. void init (unsigned long int seed);
  112. //----------------------------------------------------------
  113. // Get the next value in the sequence (range: [false, true])
  114. //----------------------------------------------------------
  115. bool nextb ();
  116. //---------------------------------------------------------------
  117. // Get the next value in the sequence (range: [0 ... 0x7fffffff])
  118. //---------------------------------------------------------------
  119. long int nexti ();
  120. //------------------------------------------------------
  121. // Get the next value in the sequence (range: [0 ... 1[)
  122. //------------------------------------------------------
  123. double nextf ();
  124. //-------------------------------------------------------------------
  125. // Get the next value in the sequence (range [rangeMin ... rangeMax[)
  126. //-------------------------------------------------------------------
  127. double nextf (double rangeMin, double rangeMax);
  128. private:
  129. unsigned short int _state[3];
  130. };
  131. //------------------------------------------------------------
  132. // Return random points uniformly distributed in a sphere with
  133. // radius 1 around the origin (distance from origin <= 1).
  134. //------------------------------------------------------------
  135. template <class Vec, class Rand>
  136. Vec
  137. solidSphereRand (Rand &rand);
  138. //-------------------------------------------------------------
  139. // Return random points uniformly distributed on the surface of
  140. // a sphere with radius 1 around the origin.
  141. //-------------------------------------------------------------
  142. template <class Vec, class Rand>
  143. Vec
  144. hollowSphereRand (Rand &rand);
  145. //-----------------------------------------------
  146. // Return random numbers with a normal (Gaussian)
  147. // distribution with zero mean and unit variance.
  148. //-----------------------------------------------
  149. template <class Rand>
  150. float
  151. gaussRand (Rand &rand);
  152. //----------------------------------------------------
  153. // Return random points whose distance from the origin
  154. // has a normal (Gaussian) distribution with zero mean
  155. // and unit variance.
  156. //----------------------------------------------------
  157. template <class Vec, class Rand>
  158. Vec
  159. gaussSphereRand (Rand &rand);
  160. //---------------------------------
  161. // erand48(), nrand48() and friends
  162. //---------------------------------
  163. IMATH_EXPORT double erand48 (unsigned short state[3]);
  164. IMATH_EXPORT double drand48 ();
  165. IMATH_EXPORT long int nrand48 (unsigned short state[3]);
  166. IMATH_EXPORT long int lrand48 ();
  167. IMATH_EXPORT void srand48 (long int seed);
  168. //---------------
  169. // Implementation
  170. //---------------
  171. inline void
  172. Rand32::init (unsigned long int seed)
  173. {
  174. _state = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL;
  175. }
  176. inline
  177. Rand32::Rand32 (unsigned long int seed)
  178. {
  179. init (seed);
  180. }
  181. inline void
  182. Rand32::next ()
  183. {
  184. _state = 1664525L * _state + 1013904223L;
  185. }
  186. inline bool
  187. Rand32::nextb ()
  188. {
  189. next ();
  190. // Return the 31st (most significant) bit, by and-ing with 2 ^ 31.
  191. return !!(_state & 2147483648UL);
  192. }
  193. inline unsigned long int
  194. Rand32::nexti ()
  195. {
  196. next ();
  197. return _state & 0xffffffff;
  198. }
  199. inline float
  200. Rand32::nextf (float rangeMin, float rangeMax)
  201. {
  202. float f = nextf();
  203. return rangeMin * (1 - f) + rangeMax * f;
  204. }
  205. inline void
  206. Rand48::init (unsigned long int seed)
  207. {
  208. seed = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL;
  209. _state[0] = (unsigned short int) (seed & 0xFFFF);
  210. _state[1] = (unsigned short int) ((seed >> 16) & 0xFFFF);
  211. _state[2] = (unsigned short int) (seed & 0xFFFF);
  212. }
  213. inline
  214. Rand48::Rand48 (unsigned long int seed)
  215. {
  216. init (seed);
  217. }
  218. inline bool
  219. Rand48::nextb ()
  220. {
  221. return nrand48 (_state) & 1;
  222. }
  223. inline long int
  224. Rand48::nexti ()
  225. {
  226. return nrand48 (_state);
  227. }
  228. inline double
  229. Rand48::nextf ()
  230. {
  231. return erand48 (_state);
  232. }
  233. inline double
  234. Rand48::nextf (double rangeMin, double rangeMax)
  235. {
  236. double f = nextf();
  237. return rangeMin * (1 - f) + rangeMax * f;
  238. }
  239. template <class Vec, class Rand>
  240. Vec
  241. solidSphereRand (Rand &rand)
  242. {
  243. Vec v;
  244. do
  245. {
  246. for (unsigned int i = 0; i < Vec::dimensions(); i++)
  247. v[i] = (typename Vec::BaseType) rand.nextf (-1, 1);
  248. }
  249. while (v.length2() > 1);
  250. return v;
  251. }
  252. template <class Vec, class Rand>
  253. Vec
  254. hollowSphereRand (Rand &rand)
  255. {
  256. Vec v;
  257. typename Vec::BaseType length;
  258. do
  259. {
  260. for (unsigned int i = 0; i < Vec::dimensions(); i++)
  261. v[i] = (typename Vec::BaseType) rand.nextf (-1, 1);
  262. length = v.length();
  263. }
  264. while (length > 1 || length == 0);
  265. return v / length;
  266. }
  267. template <class Rand>
  268. float
  269. gaussRand (Rand &rand)
  270. {
  271. float x; // Note: to avoid numerical problems with very small
  272. float y; // numbers, we make these variables singe-precision
  273. float length2; // floats, but later we call the double-precision log()
  274. // and sqrt() functions instead of logf() and sqrtf().
  275. do
  276. {
  277. x = float (rand.nextf (-1, 1));
  278. y = float (rand.nextf (-1, 1));
  279. length2 = x * x + y * y;
  280. }
  281. while (length2 >= 1 || length2 == 0);
  282. return x * sqrt (-2 * log (double (length2)) / length2);
  283. }
  284. template <class Vec, class Rand>
  285. Vec
  286. gaussSphereRand (Rand &rand)
  287. {
  288. return hollowSphereRand <Vec> (rand) * gaussRand (rand);
  289. }
  290. IMATH_INTERNAL_NAMESPACE_HEADER_EXIT
  291. #endif // INCLUDED_IMATHRANDOM_H