mathutil.hpp

#ifndef MATHUTIL_HPP
#define MATHUTIL_HPP

#ifdef HAVE_LIBQUADMATH
#include <boost/multiprecision/float128.hpp>
#endif

namespace openfpm
{
    namespace math
    {
        static inline size_t factorial(size_t f)
        {
            size_t fc = 1;

            for (size_t s = 2 ; s <= f ; s++)
            {
                fc *= s;
            }

            return fc;
        }

        static inline size_t C(size_t n, size_t k)
        {
            return factorial(n)/(factorial(k)*factorial(n-k));
        }

        static inline size_t round_big_2(size_t n)
        {
            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++;

            return n;
        }


        template<class T>
        inline constexpr size_t pow(const T base, unsigned const exponent)
        {
            // (parentheses not required in next line)
            return (exponent == 0) ? 1 : (base * pow(base, exponent-1));
        }

        /* \brief Return the positive modulo of a number
         *
         * # Example
         * \snippet mathutil_unit_test.hpp positive modulo
         *
         * \param i number
         * \param n modulo
         *
         */
        static inline long int positive_modulo(long int i, long int n)
        {
            return (i % n + n) % n;
        }

        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;
        }

        template<typename T> static inline T periodic_l(const T & pos, const T & p2, const T & p1)
        {
            T pos_tmp = pos;

            if (pos >= p2)
            {
                pos_tmp = p1 + (pos - p2);
            }
            else if (pos < p1)
            {
                pos_tmp = p2 - (p1 - pos);

                // 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;
            }

            return pos_tmp;
        }

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

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

        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};

        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];
        }

#ifdef HAVE_LIBQUADMATH

        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 */
        }

#endif
    }
}

#endif