src/utils/cstd.hh at main · doomsdayone/openmsx

◆ main openmsx / src / utils / cstd.hh

6.0 KiB · C++ · last changed chore(): initial commit 1 month ago source blame history raw

#ifndef CSTD_HH
#define CSTD_HH

#include "Math.hh"
#include "narrow.hh"
#include "xrange.hh"
#include <cassert>
#include <cmath>
#include <cstddef>
#include <functional>
#include <string_view>
#include <utility>

namespace cstd {

template<typename T>
[[nodiscard]] constexpr T abs(T t)
{
    return (t >= 0) ? t : -t;
}

// Reimplementation of various mathematical functions. You must specify an
// iteration count, this controls how accurate the result will be.
#if (defined(__GNUC__) && !defined(__clang__))

// Gcc has constexpr versions of most mathematical functions (this is a
// non-standard extension). Prefer those over our implementations.
template<int>      [[nodiscard]] constexpr double sin  (double x) { return std::sin  (x); }
template<int>      [[nodiscard]] constexpr double cos  (double x) { return std::cos  (x); }
template<int, int> [[nodiscard]] constexpr double log  (double x) { return std::log  (x); }
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?)
template<int, int> [[nodiscard]] constexpr double log10(double x) { return std::log10(x); }
template<int>      [[nodiscard]] constexpr double exp  (double x) { return std::exp  (x); }
template<int>      [[nodiscard]] constexpr double exp2 (double x) { return    ::exp2 (x); } // see log2, but apparently no need to use exp(log(2) * x) here?!
template<int, int> [[nodiscard]] constexpr double pow(double x, double y) { return std::pow(x, y); }
[[nodiscard]] constexpr double round(double x) { return ::round(x); } // should be std::round(), see above
[[nodiscard]] constexpr float  round(float  x) { return ::roundf(x); }
[[nodiscard]] constexpr double sqrt (double x) { return ::sqrt (x); }

#else

[[nodiscard]] constexpr double upow(double x, unsigned u)
{
	double y = 1.0;
	while (u) {
		if (u & 1) y *= x;
		x *= x;
		u >>= 1;
	}
	return y;
}

[[nodiscard]] constexpr double ipow(double x, int i)
{
	return (i >= 0) ? upow(x, i) : upow(x, -i);
}

template<int ITERATIONS>
[[nodiscard]] constexpr double exp(double x)
{
	// Split x into integral and fractional part:
	//   exp(x) = exp(i + f) = exp(i) * exp(f)
	//   with: i an int     (undefined if out of range)
	//         -1 < f < 1
	auto i = int(x);
	auto f = x - i;

	// Approximate exp(f) with Taylor series.
	double y = 1.0;
	double t = f;
	double n = 1.0;
	for (auto k : xrange(ITERATIONS)) {
		y += t / n;
		t *= f;
		n *= k + 2;
	}

	// Approximate exp(i) by squaring.
	int p = (i >= 0) ? i : -i; // abs(i);
	double s = upow(Math::e, p);

	// Combine the results.
	if (i >= 0) {
		return y * s;
	} else {
		return y / s;
	}
}

[[nodiscard]] constexpr double simple_fmod(double x, double y)
{
	assert(y > 0.0);
	return x - int(x / y) * y;
}

template<int ITERATIONS>
[[nodiscard]] constexpr double sin_iter(double x)
{
	double x2 = x * x;
	double y = 0.0;
	double t = x;
	double n = 1.0;
	for (int k = 1; k < (1 + 4 * ITERATIONS); /**/) {
		y += t / n;
		t *= x2;
		n *= ++k;
		n *= ++k;

		y -= t / n;
		t *= x2;
		n *= ++k;
		n *= ++k;
	}
	return y;
}

template<int ITERATIONS>
[[nodiscard]] constexpr double cos_iter(double x)
{
	double x2 = x * x;
	double y = 1.0;
	double t = x2;
	double n = 2.0;
	for (int k = 2; k < (2 + 4 * ITERATIONS); /**/) {
		y -= t / n;
		t *= x2;
		n *= ++k;
		n *= ++k;

		y += t / n;
		t *= x2;
		n *= ++k;
		n *= ++k;
	}
	return y;
}

