#include <stdint.h>

static inline unsigned int av_log2(unsigned int x) {
    unsigned int ret;
    __asm__ volatile("or $1, %1 \n\t"
             "bsr %1, %0 \n\t"
             :"=a"(ret) : "D"(x));
    return ret;
}

static const unsigned char  binvert_limb_table[128] = {
  0x01, 0xAB, 0xCD, 0xB7, 0x39, 0xA3, 0xC5, 0xEF,
  0xF1, 0x1B, 0x3D, 0xA7, 0x29, 0x13, 0x35, 0xDF,
  0xE1, 0x8B, 0xAD, 0x97, 0x19, 0x83, 0xA5, 0xCF,
  0xD1, 0xFB, 0x1D, 0x87, 0x09, 0xF3, 0x15, 0xBF,
  0xC1, 0x6B, 0x8D, 0x77, 0xF9, 0x63, 0x85, 0xAF,
  0xB1, 0xDB, 0xFD, 0x67, 0xE9, 0xD3, 0xF5, 0x9F,
  0xA1, 0x4B, 0x6D, 0x57, 0xD9, 0x43, 0x65, 0x8F,
  0x91, 0xBB, 0xDD, 0x47, 0xC9, 0xB3, 0xD5, 0x7F,
  0x81, 0x2B, 0x4D, 0x37, 0xB9, 0x23, 0x45, 0x6F,
  0x71, 0x9B, 0xBD, 0x27, 0xA9, 0x93, 0xB5, 0x5F,
  0x61, 0x0B, 0x2D, 0x17, 0x99, 0x03, 0x25, 0x4F,
  0x51, 0x7B, 0x9D, 0x07, 0x89, 0x73, 0x95, 0x3F,
  0x41, 0xEB, 0x0D, 0xF7, 0x79, 0xE3, 0x05, 0x2F,
  0x31, 0x5B, 0x7D, 0xE7, 0x69, 0x53, 0x75, 0x1F,
  0x21, 0xCB, 0xED, 0xD7, 0x59, 0xC3, 0xE5, 0x0F,
  0x11, 0x3B, 0x5D, 0xC7, 0x49, 0x33, 0x55, 0xFF
};

/*
 * Exact division. Returns the value of n/k, *provided that*
 * n%k is zero, i.e. n exactly divides k.
 * Yes, expanding this to 64 bits would just require replacing
 * "unsigned int" w/"uint64_t" and adding one more Newton step.
 * Yes, it's really that elegant.
 */
unsigned int divexact32(unsigned int n, unsigned int k)
{
    unsigned int __inv, s = av_log2(k&-k);
    n >>= s;
    k >>= s;
    __inv = binvert_limb_table[(k>>1) & 0x7F]; /*  8 bits */
    __inv = 2 * __inv - __inv * __inv * k;     /* 16 bits */
    __inv = 2 * __inv - __inv * __inv * k;     /* 32 bits */
    return __inv*n;
}

