Listing 1: tsqrt function and test program

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
static unsigned short linearInterpolation(unsigned short, unsigned short,
	                                  unsigned short, int);

/* Lookup table:  table[i] = 128*16*SquareRoot(i) */
static unsigned short table[129];
static const unsigned short *pTable = table;

/**********
* Initialize the lookup table
***********/
static void initialize_table(void) {
    int i;
    for (i=0; i<129; i++)
        table[i] = (unsigned short) (sqrt((double)i) * (128*16) + 0.5);
    }

/**********
* Calculate scaled square root
* result = 128 * SquareRoot(input)
***********/
unsigned short tsqrt(
	unsigned short x	/* Integer whose square root is desired */
	) {
    unsigned short result;
    if (x <= 127)
        result = (pTable[(int)x] + 8) / 16;
    else if (x <= 511)
        result = (linearInterpolation(x, 128, 4, 32) + 4) / 8;
    else if (x <= 2047)
        result = (linearInterpolation(x, 512, 16, 32) + 2) / 4;
    else if (x <= 8191)
        result = (linearInterpolation(x, 2048, 64, 32) + 1) / 2;
    else if (x <= 32767)
        result = linearInterpolation(x, 8192, 256, 32);
    else if (x <= 65535)
        result = (unsigned short) (linearInterpolation(x, 32768, 1024, 32) * 2);
    else
        result = 128*256;	/* How did this happen? */
    return result;
    }

/**********
* Perform linear interpolation:
* result has same scaling as the table entries.
***********/
static unsigned short linearInterpolation(
	unsigned short x,	/* Integer whose square root is desired */
	unsigned short min,	/* Minimum end of this range */
	unsigned short dx,	/* x1-x0 for this range */
	int base		/* base index for this range */
	) {
    int index0;
    unsigned short x0, y0, y1, y;
    index0 = (int) (((x - min) / dx) + base);
    x0 = (unsigned short) ((index0 - base) * dx + min);
    y0 = pTable[index0];
    y1 = pTable[index0+1];
    y = (unsigned short) (((long)(x - x0) * (y1 - y0) + (long)(y0*dx) +
                          (dx/2)) / dx);
    return y;
    }

/**********
* Initialize the lookup table,
* Compute error and percent error over entire input domain.
***********/
void main(void) {
    long maxErrorInteger[6];
    double maxPercentErrorInteger[6];
    double maxErrorValue[6];
    double maxPercentErrorValue[6];
    int subdomain;
    long x;

    /* Initialize the lookup table */
    initialize_table();

    /* Clear the maximum error arrays */
    for (subdomain=0; subdomain<6; subdomain++) {
        maxErrorInteger[subdomain] = 0;
        maxPercentErrorInteger[subdomain] = 0;
        maxErrorValue[subdomain] = 0;
        maxPercentErrorValue[subdomain] = 0;
        }

    /* Test tsqrt() results, saving maximum error in each subdomain */
    for (x=0,subdomain=0; x<65536; x++) {
        long errorInteger;
        double percentErrorInteger, errorValue, percentErrorValue;
        double tableValue, goldenValue;
        unsigned short tableInteger, goldenInteger;

        tableInteger = tsqrt((unsigned short)x);
        goldenValue = sqrt((double)x);

        tableValue = (double)tableInteger / 128.;
        goldenInteger = ((unsigned long) (goldenValue * (128*16) + 0.5) + 8) / 16;

        errorInteger = (long)goldenInteger - (long)tableInteger;
        if (errorInteger < 0)   errorInteger = -errorInteger;
        if (goldenInteger == 0)
            percentErrorInteger = 0;
        else
            percentErrorInteger = (double)errorInteger /
                                  (double)goldenInteger * 100;

        errorValue = goldenValue - tableValue;
        if (errorValue < 0)   errorValue = -errorValue;
        if (goldenValue == 0)
            percentErrorValue = 0;
        else
            percentErrorValue = errorValue / goldenValue * 100;

        if (errorInteger > maxErrorInteger[subdomain])
            maxErrorInteger[subdomain] = errorInteger;
        if (percentErrorInteger > maxPercentErrorInteger[subdomain])
            maxPercentErrorInteger[subdomain] = percentErrorInteger;
        if (errorValue > maxErrorValue[subdomain])
            maxErrorValue[subdomain] = errorValue;
        if (percentErrorValue > maxPercentErrorValue[subdomain])
            maxPercentErrorValue[subdomain] = percentErrorValue;

        if ((x == 127) || (x == 511) || (x == 2047) || (x == 8191) ||
            (x == 32767))
            subdomain++;
        }

    (void) printf("              Integer                    Floating Point\n");
    (void) printf("Range   maxError  maxPercentError   maxError
                  maxPercentError\n");
    for (subdomain=0; subdomain<6; subdomain++) {
        (void) printf("%3d         %d        %8f%%      %8f     %8f%%\n",
		subdomain, maxErrorInteger[subdomain], maxPercentErrorInteger[subdomain],
		maxErrorValue[subdomain], maxPercentErrorValue[subdomain]);
        }
    }