Нахождение обратного числа

Одна из способностей, которая проявилась у меня с детства...

Пишу программы под всякую странную аппаратуру. В странной аппаратуре плохо считается операция деления. Например, берем самую совершенную, казалось бы, архитектуру Intel, в ней операция деления выполняется [AF1] очень медленно. Все арифметические инструкции, включая умножение, запускаются по 4 шт за такт, а вот деление только одно и продолжается целых 6 тактов.
Core i7 Coffee Lake: Умеет паралелить множество инструкций, а операция умножения делается только на одном из 8 конвейеров, имеет задержку 4 такта, на прохождение (Latency) и может подряд выполняться с производительностью 1 такт на умножение (Throughput). Я не беру векторные инструкции, беру простую задачу поделить одно число на другое. Деление вызывает настоящий ступор в системе, 26 тактов на загрузку инструкции (Latency) и 6 тактов зависание на конвейере.
Берем менее совершенную архитектуру, ARM Cortex, в документации указано число циклов на операцию умножения и деления, 2 и 12 соответственно. В 6 раз операция деления медленнее операции умножения. А еще есть контроллеры, которые вообще не имеют аппаратного деления. Эмуляция деления это сотни тактов процессора.

Вот, к примеру, мне понадобилось штамп времени разложить на число, месяц и год. Потом, часы, минуты и секунды. Начинаем считать: остаток от деления штампа времени на (24*60*60) - это секунды от начала дня. Потом делим на (60*60), потом на 60 ... сплошные операции деления на константу. Или берем задачу форматирования числа, надо на каждом десятичном разряде выполнить деление на 10 с остатком. Вот тут и приходит здравая мысль - это оптимизация, которая может выполняться на этапе компиляции. Целочисленное деление на константу можно заменить на умножение со сдвигом. Фактически мы заменяем деление на умножение на дробь. Не все числа "рациональные", не все обратные числа можно представить в виде дроби с достаточной точностью. Приближение первое: мы хотим представить константу в виде дроби 1/B ≈ C/2^n. Первое что пришло в голову n=32, брать только старшие 32 бита от целочисленного умножения с разрядность 32х32бит, но это так не работает. Надо умножать на самое большое число, для которого не возникает переполнение, использовать максимально возможное n, чтобы константа С годилась для всех целых чисел. Число С должно содержать 1(единицу) в старшем разряде, не меньше 0x80000000UL.
Оптимизация деления на константу на языке Си запишется так:
uint32_t Q = (A*C)>>n,
где С и n -- константы.
Для нахождения констант предлагаю алгоритм:

uint32_t div_c(uint32_t b) {// находим число С
    uint32_t C;
    int k = __builtin_ctz(b);// число ноликов сзади
    b>>=k;
    n = 63 - __builtin_clz(b);
    C = ((1ULL<<(n))/b)+1;// сдвигаем на число ноликов спереди
    n+=k;
    return C;
}

Проверить справедливость замены можно методом перебора чисел, это занимает не много времени.

uint32_t verify(uint32_t B, uint32_t C, int n) {
    uint32_t A;
    for (A=(~0UL); A!=0; --A){// все числа 32 бит
        uint32_t Q = ((uint64_t)A*C)>>32>>(n-32);
        if (Q != (A/B)) {
            break;
        }
    }
    return A;
}

Примеры таких констант Q=A/B = (A*C)>>n:
 

Q = (A* 0xCCCСCСCDULL)>>35; // Q=A/10 кратно /5 работает для всех A
Q = (A* 0xCCCСCСCDULL)>>34; // Q=A/5 работает для всех A
Q = (A* 0xAAAAAAABULL)>>33; // Q=A/3 работает для всех A

B=  3: 0xAAAAAAAB n=33
B=  5: 0xCCCСCСCD n=34
B=  7: 0x92492493 n=34 Q=A/7 для 31 бит.
B= 60: 0x88888889 n=37
B= 10: 0xCCCСCСCD n=35
B=100: 0xA3D70A3E n=38
B=1000: 0x83126E98 n=41
B=3600: 0x91A2B3C5 n=43
B=86400: 0xC22E4507 n=48
B=1000000: 0x8637BD06 n=51
B=365: 0xB38CF9B1 n=40 Q=A/365 для 31 бит.
B=153: 0xD62B80D7 n=39
B=255: 0x80808081 n=39

Некоторые константы нельзя представить в виде 32 битного умножения, точности не хватает. Для 31 битных целых чисел всегда хватает точности.

-- Магия чисел!

[AF1] https://agner.org/optimize/instruction_tables.pdf