Быстрое сложение BCD

Этот вопрос был вдохновлен вопросом, который я недавно встретил на Stackoverflow. Это напомнило мне, что на ранних этапах разработки x86-64 ISA я написал 32-битный код x86 для сложения BCD без использования инструкций поддержки BCD, доступных в x86 ISA. Это относится к тому, что, судя по комментариям, на некоторых платформах может быть известно как упакованный BCD, то есть одна десятичная цифра, хранящаяся в тетраде битов, обычно называемая полубайтом.

Хотя было совершенно очевидно, что удаление инструкций DAA и DAS (коррекция десятичной дроби после сложения/вычитания) безопасно с точки зрения производительности, лучше всего было иметь под рукой конкретные доказательства того, как можно реализовать, скажем, сложение BCD без этих инструкций. инструкции на современном процессоре Athlon. Мой код того времени, откопанный в глубинах Интернета, слегка скорректированный и исправленный, приведен ниже. Обратите внимание, что RCR тогда не был таким большим убийцей производительности, как сейчас. Я вновь рассмотрел этот вопрос с сегодняшней точки зрения, решив остаться с обработкой 32-битными фрагментами для упрощения представления данных на различных архитектурах. Как распространить обработку 64-битных фрагментов на 64-битные ISA, очевидно.

Концептуально сложение BCD является простым. Чтобы облегчить распространение переноса между полубайтами, к каждому полубайту двоичной суммы двух исходных операндов временно добавляют 6. Для тех полубайтов, которые не производят перенос, либо вычитают 6 из полубайта (метод отмены), либо повторяют исходное двоичное сложение, но без добавления 6 к этому полубайту (метод повтора). Работая над своим исходным кодом в 1999 году, я обнаружил, что использование метода повторного выполнения приводит к несколько более высокой производительности, возможно, из-за более короткой цепочки зависимостей. Основная задача при сложении BCD заключается в определении того, где происходят переносы между полубайтами. Мой подход заключался в том, чтобы проверить младший бит следующего старшего полубайта: если его значение отличается от предварительного значения бита суммы (исключающее ИЛИ двух исходных младших битов), то должен был произойти перенос, поэтому следующий младший полубайт создал ошибку. выполнять. Несколько неприятный частный случай — обработка выноса из наиболее значимого полубайта.

Когда Кнут опубликовал том 4a TAOCP в 2011 году, я обнаружил, что он содержит алгоритм сложения BCD с помощью метода отмены и обнаружения переносов между полубайтами путем проверки старших битов соответствующих полубайтов в исходных операндах и скорректированной суммы, чтобы проверить, есть ли вынос из этого кусочка. Это делает код чрезвычайно элегантным, когда последняя функциональность выражается через вычисление большинства из 3, которое Кнут предпочитает называть более общим именем медиана.

Экспериментально я обнаружил, что метод повтора обычно превосходит метод отмены в широком диапазоне архитектур: x86, ARM, RISC-V, Sparc, Power. Что касается обнаружения переносов полубайтов, метод LSB и метод MSB обычно кажутся равными с точки зрения производительности, но метод Кнута MSB превосходит архитектуры, которые могут более эффективно реализовать медианную операцию за счет использования преимуществ дополнительных логических операций. (например, ORN, ANDN, XNOR, NOR, NAND) за пределами минимального набора логических операций (OR, AND, XOR, NOT). Это особенно актуально для современных графических процессоров NVIDIA, которые предоставляют инструкцию LOP3, способную вычислить любую логическую операцию с тремя входами. Это сопоставляет медианную операцию с любым из входных данных, возможно, отрицательным, в один экземпляр LOP3.

Хотя я считаю свои навыки работы с битами довольно солидными, прошлый опыт показал, что я часто нахожусь на отметке 85% по сравнению с лучшим решением, предложенным в ответах здесь. В частности, я хочу выяснить: существуют ли более эффективные способы решения проблемы переноса полубайтов? Использование простой ISA, такой как 32-битный RISC-V, может дать лучшее представление об этом. Любое превосходное решение, применимое к этой ISA, вряд ли повредит ISA с более широким набором доступных инструкций. Помимо этого, приветствуется любая умная идея, которая может уменьшить накладные расходы на обнаружение переноса полубайтов, о чем свидетельствует уменьшенное количество инструкций или укороченная цепочка зависимостей.

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

