// libdivide.h - Optimized integer division

// https://libdivide.com

//

// Copyright (C) 2010 - 2022 ridiculous_fish, <libdivide@ridiculousfish.com>

// Copyright (C) 2016 - 2022 Kim Walisch, <kim.walisch@gmail.com>

//

// libdivide is dual-licensed under the Boost or zlib licenses.

// You may use libdivide under the terms of either of these.

// See LICENSE.txt for more details.

#ifndef LIBDIVIDE_H

#define LIBDIVIDE_H

// *** Version numbers are auto generated - do not edit ***

#define LIBDIVIDE_VERSION "5.2.0"

#define LIBDIVIDE_VERSION_MAJOR 5

#define LIBDIVIDE_VERSION_MINOR 2

#define LIBDIVIDE_VERSION_PATCH 0

#include <stdint.h>

#if !defined(__AVR__) && __STDC_HOSTED__ != 0

#include <stdio.h>

#include <stdlib.h>

#endif

#if defined(_MSC_VER) && (defined(__cplusplus) && (__cplusplus >= 202002L)) || \

    (defined(_MSVC_LANG) && (_MSVC_LANG >= 202002L))

#include <limits.h>

#include <type_traits>

#define LIBDIVIDE_VC_CXX20

#endif

#if defined(LIBDIVIDE_SSE2)

#include <emmintrin.h>

#endif

#if defined(LIBDIVIDE_AVX2) || defined(LIBDIVIDE_AVX512)

#include <immintrin.h>

#endif

#if defined(LIBDIVIDE_NEON)

#include <arm_neon.h>

#endif

// Clang-cl prior to Visual Studio 2022 doesn't include __umulh/__mulh intrinsics

#if defined(_MSC_VER) && (!defined(__clang__) || _MSC_VER > 1930) && \

    (defined(_M_X64) || defined(_M_ARM64) || defined(_M_HYBRID_X86_ARM64) || defined(_M_ARM64EC))

#define LIBDIVIDE_MULH_INTRINSICS

#endif

#if defined(_MSC_VER)

#if defined(LIBDIVIDE_MULH_INTRINSICS) || !defined(__clang__)

#include <intrin.h>

#endif

#ifndef __clang__

#pragma warning(push)

// 4146: unary minus operator applied to unsigned type, result still unsigned

#pragma warning(disable : 4146)

// 4204: nonstandard extension used : non-constant aggregate initializer

#pragma warning(disable : 4204)

#endif

#define LIBDIVIDE_VC

#endif

#if !defined(__has_builtin)

#define __has_builtin(x) 0

#endif

#if defined(__SIZEOF_INT128__)

#define HAS_INT128_T

// clang-cl on Windows does not yet support 128-bit division

#if !(defined(__clang__) && defined(LIBDIVIDE_VC))

#define HAS_INT128_DIV

#endif

#endif

#if defined(__x86_64__) || defined(_M_X64)

#define LIBDIVIDE_X86_64

#endif

#if defined(__i386__)

#define LIBDIVIDE_i386

#endif

#if defined(__GNUC__) || defined(__clang__)

#define LIBDIVIDE_GCC_STYLE_ASM

#endif

#if defined(__cplusplus) || defined(LIBDIVIDE_VC)

#define LIBDIVIDE_FUNCTION __FUNCTION__

#else

#define LIBDIVIDE_FUNCTION __func__

#endif

// Set up forced inlining if possible.

// We need both the attribute and keyword to avoid "might not be inlineable" warnings.

#ifdef __has_attribute

#if __has_attribute(always_inline)

#define LIBDIVIDE_INLINE __attribute__((always_inline)) inline

#endif

#endif

#ifndef LIBDIVIDE_INLINE

#ifdef _MSC_VER

#define LIBDIVIDE_INLINE __forceinline

#else

#define LIBDIVIDE_INLINE inline

#endif

#endif

#if defined(__AVR__) || __STDC_HOSTED__ == 0

#define LIBDIVIDE_ERROR(msg)

#else

#define LIBDIVIDE_ERROR(msg)                                                                     \

    do {                                                                                         \

        fprintf(stderr, "libdivide.h:%d: %s(): Error: %s\n", __LINE__, LIBDIVIDE_FUNCTION, msg); \

        abort();                                                                                 \

    } while (0)

#endif

#if defined(LIBDIVIDE_ASSERTIONS_ON) && !defined(__AVR__) && __STDC_HOSTED__ != 0

