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