#define NUM_TESTS    (100000000)

#define TEST_PORTABLE_BCD64 (1) // 0: x86-specific test of 18-digit BCD arithm.
                                // 1: portable test of 16-digit BCD arithmetic
#define PORTABLE     (1)  // 0: x86-specific inline assembly
#define NATIVE_64BIT (0)  // applies only if PORTABLE==1
#define USE_KNUTH    (1)  // applies only if PORTABLE==1 && NATIVE_64BIT==0

#if PORTABLE && !NATIVE_64BIT
/* majority-of-3 operation that Knuth refers to more generically as median */
uint32_t median32 (uint32_t x, uint32_t y, uint32_t z)
{
    return (x | (y & z)) & (y | z);
}

/* 8-digit BCD addition with carry-out */
uint32_t BCD32_ADDcc (uint32_t x, uint32_t y, uint32_t *cy)
{
#if USE_KNUTH
    // derived from Knuth, TAOCP 4a (2011), answer to exercise 100 section 7.1.3
    uint32_t z, u, t;
    z = y + 0x66666666;
    u = x + z;
    t = median32 (x, z, ~u) & 0x88888888;
    *cy = t >> 31;
    return (x + y) + (t - (t >> 2));
#else // ! USE_KNUTH
    const uint32_t sixes = 0x66666666;
    const uint32_t eights = 0x88888888;
    uint32_t xe, s, u, carry_out;
    xe = x + sixes;         // add "nibble offset"; cannot overflow
    s = xe ^ y;             // preliminary sum bits
    u = xe + y;
    carry_out = u < xe;     // carry-out from most significant nibble
    s = (carry_out << 31) | ((s ^ u) >> 1); // sum bits with carry-in
    s = s & eights;         // isolate nibble lsbs with carry-in
    *cy = carry_out;        // record carry-out from most significant nibble
    u = (x + y) + (s - (s >> 2)); // add "nibble offset" to nibbles that need it
    return u;
#endif // USE_KNUTH
}

/* 8-digit BCD addition with carry-in */
uint32_t BCD32_ADDC (uint32_t x, uint32_t y, uint32_t cy)
{
#if USE_KNUTH
    // derived from Knuth, TAOCP 4a (2011), answer to exercise 100 section 7.1.3
    uint32_t z, u, t;
    x = x + cy;
    z = y + 0x66666666;
    u = x + z;
    t = median32 (x, z, ~u) & 0x88888888;
    return (x + y) + (t - (t >> 2));
#else // !USE_KNUTH
    const uint32_t sixes = 0x66666666;
    const uint32_t eights = 0x88888888;
    uint32_t xe, s, u, carry_out;
    y = y + cy;             // incorporate carry-in; cannot overflow
    xe = x + sixes;         // add "nibble offset"; cannot overflow
    s = xe ^ y;             // preliminary sum bits
    u = xe + y;              
    carry_out = u < xe;     // carry-out from most significant nibble
    s = (carry_out << 31) | ((s ^ u) >> 1); // sum bits with carry-in
    s = s & eights;         // isolate nibble lsbs with carry-in
    u = (x + y) + (s - (s >> 2)); // add "nibble offset" to nibbles that need it
    return u;
#endif // USE_KNUTH
}

/* 8-digit BCD addition with carry-in and carry-out */
uint32_t BCD32_ADDCcc (uint32_t x, uint32_t y, uint32_t *cy)
{
#if USE_KNUTH
    // derived from Knuth, TAOCP 4a (2011), answer to exercise 100 section 7.1.3
    uint32_t z, u, t;
    x = x + (*cy);
    z = y + 0x66666666;
    u = x + z;
    t = median32 (x, z, ~u) & 0x88888888;
    *cy = t >> 31;
    return (x + y) + (t - (t >> 2));
#else // !USE_KNUTH
    const uint32_t sixes = 0x66666666;
    const uint32_t eights = 0x88888888;
    uint32_t xe, s, u, carry_out;
    y = y + (*cy);          // incorporate carry-in; cannot overflow
    xe = x + sixes;         // add "nibble offset"; cannot overflow
    s = xe ^ y;             // preliminary sum bits
    u = xe + y;              
    carry_out = u < xe;     // carry-out from most significant nibble
    s = (carry_out << 31) | ((s ^ u) >> 1); // sum bits with carry-in
    s = s & eights;         // isolate nibble lsbs with carry-in
    *cy = carry_out;        // record carry-out from most significant nibble
    u = (x + y) + (s - (s >> 2)); // add "nibble offset" to nibbles that need it
    return u;
#endif // USE_KNUTH
}
#endif // PORTABLE && !NATIVE_64BIT

