1/////////////////////////////////////////////////////////////////////////// 2// 3// Copyright (c) 2002, 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#ifndef INCLUDED_IMATHRANDOM_H 37#define INCLUDED_IMATHRANDOM_H 38 39//----------------------------------------------------------------------------- 40// 41// Generators for uniformly distributed pseudo-random numbers and 42// functions that use those generators to generate numbers with 43// non-uniform distributions: 44// 45// class Rand32 46// class Rand48 47// solidSphereRand() 48// hollowSphereRand() 49// gaussRand() 50// gaussSphereRand() 51// 52// Note: class Rand48() calls erand48() and nrand48(), which are not 53// available on all operating systems. For compatibility we include 54// our own versions of erand48() and nrand48(). Our functions have 55// been reverse-engineered from the corresponding Unix/Linux man page. 56// 57//----------------------------------------------------------------------------- 58 59#include <stdlib.h> 60#include <math.h> 61 62namespace Imath { 63 64//----------------------------------------------- 65// Fast random-number generator that generates 66// a uniformly distributed sequence with a period 67// length of 2^32. 68//----------------------------------------------- 69 70class Rand32 71{ 72 public: 73 74 //------------ 75 // Constructor 76 //------------ 77 78 Rand32 (unsigned long int seed = 0); 79 80 81 //-------------------------------- 82 // Re-initialize with a given seed 83 //-------------------------------- 84 85 void init (unsigned long int seed); 86 87 88 //---------------------------------------------------------- 89 // Get the next value in the sequence (range: [false, true]) 90 //---------------------------------------------------------- 91 92 bool nextb (); 93 94 95 //--------------------------------------------------------------- 96 // Get the next value in the sequence (range: [0 ... 0xffffffff]) 97 //--------------------------------------------------------------- 98 99 unsigned long int nexti (); 100 101 102 //------------------------------------------------------ 103 // Get the next value in the sequence (range: [0 ... 1[) 104 //------------------------------------------------------ 105 106 float nextf (); 107 108 109 //------------------------------------------------------------------- 110 // Get the next value in the sequence (range [rangeMin ... rangeMax[) 111 //------------------------------------------------------------------- 112 113 float nextf (float rangeMin, float rangeMax); 114 115 116 private: 117 118 void next (); 119 120 unsigned long int _state; 121}; 122 123 124//-------------------------------------------------------- 125// Random-number generator based on the C Standard Library 126// functions erand48(), nrand48() & company; generates a 127// uniformly distributed sequence. 128//-------------------------------------------------------- 129 130class Rand48 131{ 132 public: 133 134 //------------ 135 // Constructor 136 //------------ 137 138 Rand48 (unsigned long int seed = 0); 139 140 141 //-------------------------------- 142 // Re-initialize with a given seed 143 //-------------------------------- 144 145 void init (unsigned long int seed); 146 147 148 //---------------------------------------------------------- 149 // Get the next value in the sequence (range: [false, true]) 150 //---------------------------------------------------------- 151 152 bool nextb (); 153 154 155 //--------------------------------------------------------------- 156 // Get the next value in the sequence (range: [0 ... 0x7fffffff]) 157 //--------------------------------------------------------------- 158 159 long int nexti (); 160 161 162 //------------------------------------------------------ 163 // Get the next value in the sequence (range: [0 ... 1[) 164 //------------------------------------------------------ 165 166 double nextf (); 167 168 169 //------------------------------------------------------------------- 170 // Get the next value in the sequence (range [rangeMin ... rangeMax[) 171 //------------------------------------------------------------------- 172 173 double nextf (double rangeMin, double rangeMax); 174 175 176 private: 177 178 unsigned short int _state[3]; 179}; 180 181 182//------------------------------------------------------------ 183// Return random points uniformly distributed in a sphere with 184// radius 1 around the origin (distance from origin <= 1). 185//------------------------------------------------------------ 186 187template <class Vec, class Rand> 188Vec 189solidSphereRand (Rand &rand); 190 191 192//------------------------------------------------------------- 193// Return random points uniformly distributed on the surface of 194// a sphere with radius 1 around the origin. 195//------------------------------------------------------------- 196 197template <class Vec, class Rand> 198Vec 199hollowSphereRand (Rand &rand); 200 201 202//----------------------------------------------- 203// Return random numbers with a normal (Gaussian) 204// distribution with zero mean and unit variance. 205//----------------------------------------------- 206 207template <class Rand> 208float 209gaussRand (Rand &rand); 210 211 212//---------------------------------------------------- 213// Return random points whose distance from the origin 214// has a normal (Gaussian) distribution with zero mean 215// and unit variance. 216//---------------------------------------------------- 217 218template <class Vec, class Rand> 219Vec 220gaussSphereRand (Rand &rand); 221 222 223//--------------------------------- 224// erand48(), nrand48() and friends 225//--------------------------------- 226 227double erand48 (unsigned short state[3]); 228double drand48 (); 229long int nrand48 (unsigned short state[3]); 230long int lrand48 (); 231void srand48 (long int seed); 232 233 234//--------------- 235// Implementation 236//--------------- 237 238 239inline void 240Rand32::init (unsigned long int seed) 241{ 242 _state = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL; 243} 244 245 246inline 247Rand32::Rand32 (unsigned long int seed) 248{ 249 init (seed); 250} 251 252 253inline void 254Rand32::next () 255{ 256 _state = 1664525L * _state + 1013904223L; 257} 258 259 260inline bool 261Rand32::nextb () 262{ 263 next (); 264 // Return the 31st (most significant) bit, by and-ing with 2 ^ 31. 265 return !!(_state & 2147483648UL); 266} 267 268 269inline unsigned long int 270Rand32::nexti () 271{ 272 next (); 273 return _state & 0xffffffff; 274} 275 276 277inline float 278Rand32::nextf (float rangeMin, float rangeMax) 279{ 280 float f = nextf(); 281 return rangeMin * (1 - f) + rangeMax * f; 282} 283 284 285inline void 286Rand48::init (unsigned long int seed) 287{ 288 seed = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL; 289 290 _state[0] = (unsigned short int) (seed & 0xFFFF); 291 _state[1] = (unsigned short int) ((seed >> 16) & 0xFFFF); 292 _state[2] = (unsigned short int) (seed & 0xFFFF); 293} 294 295 296inline 297Rand48::Rand48 (unsigned long int seed) 298{ 299 init (seed); 300} 301 302 303inline bool 304Rand48::nextb () 305{ 306 return Imath::nrand48 (_state) & 1; 307} 308 309 310inline long int 311Rand48::nexti () 312{ 313 return Imath::nrand48 (_state); 314} 315 316 317inline double 318Rand48::nextf () 319{ 320 return Imath::erand48 (_state); 321} 322 323 324inline double 325Rand48::nextf (double rangeMin, double rangeMax) 326{ 327 double f = nextf(); 328 return rangeMin * (1 - f) + rangeMax * f; 329} 330 331 332template <class Vec, class Rand> 333Vec 334solidSphereRand (Rand &rand) 335{ 336 Vec v; 337 338 do 339 { 340 for (unsigned int i = 0; i < Vec::dimensions(); i++) 341 v[i] = (typename Vec::BaseType) rand.nextf (-1, 1); 342 } 343 while (v.length2() > 1); 344 345 return v; 346} 347 348 349template <class Vec, class Rand> 350Vec 351hollowSphereRand (Rand &rand) 352{ 353 Vec v; 354 typename Vec::BaseType length; 355 356 do 357 { 358 for (unsigned int i = 0; i < Vec::dimensions(); i++) 359 v[i] = (typename Vec::BaseType) rand.nextf (-1, 1); 360 361 length = v.length(); 362 } 363 while (length > 1 || length == 0); 364 365 return v / length; 366} 367 368 369template <class Rand> 370float 371gaussRand (Rand &rand) 372{ 373 float x; // Note: to avoid numerical problems with very small 374 float y; // numbers, we make these variables singe-precision 375 float length2; // floats, but later we call the double-precision log() 376 // and sqrt() functions instead of logf() and sqrtf(). 377 do 378 { 379 x = float (rand.nextf (-1, 1)); 380 y = float (rand.nextf (-1, 1)); 381 length2 = x * x + y * y; 382 } 383 while (length2 >= 1 || length2 == 0); 384 385 return x * sqrt (-2 * log (double (length2)) / length2); 386} 387 388 389template <class Vec, class Rand> 390Vec 391gaussSphereRand (Rand &rand) 392{ 393 return hollowSphereRand <Vec> (rand) * gaussRand (rand); 394} 395 396} // namespace Imath 397 398#endif 399