6

One common way to produce an approximation of a function like the logarithm or the exponential is to precompute a table of values (a lookup table) for the output or some intermediate stage of the algorithms and sample that at runtime.

There are lots of ways to sample these lookup tables, depending on how the function itself is being computed, but I am interested in a very specific one: tables which are sampled by bits of the float or double input.

To take an example, this gist reproduces a function from 1986 which computes the inverse square root of a double by sampling a table. A more compact version is:

float KahanNgISR(double x) {
    double y;
    static int lookup[64]= {
        0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
        0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
        0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
        0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
        0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
        0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
        0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
        0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd
    };
    uint32_t xUPPER, yUPPER, k;
    uint64_t xINT, yINT;
    xINT = *( uint64_t *) &x;
    xUPPER = (xINT & 0xffffffff00000000) >> 32;
    k = 0x5fe80000 - (xUPPER >> 1);
    yUPPER = k - lookup[63 & (k >> 14)];
    yINT = ((uint64_t) yUPPER << 32);
    y = *( double* ) &yINT;
    y = (float) y;
    y = y * (1.5f - powf(2, -30) - (0.5f * x * y * y));
    return y;
}

The specific line here: yUPPER = k - lookup[63 & (k >> 14)];, is where the indexing occurs. An expert analysis of the above code is over on stackoverflow, for here I am more interested in it as an example of the pattern.

There is a similar pattern in Lalonde and Dawson's article in Graphics Gems (1990, pp. 424-426 and 756-757 for the appendix) computation of the square root (not inverse):

#include <math.h>
/* SPARC floating point format is as follows
BIT 31   30    23 22     0
    sign exponent mantissa
*/
static short sqrttab[0x100];
void build_table() {
    unsigned short i;
    float f;
    unsigned int *fi=&f;
    for (i = 0; i <= 0x7f; i++) {
        *fi = 0;
        /* Build a float with the bit pattern i as mantissa
         * and an exponent of 0, stored as 127 */
        */
       *fi = (i << 16) | (127 << 23);
       f = sqrt(f);
       /* Take the square root then strip the first 7 bits of
        * the mantissa into the table
        */
       sqrrtab[i] = (*fi & 0x7fffff) >> 16;
       /* Repeat the process, this time with an exponent of
        * 1, stored as 128
        */
       *fi = 0;
       *fi = (i << 16) | (128 << 23);
       f = sqrt(f);
       sqrrtab[i + 0x80] = (*fi & 0x7fffff) >> 16;
    }
}

/*
* fsqrt - fast square root by table lookup
*/

float fsqrt(float n) {
    unsigned int *num = &n;  /* to access the bits of a float in C
                              * we must misuse pointers*/
    short e;                 /* the exponent */
    if (n == 0) return (0);  /* check for square root of 0 */
    e = (*num >> 23) - 127   /* get the exponent - on a SPARC the 
                              * exponent is stored with 127 added*/
                              /* leave only the mantissa */
    if (e & 0x01) *num | = 0x800000;
                              /* the exponent is odd so we have to
                               look it up in the second half of
                               the lookup table, so we set the 
                               high bit */
    e >>= 1;                 /* divide the exponent by 2 */ 
                             /* note that in C the shift */                        
                             /* operators are sign preserving */   
                             /* for signed operands */
/* Do the table lookup, based on the quartenary mantissa,
   then reconstruct the result back into a float
*/
    *num = (sqrttab[*num >> 16] << 16) | ((e + 127) << 23);
    return (n);
}

I apologize for the long inclusion of code, but the C appendix is not widely reproduced for Graphics Gems and it is often easier to show intent with example code. The above code has been copied verbatim including comments and comment style. I have not tested to see if it runs.

The above examples are from 1986 and 1990, respectively. I know this is a much older pattern, so I am hoping that folks can offer examples of bit sampling pattern, preferably in C which are around these dates or earlier. If folks have knowledge of patterns like this without examples of their use I'd love to hear about those. Articles in academic and industry journals are useful as well.

Examples in other languages or in assembly, especially those which are older, are also of interest to me.

