Tabled rcpss (should match HW version for stated value range on Skylake, anyway)

rcpss.cpp

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

static uint32_t recip(uint32_t bits)
{
    uint32_t u;
    float f;
    memcpy(&f, &bits, sizeof(bits));
    __m128 x = _mm_load_ss(&f);
    __m128 recip = _mm_rcp_ss(x);
    _mm_store_ss(&f, recip);
    memcpy(&u, &f, sizeof(f));
    return u;
}

static void print_recips(int count)
{
    for (int i = 0; i < count; ++i)
    {
        uint32_t bits = 0x3f800000 + i;
        uint32_t r = recip(bits);
        printf("[%4d] in=0x%08x r=0x%08x\n", i, bits, r);
    }
}

static void test_recip_quant(int step)
{
    for (int base = 0; base < (1<<23); base += step)
    {
        uint32_t bits_base = 0x3f800000 + base;
        uint32_t r_block = recip(bits_base);
        for (int i = 0; i < step; ++i)
        {
            uint32_t r = recip(bits_base + i);
            if (r != r_block)
            {
                printf("step %d failure at in=0x%08x\n", step, bits_base + i);
                return;
            }
        }
    }

    printf("step %d works\n", step);
}

static void test_recip_table(int step_bits, int res_bits)
{
    int shift = 23 - res_bits;
    int step = 1<<step_bits;
    for (int i = 0; i < (1<<23); i += step)
    {
        uint32_t bits = 0x3f800000 + i;
        uint32_t r = recip(bits);
        uint32_t r_red = r >> shift;
        if ((r_red << shift) != r)
        {
            printf("res fail! x=0x%08x\n", bits);
            return;
        }

        printf("0x%08x 0x%08x\n", bits, r);
    }
}

static const int s_recip_table_lutbits = 11;
static const int s_recip_table_lutwidth = 12;
static const int s_recip_table_lutcount = 1 << s_recip_table_lutbits;
static uint16_t s_recip_table[s_recip_table_lutcount];

static void build_recip_table()
{
    for (int i = 0; i < s_recip_table_lutcount; i++)
    {
        uint32_t bits = 0x3f800000 + (i << (23 - s_recip_table_lutbits));
        uint32_t r = recip(bits);
        uint32_t mant_bits = (r - 0x3f000000) >> (23 - s_recip_table_lutwidth);
        s_recip_table[i] = (uint16_t)mant_bits;
    }
}

static uint32_t tabled_recip(uint32_t input)
{
    // not-too-large normalized values only!
    // not handled yet:
    //   0, subnormals
    //   inf, nan
    //   reciprocals of values with max finite exponent (biased exponent=254)

    // not handling special values yet
    uint32_t sign = input >> 31;
    uint32_t biased_exp = (input >> 23) & 0xff;
    uint32_t mant = input & 0x7fffff;

    // perform recip calc
    uint32_t rcp_mant = s_recip_table[mant >> (23 - s_recip_table_lutbits)] << (23 - s_recip_table_lutwidth);
    uint32_t rcp_biased_exp = 253 - biased_exp;

    // re-assemble float
    return (sign << 31) | (rcp_biased_exp << 23) | rcp_mant;
}

static void test_table_version()
{
    for (int i = 0x3e000000; i < 0x41000000; i++)
    {
        uint32_t ref = recip(i);
        uint32_t mine = tabled_recip(i);

        if (ref != mine)
        {
            printf("x=0x%08x ref=0x%08x mine=0x%08x\n", i, ref, mine);
            return;
        }
    }

    printf("looks OK\n");
}

int main()
{
    //print_recips(16384);
    //test_recip_quant(4096);
    //test_recip_table(12, 12);
    build_recip_table();
    test_table_version();
    return 0;
}

// vim:sw=4:sts=4:et