openMSX
Math.hh
Go to the documentation of this file.
1 #ifndef MATH_HH
2 #define MATH_HH
3 
4 #include "likely.hh"
5 #include <cmath>
6 #include <cstdint>
7 
8 #ifdef _MSC_VER
9 #include <intrin.h>
10 #pragma intrinsic(_BitScanForward)
11 #endif
12 
13 // These constants are very common extensions, but not guaranteed to be defined
14 // by <cmath> when compiling in a strict standards compliant mode. Also e.g.
15 // visual studio does not provide them.
16 #ifndef M_E
17 #define M_E 2.7182818284590452354
18 #endif
19 #ifndef M_LN2
20 #define M_LN2 0.69314718055994530942 // log_e 2
21 #endif
22 #ifndef M_LN10
23 #define M_LN10 2.30258509299404568402 // log_e 10
24 #endif
25 #ifndef M_PI
26 #define M_PI 3.14159265358979323846
27 #endif
28 
29 namespace Math {
30 
34 constexpr bool isPowerOfTwo(unsigned a)
35 {
36  return (a & (a - 1)) == 0;
37 }
38 
41 unsigned powerOfTwo(unsigned a);
42 
47 template <int LO, int HI>
48 inline int clip(int x)
49 {
50  static_assert(LO <= HI, "invalid clip range");
51  return unsigned(x - LO) <= unsigned(HI - LO) ? x : (x < HI ? LO : HI);
52 }
53 
57 inline int16_t clipIntToShort(int x)
58 {
59  static_assert((-1 >> 1) == -1, "right-shift must preserve sign");
60  return likely(int16_t(x) == x) ? x : (0x7FFF - (x >> 31));
61 }
62 
66 inline uint8_t clipIntToByte(int x)
67 {
68  static_assert((-1 >> 1) == -1, "right-shift must preserve sign");
69  return likely(uint8_t(x) == x) ? x : ~(x >> 31);
70 }
71 
81 inline unsigned gcd(unsigned a, unsigned b)
82 {
83  unsigned k = 0;
84  while (((a & 1) == 0) && ((b & 1) == 0)) {
85  a >>= 1; b >>= 1; ++k;
86  }
87 
88  // either a or b (or both) is odd
89  while ((a & 1) == 0) a >>= 1;
90  while ((b & 1) == 0) b >>= 1;
91 
92  // both a and b odd
93  while (a != b) {
94  if (a >= b) {
95  a -= b;
96  do { a >>= 1; } while ((a & 1) == 0);
97  } else {
98  b -= a;
99  do { b >>= 1; } while ((b & 1) == 0);
100  }
101  }
102  return b << k;
103 }
104 
109 inline unsigned reverseNBits(unsigned x, unsigned bits)
110 {
111  unsigned ret = 0;
112  while (bits--) {
113  ret = (ret << 1) | (x & 1);
114  x >>= 1;
115  }
116  return ret;
117 
118  /* Just for fun I tried the asm version below (the carry-flag trick
119  * cannot be described in plain C). It's correct and generates shorter
120  * code (both less instructions and less bytes). But it doesn't
121  * actually run faster on the machine I tested on, or only a tiny bit
122  * (possibly because of dependency chains and processor stalls???).
123  * However a big disadvantage of this asm version is that when called
124  * with compile-time constant arguments, this version performs exactly
125  * the same, while the version above can be further optimized by the
126  * compiler (constant-propagation, loop unrolling, ...).
127  unsigned ret = 0;
128  if (bits) {
129  asm (
130  "1: shr %[VAL]\n"
131  " adc %[RET],%[RET]\n"
132  " dec %[BITS]\n"
133  " jne 1b\n"
134  : [VAL] "+r" (val)
135  , [BITS] "+r" (bits)
136  , [RET] "+r" (ret)
137  );
138  }
139  return ret;
140  */
141 
142  /* Maarten suggested the following approach with O(lg(N)) time
143  * complexity (the version above is O(N)).
144  * - reverse full (32-bit) word: O(lg(N))
145  * - shift right over 32-N bits: O(1)
146  * Note: In some lower end CPU the shift-over-N-bits instruction itself
147  * is O(N), in that case this whole algorithm is O(N)
148  * Note2: Instead of '32' it's also possible to use a lower power of 2,
149  * as long as it's bigger than or equal to N.
150  * This algorithm may or may not be faster than the version above, I
151  * didn't try it yet. Also because this routine is _NOT_ performance
152  * critical _AT_ALL_ currently.
153  */
154 }
155 
159 inline uint8_t reverseByte(uint8_t a)
160 {
161  // Classical implementation (can be extended to 16 and 32 bits)
162  // a = ((a & 0xF0) >> 4) | ((a & 0x0F) << 4);
163  // a = ((a & 0xCC) >> 2) | ((a & 0x33) << 2);
164  // a = ((a & 0xAA) >> 1) | ((a & 0x55) << 1);
165  // return a;
166 
167  // The versions below are specific to reverse a single byte (can't
168  // easily be extended to wider types). Found these tricks on:
169  // http://graphics.stanford.edu/~seander/bithacks.html
170 #ifdef __x86_64
171  // on 64-bit systems this is slightly faster
172  return (((a * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL) >> 32;
173 #else
174  // on 32-bit systems this is faster
175  return (((a * 0x0802 & 0x22110) | (a * 0x8020 & 0x88440)) * 0x10101) >> 16;
176 #endif
177 }
178 
185 template<typename T> inline T floodRight(T x)
186 {
187  x |= x >> 1;
188  x |= x >> 2;
189  x |= x >> 4;
190  x |= x >> ((sizeof(x) >= 2) ? 8 : 0); // Written in a weird way to
191  x |= x >> ((sizeof(x) >= 4) ? 16 : 0); // suppress compiler warnings.
192  x |= x >> ((sizeof(x) >= 8) ? 32 : 0); // Generates equally efficient
193  return x; // code.
194 }
195 
199 inline unsigned countLeadingZeros(unsigned x)
200 {
201 #ifdef __GNUC__
202  // actually this only exists starting from gcc-3.4.x
203  return __builtin_clz(x); // undefined when x==0
204 #else
205  // gives incorrect result for x==0, but that doesn't matter here
206  unsigned lz = 0;
207  if (x <= 0x0000ffff) { lz += 16; x <<= 16; }
208  if (x <= 0x00ffffff) { lz += 8; x <<= 8; }
209  if (x <= 0x0fffffff) { lz += 4; x <<= 4; }
210  lz += (0x55ac >> ((x >> 27) & 0x1e)) & 0x3;
211  return lz;
212 #endif
213 }
214 
219 inline unsigned findFirstSet(unsigned x)
220 {
221 #if defined(__GNUC__)
222  return __builtin_ffs(x);
223 #elif defined(_MSC_VER)
224  unsigned long index;
225  return _BitScanForward(&index, x) ? index + 1 : 0;
226 #else
227  if (x == 0) return 0;
228  int pos = 0;
229  if ((x & 0xffff) == 0) { pos += 16; x >>= 16; }
230  if ((x & 0x00ff) == 0) { pos += 8; x >>= 8; }
231  if ((x & 0x000f) == 0) { pos += 4; x >>= 4; }
232  if ((x & 0x0003) == 0) { pos += 2; x >>= 2; }
233  if ((x & 0x0001) == 0) { pos += 1; }
234  return pos + 1;
235 #endif
236 }
237 
238 } // namespace Math
239 
240 #endif // MATH_HH
int clip(int x)
Clips x to the range [LO,HI].
Definition: Math.hh:48
unsigned reverseNBits(unsigned x, unsigned bits)
Reverse the lower N bits of a given value.
Definition: Math.hh:109
int16_t clipIntToShort(int x)
Clip x to range [-32768,32767].
Definition: Math.hh:57
unsigned gcd(unsigned a, unsigned b)
Calculate greatest common divider of two strictly positive integers.
Definition: Math.hh:81
unsigned findFirstSet(unsigned x)
Find the least significant bit that is set.
Definition: Math.hh:219
unsigned powerOfTwo(unsigned a)
Returns the smallest number that is both >=a and a power of two.
Definition: Math.cc:5
Definition: Math.cc:3
#define likely(x)
Definition: likely.hh:14
uint8_t clipIntToByte(int x)
Clip x to range [0,255].
Definition: Math.hh:66
uint8_t reverseByte(uint8_t a)
Reverse the bits in a byte.
Definition: Math.hh:159
constexpr bool isPowerOfTwo(unsigned a)
Is the given number an integer power of 2? Not correct for zero (according to this test 0 is a power ...
Definition: Math.hh:34
unsigned countLeadingZeros(unsigned x)
Count the number of leading zero-bits in the given word.
Definition: Math.hh:199
T floodRight(T x)
Returns the smallest number of the form 2^n-1 that is greater or equal to the given number...
Definition: Math.hh:185