template<int ITERATIONS>
[[nodiscard]] constexpr double sin(double x)
{
	double sign = 1.0;

	// reduce to [0, +inf)
	if (x < 0.0) {
		sign = -1.0;
		x = -x;
	}

	// reduce to [0, 2pi)
	x = simple_fmod(x, 2 * Math::pi);

	// reduce to [0, pi]
	if (x > Math::pi) {
		sign = -sign;
		x -= Math::pi;
	}

	// reduce to [0, pi/2]
	if (x > Math::pi / 2) {
		x = Math::pi - x;
	}

	// reduce to [0, pi/4]
	if (x > Math::pi / 4) {
		x = Math::pi / 2 - x;
		return sign * cos_iter<ITERATIONS>(x);
	} else {
		return sign * sin_iter<ITERATIONS>(x);
	}
}

template<int ITERATIONS>
[[nodiscard]] constexpr double cos(double x)
{
	double sign = 1.0;

	// reduce to [0, +inf)
	if (x < 0.0) {
		x = -x;
	}

	// reduce to [0, 2pi)
	x = simple_fmod(x, 2 * Math::pi);

	// reduce to [0, pi]
	if (x > Math::pi) {
		x = 2.0 * Math::pi - x;
	}

	// reduce to [0, pi/2]
	if (x > Math::pi / 2) {
		sign = -sign;
		x = Math::pi - x;
	}

	// reduce to [0, pi/4]
	if (x > Math::pi / 4) {
		x = Math::pi / 2 - x;
		return sign * sin_iter<ITERATIONS>(x);
	} else {
		return sign * cos_iter<ITERATIONS>(x);
	}
}


// https://en.wikipedia.org/wiki/Natural_logarithm#High_precision
template<int E_ITERATIONS, int L_ITERATIONS>
[[nodiscard]] constexpr double log(double x)
{
	int a = 0;
	while (x <= 0.25) {
		x *= Math::e;
		++a;
	}
	double y = 0.0;
	repeat(L_ITERATIONS, [&] {
		auto ey = cstd::exp<E_ITERATIONS>(y);
		y = y + 2.0 * (x - ey) / (x + ey);
	});
	return y - a;
}

template<int E_ITERATIONS, int L_ITERATIONS>
[[nodiscard]] constexpr double log2(double x)
{
	return cstd::log<E_ITERATIONS, L_ITERATIONS>(x) / Math::ln2;
}

template<int E_ITERATIONS, int L_ITERATIONS>
[[nodiscard]] constexpr double log10(double x)
{
	return cstd::log<E_ITERATIONS, L_ITERATIONS>(x) / Math::ln10;
}

template<int E_ITERATIONS, int L_ITERATIONS>
[[nodiscard]] constexpr double pow(double x, double y)
{
	return cstd::exp<E_ITERATIONS>(cstd::log<E_ITERATIONS, L_ITERATIONS>(x) * y);
}

template<int ITERATIONS>
[[nodiscard]] constexpr double exp2(double x)
{
	return cstd::exp<ITERATIONS>(Math::ln2 * x);
}

[[nodiscard]] constexpr double round(double x)
{
	return narrow_cast<double>(
		(x >= 0) ?  narrow_cast<int>( x + 0.5)
		         : -narrow_cast<int>(-x + 0.5));
}

[[nodiscard]] constexpr float round(float x)
{
	return narrow_cast<float>(
		(x >= 0) ?  narrow_cast<int>( x + 0.5f)
		         : -narrow_cast<int>(-x + 0.5f));
}

[[nodiscard]] constexpr double sqrt(double x)
{
    assert(x >= 0.0);
    double curr = x;
    double prev = 0.0;
    while (curr != prev) {
        prev = curr;
        curr = 0.5 * (curr + x / curr);
    }
    return curr;
}

#endif

} // namespace cstd

#endif

`/`	Search
`g` `c`	Code
`g` `h`	History
`g` `b`	Bookmarks
`g` `t`	Trust
`y`	Permalink to this state (files and trees)
`?`	This sheet

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
14	namespace cstd {
15
16	template<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.
28	template<int> [[nodiscard]] constexpr double sin (double x) { return std::sin (x); }
29	template<int> [[nodiscard]] constexpr double cos (double x) { return std::cos (x); }
30	template<int, int> [[nodiscard]] constexpr double log (double x) { return std::log (x); }
31	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?)
32	template<int, int> [[nodiscard]] constexpr double log10(double x) { return std::log10(x); }
33	template<int> [[nodiscard]] constexpr double exp (double x) { return std::exp (x); }
34	template<int> [[nodiscard]] constexpr double exp2 (double x) { return ::exp2 (x); } // see log2, but apparently no need to use exp(log(2) * x) here?!
35	template<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 ::roundf(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
58	template<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	auto i = int(x);
66	auto 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
96	template<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
117	template<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
138	template<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
172	template<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
207	template<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
223	template<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
229	template<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
235	template<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
241	template<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

Keyboard