Subject: New square root algorithm in molecular dynamics

From: John Marelius <[email-masked]>
Date: Sat, 7 Sep 1996 11:48:13 +0200 (MDT)

Here is a copy of Terje Mathisen's news posting describing his square root
algorithm to which I referred in my recent posting to the list.

John Marelius
---------------------------------------------------------------------------
Terje Mathisen's article:

John Marelius wrote:
> 
> This posting summarises the results of the discussion on square-root
> calculations in molecular dynamics (MD) applications. Using a novel
> algorithm by Terje Mathisen and general optimisations of the code,
> performance has increased by up to 55%, depending on the simulated system.

[snip]
 
> One person who became interested in the problem was Terje Mathisen, he
> started thinking and coding and has overwhelmed me with a stream of
> ever-improving codes for fast and accurate functions to calculate
> sqrt(1/x).
> 
> Please refer to his posting, describing the details of the algorithm!

OK, here's the story John promised I'd write:

When John made his post with a the inner loop of his simulation code,
one thing was obvious: He had two expensive operations in his inner
loop: One square root, and one division. The rest was all floating point
muls and adds, with no compares or branches.

After a while I realized that if the division and square root steps
could be combined, it would be much easier to speed his code up.

My first attempt separated the exponent and the mantissa part of the
input value, keeping the least significant bit of the exponent together
with the leading bits of the mantissa, and then used linear
interpolation with a fairly large lookup table (32KB, 4K double prec
entries).

This version was very fast, clocking in at about 20 cycles, but gave
only about 32 significant bits in the result.

I also tried higher-order interpolation functions:

A quadratic curve gave another 7-8 bits, but this doubled the processing
time.

The final solution was to use a lookup table for range reduction,
looking up 1/x, and multiplying by it, leaving a number that's very
close to 1.0.

I then used a third-order polynomial on the reduced-range value, before
multiplying the result with the product of two other table lookups: One
for the exponent and another for sqrt(1/x), i.e. the square root of the
values in the first table.

After John tested the code and found that the table-driven version gave
_exactly_ the same 8-digit results for all atoms during the first 50
simulation steps, I decided that we could probably reduce the table
sizes a bit.

At this time I finally figured out how to convert my (truncated) Taylor
polynomial into a Chebychev version, which effectively recovered all the
precision lost by reducing the lookup tables.

I also decided that by making the sqrt(1/x) table single precision, the
double precision (1/x) table would contain the exact squares of the
first table, which would both reduce the table size and increase the
accuracy.

After inspecting the final code, it was obvious that it suffered from
several pipeline stalls, due to direct linear dependencies. By writing a
version that calculated two inverse square roots in parallel,
inv_sqrt2(x, y), most of these losses could be recovered.

The final version needs a good scheduling compiler, but on an Alpha it
runs about 2.5 times faster than the compiler-generated code for
sqrt(1.0/x), leading to the 55% faster results for John's simulations.

Here's the working code:

/* Table-driven sqrt(1/x) functions
   (c) Terje Mathisen 1996

   Use freely as long as my copyright is retained

   Gives approximately 49 significant bits (15 decimal digits).
*/

#include <math.h>

#ifdef BIG_ENDIAN
  #define HIGH_WORD 0
#else
  #define HIGH_WORD 1
#endif

#define MANTISSA_BITS 20
#define MANTISSA_MASK ((1 << MANTISSA_BITS) - 1)
#define EXP_BITS  11
#define EXP_SIZE  (1 << EXP_BITS)
#define EXP_BIAS  1023
#define EXP_ZERO  (EXP_BIAS << MANTISSA_BITS)

/* Here I define the size of the mantissa lookup tables: */
#define MANT_BITS 10
#define MANT_SIZE (1 << MANT_BITS)
#define MANT_SHIFT (MANTISSA_BITS - MANT_BITS)
#define MANT_ROUND (1 << (MANT_SHIFT - 1))

static double exp_root[EXP_SIZE];   /* sqrt(1.0/(2^exp)) */
static float root_x[MANT_SIZE+1];   /* sqrt(1/x)         */
static double inv_x[MANT_SIZE+1];   /* root_x^2          */

/* Taylor terms: */
#define t0 (1.0)
#define t1 (-1.0/2.0)
#define t2 (3.0/8.0)
#define t3 (-5.0/16.0)
#define t4 (35.0/128.0)
#define t5 (-63.0/256.0)

#ifdef TAYLOR_SERIES
  #define c0 t0
  #define c1 t1
  #define c2 t2
  #define c3 t3
