initial commit of the fluid simulator.
[blender.git] / intern / elbeem / intern / ntl_rndstream.h
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2003,2004 Nils Thuerey
5  *
6  * A seperate random number stream (e.g. for photon tracing)
7  * (cf. Numerical recipes in C, sec. ed., p 283, Knuth method)
8  *
9  *****************************************************************************/
10 #ifndef NTL_RANDOMSTREAM_H
11
12
13 //! some big number
14 #define MBIG 1000000000
15
16 //! modify initial seed
17 #define MSEED 161803398
18
19 //! minimum value - no idea why this is a define?
20 #define MZ 0
21
22 //! for normalization to 0..1
23 #define FAC (1.0/MBIG)
24
25
26 /*! a stream of random numbers using Knuth's portable method */
27 class ntlRandomStream
28 {
29 public:
30   /*! Default constructor */
31   inline ntlRandomStream(long seed);
32   /*! Destructor */
33   ~ntlRandomStream() {}
34
35         /*! get a random number from the stream */
36         inline double getDouble( void );
37
38 #ifdef HAVE_GFXTYPES
39         //! gfx random functions
40         
41         /*! get a random number from the stream */
42         inline gfxReal getGfxReal( void );
43 #endif
44
45 private:
46
47         /*! random number state */
48         long idnum;
49
50         /*! pointers into number table */ 
51         int inext, inextp;
52         /*! store seed and number for subtraction */
53         long ma[56];
54
55 };
56
57
58 /* init random stream tables */
59 inline ntlRandomStream::ntlRandomStream(long seed)
60 {
61         idnum = seed;
62
63         long mj = MSEED - (idnum < 0 ? -idnum : idnum);
64         mj %= MBIG;
65         ma[55] = mj;
66         long mk = 1;
67
68         // init table once, otherwise strange results...
69         for(int i=0;i<=55;i++) ma[i] = (i*i+seed); 
70
71         // init table in random order
72         for(int i=1;i<=54;i++) {
73                 int ii = (21*i) % 56;
74                 ma[ii] = mk;
75                 mk = mj - mk;
76                 if(mk < MZ) mk += MBIG;
77                 mj = ma[ii];
78         }
79
80         // "warm up" generator
81         for(int k=1;k<=4;k++) 
82                 for(int i=1;i<=55;i++) {
83                         ma[i] -= ma[1+ (i+30) % 55];
84                         if(ma[i] < MZ) ma[i] += MBIG;
85                 }
86
87         inext = 0;
88         inextp = 31; // the special "31"
89         idnum = 1;
90 }
91
92
93 /* return one random value */
94 inline double ntlRandomStream::getDouble( void )
95 {
96         if( ++inext == 56) inext = 1;
97         if( ++inextp == 56) inextp = 1;
98
99         // generate by subtaction
100         long mj = ma[inext] - ma[inextp];
101
102         // check range
103         if(mj < MZ) mj += MBIG;
104         ma[ inext ] = mj;
105         return (double)(mj * FAC);
106 }
107
108 #ifdef HAVE_GFXTYPES
109 /* return one random value */
110 inline gfxReal ntlRandomStream::getGfxReal( void )
111 {
112         if( ++inext == 56) inext = 1;
113         if( ++inextp == 56) inextp = 1;
114
115         // generate by subtaction
116         long mj = ma[inext] - ma[inextp];
117
118         // check range
119         if(mj < MZ) mj += MBIG;
120         ma[ inext ] = mj;
121         return (gfxReal)(mj * FAC);
122 }
123 #endif
124
125 #define NTL_RANDOMSTREAM_H
126 #endif
127