1 /*
2 This file is part of Warzone 2100.
3 Copyright (C) 1999-2004 Eidos Interactive
4 Copyright (C) 2005-2020 Warzone 2100 Project
5
6 Warzone 2100 is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 Warzone 2100 is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with Warzone 2100; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20 /*!
21 * \file
22 * \brief Routines to provide simple math helper functions
23 */
24
25 #ifndef MATH_EXT_H
26 #define MATH_EXT_H
27
28 #include "wzglobal.h"
29 #include <cmath>
30 #include <complex>
31 #include <cstdlib>
32 #include <math.h>
33
34 // Also PERCENT(int,int); // returns a int value 0->100 of the percentage of the first param over the second
35 #define PERCENT(a,b) (((a)*100)/(b))
36 #define PERNUM(range,a,b) (((a)*range)/(b))
37
38 #ifndef M_PI
39 # define M_PI 3.14159265358979323846
40 #endif
41
42 #if (!defined(WZ_C99) && !defined(__cplusplus) && !defined(WZ_CC_GNU)) || (defined _MSC_VER)
43 #if (_MSC_VER <= 1600)
44 # include <float.h>
45
46
roundf(float x)47 static inline int roundf(float x)
48 {
49 // Ensure that float truncation results in a proper rounding
50 if (x < 0.0f)
51 {
52 return x - 0.5f;
53 }
54 else
55 {
56 return x + 0.5f;
57 }
58 }
59
60
61 /**
62 * nearbyint(3) implementation because that function is only available on C99
63 * compatible C libraries.
64 *
65 * This function rounds its argument to an integer value in floating point
66 * format, using the current rounding direction and without raising an
67 * @c inexact exception.
68 *
69 * @return The rounded integer value. If @c x is integral or infinite, @c x
70 * itself is returned.
71 */
nearbyint(double x)72 static inline double nearbyint(double x)
73 {
74 if (ceil(x + 0.5) == floor(x + 0.5))
75 {
76 if ((int)ceil(x) % 2 == 0)
77 {
78 return ceil(x);
79 }
80 else
81 {
82 return floor(x);
83 }
84 }
85 else
86 {
87 return floor(x + 0.5);
88 }
89 }
90 // this is already included for MSVC 2010
91 #if _MSC_VER < 1600
hypotf(float x,float y)92 static inline WZ_DECL_CONST float hypotf(float x, float y)
93 {
94 return sqrtf(x * x + y * y);
95 }
96 #endif
97 #endif
98 #endif
99
100
101 /*!
102 * Clips x to boundaries
103 * \param x Value to clip
104 * \param min Lower bound
105 * \param max Upper bound
106 */
107 template<typename T>
clip(T x,T min,T max)108 static inline WZ_DECL_CONST T clip(T x, T min, T max)
109 {
110 // std::min and std::max not constexpr until C++14.
111 return x < min ? min : x > max ? max : x;
112 }
113
114
115 /*!
116 * Clips x to boundaries
117 * \param x Value to clip
118 * \param min Lower bound
119 * \param max Upper bound
120 */
clipf(float x,float min,float max)121 static inline WZ_DECL_CONST float clipf(float x, float min, float max)
122 {
123 return x < min ? min : x > max ? max : x;
124 }
125
126 /// Finds change in y and y' after time dt, according to the differential equation y'' = -ay - fy'
solveDifferential2ndOrder(float * y_,float * dydt_,double acceleration,double friction,double dt)127 static inline void solveDifferential2ndOrder(float *y_, float *dydt_, double acceleration, double friction, double dt)
128 {
129 double y = *y_, dydt = *dydt_;
130 // Solution is y = g_1 exp(h_1 t) + g_2 exp(h_2 t), where h_i are the solutions to h^2 + f h + a = 0, which are -f/2 +- sqrt(f^2/4 - a).
131 // At t = 0, g_1 = (y' - h_2 y) / (h_1 - h_2) and g_2 = (y' - h_1 y) / (h_2 - h_1).
132
133 std::complex<double> d = friction * friction / 4 - acceleration;
134
135 std::complex<double> sqd = std::sqrt(d);
136 std::complex<double> h1 = -friction / 2 + sqd;
137 std::complex<double> h2 = -friction / 2 - sqd;
138 std::complex<double> e1 = std::exp(h1 * dt);
139 std::complex<double> e2 = std::exp(h2 * dt);
140 std::complex<double> g1 = (dydt - h2 * y) / (h1 - h2);
141 std::complex<double> g2 = (dydt - h1 * y) / (h2 - h1);
142 *y_ = (g1 * e1 + g2 * e2).real(); // .imag() should be 0.
143 *dydt_ = (g1 * h1 * e1 + g2 * h2 * e2).real(); // .imag() should be 0.
144 }
145
146 // Windows unfortunately appears to do this, so do this too for compatibility...
147 using std::abs;
148 using std::sqrt;
149 using std::pow;
150
151 #endif // MATH_EXT_H
152