#if PORTABLE
#if NATIVE_64BIT
/* majority-of-3 operation that Knuth refers to more generically as median */
uint64_t median64 (uint64_t x, uint64_t y, uint64_t z)
{
    return (x | (y & z)) & (y | z);
}

/* 16-digit BCD addition */
uint64_t BCD64_ADD (uint64_t x, uint64_t y)
{
    // derived from Knuth, TAOCP 4a (2011), answer to exercise 100 section 7.1.3
    uint64_t z, u, t;
    z = y + 0x6666666666666666ULL;
    u = x + z;
    t = median64 (x, z, ~u) & 0x8888888888888888ULL;
    return (x + y) + (t - (t >> 2));
}
#else // !NATIVE_64BIT
/* 16-digit BCD addition */
uint64_t BCD64_ADD (uint64_t x, uint64_t y)
{
    uint32_t x_lo = (uint32_t)x;
    uint32_t x_hi = (uint32_t)(x >> 32);
    uint32_t y_lo = (uint32_t)y;
    uint32_t y_hi = (uint32_t)(y >> 32);
    uint32_t cy;
    uint32_t r_lo = BCD32_ADDcc (x_lo, y_lo, &cy);
    uint32_t r_hi = BCD32_ADDC  (x_hi, y_hi,  cy);
    return ((uint64_t)r_hi << 32) | r_lo;
}
#endif // NATIVE_64BIT
#else // !PORTABLE
/*
  16-digit BCD addition; N. Juffa, about 1999. Code retrieved on 3/24/2024 from 
  https://web.archive.org/web/20001211131100/http://www.azillionmonkeys.com/qed/asmexample.html
  Bug fixed and adjusted for Intel compiler's inline asm (may need -fasm-blocks)
*/
uint64_t BCD64_ADD (uint64_t x, uint64_t y)
{
    uint64_t r;

    __asm mov   eax, dword ptr [x]   ; // x (lo)
    __asm mov   ebx, dword ptr [y]   ; // y (lo)
    __asm mov   edx, dword ptr [x+4] ; // x (hi)
    __asm mov   ecx, dword ptr [y+4] ; // y (hi)
    // x in EDX:EAX, y in EDX:EBX. Add least significant 8 BCD digits
    __asm mov   esi, eax             ; // x
    __asm lea   edi, [eax+66666666h] ; // x + 0x66666666
    __asm xor   esi, ebx             ; // x ^ y
    __asm add   eax, ebx             ; // x + y
    __asm shr   esi, 1               ; // t1 = (x ^ y) >> 1
    __asm add   edi, ebx             ; // x + y + 0x66666666
    __asm sbb   ebx, ebx             ; // capture carry: carry ? (-1) : 0
    __asm rcr   edi, 1               ; // t2 = (x + y + 0x66666666) >> 1
    __asm xor   esi, edi             ; // t1 ^ t2
    __asm and   esi, 88888888h       ; // t3 = (t1 ^ t2) & 0x88888888
    __asm add   eax, esi             ; // x + y + t3
    __asm shr   esi, 2               ; // t3 >> 2
    __asm sub   eax, esi             ; // x + y + t3 - (t3 >> 2)
    // Add most significant 8 BCD digits
    __asm sub   ecx, ebx             ; // y - carry // propagate carry
    __asm mov   esi, edx             ; // x
    __asm lea   edi, [edx+66666666h] ; // x + 0x66666666
    __asm xor   esi, ecx             ; // x ^ y
    __asm add   edx, ecx             ; // x + y
    __asm shr   esi, 1               ; // t1 = (x ^ y) >> 1
    __asm add   edi, ecx             ; // x + y + 0x66666666
    // __asm sbb   ebx, ebx             ; capture carry: carry ? (-1) : 0
    __asm rcr   edi, 1               ; // t2 = (x + y + 0x66666666) >> 1
    __asm xor   esi, edi             ; // t1 ^ t2
    __asm and   esi, 88888888h       ; // t3 = (t1 ^ t2) & 0x88888888
    __asm add   edx, esi             ; // x + y + t3
    __asm shr   esi, 2               ; // t3 >> 2
    __asm sub   edx, esi             ; // x + y + t3 - (t3 >> 2)
    // Now EDX:EAX = x + y
    __asm mov   dword ptr [r], eax   ; // save result (lo)
    __asm mov   dword ptr [r+4], edx ; // save result (hi)

    return r;
}
#endif // PORTABLE

/* convert binary operand x in [0, 9999] to BCD */
uint32_t cvt_bin2bcd_4dig (uint32_t x)
{
    uint32_t num, res = 0;
    /* convert to 21:11 fixed point */
    num = ((0x8312 * x) + 0x6000) >> 14; // (2**11 / 10**3) * 2**14
    /* pull out decimal digits one at a time */
    res =              (num >> 11);
    num = (num & ((1 << 11) - 1)) * 5;   // 22:10 fixed point
    res = (res << 4) | (num >> 10);
    num = (num & ((1 << 10) - 1)) * 5;   // 23:9 fixed point
    res = (res << 4) | (num >>  9);
    num = (num & ((1 <<  9) - 1)) * 5;   // 24:8 fixed point
    res = (res << 4) | (num >>  8);
    return res;
}

/* convert binary operand x in [0, 99999999] to BCD */
uint32_t cvt_bin2bcd_8dig (uint32_t x)
{
    uint32_t hi = x / 10000;
    uint32_t lo = x - hi * 10000;
    return (cvt_bin2bcd_4dig (hi) << 16) | cvt_bin2bcd_4dig (lo);
}

/* convert binary operand a in [0, 9999999999999999] to BCD */
uint64_t uint64_to_bcd64 (uint64_t x)
{
    uint32_t hi = (uint32_t)(x / 100000000ULL);
    uint32_t lo = (uint32_t)(x - hi * 100000000ULL);
    return (((uint64_t)cvt_bin2bcd_8dig (hi)) << 32) | cvt_bin2bcd_8dig (lo);   
}

/* return most significant 32 bits of the full product of two 32-bit operands */
uint32_t umul32hi (uint32_t a, uint32_t b)
{
    uint64_t prod = (uint64_t)a * b;
    return (uint32_t)(prod >> 32);
}

/* convert BCD operand x in [0, 99999999] to binary
   based on Knuth TAOCP 2, answer on p. 638 to exercise 19 in section 4.4
*/
uint32_t cvt_bcd2bin_8dig (uint32_t x)
{
    const uint32_t m1 = 0xf0f0f0f0;
    const uint32_t m2 = 0xff00ff00;
    const uint32_t m3 = 0xffff0000;
    const uint32_t c1 = 0x60000000; // 0.375000000000
    const uint32_t c2 = 0x9c000000; // 0.609375000000
    const uint32_t c3 = 0xd8f00000; // 0.847412109375
    x = x - umul32hi (c1, x & m1);
    x = x - umul32hi (c2, x & m2);
    x = x - umul32hi (c3, x & m3);
    return x;
}

/* convert BCD operand x in [0, 9999999999999999] to binary */
uint64_t bcd64_to_uint64 (uint64_t a)
{
    uint32_t hi = (uint32_t)(a >> 32);
    uint32_t lo = (uint32_t)a;
    return 100000000ULL * cvt_bcd2bin_8dig (hi) + cvt_bcd2bin_8dig (lo);
}