#else
  /* These Chebychev terms have been calculated while assuming
     that the 5th and higher orders of the Taylor series can be
     disregarded.
     They must be re-calculated if MANT_BITS is reduced below 10!
  */
  #define c0 (t0 -
t4/(8.0*2.0*16.0*MANT_SIZE*MANT_SIZE*MANT_SIZE*MANT_SIZE))
  #define c1 t1
  #define c2 (t2 + t4/(4.0*MANT_SIZE*MANT_SIZE))
  #define c3 t3
#endif

double inv_sqrt(double x)
{
    double dx, dx2, scale, poly;
    long *plx = (long *) &x;
    long exp, mant;

    exp = plx[HIGH_WORD];
    mant = exp & MANTISSA_MASK;
    plx[HIGH_WORD] = mant + EXP_ZERO;          /* Force exponent == 0 */
    mant = (mant + MANT_ROUND) >> MANT_SHIFT;  /* 11-bit rounded
mantissa */
    exp >>= MANTISSA_BITS;                     /* 11-bit exponent */

    dx = x * inv_x[mant] - 1.0;
    scale = exp_root[exp] * root_x[mant];     /* scale = product of
exponent and mantissa term */

    dx2 = dx*dx;
    poly = c0 + dx2*c2 + dx*(c1 + dx2*c3);
    return poly * scale;
}

double inv_sqrt2(double x, double *y)
{
    double dx, dx2, scalex, polyx;
    double dy, dy2, scaley, polyy;
    long *plx = (long *) &x;
    long *ply = (long *) y;
    long expx, mantx;
    long expy, manty;

    expx = plx[HIGH_WORD];
    expy = ply[HIGH_WORD];
    
    mantx = expx & MANTISSA_MASK;
    manty = expy & MANTISSA_MASK;
    
    plx[HIGH_WORD] = mantx + EXP_ZERO;   /* Force exponent == 0 */
    ply[HIGH_WORD] = manty + EXP_ZERO;   /* Force exponent == 0 */
    
    mantx = (mantx + MANT_ROUND) >> MANT_SHIFT;  /* rounded mantissa */
    manty = (manty + MANT_ROUND) >> MANT_SHIFT;  /* rounded mantissa */
    
    expx >>= MANTISSA_BITS;                      /* 11-bit exponent */
    expy >>= MANTISSA_BITS;                      /* 11-bit exponent */

    dx = x * inv_x[mantx] - 1.0;
    dy = *y * inv_x[manty] - 1.0;
    
    scalex = exp_root[expx] * root_x[mantx]; /* scale = product of
exponent and mantissa term */
    scaley = exp_root[expy] * root_x[manty]; /* scale = product of
exponent and mantissa term */

    dx2 = dx*dx;
    dy2 = dy*dy;
    
    polyx = c0 + dx2*c2 + dx*(c1 + dx2*c3);
    polyy = c0 + dy2*c2 + dy*(c1 + dy2*c3);

    polyx *= scalex;
    polyy *= scaley;

    *y = polyy;
    return polyx;
}

void init_tables(void)
{
    double x, root, dx;
    int i;

    /* First initialize the exp_root table: */
    root = sqrt(0.5); x = 1.0;
    for (i = 1023; i < 2047; i+=2) {
        exp_root[i] = x;
        exp_root[i+1] = x*root;
        x *= 0.5;
    }
    exp_root[2047] = 0.0;               /* sqrt(1.0/Inf) = 0.0 */
    
    x = 2.0;
    for (i = 1022; i > 0; i-=2) {
        exp_root[i] = x*root;
        exp_root[i-1] = x;
        x *= 2.0;
    }
    exp_root[0] = 0.0;
    ((long *) &exp_root[0])[HIGH_WORD] = 0x7ff00000;  /* Infinity */

    /* Next, do the two mantissa-based tables:
       Bug workaround: I must first generate the single prec table,
       and then the double prec inv_x[] table!
       Otherwise, the compiler (Watcom) refuses to use the rounded
       (f) value.
    */
    dx = 1.0 / MANT_SIZE;
    for (i = 0; i <= MANT_SIZE; i++) {
        x = 1.0 + dx * i;
        root_x[i] = sqrt(1.0/x);
    }
    for (i = 0; i <= MANT_SIZE; i++) {
        inv_x[i] = root_x[i]*root_x[i];
    }
}

-- 
- <[email-masked]>
Using self-discipline, see http://www.eiffel.com/discipline
"almost all programming can be viewed as an exercise in caching"