mathutil.hpp

#ifndef MATHUTIL_HPP
#define MATHUTIL_HPP


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 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 -= (p2 - p1);
            else if (pos < p1)
                pos_tmp += (p2 - p1);


            return pos_tmp;
        }
    }
}

#endif