1 /* Copyright 2011 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 MathFunc.cpp
16         @author Jukka Jyl�nki
17         @brief Common mathematical functions. */
18 #ifdef MATH_ENABLE_STL_SUPPORT
19 #include <utility>
20 #include <algorithm>
21 #endif
22
23 #include "Math/MathFunc.h"
24 #include "Math/float2.h"
25
26 MATH_BEGIN_NAMESPACE
27
28 bool mathBreakOnAssume = false;
29
30 void SetMathBreakOnAssume(bool isEnabled)
31 {
32         mathBreakOnAssume = isEnabled;
33 }
34
35 /// Returns the current state of the math break-on-assume flag.
36 bool MathBreakOnAssume()
37 {
38         return mathBreakOnAssume;
39 }
40
41 float Sin(float angleRadians)
42 {
43         return sin(angleRadians);
44 }
45
46 float Cos(float angleRadians)
47 {
48         return cos(angleRadians);
49 }
50
51 float Tan(float angleRadians)
52 {
53         return tan(angleRadians);
54 }
55
56 float2 SinCos(float angleRadians)
57 {
58         return float2(sin(angleRadians), cos(angleRadians));
59 }
60
61 float Asin(float x)
62 {
63         return asin(x);
64 }
65
66 float Acos(float x)
67 {
68         return acos(x);
69 }
70
71 float Atan(float x)
72 {
73         return atan(x);
74 }
75
76 float Atan2(float y, float x)
77 {
78         return atan2(y, x);
79 }
80
81 float Sinh(float x)
82 {
83         return sinh(x);
84 }
85
86 float Cosh(float x)
87 {
88         return cosh(x);
89 }
90
91 float Tanh(float x)
92 {
93         return tanh(x);
94 }
95
96 bool IsPow2(unsigned int number)
97 {
98         return (number & (number-1)) == 0;
99 }
100
101 unsigned int RoundUpPow2(unsigned int x)
102 {
103         assert(sizeof(unsigned int) <= 4);
104         --x;
105         x |= x >> 1;
106         x |= x >> 2;
107         x |= x >> 4;
108         x |= x >> 8;
109         x |= x >> 16;
110         ++x;
111
112         return x;
113 }
114
115 unsigned int RoundDownPow2(unsigned int x)
116 {
117         assert(sizeof(unsigned int) <= 4);
118         x |= x >> 1;
119         x |= x >> 2;
120         x |= x >> 4;
121         x |= x >> 8;
122         x |= x >> 16;
123         return x - (x >> 1);
124 }
125
126 float Pow(float base, float exponent)
127 {
128         return pow(base, exponent);
129 }
130
131 float Exp(float exponent)
132 {
133         return exp(exponent);
134 }
135
136 float Log(float base, float value)
137 {
138         return log(value) / log(base);
139 }
140
141 float Log2(float value)
142 {
143         return Log(2.f, value);
144 }
145
146 float Ln(float value)
147 {
148         return log(value);
149 }
150
151 float Log10(float value)
152 {
153         return Log(10.f, value);
154 }
155
156 float Ceil(float x)
157 {
158         return ceilf(x);
159 }
160
161 int CeilInt(float x)
162 {
163         return (int)ceilf(x);
164 }
165
166 float Floor(float x)
167 {
168         return floorf(x);
169 }
170
171 int FloorInt(float x)
172 {
173         return (int)floorf(x);
174 }
175
176 float Round(float x)
177 {
178         return Floor(x+0.5f);
179 }
180
181 int RoundInt(float x)
182 {
183         return (int)Round(x);
184 }
185
186 float Sign(float x)
187 {
188         return x >= 0.f ? 1.f : -1.f;
189 }
190
191 float SignOrZero(float x, float epsilon)
192 {
193         return Abs(x) <= epsilon ? 0.f : Sign(x);
194 }
195
196 float Lerp(float a, float b, float t)
197 {
198         return a + t * (b-a);
199 }
200
201 float LerpMod(float a, float b, float mod, float t)
202 {
203         a = ModPos(a, mod);
204         b = ModPos(b, mod);
205         if (Abs(b-a) * 2.f <= mod)
206                 return Lerp(a, b, t);
207         else
208         {
209                 if (a < b)
210                         return ModPos(Lerp(a + mod, b, t), mod);
211                 else
212                         return ModPos(Lerp(a, b + mod, t), mod);
213         }
214 }
215
216 float InvLerp(float a, float b, float x)
217 {
218         assume(Abs(b-a) > eps);
219         return (x - a) / (b - a);
220 }
221
222 float Step(float y, float x)
223 {
224         return (x >= y) ? 1.f : 0.f;
225 }
226
227 float SmoothStep(float min, float max, float x)
228 {
229         return x <= min ? 0.f : (x >= max ? 1.f : (x - min) / (max - min));
230 }
231
232 float PingPongMod(float x, float mod)
233 {
234         x = Mod(x, mod * 2.f);
235         return x >= mod ? (2.f * mod - x) : x;
236 }
237
238 float Mod(float x, float mod)
239 {
240         return fmod(x, mod);
241 }
242
243 float Mod(float x, int mod)
244 {
245         ///@todo Optimize.
246         return fmod(x, (float)mod);
247 }
248
249 float ModPos(float x, float mod)
250 {
251         float m = fmod(x, mod);
252         return m >= 0.f ? m : (m + mod);
253 }
254
255 float ModPos(float x, int mod)
256 {
257         ///@todo Optimize.
258         return ModPos(x, (float)mod);
259 }
260
261 float Frac(float x)
262 {
263         return x - Floor(x);
264 }
265
266 float Sqrt(float x)
267 {
268         return sqrtf(x);
269 }
270
271 float RSqrt(float x)
272 {
273         return 1.f / sqrtf(x);
274 }
275
276 /** Uses a recursive approach, not the fastest/brightest method.
277         Note that 13! = 6227020800 overflows already.
278         @return n! = n * (n-1) * (n-2) * ... * 1. */
279 int Factorial(int n)
280 {
281         int result = 1;
282         for(int i = 2; i <= n; i++)
283                 result *= i;
284         return result;
285 }
286
287 /** @return Binomial coefficients with recursion, i.e. n choose k, C(n,k) or nCk. */
288 int CombinatorialRec(int n, int k)
289 {
290         /* We could do: 
291                         return factorial(n)/(factorial(n-k)*factorial(k));
292                 But prefer the recursive approach instead, because it's not so prone
293                 to numerical overflow. This approach uses the idea of the Pascal triangle. */
294
295         if (k <= 0 || k >= n)
296                 return 1;
297         else
298                 return CombinatorialRec(n-1,k-1) + CombinatorialRec(n-1,k);
299 }
300
301 /** @return Binomial coefficients by tabulation, i.e. n choose k, C(n,k) or nCk. */
302 int CombinatorialTab(int n, int k)
303 {
304         if (k == 0 || k == n)
305                 return 1;
306         if (k < 0 || k > n)
307                 return 0;
308         // We use two auxiliary tables, one of size n-2 and one of size n-3.
309         int *table = new int[2*(k+1)]; ///@todo We can lower this size.
310         for(int i = 0; i < 2*(k+1); ++i)
311                 table[i] = 1;
312         int *t1 = &table[0];
313         int *t2 = &table[k+1];
314         // Iteratively fill the tables.
315         for(int i = 2; i <= n; ++i)
316         {
317                 for(int j = Max(1, i-n+k); j <= Min(k,i-1); ++j)
318                         t1[j] = t2[j] + t2[j-1];
319                 Swap(t1, t2);
320         }
321         int c = t2[k];
322         delete[] table;
323         return c;
324 }
325
326 float PowUInt(float base, u32 exponent)
327 {
328         // 'Fast Exponentiation': We interpret exponent in base two and calculate the power by
329         // squaring and multiplying by base.
330
331         // Find the highest bit that is set.
332         u32 e = 0x80000000;
333         while((exponent & e) == 0 && e > 0)
334                 e >>= 1;
335
336         float val = 1.f;
337         do
338         {
339                 val *= val; // Shifts the exponent one place left
340                 val *= (exponent & e) != 0 ? base : 1.f; // Adds a 1 as the LSB of the exponent
341                 e >>= 1;
342         } while(e > 0);
343
344         return val;
345 }
346
347 /** @param base Exponent base value.
348         @param exponent Integer exponent to raise base to.
349         @return pow(base,exponent) but optimized because we only use integer exponent. */
350 float PowInt(float base, int exponent)
351 {
352         if (exponent == 0)
353                 return 1.f;
354         else if (exponent < 0)
355                 return 1.f / PowUInt(base, (u32)-exponent);
356         else
357                 return PowUInt(base, (u32)exponent);
358 }
359
360 MATH_END_NAMESPACE