// credits to Daniel Lemire: https://github.com/lemire/fastmod

#ifndef FASTMOD_H
#define FASTMOD_H

#ifndef __cplusplus
#include <stdbool.h>
#include <stdint.h>
#else
// In C++ <cstdbool>/<stdbool.h> are irelevant as bool is already a type
#include <cstdint>
#endif

#ifndef __cplusplus
#define FASTMOD_API static inline
#else
// In C++ we mark all the functions inline.
// If C++14 relaxed constexpr is supported we use constexpr so functions
// can be used at compile-time.
#if __cpp_constexpr >= 201304 && !defined(_MSC_VER)
// visual studio does not like constexpr
#define FASTMOD_API constexpr
#define FASTMOD_CONSTEXPR constexpr
#else
#define FASTMOD_API inline
#define FASTMOD_CONSTEXPR
#endif
#endif

#ifdef _MSC_VER
#include <intrin.h>
#endif

#ifdef __cplusplus
namespace fastmod {
#endif

#ifdef _MSC_VER

// __umulh is only available in x64 mode under Visual Studio: don't compile to 32-bit!
FASTMOD_API uint64_t mul128_u32(uint64_t lowbits, uint32_t d) {
    return __umulh(lowbits, d);
}

#else  // _MSC_VER NOT defined

FASTMOD_API uint64_t mul128_u32(uint64_t lowbits, uint32_t d) {
    return ((__uint128_t)lowbits * d) >> 64;
}

FASTMOD_API uint64_t mul128_s32(uint64_t lowbits, int32_t d) {
    return ((__int128_t)lowbits * d) >> 64;
}

// This is for the 64-bit functions.
FASTMOD_API uint64_t mul128_u64(__uint128_t lowbits, uint64_t d) {
    __uint128_t bottom_half = (lowbits & UINT64_C(0xFFFFFFFFFFFFFFFF)) * d;  // Won't overflow
    bottom_half >>= 64;  // Only need the top 64 bits, as we'll shift the lower half away;
    __uint128_t top_half = (lowbits >> 64) * d;
    __uint128_t both_halves = bottom_half + top_half;  // Both halves are already shifted down by 64
    both_halves >>= 64;                                // Get top half of both_halves
    return (uint64_t)both_halves;
}

#endif  // _MSC_VER

/**
 * Unsigned integers.
 * Usage:
 *  uint32_t d = ... ; // divisor, should be non-zero
 *  uint64_t M = computeM_u32(d); // do once
 *  fastmod_u32(a,M,d) is a % d for all 32-bit a.
 *
 **/

// M = ceil( (1<<64) / d ), d > 0
FASTMOD_API uint64_t computeM_u32(uint32_t d) {
    return UINT64_C(0xFFFFFFFFFFFFFFFF) / d + 1;
}

// fastmod computes (a % d) given precomputed M
FASTMOD_API uint32_t fastmod_u32(uint32_t a, uint64_t M, uint32_t d) {
    uint64_t lowbits = M * a;
    return (uint32_t)(mul128_u32(lowbits, d));
}

// fastmod computes (a / d) given precomputed M for d>1
FASTMOD_API uint32_t fastdiv_u32(uint32_t a, uint64_t M) {
    return (uint32_t)(mul128_u32(M, a));
}

// given precomputed M, checks whether n % d == 0
FASTMOD_API bool is_divisible(uint32_t n, uint64_t M) {
    return n * M <= M - 1;
}

/**
 * signed integers
 * Usage:
 *  int32_t d = ... ; // should be non-zero and between [-2147483647,2147483647]
 *  int32_t positive_d = d < 0 ? -d : d; // absolute value
 *  uint64_t M = computeM_s32(d); // do once
 *  fastmod_s32(a,M,positive_d) is a % d for all 32-bit a.
 **/

// M = floor( (1<<64) / d ) + 1
// you must have that d is different from 0 and -2147483648
// if d = -1 and a = -2147483648, the result is undefined
FASTMOD_API uint64_t computeM_s32(int32_t d) {
    if (d < 0) d = -d;
    return UINT64_C(0xFFFFFFFFFFFFFFFF) / d + 1 + ((d & (d - 1)) == 0 ? 1 : 0);
}

// fastmod computes (a % d) given precomputed M,
// you should pass the absolute value of d
FASTMOD_API int32_t fastmod_s32(int32_t a, uint64_t M, int32_t positive_d) {
    uint64_t lowbits = M * a;
    int32_t highbits = mul128_u32(lowbits, positive_d);
    return highbits - ((positive_d - 1) & (a >> 31));
}

#ifndef _MSC_VER
// fastmod computes (a / d) given precomputed M, assumes that d must not
// be one of -1, 1, or -2147483648
FASTMOD_API int32_t fastdiv_s32(int32_t a, uint64_t M, int32_t d) {
    uint64_t highbits = mul128_s32(M, a);
    highbits += (a < 0 ? 1 : 0);
    if (d < 0) return -(int32_t)(highbits);
    return (int32_t)(highbits);
}

// What follows is the 64-bit functions.
// They are currently not supported on Visual Studio
// due to the lack of a mul128_u64 function.
// They may not be faster than what the compiler
// can produce.

FASTMOD_API __uint128_t computeM_u64(uint64_t d) {
    // what follows is just ((__uint128_t)0 - 1) / d) + 1 spelled out
    __uint128_t M = UINT64_C(0xFFFFFFFFFFFFFFFF);
    M <<= 64;
    M |= UINT64_C(0xFFFFFFFFFFFFFFFF);
    M /= d;
    M += 1;
    return M;
}

FASTMOD_API __uint128_t computeM_s64(int64_t d) {
    if (d < 0) d = -d;
    __uint128_t M = UINT64_C(0xFFFFFFFFFFFFFFFF);
    M <<= 64;
    M |= UINT64_C(0xFFFFFFFFFFFFFFFF);
    M /= d;
    M += 1;
    M += ((d & (d - 1)) == 0 ? 1 : 0);
    return M;
}

FASTMOD_API uint64_t fastmod_u64(uint64_t a, __uint128_t M, uint64_t d) {
    __uint128_t lowbits = M * a;
    return mul128_u64(lowbits, d);
}

FASTMOD_API uint64_t fastdiv_u64(uint64_t a, __uint128_t M) {
    return mul128_u64(M, a);
}

// End of the 64-bit functions

#endif  // #ifndef _MSC_VER

#ifdef __cplusplus

template <uint32_t d>
FASTMOD_API uint32_t fastmod(uint32_t x) {
    FASTMOD_CONSTEXPR uint64_t v = computeM_u32(d);
    return fastmod_u32(x, v, d);
}
template <uint32_t d>
FASTMOD_API uint32_t fastdiv(uint32_t x) {
    FASTMOD_CONSTEXPR uint64_t v = computeM_u32(d);
    return fastdiv_u32(x, v);
}
template <int32_t d>
FASTMOD_API int32_t fastmod(int32_t x) {
    FASTMOD_CONSTEXPR uint64_t v = computeM_s32(d);
    return fastmod_s32(x, v, d);
}
template <int32_t d>
FASTMOD_API int32_t fastdiv(int32_t x) {
    FASTMOD_CONSTEXPR uint64_t v = computeM_s32(d);
    return fastdiv_s32(x, v, d);
}

}  // fastmod
#endif

// There's no reason to polute the global scope with this macro once its use ends
// This won't create any problems as the preprocessor will have done its thing once
// it reaches this point
#undef FASTMOD_API
#undef FASTMOD_CONSTEXPR

#endif  // FASTMOD_H