1 /* Copyright Jukka Jyl�nki
2
3    Licensed under the Apache License, Version 2.0 (the "License");
4    you may not use this file except in compliance with the License.
5    You may obtain a copy of the License at
6
7        http://www.apache.org/licenses/LICENSE-2.0
8
9    Unless required by applicable law or agreed to in writing, software
10    distributed under the License is distributed on an "AS IS" BASIS,
11    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12    See the License for the specific language governing permissions and
13    limitations under the License. */
14
15 /** @file LCG.cpp
16         @author Jukka Jyl�nki
17         @brief Implementation of the linear congruential random number generator. */
18
19 #include "LCG.h"
20 #include "../../Math/MathFunc.h"
21 #include "../../Time/Clock.h"
22 #include "../../Math/MathTypes.h"
23 #include "../../Math/myassert.h"
24
25 MATH_BEGIN_NAMESPACE
26
27 LCG::LCG()
28 {
29         Seed(Clock::TickU32());
30 }
31
32 /** If you want to give different parameters for the generator, you should remember that:
33         - modulus should be prime
34         - modulus, increment and the multiplier should all be relative primes (increment can be 0)
35         - modulus should be greater than multiplier and increment
36
37         Most often you can leave increment = 0.
38
39         A list of widely used values:
40         - Park and Miller (Minimal standard): mul = 16807 (7^5)             mod = 2^31 - 1 (2147483647 == 0x7FFFFFFF)
41         - Park and Miller #2:                 mul = 48271                   mod = 2^31 - 1
42         - Park and Miller #3:                 mul = 69621                   mod = 2^31 - 1
43         - SIMSCRIPT:                          mul = 630360016               mod = 2^31 - 1
44         - URN12:                              mul = 452807053               mod = 2^31
45
46         Infamous examples (Don't use!):
47         - Classical ANSI C                    mul = 1103515245  inc = 12345 mod = 2^31
48         - RANDU                               mul = 65539                   mod = 2^31  */
49 void LCG::Seed(u32 seed, u32 mul, u32 inc, u32 mod)
50 {
51         assume((seed != 0 || inc != 0) && "Initializing LCG with seed=0 && inc=0 results in an infinite series of 0s!");
52 #ifndef MATH_SILENT_ASSUME
53         if (inc == 0 && (mul % mod == 0 || mod % mul == 0))
54                 LOGW("Warning: Multiplier %u and modulus %u are not compatible since one is a multiple of the other and the increment == 0!", mul, mod);
55 #endif
56         assume(mul != 0);
57         assume(mod > 1);
58
59         lastNumber = seed;
60         multiplier = mul;
61         increment = inc;
62         modulus = mod;
63 }
64
65 u32 LCG::IntFast()
66 {
67         assert(increment == 0);
68         assert(multiplier % 2 == 1 && "Multiplier should be odd for LCG::IntFast(), since modulus==2^32 is even!");
69 // The configurable modulus and increment are not used by this function.
70         u32 mul = lastNumber * multiplier;
71         lastNumber = mul + (mul <= lastNumber?1:0); // Whenever we overflow, flip by one to avoid even multiplier always producing even results, since modulus is even.
72         assert(lastNumber != 0); // We don't use an adder in IntFast(), so must never degenerate to zero.
73         return lastNumber;
74 }
75
76 u32 LCG::Int()
77 {
78         assert(modulus != 0);
79         /// \todo Convert to using Schrage's method for approximate factorization. (Numerical Recipes in C)
80
81         // Currently we cast everything to 64-bit to avoid overflow, which is quite dumb.
82
83         // Create the new random number
84 //#ifdef WIN32
85         u64 newNum = ((u64)lastNumber * (u64)multiplier + (u64)increment) % (u64)modulus;
86 //      u32 m = lastNumber * multiplier;
87 //      u32 i = m + increment;
88 //      u32 f = i & 0x7FFFFFFF;
89 //      u32 m = (lastNumber * 214013 + 2531011) & 0x7FFFFFFF;
90 //      unsigned __int64 newNum = (lastNumber * multiplier + increment) & 0x7FFFFFFF;
91 //#else
92         // On console platform, we rely on using smaller sequences.
93 //      unsigned long newNum = ((unsigned long)lastNumber * (unsigned long)multiplier + (unsigned long)increment) % (unsigned long)modulus;
94 //#endif
95         // Save the newly generated random number to use as seed for the next one.
96 //      lastNumber = m;//(u32)newNum;
97         assert4((((u32)newNum) != 0 || increment != 0) && "LCG degenerated to producing a stream of zeroes!", lastNumber, multiplier, increment, modulus);
98         lastNumber = (u32)newNum;
99         return lastNumber;
100 }
101
102 int LCG::Int(int a, int b)
103 {
104         assert(a <= b && "Error in range!");
105
106 //      return a + (int)(Int() * Max()/(b-a));
107         int num = a + (int)(Float() * (b-a+1));
108 //      assert(num >= a);
109 //      assert(num <= b);
110         ///\todo Some bug here - the result is not necessarily in the proper range.
111         if (num < a)
112                 num = a;
113         if (num > b)
114                 num = b;
115         return num;
116 }
117
118 float LCG::Float()
119 {
120         u32 i = ((u32)Int() & 0x007FFFFF /* random mantissa */) | 0x3F800000 /* fixed exponent */;
121         float f = ReinterpretAsFloat(i); // f is now in range [1, 2[
122         f -= 1.f; // Map to range [0, 1[
123         assert1(f >= 0.f, f);
124         assert1(f < 1.f, f);
125         return f;
126 }
127
128 float LCG::Float01Incl()
129 {
130         for(int i = 0; i < 100; ++i)
131         {
132                 u32 val = (u32)Int() & 0x00FFFFFF;
133                 if (val > 0x800000)
134                         continue;
135                 else if (val == 0x800000)
136                         return 1.0f;
137                 else
138                 {
139                         val |= 0x3F800000;
140                         float f = ReinterpretAsFloat(val) - 1.f;
141                         assert1(f >= 0.f, f);
142                         assert1(f <= 1.f, f);
143                         return f;
144                 }
145         }
146         return Float();
147 }
148
149 float LCG::FloatNeg1_1()
150 {
151         u32 i = (u32)Int();
152         u32 one = ((i & 0x00800000) << 8) /* random sign bit */ | 0x3F800000 /* fixed exponent */;
153         i = one | (i & 0x007FFFFF) /* random mantissa */;
154         float f = ReinterpretAsFloat(i); // f is now in range ]-2, -1[ union [1, 2].
155         float fone = ReinterpretAsFloat(one); // +/- 1, of same sign as f.
156         f -= fone;
157         assert1(f > -1.f, f);
158         assert1(f < 1.f, f);
159         return f;
160 }
161
162 float LCG::Float(float a, float b)
163 {
164         assume(a <= b && "LCG::Float(a,b): Error in range: b < a!");
165
166         if (a == b)
167                 return a;
168
169         for(int i = 0; i < 10; ++i)
170         {
171                 float f = a + Float() * (b-a);
172                 if (f != b)
173                 {
174                         assume2(a <= f, a, b);
175                         assume2(f < b || a == b, f, b);
176                         return f;
177                 }
178         }
179         return a;
180 }
181
182 float LCG::FloatIncl(float a, float b)
183 {
184         assume(a <= b && "LCG::Float(a,b): Error in range: b < a!");
185
186         float f = a + Float() * (b-a);
187         assume2(a <= f, a, b);
188         assume2(f <= b, f, b);
189         return f;
190 }
191
192 MATH_END_NAMESPACE