1// tRandom.cpp 
2// 
3// Functions necessary to generate pseudo-random numbers using a number of different algorithms. Random bits, integers, 
4// and floats may all be generated. The generators are objects and state isn't shared. This allows multiple generators 
5// that can optionally be run on different threads. Sharing a generator between threads is not supported. One of the 
6// generators implements the Mersenne Twistor algorithm by M. Matsumoto & T. Nishimura. 
7// 
8// Copyright (c) 2005, 2017 Tristan Grimmer. 
9// Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby 
10// granted, provided that the above copyright notice and this permission notice appear in all copies. 
11// 
12// THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL 
13// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, 
14// INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
15// AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR 
16// PERFORMANCE OF THIS SOFTWARE. 
17 
18#include "Math/tRandom.h" 
19using namespace tMath; 
20 
21 
22tRandom::tDefaultGeneratorType tRandom::DefaultGenerator; 
23 
24 
25void tRandom::tGeneratorMersenneTwister::SetSeed(uint32 seed) 
26{ 
27 StateVector[0] = seed; 
28 for (StateIndex = 1; StateIndex < Dimensions; StateIndex++) 
29 StateVector[StateIndex] = 0x6C078965 * (StateVector[StateIndex-1] ^ (StateVector[StateIndex-1] >> 30)) + StateIndex; 
30} 
31 
32 
33void tRandom::tGeneratorMersenneTwister::SetSeed(const uint32* seeds, int numSeeds) 
34{ 
35 tAssert(numSeeds > 0); 
36 SetSeed(uint32(0x54FFC97A)); 
37 
38 int i = 1; int j = 0; 
39 int countMax = tMax(int(Dimensions), numSeeds); 
40 for (int count = 0; count < countMax; count++) 
41 { 
42 StateVector[i] = (StateVector[i] ^ ((StateVector[i-1] ^ (StateVector[i-1] >> 30)) * 0x0019660D)) + seeds[j] + j; 
43 i++; 
44 j++; 
45 if (i >= Dimensions) 
46 { 
47 StateVector[0] = StateVector[Dimensions-1]; 
48 i = 1; 
49 } 
50 
51 if (j >= numSeeds) 
52 j = 0; 
53 } 
54 
55 for (int count = 0; count < Dimensions-1; count++) 
56 { 
57 StateVector[i] = (StateVector[i] ^ ((StateVector[i-1] ^ (StateVector[i-1] >> 30)) * 0x5D588B65)) - i; 
58 if (++i >= Dimensions) 
59 { 
60 StateVector[0] = StateVector[Dimensions-1]; 
61 i = 1; 
62 } 
63 } 
64 
65 StateVector[0] = 0x80000000; 
66} 
67 
68 
69uint32 tRandom::tGeneratorMersenneTwister::GetBits() const 
70{ 
71 if (StateIndex >= Dimensions) 
72 { 
73 uint32 lowerBits = (1LU << LeftShift) - 1; 
74 uint32 upperBits = 0xFFFFFFFF << LeftShift; 
75 
76 const uint32 toggle[2] = { 0, Toggle }; 
77 for (int i = 0; i < Dimensions-PartialDims; i++) 
78 { 
79 uint32 y = (StateVector[i] & upperBits) | (StateVector[i+1] & lowerBits); 
80 StateVector[i] = StateVector[i+PartialDims] ^ (y>>1) ^ toggle[ y&1 ]; 
81 } 
82 
83 for (int i = Dimensions-PartialDims; i < Dimensions-1; i++) 
84 { 
85 uint32 y = (StateVector[i] & upperBits) | (StateVector[i+1] & lowerBits); 
86 StateVector[i] = StateVector[i+PartialDims-Dimensions] ^ (y>>1) ^ toggle[ y&1 ]; 
87 } 
88 
89 uint32 y = (StateVector[Dimensions-1] & upperBits) | (StateVector[0] & lowerBits); 
90 StateVector[Dimensions-1] = StateVector[PartialDims-1] ^ (y>>1) ^ toggle[ y&1 ]; 
91 StateIndex = 0; 
92 } 
93 
94 uint32 bits = StateVector[StateIndex++]; 
95 
96 // Tempering. 
97 bits ^= bits >> TemperShiftA; 
98 bits ^= (bits << TemperShiftB) & TemperB; 
99 bits ^= (bits << TemperShiftC) & TemperC; 
100 bits ^= bits >> TemperShiftD; 
101 
102 return bits; 
103} 
104 
105 
106double tRandom::tGetDouble(const tGenerator& gen) 
107{ 
108 union 
109 { 
110 double real; 
111 uint32 pattern[2]; 
112 } convert; 
113 
114 uint32 bits = gen.GetBits(); 
115 
116 // This assumes knowledge of the IEEE double precision floating point representation. 
117 #ifdef ENDIAN_LITTLE 
118 convert.pattern[0] = bits << 20; 
119 convert.pattern[1] = (bits >> 12) | 0x3FF00000; 
120 #else 
121 convert.pattern[1] = bits << 20; 
122 convert.pattern[0] = (bits >> 12) | 0x3FF00000; 
123 #endif 
124 
125 return convert.real - 1.0; 
126} 
127 
128 
129tVector3 tRandom::tGetDir(const tVector3& centerDir, float extentAngle, const tGenerator& gen) 
130{ 
131 // @todo This doesn't work for arbitrary central directions. For example, centerDir.z should be allowed to be zero. 
132 float theta = tArcTan( tSqrt(centerDir.x*centerDir.x + centerDir.y*centerDir.y), centerDir.z ); 
133 float phi = tArcTan( centerDir.y, centerDir.x ); 
134 float normSpread = (Pi - (tGetFloat(gen)*extentAngle)) / Pi; 
135 
136 // Rotate around the axis. 
137 theta += -((extentAngle / 2.0f)) + (tGetFloat(gen) * extentAngle); 
138 
139 // Rotate away from the axis phi radians at angle theta. 
140 phi += (-extentAngle/2.0f) + (tGetFloat(gen) * normSpread * extentAngle); 
141 
142 float sinTheta = tSin(theta); 
143 float cosTheta = tCos(theta); 
144 float sinPhi = tSin(phi); 
145 float cosPhi = tCos(phi); 
146 return tVector3(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta); 
147} 
148 
149 
150// @todo Make these work. 
151 
152 
153tVector3 tRandom::tGetPointOnTriangle(const tTriangle& tri, const tGenerator& gen) 
154{ 
155 // P = a*A + b*B + c*C with a + b + c = 1. We generate a and b randomly E [0, 1/2]. Choosing 1/2 for the upper 
156 // bound ensures we don't go over 1. http://www.cgafaq.info/wiki/Random_Point_In_Triangle 
157 float a = tGetFloat(gen) / 2.0f; 
158 float b = tGetFloat(gen) / 2.0f; 
159 float c = 1.0f - a - b; 
160 return tri.A*a + tri.B*b + tri.C*c; 
161} 
162 
163 
164tVector3 tRandom::tGetPointOnUnitSquare(const tPlane& plane, const tGenerator& gen) 
165{ 
166 // Requires normalized plane. http://mathworld.wolfram.com/Plane.html 
167 tVector3 u, v; 
168 plane.ComputeOrthogonalBasisVectors(u, v); 
169 
170 float randomU = tGetFloat(gen) - 0.5f; 
171 float randomV = tGetFloat(gen) - 0.5f; 
172 tVector3 point = plane.Normal*plane.Distance + u*randomU + v*randomV; 
173 return point; 
174} 
175 
176 
177tVector3 tRandom::tGetPointOnBox(const tBox& box, const tGenerator& gen) 
178{ 
179 tVector3 extent = box.ComputeExtents(); 
180 tPlane plane; 
181 tVector3 u, v; 
182 
183 int face = tGetBounded(1, 6, gen); 
184 switch (face) 
185 { 
186 case 1: 
187 case 2: 
188 plane.Distance = extent.x; 
189 plane.Normal = tVector3::i; 
190 u = extent.y * tVector3::j; 
191 v = extent.z * tVector3::k; 
192 break; 
193 
194 case 3: 
195 case 4: 
196 plane.Distance = extent.y; 
197 plane.Normal = tVector3::j; 
198 u = extent.x * tVector3::i; 
199 v = extent.z * tVector3::k; 
200 break; 
201 
202 case 5: 
203 case 6: 
204 plane.Distance = extent.z; 
205 plane.Normal = tVector3::k; 
206 u = extent.x * tVector3::i; 
207 v = extent.y * tVector3::j; 
208 break; 
209 } 
210 
211 if (face % 2) 
212 plane.Normal = -plane.Normal; 
213  
214 float randomU = tGetBounded(-1.0f, 1.0f, gen); 
215 float randomV = tGetBounded(-1.0f, 1.0f, gen); 
216 
217 return plane.Normal*plane.Distance + u*randomU + v*randomV; 
218} 
219 
220 
221tVector3 tRandom::tGetPointOnSphere(const tSphere& sphere, const tGenerator& gen) 
222{ 
223 // See http://mathworld.wolfram.com/SphericalCoordinates.html 
224 float theta = tGetBounded(0.0f, TwoPi, gen); 
225 float phi = tGetBounded(0.0f, Pi, gen); 
226  
227 // 's' for Sin, 'c' for Cos. 
228 float ct, st, cp, sp; 
229 tCosSin(ct, st, theta); 
230 tCosSin(cp, sp, phi); 
231 
232 tVector3 point(ct*sp, st*sp, cp); 
233 point *= sphere.Radius; 
234 point += sphere.Center; 
235 return point; 
236} 
237 
238 
239tVector3 tRandom::tGetPointOnCylinder(const tCylinder& cylinder, const tGenerator& gen) 
240{ 
241 // See http://mathworld.wolfram.com/CylindricalCoordinates.html 
242 tVector3 axis = cylinder.B - cylinder.A; 
243 axis.Normalize(); 
244 
245 // This function works by first computing a plane perpendicular to the major cylinder axis 
246 // and placing a circle there. We simply randomly choose how far up the cylinder the plane 
247 // is, and then choose a random angle around the circle. 
248 tPlane plane(axis); 
249 tVector3 u, v; 
250 plane.ComputeOrthogonalBasisVectors(u, v); 
251 
252 axis *= tGetFloat(gen); // Random distance along axis. 
253 float theta = tGetBounded(0.0f, TwoPi, gen); // Random angle around circle. 
254 float st = tSin(theta); 
255 float ct = tCos(theta); 
256 
257 tVector3 point = cylinder.A + axis; // Point in middle of cylinder. 
258 point += cylinder.Radius * (u*st + v*ct); // Point on circle. 
259 return point; 
260} 
261 
262 
263tVector3 tRandom::tGetPointInBox(const tBox& box, const tGenerator& gen) 
264{ 
265 float xr = tGetFloat(gen); 
266 float yr = tGetFloat(gen); 
267 float zr = tGetFloat(gen); 
268 tVector3 offset 
269 ( 
270 xr * (box.Max.x - box.Min.x), 
271 yr * (box.Max.y - box.Min.y), 
272 zr * (box.Max.z - box.Min.z) 
273 ); 
274 return box.Min + offset; 
275} 
276 
277 
278tVector3 tRandom::tGetPointInSphere(const tSphere& sphere, const tGenerator& gen) 
279{ 
280 // See http://mathworld.wolfram.com/SphericalCoordinates.html 
281 float theta = tGetBounded(0.0f, TwoPi, gen); 
282 float phi = tGetBounded(0.0f, Pi, gen); 
283 float radius = tGetBounded(0.0f, sphere.Radius, gen); 
284  
285 // 's' for Sin, 'c' for Cos. 
286 float ct, st, cp, sp; 
287 tCosSin(ct, st, theta); 
288 tCosSin(cp, sp, phi); 
289 
290 tVector3 point(ct*sp, st*sp, cp); 
291 point *= radius; 
292 point += sphere.Center; 
293 return point; 
294} 
295 
296 
297tVector3 tRandom::tGetPointInCylinder(const tCylinder& cylinder, const tGenerator& gen) 
298{ 
299 // See http://mathworld.wolfram.com/CylindricalCoordinates.html 
300 tVector3 axis = cylinder.B - cylinder.A; 
301 axis.Normalize(); 
302 
303 // This function works by first computing a plane perpendicular to the major cylinder axis and placing a circle 
304 // there. We simply randomly choose how far up the cylinder the plane is, and then choose a random angle around 
305 // the circle. 
306 tPlane plane(axis); 
307 tVector3 u, v; 
308 plane.ComputeOrthogonalBasisVectors(u, v); 
309 
310 axis *= tGetFloat(gen); // Random distance along axis. 
311 float radius = tGetBounded(0.0f, cylinder.Radius, gen); // Random radius out. 
312 float theta = tGetBounded(0.0f, TwoPi, gen); // Random angle around circle. 
313 float st = tSin(theta); 
314 float ct = tCos(theta); 
315 
316 tVector3 point = cylinder.A + axis; // Point in middle of cylinder. 
317 point += radius * (u*st + v*ct); // Point on circle. 
318 return point; 
319} 
320