source: NonGTP/OpenEXR/include/Imath/ImathRandom.h @ 855

Revision 855, 11.1 KB checked in by igarcia, 19 years ago (diff)
RevLine 
[855]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//      different distributions:
44//
45//      class Rand32
46//      class Rand48
47//      solidSphereRand()
48//      hollowSphereRand()
49//      gaussRand()
50//      gaussSphereRand()
51//
52//-----------------------------------------------------------------------------
53
54//
55// Here is the copyright for the *rand48() functions implemented for
56// PLATFORM_DARWIN_PPC and PLATFORM_WINDOWS
57//
58
59//
60// Copyright (c) 1993 Martin Birgmeier
61// All rights reserved.
62//
63// You may redistribute unmodified or modified versions of this source
64// code provided that the above copyright notice and this and the
65// following conditions are retained.
66//
67// This software is provided ``as is'', and comes with no warranties
68// of any kind. I shall in no event be liable for anything that happens
69// to anyone/anything when using this software.
70//
71
72#include <stdlib.h>
73#include <math.h>
74
75namespace Imath {
76
77
78//-----------------------------------------------
79// Fast random-number generator that generates
80// a uniformly distributed sequence with a period
81// length of 2^32.
82//-----------------------------------------------
83
84class Rand32
85{
86  public:
87
88    //------------
89    // Constructor
90    //------------
91
92    Rand32 (unsigned long int seed = 0);
93   
94
95    //--------------------------------
96    // Re-initialize with a given seed
97    //--------------------------------
98
99    void                init (unsigned long int seed);
100
101
102    //----------------------------------------------------------
103    // Get the next value in the sequence (range: [false, true])
104    //----------------------------------------------------------
105
106    bool                nextb ();
107
108
109    //---------------------------------------------------------------
110    // Get the next value in the sequence (range: [0 ... 0xffffffff])
111    //---------------------------------------------------------------
112
113    unsigned long int   nexti ();
114
115
116    //------------------------------------------------------
117    // Get the next value in the sequence (range: [0 ... 1[)
118    //------------------------------------------------------
119
120    float               nextf ();
121
122
123    //-------------------------------------------------------------------
124    // Get the next value in the sequence (range [rangeMin ... rangeMax[)
125    //-------------------------------------------------------------------
126
127    float               nextf (float rangeMin, float rangeMax);
128
129
130  private:
131
132    void                next ();
133
134    unsigned long int   _state;
135};
136
137
138//--------------------------------------------------------
139// Random-number generator based on the C Standard Library
140// functions drand48(), lrand48() & company; generates a
141// uniformly distributed sequence.
142//--------------------------------------------------------
143
144class Rand48
145{
146  public:
147
148    //------------
149    // Constructor
150    //------------
151
152    Rand48 (unsigned long int seed = 0);
153   
154
155    //--------------------------------
156    // Re-initialize with a given seed
157    //--------------------------------
158
159    void                init (unsigned long int seed);
160
161
162    //----------------------------------------------------------
163    // Get the next value in the sequence (range: [false, true])
164    //----------------------------------------------------------
165
166    bool                nextb ();
167
168
169    //---------------------------------------------------------------
170    // Get the next value in the sequence (range: [0 ... 0x7fffffff])
171    //---------------------------------------------------------------
172
173    long int            nexti ();
174
175
176    //------------------------------------------------------
177    // Get the next value in the sequence (range: [0 ... 1[)
178    //------------------------------------------------------
179
180    double              nextf ();
181
182
183    //-------------------------------------------------------------------
184    // Get the next value in the sequence (range [rangeMin ... rangeMax[)
185    //-------------------------------------------------------------------
186
187    double              nextf (double rangeMin, double rangeMax);
188
189
190  private:
191
192    unsigned short int  _state[3];
193   
194#if defined( PLATFORM_DARWIN_PPC ) || defined ( PLATFORM_WINDOWS )
195    void                shiftState();
196#endif
197};
198
199
200//------------------------------------------------------------
201// Return random points uniformly distributed in a sphere with
202// radius 1 around the origin (distance from origin <= 1).
203//------------------------------------------------------------
204
205template <class Vec, class Rand>
206Vec             
207solidSphereRand (Rand &rand);
208
209
210//-------------------------------------------------------------
211// Return random points uniformly distributed on the surface of
212// a sphere with radius 1 around the origin.
213//-------------------------------------------------------------
214
215template <class Vec, class Rand>
216Vec             
217hollowSphereRand (Rand &rand);
218
219
220//-----------------------------------------------
221// Return random numbers with a normal (Gaussian)
222// distribution with zero mean and unit variance.
223//-----------------------------------------------
224
225template <class Rand>
226float
227gaussRand (Rand &rand);
228
229
230//----------------------------------------------------
231// Return random points whose distance from the origin
232// has a normal (Gaussian) distribution with zero mean
233// and unit variance.
234//----------------------------------------------------
235
236template <class Vec, class Rand>
237Vec
238gaussSphereRand (Rand &rand);
239
240
241//---------------
242// Implementation
243//---------------
244
245
246inline void
247Rand32::init (unsigned long int seed)
248{
249    _state = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL;
250}
251
252
253inline
254Rand32::Rand32 (unsigned long int seed)
255{
256    init (seed);
257}
258
259
260inline void
261Rand32::next ()
262{
263    _state = 1664525L * _state + 1013904223L;
264}
265
266
267inline bool
268Rand32::nextb ()
269{
270    next ();
271    // Return the 31st (most significant) bit, by and-ing with 2 ^ 31.
272    return !!(_state & 2147483648UL);
273}
274
275
276inline unsigned long int
277Rand32::nexti ()
278{
279    next ();
280    return _state & 0xffffffff;
281}
282
283
284inline float
285Rand32::nextf ()
286{
287    next ();
288    return ((int) (_state & 0xffffff)) * ((float) (1.0F / 0x1000000));
289}
290
291
292inline float
293Rand32::nextf (float rangeMin, float rangeMax)
294{
295    return rangeMin + nextf() * (rangeMax - rangeMin);
296}
297
298
299inline void
300Rand48::init (unsigned long int seed)
301{
302    seed = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL;
303
304    _state[0] = (unsigned short int) (seed);
305    _state[1] = (unsigned short int) (seed >> 16);
306    _state[2] = (unsigned short int) (seed);
307}
308
309
310inline
311Rand48::Rand48 (unsigned long int seed)
312{
313    init (seed);
314}
315
316
317#if defined( PLATFORM_DARWIN_PPC ) || defined ( PLATFORM_WINDOWS )
318
319inline void
320Rand48::shiftState()
321{
322    unsigned long   accu;
323    unsigned short  temp[2];
324
325    accu = 0xe66dUL * ( unsigned long )_state[0] + 0x000bUL;
326
327    temp[0] = ( unsigned short )accu;   /* lower 16 bits */
328    accu >>= sizeof( unsigned short ) * 8;
329
330    accu += 0xe66dUL * ( unsigned long )_state[1] +
331            0xdeecUL * ( unsigned long )_state[0];
332
333    temp[1] = ( unsigned short )accu;   /* middle 16 bits */
334    accu >>= sizeof( unsigned short ) * 8;
335
336    accu += 0xe66dUL * _state[2] +
337            0xdeecUL * _state[1] +
338            0x0005UL * _state[0];
339
340    _state[0] = temp[0];
341    _state[1] = temp[1];
342    _state[2] = ( unsigned short )accu;
343}
344
345#endif
346
347inline bool
348Rand48::nextb ()
349{
350#if defined( PLATFORM_DARWIN_PPC ) || defined ( PLATFORM_WINDOWS )
351    shiftState();
352    return ( ( long( _state[2] ) << 15 ) + ( long( _state[1] ) >> 1 ) ) & 0x1;
353#else
354    return nrand48 (_state) & 1;
355#endif
356}
357
358
359inline long int
360Rand48::nexti ()
361{
362#if defined( PLATFORM_DARWIN_PPC ) || defined ( PLATFORM_WINDOWS )
363    shiftState();
364    return ( long( _state[2] ) << 15 ) + ( long( _state[1] ) >> 1 );
365#else
366    return nrand48 (_state);
367#endif
368}
369
370
371inline double
372Rand48::nextf ()
373{
374#if defined( PLATFORM_DARWIN_PPC ) || defined ( PLATFORM_WINDOWS )
375    shiftState();
376    return ldexp( double( _state[0] ), -48 ) +
377           ldexp( double( _state[1] ), -32 ) +
378           ldexp( double( _state[2] ), -16 );
379#else
380    return erand48 (_state);
381#endif
382}
383
384
385inline double
386Rand48::nextf (double rangeMin, double rangeMax)
387{
388    return rangeMin + nextf() * (rangeMax - rangeMin);
389}
390
391
392template <class Vec, class Rand>
393Vec
394solidSphereRand (Rand &rand)
395{
396    Vec v;
397
398    do
399    {
400        for (unsigned int i = 0; i < Vec::dimensions(); i++)
401            v[i] = (typename Vec::BaseType) rand.nextf (-1, 1);
402    }
403    while (v.length2() > 1);
404
405    return v;
406}
407
408
409template <class Vec, class Rand>
410Vec
411hollowSphereRand (Rand &rand)
412{
413    Vec v;
414    typename Vec::BaseType length;
415
416    do
417    {
418        for (unsigned int i = 0; i < Vec::dimensions(); i++)
419            v[i] = (typename Vec::BaseType) rand.nextf (-1, 1);
420
421        length = v.length();
422    }
423    while (length > 1 || length == 0);
424
425    return v / length;
426}
427
428
429template <class Rand>
430float
431gaussRand (Rand &rand)
432{
433    float x;            // Note: to avoid numerical problems with very small
434    float y;            // numbers, we make these variables singe-precision
435    float length2;      // floats, but later we call the double-precision log()
436                        // and sqrt() functions instead of logf() and sqrtf().
437    do
438    {
439        x = float (rand.nextf (-1, 1));
440        y = float (rand.nextf (-1, 1));
441        length2 = x * x + y * y;
442    }
443    while (length2 >= 1 || length2 == 0);
444
445    return x * sqrt (-2 * log (length2) / length2);
446}
447
448
449template <class Vec, class Rand>
450Vec
451gaussSphereRand (Rand &rand)
452{
453    return hollowSphereRand <Vec> (rand) * gaussRand (rand);
454}
455
456double drand48();
457long int lrand48();
458
459} // namespace Imath
460
461#endif
Note: See TracBrowser for help on using the repository browser.