Scene Storm/Scene Storm - Volume 1.iso/coding/c/pdc/libsrc/math/sqrt.c

C/C++ Source or Header | 1990-04-07 | 5.5 KB | 179 lines

/************************************************************************
 *                                                                      *
 *                              N O T I C E                             *
 *                                                                      *
 *                      Copyright Abandoned, 1987, Fred Fish            *
 *                                                                      *
 *      This previously copyrighted work has been placed into the       *
 *      public domain by the author (Fred Fish) and may be freely used  *
 *      for any purpose, private or commercial.  I would appreciate     *
 *      it, as a courtesy, if this notice is left in all copies and     *
 *      derivative works.  Thank you, and enjoy...                      *
 *                                                                      *
 *      The author makes no warranty of any kind with respect to this   *
 *      product and explicitly disclaims any implied warranties of      *
 *      merchantability or fitness for any particular purpose.          *
 *                                                                      *
 ************************************************************************
 */

/*
 *  FUNCTION
 *
 *      sqrt   double precision square root
 *
 *  KEY WORDS
 *
 *      sqrt
 *      machine independent routines
 *      math libraries
 *
 *  DESCRIPTION
 *
 *      Returns double precision square root of double precision
 *      floating point argument.
 *
 *  USAGE
 *
 *      double sqrt (x)
 *      double x;
 *
 *  REFERENCES
 *
 *      Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
 *
 *      Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *      1968, pp. 89-96.
 *
 *  RESTRICTIONS
 *
 *      The relative error is 10**(-30.1) after three applications
 *      of Heron's iteration for the square root.
 *
 *      However, this assumes exact arithmetic in the iterations
 *      and initial approximation.  Additional errors may occur
 *      due to truncation, rounding, or machine precision limits.
 *
 *  PROGRAMMER
 *
 *      Fred Fish
 *
 *  INTERNALS
 *
 *      Computes square root by:
 *
 *        (1)   Range reduction of argument to [0.5,1.0]
 *              by application of identity:
 *
 *              sqrt(x)  =  2**(k/2) * sqrt(x * 2**(-k))
 *
 *              k is the exponent when x is written as
 *              a mantissa times a power of 2 (m * 2**k).
 *              It is assumed that the mantissa is
 *              already normalized (0.5 =< m < 1.0).
 *
 *        (2)   An approximation to sqrt(m) is obtained
 *              from:
 *
 *              u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
 *
 *              P0 = 0.594604482
 *              P1 = 2.54164041
 *              Q0 = 2.13725758
 *              Q1 = 1.0
 *
 *              (coefficients from HART table #350 pg 193)
 *
 *        (3)   Three applications of Heron's iteration are
 *              performed using:
 *
 *              y[n+1] = 0.5 * (y[n] + (m/y[n]))
 *
 *              where y[0] = u = approx sqrt(m)
 *
 *        (4)   If the value of k was odd then y is either
 *              multiplied by the square root of two or
 *              divided by the square root of two for k positive
 *              or negative respectively.  This rescales y
 *              by multiplying by 2**frac(k/2).
 *
 *        (5)   Finally, y is rescaled by int(k/2) which
 *              is equivalent to multiplication by 2**int(k/2).
 *
 *              The result of steps 4 and 5 is that the value
 *              of y between 0.5 and 1.0 has been rescaled by
 *              2**(k/2) which removes the original rescaling
 *              done prior to finding the mantissa square root.
 *
 *  NOTES
 *
 *      The Convergent Technologies compiler optimizes division
 *      by powers of two to "arithmetic shift right" instructions.
 *      This is ok for positive integers but gives different
 *      results than other C compilers for negative integers.
 *      For example, (-1)/2 is -1 on the CT box, 0 on every other
 *      machine I tried.
 *
 */

#include <stdio.h>
#include "pml.h"

#define P0 0.594604482                  /* Approximation coeff  */
#define P1 2.54164041                   /* Approximation coeff  */
#define Q0 2.13725758                   /* Approximation coeff  */
#define Q1 1.0                          /* Approximation coeff  */

#define ITERATIONS 3                    /* Number of iterations */

static char funcname[] = "sqrt";

extern double frexp ();
extern double ldexp ();

double
sqrt (x)
        double x;
{
    int k;
    int bugfix;
    int kmod2;
    int count;
    int exponent;
    double m;
    double u;
    double y;
    double rtnval;
    struct exception xcpt;

    DBUG_ENTER ("sqrt");
    DBUG_3 ("sqrtin", "arg %le", x);
    if (x == 0.0) {
        rtnval = 0.0;
    } else if (x < 0.0) {
        xcpt.type = DOMAIN;
        xcpt.name = funcname;
        xcpt.arg1 = x;
        if (!matherr (&xcpt)) {
            fprintf (stderr, "%s: DOMAIN error\n", funcname);
            errno = EDOM;
            xcpt.retval = 0.0;
        }
    } else {
        m = frexp (x, &k);
        u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
        for (count = 0, y = u; count < ITERATIONS; count++) {
            y = 0.5 * (y + (m / y));
        }
        if ((kmod2 = (k % 2)) < 0) {
            y /= SQRT2;
        } else if (kmod2 > 0) {
            y *= SQRT2;
        }
        bugfix = 2;
        xcpt.retval = ldexp (y, k/bugfix);
    }
    DBUG_3 ("sqrtout", "result %le", xcpt.retval);
    DBUG_RETURN (xcpt.retval);
}