Сформировать таблицу поиска, чтобы разделить на 10 битное целое посредством умножения

Question

У меня есть целое число, это кратно 256 между 256 и 131072 , и я хочу, чтобы разделить его на целое число между 1 и 1024

Это hotspot во внутреннем цикле моего кода и комментариях, это значительно ускоряет мое приложение.

Могу ли я сделать таблицу поиска 1024 размера, которая поможет преобразовать деление в «умножение плюс сдвиг» за меньшее время, чем фактическое деление на процессор x86_64?

Может кто-нибудь помочь придумать код для генерации таблицы поиска, которая позволяет эффективно делить один из этих 1024 возможных делителей?

Мне очень хотелось бы увидеть шаблон метапрограммирования способа создания соответствующей таблицы как constexpr.

Answer 1

Поскольку выбор дивидендов и дивизоров очень ограничен, можно использовать более простой подход, чем использовать в paper Торбьёрном Гранлундом и Питером Монтгомери, которые вы обнаружили. Я не разбираюсь в мета-программировании шаблонов C++, но я могу продемонстрировать подход к созданию и использованию таблицы правильно масштабированных взаимных ссылок.

Прежде всего отметим, что дивиденды все умножены на 256, поэтому их можно свести к 1 ... 0x200 простой предварительной сменой справа от 8 бит. Поскольку мы не хотим переполнять 32-битное целое число без знака во время умножения уменьшенного дивиденда на масштабное обратное, обратное идеально масштабируется в диапазоне 0x00200000 < rcp <= 0x00400000.

Если быстрый подсчет ведущего-ноль инструкция доступна, она может быть использована для масштабирования до обратных в этот диапазон в течение таблицы предварительного вычисления, на основе логарифма базового-2 делителя, то при беге используют тот же динамически рассчитанный масштабный коэффициент для шкалы down продукт уменьшенного дивиденда и масштабированный обратный на тот же коэффициент. При масштабировании обратного нам нужно округлить до результат к следующему целому числу, чтобы компенсировать усечение природы масштабирования с помощью сдвига вправо. Вариант 0 в коде ниже использует этот подход.

Что нам делать, если нет инструкции по быстрому счету-ведущему нулю? Нам нужно сохранить масштабированное взаимное число с достаточным количеством бит для поддержания точности вычисления. Оказывается, нам повезло из-за ограниченного ограничения на диапазон дивизоров и может обойтись только с двумя разными масштабными коэффициентами, которые можно легко вычислить из делителя во время выполнения: один фактор для делителей < = 32, другой для делителей в (32, 1024). Это используется в варианте 1 кода, где два масштабных коэффициента составлены до 2 и 2 .

Наконец, мы можем не захотеть вычислить масштабный коэффициент «на лету», а сохранить его вместе с масштабированным обратным, используя самые значащие биты каждой записи таблицы для их хранения, в то время как менее значимые биты используются для наоборот. Одним из недостатков является необходимость дополнительных операций по извлечению масштабированного взаимного и масштабного коэффициента из записи таблицы, что делает этот подход менее подходящим, чем два других. Это показано в кодовом варианте 2 ниже.

#include <stdio.h> 
#include <stdlib.h> 
#include <stdint.h> 
#include <math.h> 

#define VARIANT 0 // variants available: 0 ... 2 

extern int __clz (uint32_t x); // intrinsic for access to CLZ instruction 

uint32_t rcp [1025]; // table of scaled reciprocals (& possibly scale factors) 

#define PRE_SCALE  (8) // downscaling for dividends, since multiples of 256 
#define POST_SCALE_1 (14) // upscale factor #1 for reciprocals 
#define POST_SCALE_2 (19) // upscale factor #2 for reciprocals 
#define RECIP_BITS (24) // bits used in table entry for scaled reciprocal 

// helper function: logarithm base-2 of an 32-bit integer  
int ilog2 (uint32_t x) 
{ 
    return 31 - __clz (x); 
} 

// division for dividends n*256, n in [1,512], divisor in [1,1024]  
uint32_t fast_div (uint32_t dividend, uint32_t divisor) 
{ 
#if VARIANT == 0 
    uint32_t scale = POST_SCALE_1 + ilog2 (divisor); 
    return ((dividend >> PRE_SCALE) * rcp[divisor]) >> scale; 
#elif VARIANT == 1 
    uint32_t scale = (divisor > 0x20) ? POST_SCALE_2 : POST_SCALE_1; 
    return ((dividend >> PRE_SCALE) * rcp[divisor]) >> scale; 
#elif VARIANT == 2 
    uint32_t scale = rcp[divisor] >> RECIP_BITS; 
    return ((dividend >> PRE_SCALE) * (rcp[divisor] & ((1 << RECIP_BITS) - 1))) >> scale; 
#else 
#error non-existing VARIANT 
#endif 
} 

int main (void) 
{ 
    uint32_t dividend, divisor, res, ref; 
    int i; 

    // precompute table of recprocals 
    for (i = 1; i < 1025; i++) { 
#if VARIANT == 0 
     uint32_t scale = POST_SCALE_1 + ilog2 (i); 
     rcp[i] = ((uint32_t)(pow (2.0, PRE_SCALE + scale)/i + 0.99999)); 
#elif VARIANT == 1 
     uint32_t scale = (i > 0x20) ? POST_SCALE_2 : POST_SCALE_1; 
     rcp[i] = ((uint32_t)(pow (2.0, PRE_SCALE + scale)/i + 0.99999)); 
#elif VARIANT == 2 
     uint32_t scale = (i > 0x20) ? POST_SCALE_2 : POST_SCALE_1; 
     rcp[i] = ((uint32_t)(pow (2.0, PRE_SCALE + scale)/i + 0.99999) + 
        (scale << RECIP_BITS)); 
#else 
#error non-existing VARIANT 
#endif 
    } 

    // test all supported dividens and divisors exhaustively 
    divisor = 1; 
    while (divisor <= 1024) { 
     dividend = 256; 
     while (dividend <= 131072) { 
      res = fast_div (dividend, divisor); 
      ref = dividend/divisor; 
      if (res != ref) { 
       printf ("n=%08x d=%08x res=%08x ref=%08x rcp=%08x\n", 
         dividend, divisor, res, ref, rcp[divisor]); 
       return EXIT_FAILURE; 
      } 
      dividend += 256; 
     } 
     divisor++; 
    } 
    printf ("division test passed\n"); 
    return EXIT_SUCCESS; 
}