source: Revenant/marslib/include/mars_calc.h@ 25c4774

port/mars-tycoon
Last change on this file since 25c4774 was 80a6a52, checked in by Jonathan Neufeld <support@…>, 3 years ago

Get to a compile state for terrain procedural generation
  • Property mode set to 100755
File size: 16.7 KB
RevLine 
[80a6a52]1#pragma once
2
3#include <math.h>
4#include <stdlib.h>
5#include <exception>
6#include <assert.h>
7
8#include "mars_meta.h"
9
10namespace mars
11{
12 const double PI = 3.14159265358979323846264338327950;
13
14 template <typename T>
15 inline T SQ(const T & a) { return a * a; }
16 template <typename T>
17 inline T CUBE(const T & a) { return a * a * a; }
18 inline int RAND(const int n) { return (rand() * (n - 1) + (RAND_MAX / 2)) / RAND_MAX; }
19 inline double toRadians(int deg) { return static_cast <double> (deg) / 180.0 * PI; }
20 inline short toDegrees(double rad) { return static_cast <short> (rad / PI * 180.0); }
21 template< typename T >
22 inline T RANDf(const T a = 1.0f) { return static_cast <T> (rand()) / static_cast <T> (RAND_MAX) * a; }
23 inline int RNDi(double f) { return static_cast <int> (f < 0 ? f - 0.5f : f + 0.5f); } // TODO: Optimize tertiary operator
24 template < typename T >
25 inline T ALIGN(const T x, const T unit)
26 {
27 assert(unit > 0);
28 return x - (x % unit) - (x < 0 ? unit : 0);
29 }
30 template <typename T>
31 inline T AVG(const T a, const T b) { return (a + b) / static_cast< T > (2); }
32 template< typename T, typename UT >
33 inline UT WRAP(const T x, const UT w) { return static_cast< UT > ((x % static_cast< T > (w) + static_cast< T > (w)) % w); }
34
35 static const double ABSCISSAE_5[5][2] = {
36 -0.90617984593866399280, 0.23692688505618908751,
37 -0.53846931010568309104, 0.47862867049936646804,
38 0.0, 0.56888888888888888889,
39 +0.53846931010568309104, 0.47862867049936646804,
40 +0.90617984593866399280, 0.23692688505618908751
41 };
42
43 static const char LOG2_TABLE_256[256] =
44 {
45 #define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
46 -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
47 LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
48 LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
49 };
50
51 inline unsigned int LOG2(const unsigned int v)
52 {
53 unsigned int t, tt;
54
55 if (tt = v >> 16)
56 return (t = tt >> 8) ? 24 + LOG2_TABLE_256[t] : 16 + LOG2_TABLE_256[tt];
57 else
58 return (t = v >> 8) ? 8 + LOG2_TABLE_256[t] : LOG2_TABLE_256[v];
59 }
60
61 class CalcEx : public std::exception {};
62
63 template <typename T, typename Derived, unsigned NN>
64 class GaussLegendreInt
65 {
66 private:
67 template<unsigned N>
68 T abscissa (const unsigned int i, const unsigned int j, quantity<N>) const
69 {
70 throw CalcEx ();
71 }
72 T abscissa (const unsigned int i, const unsigned int j, quantity<5>) const
73 {
74 return static_cast <T> (ABSCISSAE_5[i][j]);
75 }
76
77 template <unsigned N>
78 T abscissa (const unsigned int i, const unsigned int j) const
79 {
80 return abscissa(i, j, quantity<N>());
81 }
82
83 template <unsigned N>
84 T summate (const T val, const T a, const T b, quantity<N>) const
85 {
86 const Derived * derived = static_cast <const Derived *> (this);
87 return val + summate <N - 1> (
88 abscissa <NN> (N - 1, 1) *
89 derived->f(
90 (b + a) / 2 +
91 ((b - a) / 2) * abscissa <NN> (N - 1, 0)
92 ),
93 a, b
94 );
95 }
96 T summate (const T val, const T a, const T b, quantity<0>) const { return val; }
97
98 template <unsigned N>
99 T summate (const T val, const T a, const T b) const
100 {
101 return summate(val, a, b, quantity<N>());
102 }
103
104 protected:
105 inline T compute (const T a, const T b) const
106 {
107 return ((b - a) / 2) * summate <NN> (0, a, b);
108 }
109 };
110
111 template <typename T, class EqLeft, class EqRight>
112 class SolutionEquation
113 {
114 private:
115 const EqLeft & _left;
116 const EqRight & _right;
117
118 public:
119 inline SolutionEquation (const EqLeft & left, const EqRight & right)
120 : _left(left), _right(right) {}
121
122 inline const T f(const T x) const { return _left.f(x) - _right.f(x); }
123 inline const T df(const T x) const { return _left.df(x) - _right.df(x); }
124 inline const T f2(const T x) const { return _left.f2(x) - _right.f2(x); }
125 };
126
127 template <typename T, class Eq, unsigned NN = 5>
128 class NewtonsMethod
129 {
130 private:
131 const Eq & _eq;
132
133 template <unsigned N> T iterate (const T x0, quantity<N>) const
134 {
135 return iterate <N-1> (x0 - _eq.f(x0) / _eq.df(x0));
136 }
137 T iterate (const T x0, quantity<0>) const { return x0; }
138
139 template <unsigned N> T iterate (const T x0) const
140 {
141 return iterate(x0, quantity<N>());
142 }
143
144 inline NewtonsMethod (const Eq & eq)
145 : _eq(eq) {}
146
147 inline T solve (const T x0) const { return iterate <NN> (x0); }
148
149 public:
150 // TODO: Cache the result!
151 static inline T solve (const Eq & eq, const T x0) { return NewtonsMethod (eq).solve(x0); }
152 static inline T solve (const Eq & eq) { return NewtonsMethod (eq).solve(0); }
153 };
154
155 template <typename T, class Eq, unsigned NN = 5>
156 class BisectMethod
157 {
158 private:
159 const Eq & _eq;
160
161 template <unsigned N> T iterate (const T l, const T r, quantity<N>) const
162 {
163 T m = (l + r) / 2;
164
165 if (_eq.f(l) * _eq.f(m) > 0)
166 return iterate <N-1> (m, r);
167 else
168 return iterate <N-1> (l, m);
169 }
170 T iterate (const T l, const T r, quantity<0>) const
171 {
172 return (r + l) / 2;
173 }
174
175 template <unsigned N> T iterate (const T l, const T r) const
176 {
177 return iterate(l, r, quantity<N>());
178 }
179
180 inline BisectMethod (const Eq & eq)
181 : _eq(eq) {}
182
183 public:
184 static inline T solve (const Eq & eq, const T l, const T r) { return BisectMethod(eq).iterate <NN> (l, r); }
185 };
186
187 template <typename T, int DIR = +1>
188 class ParabolicFn
189 {
190 public:
191 T scale, broadness, depth;
192
193 ParabolicFn () : scale(0), broadness(0), depth(0) {}
194 ParabolicFn (const T scale, const T broadness, const T depth)
195 : scale(scale), broadness(broadness), depth(depth) {}
196
197 inline const T f(const T x) const { return DIR * (mars::SQ(x / scale) / broadness + depth) * scale; }
198 inline const T f(const T x, const T y) const
199 {
200 return DIR * ((mars::SQ(x / scale) + mars::SQ(y / scale)) / broadness + depth) * scale;
201 }
202 inline const T invf (const T y) const { return sqrt((y / scale * DIR - depth) * broadness) * scale; }
203 inline const T df(const T x) const { return DIR * (2 * x / mars::SQ(scale)) / broadness * scale; }
204 inline const T f2(const T x) const { return DIR * (mars::CUBE(x) / 3 / mars::SQ(scale) / broadness + x * depth) * scale; }
205 };
206
207 template <typename T, unsigned NN = 5>
208 class ErrorFn
209 {
210 private:
211 T _sqrtPI;
212
213 template <unsigned L> static T computeEpsilon () { return computeEpsilon(quantity<L>()); }
214 template <unsigned M> static T s (const T j) {
215 return s(j, quantity<M>());
216 }
217 template <unsigned L> static T computeEpsilon (quantity<L>) { return s<L>(0); }
218 inline static T computeEpsilon (quantity<0>) { return 1; }
219 template <unsigned L, unsigned M> static T s (const T j, quantity<L>, quantity<M>)
220 { return j + s<M - 1> (computeEpsilon<L - 1 - (M-1)>() * computeEpsilon<M-1>() / (((M-1) + 1)*(2 * (M-1) + 1))); }
221 template <unsigned M > static T s (const T j, quantity<0>, quantity<M>) { return j; }
222
223 template <unsigned N> T prod (const T j, const T z, quantity<N>) const
224 {
225 const int k = (NN - N + 1);
226
227 return j * prod <N - 1> (-(2 * k - 1) * z * z / (k * (2 * k + 1)), z);
228 }
229 inline T prod (const T j, const T z, quantity<0>) const { return j; }
230 template <unsigned N> T prod (const T j, const T z) const { return prod(j, z, quantity<N>()); }
231
232 template <unsigned N> T taylor (const T j, const T z, quantity<N>) const
233 { return j + taylor<N - 1> (prod <N>(1, z), z); }
234 inline T taylor (const T j, const T z, quantity<0>) const { return j; }
235 template <unsigned N> T taylor (const T j, const T z) const { return taylor(j, z, quantity<N>()); }
236
237 template <unsigned N> T invfsum (const T j, const T z, quantity<N>) const
238 { return j + invfsum <N - 1> (computeEpsilon<N-1>() / (2*(N-1) + 1) * pow(static_cast <T> (z * _sqrtPI/2), static_cast <int> (2*(N-1) + 1)), z); }
239 inline T invfsum (const T j, const T z, quantity<0>) const { return j; }
240 template <unsigned N> T invfsum (const T j, const T z) const { return invfsum(j, z, quantity<N>()); }
241
242 static T createRootPI () { return static_cast< T > (sqrt(mars::PI)); }
243
244 public:
245 T scale, depth, shift;
246
247 ErrorFn ()
248 : scale(0), depth(0), shift(0), _sqrtPI(createRootPI()) {}
249 ErrorFn (const T frScale, const T fDepth, const T fShift)
250 : scale(frScale), depth (fDepth), shift (fShift), _sqrtPI(createRootPI()) {}
251
252 inline ErrorFn & operator = (const ErrorFn & copy)
253 {
254 const_cast< T & > (scale) = copy.scale;
255 const_cast< T & > (depth) = copy.depth;
256 const_cast< T & > (shift) = copy.shift;
257 return *this;
258 }
259
260 inline const T f (const T x) const { return (static_cast <T> (2.0) / _sqrtPI * taylor <NN> (0, (x - shift) / scale) + depth) * scale; } // TODO: Redundant factors for readability
261 inline const T invf (const T y) const { return invfsum <NN + 1> (0, y / scale - depth) * scale + shift; }
262 inline const T df (const T x) const { return static_cast <T> (2.0) / _sqrtPI * exp(-SQ((x - shift) / scale)) * scale; }
263 inline const T f2 (const T x) const { return (exp(-SQ((x - shift) / scale)) / _sqrtPI + depth * ((x - shift) / scale)) * scale; }
264 };
265
266 template <typename T, int DIR = +1, unsigned NN = 5>
267 class GaussianFn
268 {
269 private:
270 ErrorFn <T, NN> _errorfn;
271
272 protected:
273 inline const T f_base (const T x) const
274 {
275 return f_base(x, broadness, amplitude);
276 }
277
278 inline static const T f_base(const T x, const T broadness, const T amplitude)
279 {
280 return amplitude * exp(
281 -mars::SQ(x) /
282 (2.0f * mars::SQ(broadness))
283 );
284 }
285 inline static const T f_base(const T x, const T y, const T bx, const T by, const T amplitude)
286 {
287 return amplitude * exp(
288 -(
289 mars::SQ(x) / (2.0f * mars::SQ(bx)) +
290 mars::SQ(y) / (2.0f * mars::SQ(by))
291 )
292 );
293 }
294
295 public:
296 const T scale, broadness, depth, amplitude;
297
298 GaussianFn ()
299 : scale(0), amplitude(0), broadness(0), depth(0) {}
300 GaussianFn (const T scale, const T amplitude = 1, const T broadness = 1, const T depth = 0)
301 : scale(scale), amplitude(amplitude), broadness(broadness), depth(depth),
302 _errorfn (scale, depth, 0) {}
303
304 inline GaussianFn & operator = (const GaussianFn & copy)
305 {
306 _errorfn = copy._errorfn;
307 const_cast< T & > (scale) = copy.scale;
308 const_cast< T & > (broadness) = copy.broadness;
309 const_cast< T & > (depth) = copy.depth;
310 const_cast< T & > (amplitude) = copy.amplitude;
311 return *this;
312 }
313
314 inline const T f(const T x) const
315 {
316 return DIR * (scale * (f_base(x / scale, broadness, amplitude) + depth));
317 }
318 inline static const T f (const T x, const T scale, const T broadness, const T amplitude, const T depth)
319 {
320 return DIR * scale * (f_base(x / scale, broadness, amplitude) + depth);
321 }
322
323 inline const T f(const T x, const T y) const
324 {
325 return DIR * scale * (f_base(x / scale, y / scale, broadness, broadness, amplitude) + depth);
326 }
327
328 // Inverse function, warning this is computationally intensive
329 // Valid for all y > 0
330 inline const T invf (const T y) const
331 {
332 return sqrt(log(pow((y / scale * DIR - depth) / amplitude, -2*mars::SQ(broadness)))) * scale;
333 }
334
335 // The integral of the Gaussian is the error-function, but we will represent it using the iterative Taylor series,
336 // which is the product of a sequence. We will implement this using a recursive template algorithm.
337 inline const T f2(const T x) const
338 {
339 return _errorfn.f(x);
340 }
341
342 // First-derivative of the Gaussian is the Gaussian multiplied by the first Hermite polynomial (or x)
343 inline const T df (const T x) const
344 {
345 return DIR * (scale * (f_base(x / scale) * (x / scale)));
346 }
347 };
348
349 template <typename T, int DIR>
350 class InvertedSquareFn
351 {
352 public:
353 T scale, broadness, depth;
354
355 inline InvertedSquareFn (const T scale, const T broadness, const T depth)
356 : scale(scale), broadness (broadness), depth(depth) {}
357
358 inline const T f(const T x) const { return DIR * ((1.0f / mars::SQ(x / scale)) * broadness + depth) * scale; }
359 inline const T f(const T x, const T y) const { return DIR * ((1.0f / (mars::SQ(x / scale) + mars::SQ(y / scale))) * broadness + depth) * scale; }
360 inline const T invf(const T y) const { return sqrt(broadness / ((DIR * y / scale) - depth)) * scale; }
361 inline const T df(const T x) const { return DIR * (-2.0f / (mars::CUBE(x) / mars::SQ(scale))) * broadness * scale; }
362 inline const T f2(const T x) const { return DIR * (-1.0f / (x / scale) * broadness + (x / scale) * depth) * scale; }
363 };
364
365 template <typename T, int DIR>
366 class CosineFn
367 {
368 public:
369 T scale, period, amplitude, depth;
370
371 inline CosineFn (const T scale, const T period, const T amplitude, const T depth)
372 : scale(scale), period (period), amplitude(amplitude), depth(depth) {}
373
374 inline const T f(const T x) const { return DIR * (cos(x * period / scale) / amplitude + depth) * scale; }
375 inline const T invf(const T y) const { return acos(((DIR * y / scale) - depth) * amplitude) * scale / period; }
376 inline const T df(const T x) const { return DIR * (-sin(x * period / scale) * period / scale) / amplitude * scale; }
377 inline const T f2(const T x) const { return DIR * (sin(x * period / scale) / amplitude / (period / scale) + (x / scale)) * scale; }
378 };
379
380 namespace FnVar
381 {
382 enum Type
383 {
384 Normal = 0,
385 Derivative = 1,
386 Inverse = 2,
387 Integral = 3
388 };
389
390 enum Flag
391 {
392 Flag_Normal = 1 << Normal,
393 Flag_Derivative = 1 << Derivative,
394 Flag_Inverse = 1 << Inverse,
395 Flag_Integral = 1 << Integral
396 };
397 };
398
399 template <typename T, typename Fn, FnVar::Flag FLAGS = FnVar::Flag_Normal>
400 class CacheFnProxy
401 {
402 public:
403 CacheFnProxy (Fn fn, const T fMin, const T fMax, const unsigned short nResolution)
404 : _fn(fn), _fMin(fMin), _fMax(fMax),
405 _fStep ((fMax - fMin) / static_cast< T > (nResolution)),
406 _nResolution(nResolution + 1)
407 {
408 assert(fMax > fMin);
409 assert(_fStep > 0);
410
411 _ffTable[FnVar::Normal] =
412 _ffTable[FnVar::Derivative] =
413 _ffTable[FnVar::Inverse] =
414 _ffTable[FnVar::Integral] = NULL;
415
416 processFlag< FLAGS & FnVar::Flag_Normal > ();
417 processFlag< FLAGS & FnVar::Flag_Derivative > ();
418 processFlag< FLAGS & FnVar::Flag_Inverse > ();
419 processFlag< FLAGS & FnVar::Flag_Integral > ();
420 }
421
422 inline const unsigned short getResolution () const { return _nResolution - 1; }
423
424 inline const T f(const T x) const { return interpolate<FnVar::Normal>(x); }
425 inline const T invf(const T y) const { return interpolate<FnVar::Inverse>(y); }
426 inline const T df(const T x) const { return interpolate<FnVar::Derivative>(x); }
427 inline const T f2(const T x) const { return interpolate<FnVar::Integral>(x); }
428
429 ~CacheFnProxy()
430 {
431 using namespace FnVar;
432
433 delete [] _ffTable[FnVar::Normal];
434 delete [] _ffTable[FnVar::Derivative];
435 delete [] _ffTable[FnVar::Inverse];
436 delete [] _ffTable[FnVar::Integral];
437 }
438
439 private:
440 Fn _fn;
441
442 T _fMin, _fMax, _fStep;
443 unsigned short _nResolution;
444
445 T * _ffTable[4];
446
447 template< int I >
448 inline const T interpolate(const T x) const
449 {
450 assert(_ffTable[I] != NULL);
451
452 const T
453 mx = (x - _fMin) / _fStep,
454 x0 = floor(mx);
455 const unsigned int
456 nx0 = static_cast< unsigned int > (mx);
457
458 assert(nx0 >= 0 && nx0 < _nResolution);
459
460 return
461 _ffTable[I][nx0] + (_ffTable[I][nx0 + 1] - _ffTable[I][nx0]) * (mx - x0);
462 }
463
464 template< int FLAG > void processFlag () { processFlag(quantity<FLAG>()); }
465 void processFlag (quantity<0>) {}
466
467 void processFlag (quantity<FnVar::Flag_Normal>)
468 { computeTable< &Fn::f, FnVar::Normal > (); }
469 void processFlag (quantity<FnVar::Flag_Derivative>)
470 { computeTable< &Fn::df, FnVar::Derivative > (); }
471 void processFlag (quantity<FnVar::Flag_Integral>)
472 { computeTable< &Fn::f2, FnVar::Integral > (); }
473 void processFlag (quantity<FnVar::Flag_Inverse>)
474 { computeTable< &Fn::invf, FnVar::Inverse > (); }
475
476 template< const T (Fn::* FnFn)(const T) const, int I >
477 void computeTable ()
478 {
479 if (FLAGS & (1 << I))
480 {
481 _ffTable[I] = new T[_nResolution];
482 for (unsigned int i = 0; i < _nResolution; ++i)
483 _ffTable[I][i] = (_fn.*FnFn)(_fMin + _fStep * static_cast< T > (i));
484 } else
485 _ffTable[I] = NULL;
486 }
487 };
488
489 template< typename T >
490 class ParabolicScale : protected ParabolicFn< float >
491 {
492 private:
493 const size_t _nCount;
494 const float _fScale;
495 float * _fIndex;
496
497 public:
498 ParabolicScale(const size_t nCount, const float fFactor = 0.5f, const float fScale = 1.0f)
499 : _nCount (nCount), _fScale(fScale), _fIndex(new float[nCount]), ParabolicFn(static_cast< float > (nCount - 1), 1.0f / fFactor, 0.0f)
500 {
501 const float fTotalArea = f2(static_cast< float > (nCount));
502 for (size_t c = 0; c < nCount; ++c)
503 _fIndex[c] = f2(static_cast< float > (c)) / fTotalArea * fScale;
504 }
505 ~ParabolicScale()
506 { delete [] _fIndex; }
507
508 inline T operator () (const float f) const
509 {
510 assert (f >= 0.0f && f <= _fScale);
511 for (size_t c = 0; c < _nCount; ++c)
512 if (f < _fIndex[c])
513 return static_cast< T > (c);
514
515 return static_cast< T > (_nCount - 1);
516 }
517 inline float operator [] (const T enps) const
518 { return _fIndex[enps]; }
519 };
520
521 template< typename T, T COUNT >
522 class StaticParabolicScale : public ParabolicScale< T >
523 {
524 public:
525 StaticParabolicScale (const float fFactor = 0.5f, const float fScale = 1.0f)
526 : ParabolicScale<T>(COUNT, fFactor, fScale) {}
527 };
528}
Note: See TracBrowser for help on using the repository browser.