openMSX
Math.hh
Go to the documentation of this file.
1#ifndef MATH_HH
2#define MATH_HH
3
4#include <bit>
5#include <cassert>
6#include <cmath>
7#include <concepts>
8#include <cstdint>
9#include <numbers>
10
11#ifdef _MSC_VER
12#include <intrin.h>
13#pragma intrinsic(_BitScanForward)
14#endif
15
16namespace Math {
17
18constexpr inline double e = std::numbers::e_v <double>;
19constexpr inline double ln2 = std::numbers::ln2_v <double>;
20constexpr inline double ln10 = std::numbers::ln10_v<double>;
21constexpr inline double pi = std::numbers::pi_v <double>;
22
29[[nodiscard]] constexpr auto floodRight(std::unsigned_integral auto x) noexcept
30{
31 x |= x >> 1;
32 x |= x >> 2;
33 x |= x >> 4;
34 x |= x >> ((sizeof(x) >= 2) ? 8 : 0); // Written in a weird way to
35 x |= x >> ((sizeof(x) >= 4) ? 16 : 0); // suppress compiler warnings.
36 x |= x >> ((sizeof(x) >= 8) ? 32 : 0); // Generates equally efficient
37 return x; // code.
38}
39
43[[nodiscard]] inline int16_t clipIntToShort(int x)
44{
45 static_assert((-1 >> 1) == -1, "right-shift must preserve sign");
46 if (int16_t(x) == x) [[likely]] {
47 return x;
48 } else {
49 return 0x7FFF - (x >> 31);
50 }
51}
52
56[[nodiscard]] inline uint8_t clipIntToByte(int x)
57{
58 static_assert((-1 >> 1) == -1, "right-shift must preserve sign");
59 if (uint8_t(x) == x) [[likely]] {
60 return x;
61 } else {
62 return ~(x >> 31);
63 }
64}
65
70[[nodiscard]] constexpr unsigned reverseNBits(unsigned x, unsigned bits)
71{
72 unsigned ret = 0;
73 while (bits--) {
74 ret = (ret << 1) | (x & 1);
75 x >>= 1;
76 }
77 return ret;
78
79 /* Just for fun I tried the asm version below (the carry-flag trick
80 * cannot be described in plain C). It's correct and generates shorter
81 * code (both less instructions and less bytes). But it doesn't
82 * actually run faster on the machine I tested on, or only a tiny bit
83 * (possibly because of dependency chains and processor stalls???).
84 * However a big disadvantage of this asm version is that when called
85 * with compile-time constant arguments, this version performs exactly
86 * the same, while the version above can be further optimized by the
87 * compiler (constant-propagation, loop unrolling, ...).
88 unsigned ret = 0;
89 if (bits) {
90 asm (
91 "1: shr %[VAL]\n"
92 " adc %[RET],%[RET]\n"
93 " dec %[BITS]\n"
94 " jne 1b\n"
95 : [VAL] "+r" (val)
96 , [BITS] "+r" (bits)
97 , [RET] "+r" (ret)
98 );
99 }
100 return ret;
101 */
102
103 /* Maarten suggested the following approach with O(lg(N)) time
104 * complexity (the version above is O(N)).
105 * - reverse full (32-bit) word: O(lg(N))
106 * - shift right over 32-N bits: O(1)
107 * Note: In some lower end CPU the shift-over-N-bits instruction itself
108 * is O(N), in that case this whole algorithm is O(N)
109 * Note2: Instead of '32' it's also possible to use a lower power of 2,
110 * as long as it's bigger than or equal to N.
111 * This algorithm may or may not be faster than the version above, I
112 * didn't try it yet. Also because this routine is _NOT_ performance
113 * critical _AT_ALL_ currently.
114 */
115}
116
120[[nodiscard]] constexpr uint8_t reverseByte(uint8_t a)
121{
122 // Classical implementation (can be extended to 16 and 32 bits)
123 // a = ((a & 0xF0) >> 4) | ((a & 0x0F) << 4);
124 // a = ((a & 0xCC) >> 2) | ((a & 0x33) << 2);
125 // a = ((a & 0xAA) >> 1) | ((a & 0x55) << 1);
126 // return a;
127
128 // The versions below are specific to reverse a single byte (can't
129 // easily be extended to wider types). Found these tricks on:
130 // http://graphics.stanford.edu/~seander/bithacks.html
131#ifdef __x86_64
132 // on 64-bit systems this is slightly faster
133 return (((a * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL) >> 32;
134#else
135 // on 32-bit systems this is faster
136 return (((a * 0x0802 & 0x22110) | (a * 0x8020 & 0x88440)) * 0x10101) >> 16;
137#endif
138}
139
144[[nodiscard]] inline /*constexpr*/ unsigned findFirstSet(uint32_t x)
145{
146 return x ? std::countr_zero(x) + 1 : 0;
147}
148
149// Cubic Hermite Interpolation:
150// Given 4 points: (-1, y[-1]), (0, y[0]), (1, y[1]), (2, y[2])
151// Fit a polynomial: f(x) = a*x^3 + b*x^2 + c*x + d
152// which passes through the given points at x=0 and x=1
153// f(0) = y[0]
154// f(1) = y[1]
155// and which has specific derivatives at x=0 and x=1
156// f'(0) = (y[1] - y[-1]) / 2
157// f'(1) = (y[2] - y[ 0]) / 2
158// Then evaluate this polynomial at the given x-position (x in [0, 1]).
159// For more details see:
160// https://en.wikipedia.org/wiki/Cubic_Hermite_spline
161// https://www.paulinternet.nl/?page=bicubic
162[[nodiscard]] constexpr float cubicHermite(const float* y, float x)
163{
164 assert(0.0f <= x); assert(x <= 1.0f);
165 float a = -0.5f*y[-1] + 1.5f*y[0] - 1.5f*y[1] + 0.5f*y[2];
166 float b = y[-1] - 2.5f*y[0] + 2.0f*y[1] - 0.5f*y[2];
167 float c = -0.5f*y[-1] + 0.5f*y[1];
168 float d = y[0];
169 float x2 = x * x;
170 float x3 = x * x2;
171 return a*x3 + b*x2 + c*x + d;
172}
173
182};
183constexpr QuotientRemainder div_mod_floor(int dividend, int divisor) {
184 int q = dividend / divisor;
185 int r = dividend % divisor;
186 if ((r != 0) && ((r < 0) != (divisor < 0))) {
187 --q;
188 r += divisor;
189 }
190 return {q, r};
191}
192constexpr int div_floor(int dividend, int divisor) {
193 return div_mod_floor(dividend, divisor).quotient;
194}
195constexpr int mod_floor(int dividend, int divisor) {
196 return div_mod_floor(dividend, divisor).remainder;
197}
198
199} // namespace Math
200
201#endif // MATH_HH
Definition: Math.hh:16
constexpr QuotientRemainder div_mod_floor(int dividend, int divisor)
Definition: Math.hh:183
constexpr unsigned reverseNBits(unsigned x, unsigned bits)
Reverse the lower N bits of a given value.
Definition: Math.hh:70
unsigned findFirstSet(uint32_t x)
Find the least significant bit that is set.
Definition: Math.hh:144
constexpr auto floodRight(std::unsigned_integral auto x) noexcept
Returns the smallest number of the form 2^n-1 that is greater or equal to the given number.
Definition: Math.hh:29
uint8_t clipIntToByte(int x)
Clip x to range [0,255].
Definition: Math.hh:56
constexpr double pi
Definition: Math.hh:21
constexpr int div_floor(int dividend, int divisor)
Definition: Math.hh:192
constexpr int mod_floor(int dividend, int divisor)
Definition: Math.hh:195
constexpr double ln2
Definition: Math.hh:19
constexpr double ln10
Definition: Math.hh:20
constexpr uint8_t reverseByte(uint8_t a)
Reverse the bits in a byte.
Definition: Math.hh:120
int16_t clipIntToShort(int x)
Clip x to range [-32768,32767].
Definition: Math.hh:43
constexpr float cubicHermite(const float *y, float x)
Definition: Math.hh:162
constexpr double e
Definition: Math.hh:18
constexpr KeyMatrixPosition x
Keyboard bindings.
Definition: Keyboard.cc:127
Divide one integer by another, rounding towards minus infinity.
Definition: Math.hh:179