/*
  https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
  From: geo <[email protected]>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, \
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

#if TEST_PORTABLE_BCD64
int main (void)
{
    uint64_t bx, by, x, y, ref, res, i;

    printf ("Testing 16-digit BCD arithmetic\n");
    for (i = 0; i < NUM_TESTS; i++) {
        x = KISS64 % 10000000000000000ULL;
        y = KISS64 % 10000000000000000ULL;
        bx = uint64_to_bcd64 (x);
        by = uint64_to_bcd64 (y);
        res = BCD64_ADD (bx, by);
        ref = uint64_to_bcd64 ((x + y) % 10000000000000000ULL);
        if (res != ref) {
            printf ("!!!! x=%016llx  y=%016llx  bx=%016llx  by=%016llx  res=%016llx  ref=%016llx\n", 
                    x, y, bx, by, res, ref);
            return EXIT_FAILURE;
        }
    }
    printf ("test PASSED\n");
    return EXIT_SUCCESS;
}

#else // TEST_PORTABLE_BCD64

typedef struct ten_byte {
    uint64_t low64;
    uint16_t high16;
} ten_byte;

/* 18-digit BCD addition modulo 10**18 using the x87 FPU */
ten_byte bcd18dig_add_ref (ten_byte a, ten_byte b)
{
    const double exp10_18 = 1e18;
    const uint16_t cw_extended_chop = 0xf7f;
    ten_byte aa, bb, r;
    uint16_t cw_original;

    aa = a; // workaround: bad code generated when data not local to function
    bb = b;
    __asm fstcw  [cw_original]        ;// save current CW
    __asm fldcw  [cw_extended_chop]   ;// RC=truncate PC=extended
    __asm fld    qword ptr [exp10_18] ;// 10**18
    __asm fbld   tbyte ptr [aa]       ;// a, 10**8
    __asm fbld   tbyte ptr [bb]       ;// b, a, 10**18
    __asm fadd   st, st(1)            ;// a+b, a, 10**18
    __asm fst    st(1)                ;// a+b, a+b, 10**18
    __asm fdiv   st, st(2)            ;// (a+b)/10**18, a+b, 10**18
    __asm frndint                     ;// trunc((a+b)/10**18), a+b, 10**18
    __asm fmulp  st(2), st            ;// a+b, trunc((a+b)/10**18)*10**18
    __asm fsubrp st(1),st             ;// a+b - trunc((a+b)/10**18)*10**18;
    __asm fbstp  tbyte ptr [r]        ;// <empty>
    __asm fldcw  [cw_original]        ;// restore original CW

    return r;
}

/* 18-digit BCD addition modulo 10**18, processed in 8-digit chunks */
ten_byte bcd18dig_add (ten_byte x, ten_byte y)
{
    ten_byte r;

    uint32_t x_lo = (uint32_t)(x.low64);
    uint32_t x_md = (uint32_t)(x.low64 >> 32);
    uint32_t x_hi = (uint32_t)(x.high16);
    uint32_t y_lo = (uint32_t)(y.low64);
    uint32_t y_md = (uint32_t)(y.low64 >> 32);
    uint32_t y_hi = (uint32_t)(y.high16);
    uint32_t cy;
    uint32_t r_lo = BCD32_ADDcc (x_lo, y_lo, &cy);
    uint32_t r_md = BCD32_ADDCcc(x_md, y_md, &cy);
    uint32_t r_hi = BCD32_ADDC  (x_hi, y_hi,  cy);
    r_hi &= 0xff;
    r.low64 = ((uint64_t)r_md << 32) | r_lo;
    r.high16 = (uint16_t)r_hi;

    return r;
}

int main (void)
{
    uint64_t x, y;
    ten_byte bcd_augend, bcd_addend, bcd_sum_ref, bcd_sum_res;
    
    printf ("Testing 18-digit BCD arithmetic against x87 FPU\n");
    for (int i = 0; i < NUM_TESTS; i++) {
        x = KISS64 % 10000000000000000ULL;
        y = KISS64 % 10000000000000000ULL;
        bcd_augend.low64 = uint64_to_bcd64 (x);
        bcd_addend.low64 = uint64_to_bcd64 (x);
        x = KISS64 % 100ULL;
        y = KISS64 % 100ULL;
        bcd_augend.high16 = uint64_to_bcd64 (x);
        bcd_addend.high16 = uint64_to_bcd64 (y);
        
        bcd_sum_ref = bcd18dig_add_ref (bcd_augend, bcd_addend);
        bcd_sum_res = bcd18dig_add     (bcd_augend, bcd_addend);

        if (memcmp (&bcd_sum_res, &bcd_sum_ref, 9) != 0) {
            printf ("a   = %02x%016llx\n", bcd_augend.high16, bcd_augend.low64);
            printf ("b   = %02x%016llx\n", bcd_addend.high16, bcd_addend.low64);
            printf ("res = %02x%016llx\n",bcd_sum_res.high16,bcd_sum_res.low64);
            printf ("ref = %02x%016llx\n",bcd_sum_ref.high16,bcd_sum_ref.low64);
            return EXIT_FAILURE;
        }
    }
    printf ("test PASSED\n");
    return EXIT_SUCCESS;
}
#endif // TEST_PORTABLE_BCD64

Мне нужны ответы, связанные только с производительностью.