openMSX
cstd.hh
Go to the documentation of this file.
1 #ifndef CSTD_HH
2 #define CSTD_HH
3 
4 #include <cassert>
5 #include <cmath>
6 #include <cstddef>
7 #include <functional>
8 #include <string_view>
9 #include <utility>
10 
11 namespace cstd {
12 
13 // Inspired by this post:
14 // http://tristanbrindle.com/posts/a-more-useful-compile-time-quicksort
15 
16 // Constexpr reimplementations of standard algorithms or data-structures.
17 //
18 // Everything implemented in 'cstd' will very likely become part of standard
19 // C++ in the future. So then we can remove our (re)implementation and change
20 // the callers from 'cstd::xxx()' to 'std::xxx()'.
21 
22 //
23 // Various constexpr reimplementation of STL algorithms.
24 //
25 
26 template<typename Iter1, typename Iter2>
27 constexpr void iter_swap(Iter1 a, Iter2 b)
28 {
29  auto temp = std::move(*a);
30  *a = std::move(*b);
31  *b = std::move(temp);
32 }
33 
34 template<typename InputIt, typename UnaryPredicate>
35 [[nodiscard]] constexpr InputIt find_if_not(InputIt first, InputIt last, UnaryPredicate p)
36 {
37  for (; first != last; ++first) {
38  if (!p(*first)) {
39  return first;
40  }
41  }
42  return last;
43 }
44 
45 template<typename ForwardIt, typename UnaryPredicate>
46 [[nodiscard]] constexpr ForwardIt partition(ForwardIt first, ForwardIt last, UnaryPredicate p)
47 {
48  first = cstd::find_if_not(first, last, p);
49  if (first == last) return first;
50 
51  for (ForwardIt i = first + 1; i != last; ++i) {
52  if (p(*i)) {
53  cstd::iter_swap(i, first);
54  ++first;
55  }
56  }
57  return first;
58 }
59 
60 // Textbook implementation of quick sort. Not optimized, but the intention is
61 // that it only gets executed during compile-time.
62 template<typename RAIt, typename Compare = std::less<>>
63 constexpr void sort(RAIt first, RAIt last, Compare cmp = Compare{})
64 {
65  auto const N = last - first;
66  if (N <= 1) return;
67  auto const pivot = *(first + N / 2);
68  auto const middle1 = cstd::partition(
69  first, last, [&](const auto& elem) { return cmp(elem, pivot); });
70  auto const middle2 = cstd::partition(
71  middle1, last, [&](const auto& elem) { return !cmp(pivot, elem); });
72  cstd::sort(first, middle1, cmp);
73  cstd::sort(middle2, last, cmp);
74 }
75 
76 template<typename RandomAccessRange, typename Compare = std::less<>>
77 constexpr void sort(RandomAccessRange&& range, Compare cmp = Compare{})
78 {
79  cstd::sort(std::begin(range), std::end(range), cmp);
80 }
81 
82 template<typename InputIt1, typename InputIt2>
83 [[nodiscard]] constexpr bool lexicographical_compare(InputIt1 first1, InputIt1 last1,
84  InputIt2 first2, InputIt2 last2)
85 {
86  for (; (first1 != last1) && (first2 != last2); ++first1, ++first2) {
87  if (*first1 < *first2) return true;
88  if (*first2 < *first1) return false;
89  }
90  return (first1 == last1) && (first2 != last2);
91 }
92 
93 template<typename T>
94 constexpr T abs(T t)
95 {
96  return (t >= 0) ? t : -t;
97 }
98 
99 // Reimplementation of various mathematical functions. You must specify an
100 // iteration count, this controls how accurate the result will be.
101 #if (defined(__GNUC__) && !defined(__clang__))
102 
103 // Gcc has constexpr versions of most mathematical functions (this is a
104 // non-standard extension). Prefer those over our implementations.
105 template<int> [[nodiscard]] constexpr double sin (double x) { return std::sin (x); }
106 template<int> [[nodiscard]] constexpr double cos (double x) { return std::cos (x); }
107 template<int, int> [[nodiscard]] constexpr double log (double x) { return std::log (x); }
108 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?)
109 template<int, int> [[nodiscard]] constexpr double log10(double x) { return std::log10(x); }
110 template<int> [[nodiscard]] constexpr double exp (double x) { return std::exp (x); }
111 template<int> [[nodiscard]] constexpr double exp2 (double x) { return ::exp2 (x); } // see log2, but apparently no need to use exp(log(2) * x) here?!
112 template<int, int> [[nodiscard]] constexpr double pow(double x, double y) { return std::pow(x, y); }
113 [[nodiscard]] inline constexpr double round(double x) { return ::round(x); } // should be std::round(), see above
114 [[nodiscard]] inline constexpr float round(float x) { return ::round(x); }
115 [[nodiscard]] inline constexpr double sqrt (double x) { return ::sqrt (x); }
116 
117 #else
118 
119 [[nodiscard]] constexpr double upow(double x, unsigned u)
120 {
121  double y = 1.0;
122  while (u) {
123  if (u & 1) y *= x;
124  x *= x;
125  u >>= 1;
126  }
127  return y;
128 }
129 
130 [[nodiscard]] constexpr double ipow(double x, int i)
131 {
132  return (i >= 0) ? upow(x, i) : upow(x, -i);
133 }
134 
135 template<int ITERATIONS>
136 [[nodiscard]] constexpr double exp(double x)
137 {
138  // Split x into integral and fractional part:
139  // exp(x) = exp(i + f) = exp(i) * exp(f)
140  // with: i an int (undefined if out of range)
141  // -1 < f < 1
142  int i = int(x);
143  double f = x - i;
144 
145  // Approximate exp(f) with Taylor series.
146  double y = 1.0;
147  double t = f;
148  double n = 1.0;
149  for (int k = 2; k < (2 + ITERATIONS); ++k) {
150  y += t / n;
151  t *= f;
152  n *= k;
153  }
154 
155  // Approximate exp(i) by squaring.
156  int p = (i >= 0) ? i : -i; // abs(i);
157  double s = upow(M_E, p);
158 
159  // Combine the results.
160  if (i >= 0) {
161  return y * s;
162  } else {
163  return y / s;
164  }
165 }
166 
167 [[nodiscard]] constexpr double simple_fmod(double x, double y)
168 {
169  assert(y > 0.0);
170  return x - int(x / y) * y;
171 }
172 
173 template<int ITERATIONS>
174 [[nodiscard]] constexpr double sin_iter(double x)
175 {
176  double x2 = x * x;
177  double y = 0.0;
178  double t = x;
179  double n = 1.0;
180  for (int k = 1; k < (1 + 4 * ITERATIONS); ) {
181  y += t / n;
182  t *= x2;
183  n *= ++k;
184  n *= ++k;
185 
186  y -= t / n;
187  t *= x2;
188  n *= ++k;
189  n *= ++k;
190  }
191  return y;
192 }
193 
194 template<int ITERATIONS>
195 [[nodiscard]] constexpr double cos_iter(double x)
196 {
197  double x2 = x * x;
198  double y = 1.0;
199  double t = x2;
200  double n = 2.0;
201  for (int k = 2; k < (2 + 4 * ITERATIONS); ) {
202  y -= t / n;
203  t *= x2;
204  n *= ++k;
205  n *= ++k;
206 
207  y += t / n;
208  t *= x2;
209  n *= ++k;
210  n *= ++k;
211  }
212  return y;
213 }
214 
215 template<int ITERATIONS>
216 [[nodiscard]] constexpr double sin(double x)
217 {
218  double sign = 1.0;
219 
220  // reduce to [0, +inf)
221  if (x < 0.0) {
222  sign = -1.0;
223  x = -x;
224  }
225 
226  // reduce to [0, 2pi)
227  x = simple_fmod(x, 2 * M_PI);
228 
229  // reduce to [0, pi]
230  if (x > M_PI) {
231  sign = -sign;
232  x -= M_PI;
233  }
234 
235  // reduce to [0, pi/2]
236  if (x > M_PI/2) {
237  x = M_PI - x;
238  }
239 
240  // reduce to [0, pi/4]
241  if (x > M_PI/4) {
242  x = M_PI/2 - x;
243  return sign * cos_iter<ITERATIONS>(x);
244  } else {
245  return sign * sin_iter<ITERATIONS>(x);
246  }
247 }
248 
249 template<int ITERATIONS>
250 [[nodiscard]] constexpr double cos(double x)
251 {
252  double sign = 1.0;
253 
254  // reduce to [0, +inf)
255  if (x < 0.0) {
256  x = -x;
257  }
258 
259  // reduce to [0, 2pi)
260  x = simple_fmod(x, 2 * M_PI);
261 
262  // reduce to [0, pi]
263  if (x > M_PI) {
264  x = 2.0 * M_PI - x;
265  }
266 
267  // reduce to [0, pi/2]
268  if (x > M_PI/2) {
269  sign = -sign;
270  x = M_PI - x;
271  }
272 
273  // reduce to [0, pi/4]
274  if (x > M_PI/4) {
275  x = M_PI/2 - x;
276  return sign * sin_iter<ITERATIONS>(x);
277  } else {
278  return sign * cos_iter<ITERATIONS>(x);
279  }
280 }
281 
282 
283 // https://en.wikipedia.org/wiki/Natural_logarithm#High_precision
284 template<int E_ITERATIONS, int L_ITERATIONS>
285 [[nodiscard]] constexpr double log(double x)
286 {
287  int a = 0;
288  while (x <= 0.25) {
289  x *= M_E;
290  ++a;
291  }
292  double y = 0.0;
293  for (int i = 0; i < L_ITERATIONS; ++i) {
294  auto ey = cstd::exp<E_ITERATIONS>(y);
295  y = y + 2.0 * (x - ey) / (x + ey);
296  }
297  return y - a;
298 }
299 
300 template<int E_ITERATIONS, int L_ITERATIONS>
301 [[nodiscard]] constexpr double log2(double x)
302 {
303  return cstd::log<E_ITERATIONS, L_ITERATIONS>(x) / M_LN2;
304 }
305 
306 template<int E_ITERATIONS, int L_ITERATIONS>
307 [[nodiscard]] constexpr double log10(double x)
308 {
309  return cstd::log<E_ITERATIONS, L_ITERATIONS>(x) / M_LN10;
310 }
311 
312 template<int E_ITERATIONS, int L_ITERATIONS>
313 [[nodiscard]] constexpr double pow(double x, double y)
314 {
315  return cstd::exp<E_ITERATIONS>(cstd::log<E_ITERATIONS, L_ITERATIONS>(x) * y);
316 }
317 
318 template<int ITERATIONS>
319 [[nodiscard]] constexpr double exp2(double x)
320 {
321  return cstd::exp<ITERATIONS>(M_LN2 * x);
322 }
323 
324 [[nodiscard]] constexpr double round(double x)
325 {
326  return (x >= 0) ? int( x + 0.5)
327  : -int(-x + 0.5);
328 }
329 
330 [[nodiscard]] constexpr float round(float x)
331 {
332  return (x >= 0) ? int( x + 0.5f)
333  : -int(-x + 0.5f);
334 }
335 
336 [[nodiscard]] constexpr double sqrt(double x)
337 {
338  assert(x >= 0.0);
339  double curr = x;
340  double prev = 0.0;
341  while (curr != prev) {
342  prev = curr;
343  curr = 0.5 * (curr + x / curr);
344  }
345  return curr;
346 }
347 
348 #endif
349 
350 } // namespace cstd
351 
352 #endif
cstd::sin_iter
constexpr double sin_iter(double x)
Definition: cstd.hh:174
cstd::sort
constexpr void sort(RAIt first, RAIt last, Compare cmp=Compare{})
Definition: cstd.hh:63
cstd::round
constexpr double round(double x)
Definition: cstd.hh:324
cstd::ipow
constexpr double ipow(double x, int i)
Definition: cstd.hh:130
cstd::cos_iter
constexpr double cos_iter(double x)
Definition: cstd.hh:195
M_LN10
#define M_LN10
Definition: Math.hh:24
cstd
Definition: cstd.hh:11
M_E
#define M_E
Definition: Math.hh:18
t
TclObject t
Definition: TclObject_test.cc:264
cstd::exp
constexpr double exp(double x)
Definition: cstd.hh:136
cstd::cos
constexpr double cos(double x)
Definition: cstd.hh:250
cstd::simple_fmod
constexpr double simple_fmod(double x, double y)
Definition: cstd.hh:167
M_PI
#define M_PI
Definition: Math.hh:27
cstd::upow
constexpr double upow(double x, unsigned u)
Definition: cstd.hh:119
cstd::partition
constexpr ForwardIt partition(ForwardIt first, ForwardIt last, UnaryPredicate p)
Definition: cstd.hh:46
cstd::round
constexpr float round(float x)
Definition: cstd.hh:330
cstd::sin
constexpr double sin(double x)
Definition: cstd.hh:216
cstd::lexicographical_compare
constexpr bool lexicographical_compare(InputIt1 first1, InputIt1 last1, InputIt2 first2, InputIt2 last2)
Definition: cstd.hh:83
cstd::find_if_not
constexpr InputIt find_if_not(InputIt first, InputIt last, UnaryPredicate p)
Definition: cstd.hh:35
cstd::sqrt
constexpr double sqrt(double x)
Definition: cstd.hh:336
cstd::log
constexpr double log(double x)
Definition: cstd.hh:285
openmsx::N
constexpr unsigned N
Definition: ResampleHQ.cc:225
cstd::abs
constexpr T abs(T t)
Definition: cstd.hh:94
cstd::log10
constexpr double log10(double x)
Definition: cstd.hh:307
openmsx::x
constexpr KeyMatrixPosition x
Keyboard bindings.
Definition: Keyboard.cc:1419
cstd::iter_swap
constexpr void iter_swap(Iter1 a, Iter2 b)
Definition: cstd.hh:27
M_LN2
#define M_LN2
Definition: Math.hh:21
cstd::pow
constexpr double pow(double x, double y)
Definition: cstd.hh:313
cstd::log2
constexpr double log2(double x)
Definition: cstd.hh:301
cstd::exp2
constexpr double exp2(double x)
Definition: cstd.hh:319