mathutil.hpp 4.83 KB
Newer Older
incardon's avatar
incardon committed
1
2
3
#ifndef MATHUTIL_HPP
#define MATHUTIL_HPP

incardon's avatar
incardon committed
4
#ifdef HAVE_LIBQUADMATH
incardon's avatar
incardon committed
5
#include <boost/multiprecision/float128.hpp>
incardon's avatar
incardon committed
6
#endif
incardon's avatar
incardon committed
7

8
namespace openfpm
incardon's avatar
incardon committed
9
{
10
	namespace math
incardon's avatar
incardon committed
11
	{
incardon's avatar
incardon committed
12
13
14
15
		/*! \brief calculate the factorial
		 *
		 * calculate the factorial of a number
		 *
16
17
18
		 * # Example
		 * \snippet mathutil_unit_test.hpp factorial usage
		 *
incardon's avatar
incardon committed
19
20
21
22
		 * \param f number
		 * \return the factorial
		 *
		 */
23
		static inline size_t factorial(size_t f)
24
25
		{
			size_t fc = 1;
incardon's avatar
incardon committed
26

27
28
29
30
			for (size_t s = 2 ; s <= f ; s++)
			{
				fc *= s;
			}
incardon's avatar
incardon committed
31

32
33
			return fc;
		}
incardon's avatar
incardon committed
34

35
36
37
38
39
40
		/*! \brief C(n,k) Combination of n objects taken on group of k elements
		 *
		 * C(n,k) Combination of n objects taken on group of k elements, defined as
		 *
		 * n!/(k!(n-k)!)
		 *
41
42
43
44
45
46
		 * # Example
		 * \snippet mathutil_unit_test.hpp Combination usage
		 *
		 * \param n
		 * \param k
		 *
47
		 */
48
		static inline size_t C(size_t n, size_t k)
49
50
51
		{
			return factorial(n)/(factorial(k)*factorial(n-k));
		}
incardon's avatar
incardon committed
52

53
54
		/*! \brief Round to the nearest bigger power of 2 number
		 *
incardon's avatar
incardon committed
55
		 * \param n number
56
57
		 * \return nearest bigger power of 2 number
		 *
58
59
60
61
		 * # Example
		 * \snippet mathutil_unit_test.hpp round to big pow
		 *
		 *
62
		 */
63
		static inline size_t round_big_2(size_t n)
64
65
66
67
68
69
70
71
		{
			n--;
			n |= n >> 1;   // Divide by 2^k for consecutive doublings of k up to 32,
			n |= n >> 2;   // and then or the results.
			n |= n >> 4;
			n |= n >> 8;
			n |= n >> 16;
			n++;
incardon's avatar
incardon committed
72

73
74
75
76
77
			return n;
		}


		/*! \brief It calculate at compile-time and runtime the power with recursion
78
79
80
		 *
		 * # Example
		 * \snippet mathutil_unit_test.hpp pow
81
82
83
84
85
86
87
88
89
		 *
		 * \tparam type of the pow expression
		 *
		 * \param base
		 * \param exponent
		 *
		 */

		template<class T>
incardon's avatar
incardon committed
90
		inline constexpr size_t pow(const T base, unsigned const exponent)
91
92
93
94
		{
			// (parentheses not required in next line)
			return (exponent == 0) ? 1 : (base * pow(base, exponent-1));
		}
incardon's avatar
incardon committed
95
96

		/* \brief Return the positive modulo of a number
97
98
99
		 *
		 * # Example
		 * \snippet mathutil_unit_test.hpp positive modulo
incardon's avatar
incardon committed
100
101
102
103
104
		 *
		 * \param i number
		 * \param n modulo
		 *
		 */
105
		static inline long int positive_modulo(long int i, long int n)
incardon's avatar
incardon committed
106
107
108
		{
		    return (i % n + n) % n;
		}
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

		/*! \brief Bound the position to be inside p2 and p1
		 *
		 * Given pos = 10.9, p2 = 1.0 and p1 = 0.1 and l = p2 - p1 = 1.0,
		 * it return 10.9 - (integer)( (10.9 - 0.1) / 1.0) * 1.0
		 *
		 * # Example
		 * \snippet mathutil_unit_test.hpp periodic
		 *
		 * \param pos
		 * \param l
		 * \param b
		 *
		 * \return the bound number
		 *
		 */
		template<typename T> static inline T periodic(const T & pos, const T & p2, const T & p1)
		{
			T pos_tmp;

			pos_tmp = pos - (p2 - p1) * (long int)( (pos -p1) / (p2 - p1));
			pos_tmp += (pos < p1)?(p2 - p1):0;

			return pos_tmp;
		}

		/*! \brief Bound the position to be inside p2 and p1
		 *
		 * It is like periodic but faster
		 *
		 * \warning pos should not overshoot the domain more than one time, for example
		 *          if p2 = 1.1 and p1 = 0.1, pos can be between 2.1 and -0.9
		 *
		 * # Example
		 * \snippet mathutil_unit_test.hpp periodic_l
		 *
		 * \param pos
		 * \param l
		 * \param b
		 *
		 * \return the bound number
		 *
		 */
		template<typename T> static inline T periodic_l(const T & pos, const T & p2, const T & p1)
		{
			T pos_tmp = pos;

Pietro Incardona's avatar
Pietro Incardona committed
156
			if (pos >= p2)
157
158
159
			{
				pos_tmp = p1 + (pos - p2);
			}
160
			else if (pos < p1)
161
162
			{
				pos_tmp = p2 - (p1 - pos);
163

164
165
166
167
168
				// This is a round off error fix
				// if the shift bring exactly on p2 p2 we back the particle to p1
				if (pos_tmp == p2)
					pos_tmp = p1;
			}
169
170
171

			return pos_tmp;
		}
incardon's avatar
incardon committed
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192

		/*! \brief floor math function
		 *
		 *
		 */
		inline long int size_t_floor(double x)
		{
		  size_t i = (long int)x; /* truncate */
		  return i - ( i > x ); /* convert trunc to floor */
		}

		/*! \brief floor math function
		 *
		 *
		 */
		inline long int size_t_floor(float x)
		{
		  size_t i = (long int)x; /* truncate */
		  return i - ( i > x ); /* convert trunc to floor */
		}

incardon's avatar
incardon committed
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
		const int tab64[64] = {
		    63,  0, 58,  1, 59, 47, 53,  2,
		    60, 39, 48, 27, 54, 33, 42,  3,
		    61, 51, 37, 40, 49, 18, 28, 20,
		    55, 30, 34, 11, 43, 14, 22,  4,
		    62, 57, 46, 52, 38, 26, 32, 41,
		    50, 36, 17, 19, 29, 10, 13, 21,
		    56, 45, 25, 31, 35, 16,  9, 12,
		    44, 24, 15,  8, 23,  7,  6,  5};

		/*! \brief Calculate the logarith 2 of a 64 bit integer
		 *
		 * \param value
		 *
		 */
		inline int log2_64 (uint64_t value)
		{
		    value |= value >> 1;
		    value |= value >> 2;
		    value |= value >> 4;
		    value |= value >> 8;
		    value |= value >> 16;
		    value |= value >> 32;
		    return tab64[((uint64_t)((value - (value >> 1))*0x07EDD5E59A4E28C2)) >> 58];
		}

incardon's avatar
incardon committed
219
220
#ifdef HAVE_LIBQUADMATH

incardon's avatar
incardon committed
221
222
223
224
225
226
227
228
229
		/*! \brief floor math function
		 *
		 *
		 */
		inline long int size_t_floor(boost::multiprecision::float128 x)
		{
		  size_t i = (long int)x; /* truncate */
		  return i - ( i > x ); /* convert trunc to floor */
		}
incardon's avatar
incardon committed
230
231

#endif
232
233
	}
}
incardon's avatar
incardon committed
234
235

#endif