#define LIBDIVIDE_ASSERT(x)                                                           \

    do {                                                                              \

        if (!(x)) {                                                                   \

            fprintf(stderr, "libdivide.h:%d: %s(): Assertion failed: %s\n", __LINE__, \

                LIBDIVIDE_FUNCTION, #x);                                              \

            abort();                                                                  \

        }                                                                             \

    } while (0)

#else

#define LIBDIVIDE_ASSERT(x)

#endif

#ifdef __cplusplus

// For constexpr zero initialization, c++11 might handle things ok,

// but just limit to at least c++14 to ensure we don't break anyone's code:

// Use https://en.cppreference.com/w/cpp/feature_test#cpp_constexpr

#if defined(__cpp_constexpr) && (__cpp_constexpr >= 201304L)

#define LIBDIVIDE_CONSTEXPR constexpr LIBDIVIDE_INLINE

// Supposedly, MSVC might not implement feature test macros right:

// https://stackoverflow.com/questions/49316752/feature-test-macros-not-working-properly-in-visual-c

// so check that _MSVC_LANG corresponds to at least c++14, and _MSC_VER corresponds to at least VS

// 2017 15.0 (for extended constexpr support:

// https://learn.microsoft.com/en-us/cpp/overview/visual-cpp-language-conformance?view=msvc-170)

#elif (defined(_MSC_VER) && _MSC_VER >= 1910) && (defined(_MSVC_LANG) && _MSVC_LANG >= 201402L)

#define LIBDIVIDE_CONSTEXPR constexpr LIBDIVIDE_INLINE

#else

#define LIBDIVIDE_CONSTEXPR LIBDIVIDE_INLINE

#endif

namespace libdivide {

#endif

#if defined(_MSC_VER) && !defined(__clang__)

#if defined(LIBDIVIDE_VC_CXX20)

static LIBDIVIDE_CONSTEXPR int __builtin_clz(unsigned x) {

    if (std::is_constant_evaluated()) {

        for (int i = 0; i < sizeof(x) * CHAR_BIT; ++i) {

            if (x >> (sizeof(x) * CHAR_BIT - 1 - i)) return i;

        }

        return sizeof(x) * CHAR_BIT;

    }

#else

static LIBDIVIDE_INLINE int __builtin_clz(unsigned x) {

#endif

#if defined(_M_ARM) || defined(_M_ARM64) || defined(_M_HYBRID_X86_ARM64) || defined(_M_ARM64EC)

    return (int)_CountLeadingZeros(x);

#elif defined(__AVX2__) || defined(__LZCNT__)

    return (int)_lzcnt_u32(x);

#else

    unsigned long r;

    _BitScanReverse(&r, x);

    return (int)(r ^ 31);

#endif

}

#if defined(LIBDIVIDE_VC_CXX20)

static LIBDIVIDE_CONSTEXPR int __builtin_clzll(unsigned long long x) {

    if (std::is_constant_evaluated()) {

        for (int i = 0; i < sizeof(x) * CHAR_BIT; ++i) {

            if (x >> (sizeof(x) * CHAR_BIT - 1 - i)) return i;

        }

        return sizeof(x) * CHAR_BIT;

    }

#else

static LIBDIVIDE_INLINE int __builtin_clzll(unsigned long long x) {

#endif

#if defined(_M_ARM) || defined(_M_ARM64) || defined(_M_HYBRID_X86_ARM64) || defined(_M_ARM64EC)

    return (int)_CountLeadingZeros64(x);

#elif defined(_WIN64)

#if defined(__AVX2__) || defined(__LZCNT__)

    return (int)_lzcnt_u64(x);

#else

    unsigned long r;

    _BitScanReverse64(&r, x);

    return (int)(r ^ 63);

#endif

#else

    int l = __builtin_clz((unsigned)x) + 32;

    int h = __builtin_clz((unsigned)(x >> 32));

    return !!((unsigned)(x >> 32)) ? h : l;

#endif

}

#endif // defined(_MSC_VER) && !defined(__clang__)

// pack divider structs to prevent compilers from padding.

// This reduces memory usage by up to 43% when using a large

// array of libdivide dividers and improves performance

// by up to 10% because of reduced memory bandwidth.

#pragma pack(push, 1)

struct libdivide_u16_t {

    uint16_t magic;

    uint8_t more;

};

struct libdivide_s16_t {

    int16_t magic;

    uint8_t more;

};

struct libdivide_u32_t {

    uint32_t magic;

    uint8_t more;

};

struct libdivide_s32_t {

    int32_t magic;

    uint8_t more;

};

struct libdivide_u64_t {

    uint64_t magic;

    uint8_t more;

};

struct libdivide_s64_t {

    int64_t magic;

    uint8_t more;

};

struct libdivide_u16_branchfree_t {

    uint16_t magic;

    uint8_t more;

};

struct libdivide_s16_branchfree_t {

    int16_t magic;

    uint8_t more;

};

struct libdivide_u32_branchfree_t {

    uint32_t magic;

    uint8_t more;

};

struct libdivide_s32_branchfree_t {

    int32_t magic;

    uint8_t more;

};

struct libdivide_u64_branchfree_t {

    uint64_t magic;

    uint8_t more;

};

struct libdivide_s64_branchfree_t {

    int64_t magic;

    uint8_t more;

};

#pragma pack(pop)

// Explanation of the "more" field:

//

// * Bits 0-5 is the shift value (for shift path or mult path).

// * Bit 6 is the add indicator for mult path.

// * Bit 7 is set if the divisor is negative. We use bit 7 as the negative

//   divisor indicator so that we can efficiently use sign extension to

//   create a bitmask with all bits set to 1 (if the divisor is negative)

//   or 0 (if the divisor is positive).

//

// u32: [0-4] shift value

//      [5] ignored

//      [6] add indicator

//      magic number of 0 indicates shift path

//

// s32: [0-4] shift value

//      [5] ignored

//      [6] add indicator

//      [7] indicates negative divisor

//      magic number of 0 indicates shift path

//

// u64: [0-5] shift value

//      [6] add indicator

//      magic number of 0 indicates shift path

//

// s64: [0-5] shift value

//      [6] add indicator

//      [7] indicates negative divisor

//      magic number of 0 indicates shift path

//

// In s32 and s64 branchfree modes, the magic number is negated according to

// whether the divisor is negated. In branchfree strategy, it is not negated.

enum {

    LIBDIVIDE_16_SHIFT_MASK = 0x1F,

    LIBDIVIDE_32_SHIFT_MASK = 0x1F,

    LIBDIVIDE_64_SHIFT_MASK = 0x3F,

    LIBDIVIDE_ADD_MARKER = 0x40,

    LIBDIVIDE_NEGATIVE_DIVISOR = 0x80

};

static LIBDIVIDE_INLINE struct libdivide_s16_t libdivide_s16_gen(int16_t d);

static LIBDIVIDE_INLINE struct libdivide_u16_t libdivide_u16_gen(uint16_t d);

static LIBDIVIDE_INLINE struct libdivide_s32_t libdivide_s32_gen(int32_t d);

static LIBDIVIDE_INLINE struct libdivide_u32_t libdivide_u32_gen(uint32_t d);

static LIBDIVIDE_INLINE struct libdivide_s64_t libdivide_s64_gen(int64_t d);

static LIBDIVIDE_INLINE struct libdivide_u64_t libdivide_u64_gen(uint64_t d);

static LIBDIVIDE_INLINE struct libdivide_s16_branchfree_t libdivide_s16_branchfree_gen(int16_t d);

static LIBDIVIDE_INLINE struct libdivide_u16_branchfree_t libdivide_u16_branchfree_gen(uint16_t d);

static LIBDIVIDE_INLINE struct libdivide_s32_branchfree_t libdivide_s32_branchfree_gen(int32_t d);

static LIBDIVIDE_INLINE struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t d);

static LIBDIVIDE_INLINE struct libdivide_s64_branchfree_t libdivide_s64_branchfree_gen(int64_t d);

static LIBDIVIDE_INLINE struct libdivide_u64_branchfree_t libdivide_u64_branchfree_gen(uint64_t d);

static LIBDIVIDE_INLINE int16_t libdivide_s16_do_raw(

    int16_t numer, int16_t magic, uint8_t more);

static LIBDIVIDE_INLINE int16_t libdivide_s16_do(

    int16_t numer, const struct libdivide_s16_t *denom);

static LIBDIVIDE_INLINE uint16_t libdivide_u16_do_raw(

    uint16_t numer, uint16_t magic, uint8_t more);

static LIBDIVIDE_INLINE uint16_t libdivide_u16_do(

    uint16_t numer, const struct libdivide_u16_t *denom);

static LIBDIVIDE_INLINE int32_t libdivide_s32_do_raw(

    int32_t numer, int32_t magic, uint8_t more);

static LIBDIVIDE_INLINE int32_t libdivide_s32_do(

    int32_t numer, const struct libdivide_s32_t *denom);

static LIBDIVIDE_INLINE uint32_t libdivide_u32_do_raw(

    uint32_t numer, uint32_t magic, uint8_t more);

static LIBDIVIDE_INLINE uint32_t libdivide_u32_do(

    uint32_t numer, const struct libdivide_u32_t *denom);

static LIBDIVIDE_INLINE int64_t libdivide_s64_do_raw(

    int64_t numer, int64_t magic, uint8_t more);

static LIBDIVIDE_INLINE int64_t libdivide_s64_do(

    int64_t numer, const struct libdivide_s64_t *denom);

static LIBDIVIDE_INLINE uint64_t libdivide_u64_do_raw(

    uint64_t numer, uint64_t magic, uint8_t more);

static LIBDIVIDE_INLINE uint64_t libdivide_u64_do(

    uint64_t numer, const struct libdivide_u64_t *denom);

static LIBDIVIDE_INLINE int16_t libdivide_s16_branchfree_do(

    int16_t numer, const struct libdivide_s16_branchfree_t *denom);

static LIBDIVIDE_INLINE uint16_t libdivide_u16_branchfree_do(

    uint16_t numer, const struct libdivide_u16_branchfree_t *denom);

static LIBDIVIDE_INLINE int32_t libdivide_s32_branchfree_do(

    int32_t numer, const struct libdivide_s32_branchfree_t *denom);

static LIBDIVIDE_INLINE uint32_t libdivide_u32_branchfree_do(

    uint32_t numer, const struct libdivide_u32_branchfree_t *denom);

static LIBDIVIDE_INLINE int64_t libdivide_s64_branchfree_do(

    int64_t numer, const struct libdivide_s64_branchfree_t *denom);

static LIBDIVIDE_INLINE uint64_t libdivide_u64_branchfree_do(

    uint64_t numer, const struct libdivide_u64_branchfree_t *denom);

static LIBDIVIDE_INLINE int16_t libdivide_s16_recover(const struct libdivide_s16_t *denom);

static LIBDIVIDE_INLINE uint16_t libdivide_u16_recover(const struct libdivide_u16_t *denom);

static LIBDIVIDE_INLINE int32_t libdivide_s32_recover(const struct libdivide_s32_t *denom);

static LIBDIVIDE_INLINE uint32_t libdivide_u32_recover(const struct libdivide_u32_t *denom);

static LIBDIVIDE_INLINE int64_t libdivide_s64_recover(const struct libdivide_s64_t *denom);

static LIBDIVIDE_INLINE uint64_t libdivide_u64_recover(const struct libdivide_u64_t *denom);

static LIBDIVIDE_INLINE int16_t libdivide_s16_branchfree_recover(

    const struct libdivide_s16_branchfree_t *denom);

static LIBDIVIDE_INLINE uint16_t libdivide_u16_branchfree_recover(

    const struct libdivide_u16_branchfree_t *denom);

static LIBDIVIDE_INLINE int32_t libdivide_s32_branchfree_recover(

    const struct libdivide_s32_branchfree_t *denom);

static LIBDIVIDE_INLINE uint32_t libdivide_u32_branchfree_recover(

    const struct libdivide_u32_branchfree_t *denom);

static LIBDIVIDE_INLINE int64_t libdivide_s64_branchfree_recover(

    const struct libdivide_s64_branchfree_t *denom);

static LIBDIVIDE_INLINE uint64_t libdivide_u64_branchfree_recover(

    const struct libdivide_u64_branchfree_t *denom);

//////// Internal Utility Functions

static LIBDIVIDE_INLINE uint16_t libdivide_mullhi_u16(uint16_t x, uint16_t y) {

    uint32_t xl = x, yl = y;

    uint32_t rl = xl * yl;

    return (uint16_t)(rl >> 16);

}

static LIBDIVIDE_INLINE int16_t libdivide_mullhi_s16(int16_t x, int16_t y) {

    int32_t xl = x, yl = y;

    int32_t rl = xl * yl;

    // needs to be arithmetic shift

    return (int16_t)(rl >> 16);

}

static LIBDIVIDE_INLINE uint32_t libdivide_mullhi_u32(uint32_t x, uint32_t y) {

    uint64_t xl = x, yl = y;

    uint64_t rl = xl * yl;

    return (uint32_t)(rl >> 32);

}

static LIBDIVIDE_INLINE int32_t libdivide_mullhi_s32(int32_t x, int32_t y) {

    int64_t xl = x, yl = y;

    int64_t rl = xl * yl;

    // needs to be arithmetic shift

    return (int32_t)(rl >> 32);

}

static LIBDIVIDE_INLINE uint64_t libdivide_mullhi_u64(uint64_t x, uint64_t y) {

#if defined(LIBDIVIDE_MULH_INTRINSICS)

    return __umulh(x, y);

#elif defined(HAS_INT128_T)

    __uint128_t xl = x, yl = y;

    __uint128_t rl = xl * yl;

    return (uint64_t)(rl >> 64);

#else

    // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64)

    uint32_t mask = 0xFFFFFFFF;

    uint32_t x0 = (uint32_t)(x & mask);

    uint32_t x1 = (uint32_t)(x >> 32);

    uint32_t y0 = (uint32_t)(y & mask);

    uint32_t y1 = (uint32_t)(y >> 32);

    uint32_t x0y0_hi = libdivide_mullhi_u32(x0, y0);

    uint64_t x0y1 = x0 * (uint64_t)y1;

    uint64_t x1y0 = x1 * (uint64_t)y0;

    uint64_t x1y1 = x1 * (uint64_t)y1;

    uint64_t temp = x1y0 + x0y0_hi;

    uint64_t temp_lo = temp & mask;

    uint64_t temp_hi = temp >> 32;

    return x1y1 + temp_hi + ((temp_lo + x0y1) >> 32);

#endif

}

static LIBDIVIDE_INLINE int64_t libdivide_mullhi_s64(int64_t x, int64_t y) {

#if defined(LIBDIVIDE_MULH_INTRINSICS)

    return __mulh(x, y);

#elif defined(HAS_INT128_T)

    __int128_t xl = x, yl = y;

    __int128_t rl = xl * yl;

    return (int64_t)(rl >> 64);

#else

    // full 128 bits are x0 * y0 + (x0 * y1 << 32) + (x1 * y0 << 32) + (x1 * y1 << 64)

    uint32_t mask = 0xFFFFFFFF;

    uint32_t x0 = (uint32_t)(x & mask);

    uint32_t y0 = (uint32_t)(y & mask);

    int32_t x1 = (int32_t)(x >> 32);

    int32_t y1 = (int32_t)(y >> 32);

    uint32_t x0y0_hi = libdivide_mullhi_u32(x0, y0);

    int64_t t = x1 * (int64_t)y0 + x0y0_hi;

    int64_t w1 = x0 * (int64_t)y1 + (t & mask);

    return x1 * (int64_t)y1 + (t >> 32) + (w1 >> 32);

#endif

}

static LIBDIVIDE_INLINE int16_t libdivide_count_leading_zeros16(uint16_t val) {

#if defined(__AVR__)

    // Fast way to count leading zeros

    // On the AVR 8-bit architecture __builtin_clz() works on a int16_t.

    return __builtin_clz(val);

#elif defined(__GNUC__) || __has_builtin(__builtin_clz) || defined(_MSC_VER)

    // Fast way to count leading zeros

    return (int16_t)(__builtin_clz(val) - 16);

#else

    if (val == 0) return 16;

    int16_t result = 4;

    uint16_t hi = 0xFU << 12;

    while ((val & hi) == 0) {

        hi >>= 4;

        result += 4;

    }

    while (val & hi) {

        result -= 1;

        hi <<= 1;

    }

    return result;

#endif

}

static LIBDIVIDE_INLINE int32_t libdivide_count_leading_zeros32(uint32_t val) {

#if defined(__AVR__)

    // Fast way to count leading zeros

    return __builtin_clzl(val);

#elif defined(__GNUC__) || __has_builtin(__builtin_clz) || defined(_MSC_VER)

    // Fast way to count leading zeros

    return __builtin_clz(val);

#else

    if (val == 0) return 32;

    int32_t result = 8;

    uint32_t hi = 0xFFU << 24;

    while ((val & hi) == 0) {

        hi >>= 8;

        result += 8;

    }

    while (val & hi) {

        result -= 1;

        hi <<= 1;

    }

    return result;

#endif

}

static LIBDIVIDE_INLINE int32_t libdivide_count_leading_zeros64(uint64_t val) {

#if defined(__GNUC__) || __has_builtin(__builtin_clzll) || defined(_MSC_VER)

    // Fast way to count leading zeros

    return __builtin_clzll(val);

#else

    uint32_t hi = val >> 32;

    uint32_t lo = val & 0xFFFFFFFF;

    if (hi != 0) return libdivide_count_leading_zeros32(hi);

    return 32 + libdivide_count_leading_zeros32(lo);

#endif

}

// libdivide_32_div_16_to_16: divides a 32-bit uint {u1, u0} by a 16-bit

// uint {v}. The result must fit in 16 bits.

// Returns the quotient directly and the remainder in *r

static LIBDIVIDE_INLINE uint16_t libdivide_32_div_16_to_16(

    uint16_t u1, uint16_t u0, uint16_t v, uint16_t *r) {

    uint32_t n = ((uint32_t)u1 << 16) | u0;

    uint16_t result = (uint16_t)(n / v);

    *r = (uint16_t)(n - result * (uint32_t)v);

    return result;

}

// libdivide_64_div_32_to_32: divides a 64-bit uint {u1, u0} by a 32-bit

// uint {v}. The result must fit in 32 bits.

// Returns the quotient directly and the remainder in *r

static LIBDIVIDE_INLINE uint32_t libdivide_64_div_32_to_32(

    uint32_t u1, uint32_t u0, uint32_t v, uint32_t *r) {

#if (defined(LIBDIVIDE_i386) || defined(LIBDIVIDE_X86_64)) && defined(LIBDIVIDE_GCC_STYLE_ASM)

    uint32_t result;

    __asm__("divl %[v]" : "=a"(result), "=d"(*r) : [v] "r"(v), "a"(u0), "d"(u1));

    return result;

#else

    uint64_t n = ((uint64_t)u1 << 32) | u0;

    uint32_t result = (uint32_t)(n / v);

    *r = (uint32_t)(n - result * (uint64_t)v);

    return result;

#endif

}

// libdivide_128_div_64_to_64: divides a 128-bit uint {numhi, numlo} by a 64-bit uint {den}. The

// result must fit in 64 bits. Returns the quotient directly and the remainder in *r

static LIBDIVIDE_INLINE uint64_t libdivide_128_div_64_to_64(

    uint64_t numhi, uint64_t numlo, uint64_t den, uint64_t *r) {

    // N.B. resist the temptation to use __uint128_t here.

    // In LLVM compiler-rt, it performs a 128/128 -> 128 division which is many times slower than

    // necessary. In gcc it's better but still slower than the divlu implementation, perhaps because

    // it's not LIBDIVIDE_INLINEd.

#if defined(LIBDIVIDE_X86_64) && defined(LIBDIVIDE_GCC_STYLE_ASM)

    uint64_t result;

    __asm__("div %[v]" : "=a"(result), "=d"(*r) : [v] "r"(den), "a"(numlo), "d"(numhi));

    return result;

#else

    // We work in base 2**32.

    // A uint32 holds a single digit. A uint64 holds two digits.

    // Our numerator is conceptually [num3, num2, num1, num0].

    // Our denominator is [den1, den0].

    const uint64_t b = ((uint64_t)1 << 32);

    // The high and low digits of our computed quotient.

    uint32_t q1;

    uint32_t q0;

    // The normalization shift factor.

    int shift;

    // The high and low digits of our denominator (after normalizing).

    // Also the low 2 digits of our numerator (after normalizing).

    uint32_t den1;

    uint32_t den0;

    uint32_t num1;

    uint32_t num0;

    // A partial remainder.

    uint64_t rem;

    // The estimated quotient, and its corresponding remainder (unrelated to true remainder).

    uint64_t qhat;

    uint64_t rhat;

    // Variables used to correct the estimated quotient.

    uint64_t c1;

    uint64_t c2;

    // Check for overflow and divide by 0.

    if (numhi >= den) {

        if (r) *r = ~0ull;

        return ~0ull;

    }

    // Determine the normalization factor. We multiply den by this, so that its leading digit is at

    // least half b. In binary this means just shifting left by the number of leading zeros, so that

    // there's a 1 in the MSB.

    // We also shift numer by the same amount. This cannot overflow because numhi < den.

    // The expression (-shift & 63) is the same as (64 - shift), except it avoids the UB of shifting

    // by 64. The funny bitwise 'and' ensures that numlo does not get shifted into numhi if shift is

    // 0. clang 11 has an x86 codegen bug here: see LLVM bug 50118. The sequence below avoids it.

    shift = libdivide_count_leading_zeros64(den);

    den <<= shift;

    numhi <<= shift;

    numhi |= (numlo >> (-shift & 63)) & (uint64_t)(-(int64_t)shift >> 63);

    numlo <<= shift;

    // Extract the low digits of the numerator and both digits of the denominator.

    num1 = (uint32_t)(numlo >> 32);

    num0 = (uint32_t)(numlo & 0xFFFFFFFFu);

    den1 = (uint32_t)(den >> 32);

    den0 = (uint32_t)(den & 0xFFFFFFFFu);

    // We wish to compute q1 = [n3 n2 n1] / [d1 d0].

    // Estimate q1 as [n3 n2] / [d1], and then correct it.

    // Note while qhat may be 2 digits, q1 is always 1 digit.

    qhat = numhi / den1;

    rhat = numhi % den1;

    c1 = qhat * den0;

    c2 = rhat * b + num1;

    if (c1 > c2) qhat -= (c1 - c2 > den) ? 2 : 1;

    q1 = (uint32_t)qhat;

    // Compute the true (partial) remainder.

    rem = numhi * b + num1 - q1 * den;

    // We wish to compute q0 = [rem1 rem0 n0] / [d1 d0].

    // Estimate q0 as [rem1 rem0] / [d1] and correct it.

    qhat = rem / den1;

    rhat = rem % den1;

    c1 = qhat * den0;

    c2 = rhat * b + num0;

    if (c1 > c2) qhat -= (c1 - c2 > den) ? 2 : 1;

    q0 = (uint32_t)qhat;

    // Return remainder if requested.

    if (r) *r = (rem * b + num0 - q0 * den) >> shift;

    return ((uint64_t)q1 << 32) | q0;

#endif

}

#if !(defined(HAS_INT128_T) && \

      defined(HAS_INT128_DIV))

// Bitshift a u128 in place, left (signed_shift > 0) or right (signed_shift < 0)

static LIBDIVIDE_INLINE void libdivide_u128_shift(

    uint64_t *u1, uint64_t *u0, int32_t signed_shift) {

    if (signed_shift > 0) {

        uint32_t shift = signed_shift;

        *u1 <<= shift;

        *u1 |= *u0 >> (64 - shift);

        *u0 <<= shift;

    } else if (signed_shift < 0) {

        uint32_t shift = -signed_shift;

        *u0 >>= shift;

        *u0 |= *u1 << (64 - shift);

        *u1 >>= shift;

    }

}

#endif

// Computes a 128 / 128 -> 64 bit division, with a 128 bit remainder.

static LIBDIVIDE_INLINE uint64_t libdivide_128_div_128_to_64(

    uint64_t u_hi, uint64_t u_lo, uint64_t v_hi, uint64_t v_lo, uint64_t *r_hi, uint64_t *r_lo) {

#if defined(HAS_INT128_T) && defined(HAS_INT128_DIV)

    __uint128_t ufull = u_hi;

    __uint128_t vfull = v_hi;

    ufull = (ufull << 64) | u_lo;

    vfull = (vfull << 64) | v_lo;

    uint64_t res = (uint64_t)(ufull / vfull);

    __uint128_t remainder = ufull - (vfull * res);

    *r_lo = (uint64_t)remainder;

    *r_hi = (uint64_t)(remainder >> 64);

    return res;

#else

    // Adapted from "Unsigned Doubleword Division" in Hacker's Delight

    // We want to compute u / v

    typedef struct {

        uint64_t hi;

        uint64_t lo;

    } u128_t;

    u128_t u = {u_hi, u_lo};

    u128_t v = {v_hi, v_lo};

    if (v.hi == 0) {

        // divisor v is a 64 bit value, so we just need one 128/64 division

        // Note that we are simpler than Hacker's Delight here, because we know

        // the quotient fits in 64 bits whereas Hacker's Delight demands a full

        // 128 bit quotient

        *r_hi = 0;

        return libdivide_128_div_64_to_64(u.hi, u.lo, v.lo, r_lo);

    }

    // Here v >= 2**64

    // We know that v.hi != 0, so count leading zeros is OK

    // We have 0 <= n <= 63

    uint32_t n = libdivide_count_leading_zeros64(v.hi);

    // Normalize the divisor so its MSB is 1

    u128_t v1t = v;

    libdivide_u128_shift(&v1t.hi, &v1t.lo, n);

    uint64_t v1 = v1t.hi;  // i.e. v1 = v1t >> 64

    // To ensure no overflow

    u128_t u1 = u;

    libdivide_u128_shift(&u1.hi, &u1.lo, -1);

    // Get quotient from divide unsigned insn.

    uint64_t rem_ignored;

    uint64_t q1 = libdivide_128_div_64_to_64(u1.hi, u1.lo, v1, &rem_ignored);

    // Undo normalization and division of u by 2.

    u128_t q0 = {0, q1};

    libdivide_u128_shift(&q0.hi, &q0.lo, n);

    libdivide_u128_shift(&q0.hi, &q0.lo, -63);

    // Make q0 correct or too small by 1

    // Equivalent to `if (q0 != 0) q0 = q0 - 1;`

    if (q0.hi != 0 || q0.lo != 0) {

        q0.hi -= (q0.lo == 0);  // borrow

        q0.lo -= 1;

    }

    // Now q0 is correct.

    // Compute q0 * v as q0v

    // = (q0.hi << 64 + q0.lo) * (v.hi << 64 + v.lo)

    // = (q0.hi * v.hi << 128) + (q0.hi * v.lo << 64) +

    //   (q0.lo * v.hi <<  64) + q0.lo * v.lo)

    // Each term is 128 bit

    // High half of full product (upper 128 bits!) are dropped

    u128_t q0v = {0, 0};

    q0v.hi = q0.hi * v.lo + q0.lo * v.hi + libdivide_mullhi_u64(q0.lo, v.lo);

    q0v.lo = q0.lo * v.lo;

    // Compute u - q0v as u_q0v

    // This is the remainder

    u128_t u_q0v = u;

    u_q0v.hi -= q0v.hi + (u.lo < q0v.lo);  // second term is borrow

    u_q0v.lo -= q0v.lo;

    // Check if u_q0v >= v

    // This checks if our remainder is larger than the divisor

    if ((u_q0v.hi > v.hi) || (u_q0v.hi == v.hi && u_q0v.lo >= v.lo)) {

        // Increment q0

        q0.lo += 1;

        q0.hi += (q0.lo == 0);  // carry

        // Subtract v from remainder

        u_q0v.hi -= v.hi + (u_q0v.lo < v.lo);

        u_q0v.lo -= v.lo;

    }

    *r_hi = u_q0v.hi;

    *r_lo = u_q0v.lo;

    LIBDIVIDE_ASSERT(q0.hi == 0);

    return q0.lo;

#endif

}

////////// UINT16

static LIBDIVIDE_INLINE struct libdivide_u16_t libdivide_internal_u16_gen(

    uint16_t d, int branchfree) {

    if (d == 0) {

        LIBDIVIDE_ERROR("divider must be != 0");

    }

    struct libdivide_u16_t result;

    uint8_t floor_log_2_d = (uint8_t)(15 - libdivide_count_leading_zeros16(d));

    // Power of 2

    if ((d & (d - 1)) == 0) {

        // We need to subtract 1 from the shift value in case of an unsigned

        // branchfree divider because there is a hardcoded right shift by 1

        // in its division algorithm. Because of this we also need to add back

        // 1 in its recovery algorithm.

        result.magic = 0;

        result.more = (uint8_t)(floor_log_2_d - (branchfree != 0));

    } else {

        uint8_t more;

        uint16_t rem, proposed_m;

        proposed_m = libdivide_32_div_16_to_16((uint16_t)1 << floor_log_2_d, 0, d, &rem);

        LIBDIVIDE_ASSERT(rem > 0 && rem < d);

        const uint16_t e = d - rem;

        // This power works if e < 2**floor_log_2_d.

        if (!branchfree && (e < ((uint16_t)1 << floor_log_2_d))) {

            // This power works

            more = floor_log_2_d;

        } else {

            // We have to use the general 17-bit algorithm.  We need to compute

            // (2**power) / d. However, we already have (2**(power-1))/d and

            // its remainder.  By doubling both, and then correcting the

            // remainder, we can compute the larger division.

            // don't care about overflow here - in fact, we expect it

            proposed_m += proposed_m;

            const uint16_t twice_rem = rem + rem;

            if (twice_rem >= d || twice_rem < rem) proposed_m += 1;

            more = floor_log_2_d | LIBDIVIDE_ADD_MARKER;

        }

        result.magic = 1 + proposed_m;

        result.more = more;

        // result.more's shift should in general be ceil_log_2_d. But if we

        // used the smaller power, we subtract one from the shift because we're

        // using the smaller power. If we're using the larger power, we

        // subtract one from the shift because it's taken care of by the add

        // indicator. So floor_log_2_d happens to be correct in both cases.

    }

    return result;

}

static LIBDIVIDE_INLINE struct libdivide_u16_t libdivide_u16_gen(uint16_t d) {

    return libdivide_internal_u16_gen(d, 0);

}

static LIBDIVIDE_INLINE struct libdivide_u16_branchfree_t libdivide_u16_branchfree_gen(uint16_t d) {

    if (d == 1) {

        LIBDIVIDE_ERROR("branchfree divider must be != 1");

    }

    struct libdivide_u16_t tmp = libdivide_internal_u16_gen(d, 1);

    struct libdivide_u16_branchfree_t ret = {

        tmp.magic, (uint8_t)(tmp.more & LIBDIVIDE_16_SHIFT_MASK)};

    return ret;

}

// The original libdivide_u16_do takes a const pointer. However, this cannot be used

// with a compile time constant libdivide_u16_t: it will generate a warning about

// taking the address of a temporary. Hence this overload.

static LIBDIVIDE_INLINE uint16_t libdivide_u16_do_raw(uint16_t numer, uint16_t magic, uint8_t more) {

    if (!magic) {

        return numer >> more;

    } else {

        uint16_t q = libdivide_mullhi_u16(numer, magic);

        if (more & LIBDIVIDE_ADD_MARKER) {

            uint16_t t = ((numer - q) >> 1) + q;

            return t >> (more & LIBDIVIDE_16_SHIFT_MASK);

        } else {

            // All upper bits are 0,

            // don't need to mask them off.

            return q >> more;

        }

    }

}

static LIBDIVIDE_INLINE uint16_t libdivide_u16_do(uint16_t numer, const struct libdivide_u16_t *denom) {

    return libdivide_u16_do_raw(numer, denom->magic, denom->more);

}

static LIBDIVIDE_INLINE uint16_t libdivide_u16_branchfree_do(

    uint16_t numer, const struct libdivide_u16_branchfree_t *denom) {

    uint16_t q = libdivide_mullhi_u16(numer, denom->magic);

    uint16_t t = ((numer - q) >> 1) + q;

    return t >> denom->more;

}

static LIBDIVIDE_INLINE uint16_t libdivide_u16_recover(const struct libdivide_u16_t *denom) {

    uint8_t more = denom->more;

    uint8_t shift = more & LIBDIVIDE_16_SHIFT_MASK;

    if (!denom->magic) {

        return (uint16_t)1 << shift;

    } else if (!(more & LIBDIVIDE_ADD_MARKER)) {

        // We compute q = n/d = n*m / 2^(16 + shift)

        // Therefore we have d = 2^(16 + shift) / m

        // We need to ceil it.

        // We know d is not a power of 2, so m is not a power of 2,

        // so we can just add 1 to the floor

        uint16_t hi_dividend = (uint16_t)1 << shift;

        uint16_t rem_ignored;

        return 1 + libdivide_32_div_16_to_16(hi_dividend, 0, denom->magic, &rem_ignored);

    } else {

        // Here we wish to compute d = 2^(16+shift+1)/(m+2^16).

        // Notice (m + 2^16) is a 17 bit number. Use 32 bit division for now

        // Also note that shift may be as high as 15, so shift + 1 will

        // overflow. So we have to compute it as 2^(16+shift)/(m+2^16), and

        // then double the quotient and remainder.

        uint32_t half_n = (uint32_t)1 << (16 + shift);

        uint32_t d = ((uint32_t)1 << 16) | denom->magic;

        // Note that the quotient is guaranteed <= 16 bits, but the remainder

        // may need 17!

        uint16_t half_q = (uint16_t)(half_n / d);

        uint32_t rem = half_n % d;

        // We computed 2^(16+shift)/(m+2^16)

        // Need to double it, and then add 1 to the quotient if doubling th

        // remainder would increase the quotient.

        // Note that rem<<1 cannot overflow, since rem < d and d is 17 bits

        uint16_t full_q = half_q + half_q + ((rem << 1) >= d);

        // We rounded down in gen (hence +1)

        return full_q + 1;

    }

}

static LIBDIVIDE_INLINE uint16_t libdivide_u16_branchfree_recover(const struct libdivide_u16_branchfree_t *denom) {

    uint8_t more = denom->more;

    uint8_t shift = more & LIBDIVIDE_16_SHIFT_MASK;

    if (!denom->magic) {

        return (uint16_t)1 << (shift + 1);

    } else {

        // Here we wish to compute d = 2^(16+shift+1)/(m+2^16).

        // Notice (m + 2^16) is a 17 bit number. Use 32 bit division for now

        // Also note that shift may be as high as 15, so shift + 1 will

        // overflow. So we have to compute it as 2^(16+shift)/(m+2^16), and

        // then double the quotient and remainder.

        uint32_t half_n = (uint32_t)1 << (16 + shift);

        uint32_t d = ((uint32_t)1 << 16) | denom->magic;

        // Note that the quotient is guaranteed <= 16 bits, but the remainder

        // may need 17!

        uint16_t half_q = (uint16_t)(half_n / d);

        uint32_t rem = half_n % d;

        // We computed 2^(16+shift)/(m+2^16)

        // Need to double it, and then add 1 to the quotient if doubling th

        // remainder would increase the quotient.

        // Note that rem<<1 cannot overflow, since rem < d and d is 33 bits

        uint16_t full_q = half_q + half_q + ((rem << 1) >= d);

        // We rounded down in gen (hence +1)

        return full_q + 1;

    }

}

////////// UINT32

static LIBDIVIDE_INLINE struct libdivide_u32_t libdivide_internal_u32_gen(

    uint32_t d, int branchfree) {

    if (d == 0) {

        LIBDIVIDE_ERROR("divider must be != 0");

    }

    struct libdivide_u32_t result;

    uint32_t floor_log_2_d = 31 - libdivide_count_leading_zeros32(d);

    // Power of 2

    if ((d & (d - 1)) == 0) {

        // We need to subtract 1 from the shift value in case of an unsigned

        // branchfree divider because there is a hardcoded right shift by 1

        // in its division algorithm. Because of this we also need to add back

        // 1 in its recovery algorithm.

        result.magic = 0;

        result.more = (uint8_t)(floor_log_2_d - (branchfree != 0));

    } else {

        uint8_t more;

        uint32_t rem, proposed_m;

        proposed_m = libdivide_64_div_32_to_32((uint32_t)1 << floor_log_2_d, 0, d, &rem);

        LIBDIVIDE_ASSERT(rem > 0 && rem < d);

        const uint32_t e = d - rem;

        // This power works if e < 2**floor_log_2_d.

        if (!branchfree && (e < ((uint32_t)1 << floor_log_2_d))) {

            // This power works

            more = (uint8_t)floor_log_2_d;

        } else {

            // We have to use the general 33-bit algorithm.  We need to compute

            // (2**power) / d. However, we already have (2**(power-1))/d and

            // its remainder.  By doubling both, and then correcting the

            // remainder, we can compute the larger division.

            // don't care about overflow here - in fact, we expect it

            proposed_m += proposed_m;

            const uint32_t twice_rem = rem + rem;

            if (twice_rem >= d || twice_rem < rem) proposed_m += 1;

            more = (uint8_t)(floor_log_2_d | LIBDIVIDE_ADD_MARKER);

        }

        result.magic = 1 + proposed_m;

        result.more = more;

        // result.more's shift should in general be ceil_log_2_d. But if we

        // used the smaller power, we subtract one from the shift because we're

        // using the smaller power. If we're using the larger power, we

        // subtract one from the shift because it's taken care of by the add

        // indicator. So floor_log_2_d happens to be correct in both cases.

    }

    return result;

}

static LIBDIVIDE_INLINE struct libdivide_u32_t libdivide_u32_gen(uint32_t d) {

    return libdivide_internal_u32_gen(d, 0);

}

static LIBDIVIDE_INLINE struct libdivide_u32_branchfree_t libdivide_u32_branchfree_gen(uint32_t d) {

    if (d == 1) {

        LIBDIVIDE_ERROR("branchfree divider must be != 1");

    }

    struct libdivide_u32_t tmp = libdivide_internal_u32_gen(d, 1);

    struct libdivide_u32_branchfree_t ret = {

        tmp.magic, (uint8_t)(tmp.more & LIBDIVIDE_32_SHIFT_MASK)};

    return ret;

}