What is the earliest use in C of indexing the bits of a float or double to sample a table lookup?

Question

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.

Answer 1

This is also not an exact answer because it is not in C, however...

R.W. Bemer reported in Communications of the ACM (CACM) May 1958 on A subroutine method for calculating logarithms for binary and decimal floating point machines. His method included construction of a table and indexing a table lookup by bits in the floating point number. The method of extracting the bits is not prescribed, with page 7 noting "[a]ddress formation may thus be done by extraction or some other direct method."

The table is also included on the same page, reproduced here as a scan:

The same author writes in a later (1963) CACM note that the techniques don't seem to have been taken up widely and that "The table operations in [A subroutine method for calculating logarithms] may be difficult to handle cleanly on many machines" (p. 306 of A note on range transformations for square root and logarithm) when proposing alternate methods, so it's unclear to me how widespread the influence of this was, but I this is the earliest example of a table lookup indexed by float bits I've seen in print. However, the late 1950s is when these publications started up with both CACM and Computer starting publication in 1958 so the lack of evidence for table lookups prior to this likely does not indicate lack of their use.

Answer 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.