1#ifndef MATHUTIL_HPP
2#define MATHUTIL_HPP
3
4#ifdef HAVE_LIBQUADMATH
5#include <boost/multiprecision/float128.hpp>
6#endif
7
8#include "util/cuda_util.hpp"
9
10namespace openfpm
11{
12 namespace math
13 {
14 /*! \brief return the sign of an integer
15 *
16 * \param val integer
17 *
18 * \return the sign of the function
19 *
20 */
21 template <typename T>
22 int sgn(T val)
23 {
24 return (T(0) < val) - (val < T(0));
25 }
26
27 /*! \brief calculate the factorial
28 *
29 * calculate the factorial of a number
30 *
31 * # Example
32 * \snippet mathutil_unit_test.hpp factorial usage
33 *
34 * \param f number
35 * \return the factorial
36 *
37 */
38 static inline size_t factorial(size_t f)
39 {
40 size_t fc = 1;
41
42 for (size_t s = 2 ; s <= f ; s++)
43 {
44 fc *= s;
45 }
46
47 return fc;
48 }
49
50 /*! \brief C(n,k) Combination of n objects taken on group of k elements
51 *
52 * C(n,k) Combination of n objects taken on group of k elements, defined as
53 *
54 * n!/(k!(n-k)!)
55 *
56 * # Example
57 * \snippet mathutil_unit_test.hpp Combination usage
58 *
59 * \param n
60 * \param k
61 *
62 */
63 static inline size_t C(size_t n, size_t k)
64 {
65 return factorial(n)/(factorial(k)*factorial(n-k));
66 }
67
68 /*! \brief Round to the nearest bigger power of 2 number
69 *
70 * \param n number
71 * \return nearest bigger power of 2 number
72 *
73 * # Example
74 * \snippet mathutil_unit_test.hpp round to big pow
75 *
76 *
77 */
78 static inline size_t round_big_2(size_t n)
79 {
80 n--;
81 n |= n >> 1; // Divide by 2^k for consecutive doublings of k up to 32,
82 n |= n >> 2; // and then or the results.
83 n |= n >> 4;
84 n |= n >> 8;
85 n |= n >> 16;
86 n++;
87
88 return n;
89 }
90
91
92 /*! \brief It calculate at compile-time and runtime the power with recursion
93 *
94 * # Example
95 * \snippet mathutil_unit_test.hpp pow
96 *
97 * \tparam type of the pow expression
98 *
99 * \param base
100 * \param exponent
101 *
102 */
103
104 template<class T>
105 inline constexpr size_t pow(const T base, unsigned const exponent)
106 {
107 // (parentheses not required in next line)
108 return (exponent == 0) ? 1 : (base * pow(base, exponent-1));
109 }
110
111 /* \brief Return the positive modulo of a number
112 *
113 * # Example
114 * \snippet mathutil_unit_test.hpp positive modulo
115 *
116 * \param i number
117 * \param n modulo
118 *
119 */
120 static inline long int positive_modulo(long int i, long int n)
121 {
122 return (i % n + n) % n;
123 }
124
125 /*! \brief Bound the position to be inside p2 and p1
126 *
127 * Given pos = 10.9, p2 = 1.0 and p1 = 0.1 and l = p2 - p1 = 1.0,
128 * it return 10.9 - (integer)( (10.9 - 0.1) / 1.0) * 1.0
129 *
130 * # Example
131 * \snippet mathutil_unit_test.hpp periodic
132 *
133 * \param pos
134 * \param l
135 * \param b
136 *
137 * \return the bound number
138 *
139 */
140 template<typename T> static inline T periodic(const T & pos, const T & p2, const T & p1)
141 {
142 T pos_tmp;
143
144 pos_tmp = pos - (p2 - p1) * (long int)( (pos -p1) / (p2 - p1));
145 pos_tmp += (pos < p1)?(p2 - p1):0;
146
147 return pos_tmp;
148 }
149
150 /*! \brief Bound the position to be inside p2 and p1
151 *
152 * It is like periodic but faster
153 *
154 * \warning pos should not overshoot the domain more than one time, for example
155 * if p2 = 1.1 and p1 = 0.1, pos can be between 2.1 and -0.9
156 *
157 * # Example
158 * \snippet mathutil_unit_test.hpp periodic_l
159 *
160 * \param pos
161 * \param l
162 * \param b
163 *
164 * \return the bound number
165 *
166 */
167 template<typename T> __device__ __host__ static inline T periodic_l(const T & pos, const T & p2, const T & p1)
168 {
169 T pos_tmp = pos;
170
171 if (pos >= p2)
172 {
173 pos_tmp = p1 + (pos - p2);
174 }
175 else if (pos < p1)
176 {
177 pos_tmp = p2 - (p1 - pos);
178
179 // This is a round off error fix
180 // if the shift bring exactly on p2 p2 we back the particle to p1
181 if (pos_tmp == p2)
182 pos_tmp = p1;
183 }
184
185 return pos_tmp;
186 }
187
188 /*! \brief floor math function
189 *
190 *
191 */
192 inline long int size_t_floor(double x)
193 {
194 size_t i = (long int)x; /* truncate */
195 return i - ( i > x ); /* convert trunc to floor */
196 }
197
198 /*! \brief floor math function
199 *
200 *
201 */
202 inline long int size_t_floor(float x)
203 {
204 size_t i = (long int)x; /* truncate */
205 return i - ( i > x ); /* convert trunc to floor */
206 }
207
208 /*! \brief floor math function
209 *
210 *
211 */
212 __device__ __host__ inline int uint_floor(double x)
213 {
214 unsigned int i = (int)x; /* truncate */
215 return i - ( i > x ); /* convert trunc to floor */
216 }
217
218 /*! \brief floor math function
219 *
220 *
221 */
222 __device__ __host__ inline int uint_floor(float x)
223 {
224 unsigned int i = (int)x; /* truncate */
225 return i - ( i > x ); /* convert trunc to floor */
226 }
227
228 const int tab64[64] = {
229 63, 0, 58, 1, 59, 47, 53, 2,
230 60, 39, 48, 27, 54, 33, 42, 3,
231 61, 51, 37, 40, 49, 18, 28, 20,
232 55, 30, 34, 11, 43, 14, 22, 4,
233 62, 57, 46, 52, 38, 26, 32, 41,
234 50, 36, 17, 19, 29, 10, 13, 21,
235 56, 45, 25, 31, 35, 16, 9, 12,
236 44, 24, 15, 8, 23, 7, 6, 5};
237
238 /*! \brief Calculate the logarith 2 of a 64 bit integer
239 *
240 * \param value
241 *
242 */
243 inline int log2_64 (uint64_t value)
244 {
245 value |= value >> 1;
246 value |= value >> 2;
247 value |= value >> 4;
248 value |= value >> 8;
249 value |= value >> 16;
250 value |= value >> 32;
251 return tab64[((uint64_t)((value - (value >> 1))*0x07EDD5E59A4E28C2)) >> 58];
252 }
253
254
255#ifdef HAVE_LIBQUADMATH
256
257 /*! \brief floor math function
258 *
259 *
260 */
261 inline long int size_t_floor(boost::multiprecision::float128 x)
262 {
263 size_t i = (long int)x; /* truncate */
264 return i - ( i > x ); /* convert trunc to floor */
265 }
266
267#endif
268 }
269}
270
271#endif
272