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

#c #алгоритм #binary

На 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)\n",
       d, r,
       v / 10, v % 10);
  }


Похоже, мое изменение работает.

Вопрос, собственно в следующем -- может ли кто-нибудь в самом деле разбирающийся
в двоичной арифметике подтвердить, что divu64_10() работает верно?
    

Ответы

Ответ 1

Разберем как работает исходный код.

Он вычисляет 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.