openMSX
cstd.hh
Go to the documentation of this file.
1 #ifndef CSTD_HH
2 #define CSTD_HH
3 
4 #include "Math.hh"
5 #include "xrange.hh"
6 #include <cassert>
7 #include <cmath>
8 #include <cstddef>
9 #include <functional>
10 #include <string_view>
11 #include <utility>
12 
13 namespace cstd {
14 
15 template<typename T>
16 [[nodiscard]] constexpr T abs(T t)
17 {
18  return (t >= 0) ? t : -t;
19 }
20 
21 // Reimplementation of various mathematical functions. You must specify an
22 // iteration count, this controls how accurate the result will be.
23 #if (defined(__GNUC__) && !defined(__clang__))
24 
25 // Gcc has constexpr versions of most mathematical functions (this is a
26 // non-standard extension). Prefer those over our implementations.
27 template<int> [[nodiscard]] constexpr double sin (double x) { return std::sin (x); }
28 template<int> [[nodiscard]] constexpr double cos (double x) { return std::cos (x); }
29 template<int, int> [[nodiscard]] constexpr double log (double x) { return std::log (x); }
30 template<int, int> [[nodiscard]] constexpr double log2 (double x) { return ::log (x) / ::log(2); } // should be std::log2(x) but this doesn't seem to compile in g++-4.8/g++-4.9 (bug?)
31 template<int, int> [[nodiscard]] constexpr double log10(double x) { return std::log10(x); }
32 template<int> [[nodiscard]] constexpr double exp (double x) { return std::exp (x); }
33 template<int> [[nodiscard]] constexpr double exp2 (double x) { return ::exp2 (x); } // see log2, but apparently no need to use exp(log(2) * x) here?!
34 template<int, int> [[nodiscard]] constexpr double pow(double x, double y) { return std::pow(x, y); }
35 [[nodiscard]] constexpr double round(double x) { return ::round(x); } // should be std::round(), see above
36 [[nodiscard]] constexpr float round(float x) { return ::round(x); }
37 [[nodiscard]] constexpr double sqrt (double x) { return ::sqrt (x); }
38 
39 #else
40 
41 [[nodiscard]] constexpr double upow(double x, unsigned u)
42 {
43  double y = 1.0;
44  while (u) {
45  if (u & 1) y *= x;
46  x *= x;
47  u >>= 1;
48  }
49  return y;
50 }
51 
52 [[nodiscard]] constexpr double ipow(double x, int i)
53 {
54  return (i >= 0) ? upow(x, i) : upow(x, -i);
55 }
56 
57 template<int ITERATIONS>
58 [[nodiscard]] constexpr double exp(double x)
59 {
60  // Split x into integral and fractional part:
61  // exp(x) = exp(i + f) = exp(i) * exp(f)
62  // with: i an int (undefined if out of range)
63  // -1 < f < 1
64  int i = int(x);
65  double f = x - i;
66 
67  // Approximate exp(f) with Taylor series.
68  double y = 1.0;
69  double t = f;
70  double n = 1.0;
71  for (auto k : xrange(ITERATIONS)) {
72  y += t / n;
73  t *= f;
74  n *= k + 2;
75  }
76 
77  // Approximate exp(i) by squaring.
78  int p = (i >= 0) ? i : -i; // abs(i);
79  double s = upow(Math::e, p);
80 
81  // Combine the results.
82  if (i >= 0) {
83  return y * s;
84  } else {
85  return y / s;
86  }
87 }
88 
89 [[nodiscard]] constexpr double simple_fmod(double x, double y)
90 {
91  assert(y > 0.0);
92  return x - int(x / y) * y;
93 }
94 
95 template<int ITERATIONS>
96 [[nodiscard]] constexpr double sin_iter(double x)
97 {
98  double x2 = x * x;
99  double y = 0.0;
100  double t = x;
101  double n = 1.0;
102  for (int k = 1; k < (1 + 4 * ITERATIONS); ) {
103  y += t / n;
104  t *= x2;
105  n *= ++k;
106  n *= ++k;
107 
108  y -= t / n;
109  t *= x2;
110  n *= ++k;
111  n *= ++k;
112  }
113  return y;
114 }
115 
116 template<int ITERATIONS>
117 [[nodiscard]] constexpr double cos_iter(double x)
118 {
119  double x2 = x * x;
120  double y = 1.0;
121  double t = x2;
122  double n = 2.0;
123  for (int k = 2; k < (2 + 4 * ITERATIONS); ) {
124  y -= t / n;
125  t *= x2;
126  n *= ++k;
127  n *= ++k;
128 
129  y += t / n;
130  t *= x2;
131  n *= ++k;
132  n *= ++k;
133  }
134  return y;
135 }
136 
137 template<int ITERATIONS>
138 [[nodiscard]] constexpr double sin(double x)
139 {
140  double sign = 1.0;
141 
142  // reduce to [0, +inf)
143  if (x < 0.0) {
144  sign = -1.0;
145  x = -x;
146  }
147 
148  // reduce to [0, 2pi)
149  x = simple_fmod(x, 2 * Math::pi);
150 
151  // reduce to [0, pi]
152  if (x > Math::pi) {
153  sign = -sign;
154  x -= Math::pi;
155  }
156 
157  // reduce to [0, pi/2]
158  if (x > Math::pi / 2) {
159  x = Math::pi - x;
160  }
161 
162  // reduce to [0, pi/4]
163  if (x > Math::pi / 4) {
164  x = Math::pi / 2 - x;
165  return sign * cos_iter<ITERATIONS>(x);
166  } else {
167  return sign * sin_iter<ITERATIONS>(x);
168  }
169 }
170 
171 template<int ITERATIONS>
172 [[nodiscard]] constexpr double cos(double x)
173 {
174  double sign = 1.0;
175 
176  // reduce to [0, +inf)
177  if (x < 0.0) {
178  x = -x;
179  }
180 
181  // reduce to [0, 2pi)
182  x = simple_fmod(x, 2 * Math::pi);
183 
184  // reduce to [0, pi]
185  if (x > Math::pi) {
186  x = 2.0 * Math::pi - x;
187  }
188 
189  // reduce to [0, pi/2]
190  if (x > Math::pi / 2) {
191  sign = -sign;
192  x = Math::pi - x;
193  }
194 
195  // reduce to [0, pi/4]
196  if (x > Math::pi / 4) {
197  x = Math::pi / 2 - x;
198  return sign * sin_iter<ITERATIONS>(x);
199  } else {
200  return sign * cos_iter<ITERATIONS>(x);
201  }
202 }
203 
204 
205 // https://en.wikipedia.org/wiki/Natural_logarithm#High_precision
206 template<int E_ITERATIONS, int L_ITERATIONS>
207 [[nodiscard]] constexpr double log(double x)
208 {
209  int a = 0;
210  while (x <= 0.25) {
211  x *= Math::e;
212  ++a;
213  }
214  double y = 0.0;
215  repeat(L_ITERATIONS, [&] {
216  auto ey = cstd::exp<E_ITERATIONS>(y);
217  y = y + 2.0 * (x - ey) / (x + ey);
218  });
219  return y - a;
220 }
221 
222 template<int E_ITERATIONS, int L_ITERATIONS>
223 [[nodiscard]] constexpr double log2(double x)
224 {
225  return cstd::log<E_ITERATIONS, L_ITERATIONS>(x) / Math::ln2;
226 }
227 
228 template<int E_ITERATIONS, int L_ITERATIONS>
229 [[nodiscard]] constexpr double log10(double x)
230 {
231  return cstd::log<E_ITERATIONS, L_ITERATIONS>(x) / Math::ln10;
232 }
233 
234 template<int E_ITERATIONS, int L_ITERATIONS>
235 [[nodiscard]] constexpr double pow(double x, double y)
236 {
237  return cstd::exp<E_ITERATIONS>(cstd::log<E_ITERATIONS, L_ITERATIONS>(x) * y);
238 }
239 
240 template<int ITERATIONS>
241 [[nodiscard]] constexpr double exp2(double x)
242 {
243  return cstd::exp<ITERATIONS>(Math::ln2 * x);
244 }
245 
246 [[nodiscard]] constexpr double round(double x)
247 {
248  return (x >= 0) ? int( x + 0.5)
249  : -int(-x + 0.5);
250 }
251 
252 [[nodiscard]] constexpr float round(float x)
253 {
254  return (x >= 0) ? int( x + 0.5f)
255  : -int(-x + 0.5f);
256 }
257 
258 [[nodiscard]] constexpr double sqrt(double x)
259 {
260  assert(x >= 0.0);
261  double curr = x;
262  double prev = 0.0;
263  while (curr != prev) {
264  prev = curr;
265  curr = 0.5 * (curr + x / curr);
266  }
267  return curr;
268 }
269 
270 #endif
271 
272 } // namespace cstd
273 
274 #endif
TclObject t
constexpr double pi
Definition: Math.hh:21
constexpr double ln2
Definition: Math.hh:19
constexpr double ln10
Definition: Math.hh:20
constexpr double e
Definition: Math.hh:18
Definition: cstd.hh:13
constexpr double pow(double x, double y)
Definition: cstd.hh:235
constexpr double log(double x)
Definition: cstd.hh:207
constexpr float round(float x)
Definition: cstd.hh:252
constexpr double sin_iter(double x)
Definition: cstd.hh:96
constexpr double simple_fmod(double x, double y)
Definition: cstd.hh:89
constexpr double exp2(double x)
Definition: cstd.hh:241
constexpr double sqrt(double x)
Definition: cstd.hh:258
constexpr double round(double x)
Definition: cstd.hh:246
constexpr T abs(T t)
Definition: cstd.hh:16
constexpr double exp(double x)
Definition: cstd.hh:58
constexpr double cos_iter(double x)
Definition: cstd.hh:117
constexpr double ipow(double x, int i)
Definition: cstd.hh:52
constexpr double log10(double x)
Definition: cstd.hh:229
constexpr double sin(double x)
Definition: cstd.hh:138
constexpr double cos(double x)
Definition: cstd.hh:172
constexpr double log2(double x)
Definition: cstd.hh:223
constexpr double upow(double x, unsigned u)
Definition: cstd.hh:41
constexpr KeyMatrixPosition x
Keyboard bindings.
Definition: Keyboard.cc:127
constexpr void repeat(T n, Op op)
Repeat the given operation 'op' 'n' times.
Definition: xrange.hh:148
constexpr auto xrange(T e)
Definition: xrange.hh:133