Деление с остатком uint64_t на 10

На 32-bit микропроцессоре нужна функция (эффективный алгоритм) для получения частного и остатка от деления uint64_t на 10 (как всегда для печати).
Полазив по интернету, в http://www.hackersdelight.org/divcMore.pdf нашел лишь нечто похожее (для uint32_t) с вот таким кодом
unsigned divu10(unsigned n) { unsigned q, r; q = (n >> 1) + (n >> 2); q = q + (q >> 4); q = q + (q >> 8); q = q + (q >> 16); q = q >> 3; r = n - q*10; return q + ((r + 6) >> 4); // return q + (r > 9); }
Для своей цели изменил ее вот так (фактически, добавил оператор q = q + (q >> 32); :
uint64_t divu64_10 (uint64_t n, uint32_t *rem) { uint64_t q, r; q = (n >> 1) + (n >> 2); q = q + (q >> 4); q = q + (q >> 8); q = q + (q >> 16); q = q + (q >> 32); q = q >> 3; r = n - ((q << 3) + (q << 1)); *rem = r > 9 ? r - 10 : r;
// orig: return q + ((r + 6) >> 4); return q + (r > 9); }
Проверил, вводя разные числа вот в таком коде на нормальном 64-бит компе:
unsigned long long v;
while (scanf("%llu", &v) == 1) { unsigned long long d; uint32_t r;
d = divu64_10(v, &r); printf("divu10: %llu %u (C: %llu %llu)
", d, r, v / 10, v % 10); }
Похоже, мое изменение работает.
Вопрос, собственно в следующем -- может ли кто-нибудь в самом деле разбирающийся в двоичной арифметике подтвердить, что divu64_10() работает верно?

Ответ

Разберем как работает исходный код.
Он вычисляет q = n * 0x33333333 >> 33, что примерно равно n * 0.0999999999... Затем вычисляет ошибку n - q*10, и если она меньше 10, то всё правильно, если больше 10, то к результату добавляется 1.
Расмотрим значения ошибок.
10 * 0x33333333 / 2 ** 33 = 0b00.111111111111111111111111111111110 11 * 0x33333333 / 2 ** 33 = 0b01.000110011001100110011001100110001 19 * 0x33333333 / 2 ** 33 = 0b01.111001100110011001100110011001001 20 * 0x33333333 / 2 ** 33 = 0b01.111111111111111111111111111111100 21 * 0x33333333 / 2 ** 33 = 0b10.000110011001100110011001100101111 ... 65530 * 0x33333333 / 2 ** 33 = 6552.999998474261 = 0b1100110011000.111111111111111111100110011001110 ... 4294967290 * 0x33333333 / 2 ** 33 = 429496728.9 = 0b11001100110011001100110011000.111001100110011001100110011001110
Видно что для n = 10*d + r, ошибка равна 10 если r=0.
Можно посчитать, что если бы вместо сдвига было бы деление, то ошибка была бы равна
n - n * 10 * 0x33333333 / 2**33 = n / 2 ** 32
Так как n < 2 ** 32, то ошибка результата не превышает единицу.

Для 64 разрядов всё так же, только считаем n * 0x3333333333333333 >> 65
Добавление q += q >> 32 эквивалентно умножению на (2**32+1)/2**32, при этом (2**32+1)*0x33333333 = 0x3333333333333333