/src/FreeImage/Source/OpenEXR/Imath/ImathRandom.h

  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