uint64_muldiv.h

// Through https://github.com/nordlow/justcxx/blob/master/muldiv64.h
 
#pragma once
 
#ifdef __cplusplus
extern "C" {
#endif

   static inline ::u32 uint64_msb(uint64_t value)
   {
      const ::i32 MAX_LOGLOG = 6;
      const uint64_t BIT_LL[MAX_LOGLOG] = { 0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000, 0xFFFFFFFF00000000LL };
      const ::u32 EXP_LL[MAX_LOGLOG] = { 1, 2, 4, 8, 16, 32 };
      ::u32 r = 0;
      for (::i32 i = MAX_LOGLOG - 1; i >= 0; i--) {
         if (value & BIT_LL[i]) {
            value >>= EXP_LL[i];
            r |= EXP_LL[i];
         }
      }
      return r;
   }

   static inline uint64_t uint64_muldiv(uint64_t a, uint32_t b, uint64_t c)
   {
      const uint32_t aH = ((uint32_t*)&a)[1];
      const uint32_t aL = ((uint32_t*)&a)[0];
 
      // lower 32-bit portions yield 64-bit product that can be divided
      // directly giving 64-bits of result,
      uint64_t res = (aL * b) / c;
 
      // upper 32-bit have to be shifted, so calculate modulus 2^32
      const uint64_t prod_h = aH * b;
      const uint64_t div_h = prod_h / c; // division main
      const uint64_t mod_h = prod_h - c * div_h; // modulus
 
      // upper bits
      res += div_h << 32;
      if (mod_h == 0) {
         return res;
      }
 
      // remainder of division. if msb modulus < 32 we can be quick about
      // it
      if ((mod_h >> 32) == 0) {
         res += (mod_h << 32) / c;
         return res;
      }
 
      // if we reach this point we have full 64-bit values, i.e. a 96-bit
      // dividend. calculate an approximate result by shifting according
      // to msb set
      ::i32 msb_nominator = uint64_msb(mod_h) + 32;
      ::i32 msb_c = uint64_msb(c);
      ::i32 msb = maximum(msb_nominator, msb_c);
      ::i32 shift = msb - 63;
      res += (mod_h << (32 - shift)) / (c >> shift);
 
      return res;
   }
 
#ifdef __cplusplus
}
#endif