12
  • 1
    FWIW, the 1985 BSD math library sources, also based (largely?) on the work of Kahan and Ng, appear to contain no such constructs. Obvious candidates would be implementations of sqrt, cbrt, log, exp, pow. Jerome Coonen's Ph.D. thesis from 1984 (Pascal code) does not seem to contain any instance of such a construct.
    – njuffa
    Feb 2, 2023 at 21:33
  • 2
    I’m voting to close this question because it belongs elsewhere, either Stack Overflow or Computer Science perhaps. Feb 2, 2023 at 21:58
  • It is not a computer science question. I'm happy to ask it on SO if it doesn't belong here, but I submit it does. Feb 2, 2023 at 22:09
  • 1
    Best I can tell, Hirondo Kuki's 1966 math library for the IBM 7094 (written in assembly language) does not use any such table look-ups.
    – njuffa
    Feb 2, 2023 at 22:11
  • 4
    This question is asking about the earliest use, which is likely in the retro era, so this question is on-topic. If it had been asking about details of the algorithm, that would be off-topic and better asked on SO.
    – DrSheldon
    Feb 3, 2023 at 3:10

1 Answer 1

2

The below isn't the best answer because it is newer than 1986 and 1990, but neither should it be a comment or an extension of the question.

Shane Peelar uncovered an implementation of the "fast inverse square root" (c.f. the Wikipedia article) in a game released in 1997 called Interstate 76: Fast reciprocal square root... in 1997?!

A version of the code is below (with a precomputed LUT); you can read it on their blog post or on this repo. They use the bits of a floating-point number (in this case a double) to index a 256 bit table.

static uint8_t LUT[256] = {
    106, 105, 103, 102, 101, 99, 98, 97, 95, 94, 93, 91,
    90, 89, 88, 87, 85, 84, 83, 82, 81, 80, 78, 77, 76,
    75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63,
    62, 61, 60, 59, 58, 57, 56, 55, 55, 54, 53, 52, 51,
    50, 49, 48, 48, 47, 46, 45, 44, 44, 43, 42, 41, 40,
    40, 39, 38, 37, 37, 36, 35, 34, 34, 33, 32, 31, 31,
    30, 29, 29, 28, 27, 27, 26, 25, 25, 24, 23, 23, 22,
    21, 21, 20, 20, 19, 18, 18, 17, 16, 16, 15, 15, 14,
    13, 13, 12, 12, 11, 11, 10, 10, 9, 8, 8, 7, 7, 6, 6,
    5, 5, 4, 4, 3, 3, 2, 2, 1, 1, 255, 254, 252, 250, 248,
    246, 244, 243, 241, 239, 237, 235, 234, 232, 230, 228,
    227, 225, 223, 222, 220, 219, 217, 215, 214, 212, 211,
    209, 208, 206, 205, 203, 202, 201, 199, 198, 196, 195,
    194, 192, 191, 190, 188, 187, 186, 184, 183, 182, 181,
    179, 178, 177, 176, 175, 173, 172, 171, 170, 169, 168,
    166, 165, 164, 163, 162, 161, 160, 159, 158, 157, 156,
    155, 154, 153, 152, 151, 150, 149, 148, 147, 146, 145,
    144, 143, 142, 141, 140, 139, 138, 137, 136, 135, 135,
    134, 133, 132, 131, 130, 129, 128, 128, 127, 126, 125,
    124, 123, 123, 122, 121, 120, 119, 119, 118, 117, 116,
    116, 115, 114, 113, 113, 112, 111, 110, 110, 109, 108,
    107, 107
};

double i76_rsqrt(double const scalar) {
    uint64_t const scalar_bits = std::bit_cast<uint64_t>(scalar);

    uint8_t const index = (scalar_bits >> 0x2d) & 0xff;

    //LUT[index] contains the 8 most significant bits of the mantissa, rounded up.
    //Treat all lower 44 bits as zeroed out
    uint64_t const mantissa_bits = static_cast<uint64_t>(LUT[index]) << 0x2c;
   
    //Exponent bits are calculated based on the formula in the article
    uint64_t const exponent_bits = ((0xbfcUL - (scalar_bits >> 0x34)) >> 1) << 0x34;
    
    //exponent_bits have form 0xYYY00000 00000000
    //mantissa_bits have form 0x000ZZ000 00000000
    //so combined, we have    0xYYYZZ000 00000000 -- a complete float64 for the guess
    uint64_t const combined_bits = exponent_bits | mantissa_bits;

    auto const initial_guess = std::bit_cast<double>(combined_bits);

    auto const half_initial_guess = initial_guess * 0.5;
    auto const initial_guess_squared = initial_guess * initial_guess;
    
    //One iteration of Newton-Raphson
    auto const newton_raphson = (3.0 - (scalar * initial_guess_squared)) * half_initial_guess;

    //post-hoc fixup
    auto const fixup = newton_raphson * 1.00001;

    return fixup; 
}

They sample differently from Kahan & Ng as well as LaLonde and Dawson but I believe this is the same pattern.

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .