home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / icon8_fx / sqrt.c < prev   
Encoding:
C/C++ Source or Header  |  1993-10-23  |  5.9 KB  |  232 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    sqrt   double precision square root
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    sqrt
  29.  *    machine independent routines
  30.  *    math libraries
  31.  *
  32.  *  DESCRIPTION
  33.  *
  34.  *    Returns double precision square root of double precision
  35.  *    floating point argument.
  36.  *
  37.  *  USAGE
  38.  *
  39.  *    double sqrt (x)
  40.  *    double x;
  41.  *
  42.  *  REFERENCES
  43.  *
  44.  *    Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
  45.  *
  46.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  47.  *    1968, pp. 89-96.
  48.  *
  49.  *  RESTRICTIONS
  50.  *
  51.  *    The relative error is 10**(-30.1) after three applications
  52.  *    of Heron's iteration for the square root.
  53.  *
  54.  *    However, this assumes exact arithmetic in the iterations
  55.  *    and initial approximation.  Additional errors may occur
  56.  *    due to truncation, rounding, or machine precision limits.
  57.  *    
  58.  *  PROGRAMMER
  59.  *
  60.  *    Fred Fish
  61.  *
  62.  *  INTERNALS
  63.  *
  64.  *    Computes square root by:
  65.  *
  66.  *      (1)    Range reduction of argument to [0.5,1.0]
  67.  *        by application of identity:
  68.  *
  69.  *        sqrt(x)  =  2**(k/2) * sqrt(x * 2**(-k))
  70.  *
  71.  *        k is the exponent when x is written as
  72.  *        a mantissa times a power of 2 (m * 2**k).
  73.  *        It is assumed that the mantissa is
  74.  *        already normalized (0.5 =< m < 1.0).
  75.  *
  76.  *      (2)    An approximation to sqrt(m) is obtained
  77.  *        from:
  78.  *
  79.  *        u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
  80.  *
  81.  *        P0 = 0.594604482
  82.  *        P1 = 2.54164041
  83.  *        Q0 = 2.13725758
  84.  *        Q1 = 1.0
  85.  *
  86.  *        (coefficients from HART table #350 pg 193)
  87.  *
  88.  *      (3)    Three applications of Heron's iteration are
  89.  *        performed using:
  90.  *
  91.  *        y[n+1] = 0.5 * (y[n] + (m/y[n]))
  92.  *
  93.  *        where y[0] = u = approx sqrt(m)
  94.  *
  95.  *      (4)    If the value of k was odd then y is either
  96.  *        multiplied by the square root of two or
  97.  *        divided by the square root of two for k positive
  98.  *        or negative respectively.  This rescales y
  99.  *        by multiplying by 2**frac(k/2).
  100.  *
  101.  *      (5)    Finally, y is rescaled by int(k/2) which
  102.  *        is equivalent to multiplication by 2**int(k/2).
  103.  *
  104.  *        The result of steps 4 and 5 is that the value
  105.  *        of y between 0.5 and 1.0 has been rescaled by
  106.  *        2**(k/2) which removes the original rescaling
  107.  *        done prior to finding the mantissa square root.
  108.  *
  109.  *  NOTES
  110.  *
  111.  *    The Convergent Technologies compiler optimizes division
  112.  *    by powers of two to "arithmetic shift right" instructions.
  113.  *    This is ok for positive integers but gives different
  114.  *    results than other C compilers for negative integers.
  115.  *    For example, (-1)/2 is -1 on the CT box, 0 on every other
  116.  *    machine I tried.
  117.  *
  118.  */
  119. /*
  120.  *  MODIFICATIONS
  121.  *
  122.  *    This routine modified to use polynomial, instead of rational,
  123.  *    approximation with coefficients  
  124.  *
  125.  *        P0  0.22906994529e+00
  126.  *        P1  0.13006690496e+01
  127.  *        P2 -0.90932104982e+00
  128.  *        P3  0.50104207633e+00
  129.  *        P4 -0.12146838249e+00
  130.  *
  131.  *    as given by Hart (op. cit.) in SQRT 0132.  This approximation
  132.  *    gives around 5 digits correct.  Two applications of Heron's
  133.  *    approximation will give more then practically achievable
  134.  *    13-15 decimal digits.  Multiplications by powers of 2 are
  135.  *    replaced by explicit calls to ldexp.
  136.  *
  137.  *    Michal Jaegermann, 24 October 1990
  138.  */
  139.  
  140. #include <stdio.h>
  141. #include <pmluser.h>
  142. #include "pml.h"
  143.  
  144. #ifdef OLD
  145. #define P0 0.594604482            /* Approximation coeff    */
  146. #define P1 2.54164041            /* Approximation coeff    */
  147. #define Q0 2.13725758            /* Approximation coeff    */
  148. #define Q1 1.0                /* Approximation coeff    */
  149.  
  150. #define ITERATIONS 3            /* Number of iterations    */
  151.  
  152. #endif
  153.  
  154. #define P0  0.22906994529e+00        /* Hart SQRT 0132 */
  155. #define P1  0.13006690496e+01
  156. #define P2 -0.90932104982e+00
  157. #define P3  0.50104207633e+00
  158. #define P4 -0.12146838249e+00
  159.  
  160. static char funcname[] = "sqrt";
  161.  
  162.  
  163. double sqrt (x)
  164. double x;
  165. {
  166. #ifdef OLD
  167.     auto int k;
  168.     register int bugfix;
  169.     register int kmod2;
  170.     register int count;
  171.     auto int exponent;
  172. #else
  173. #ifdef __MSHORT__
  174.     auto short k;
  175. #else
  176.     auto int k;
  177. #endif
  178. #endif
  179.     auto double m;
  180.     auto double u;
  181. #ifdef OLD
  182.     auto double y;
  183. #endif
  184.     auto struct exception xcpt;
  185.     extern double frexp ();
  186.     extern double ldexp ();
  187.  
  188.     DBUG_ENTER ("sqrt");
  189.     DBUG_3 ("sqrtin", "arg %le", x);
  190.     if (x == 0.0) {
  191.     xcpt.retval = 0.0;
  192.     } else if (x < 0.0) {
  193.     xcpt.type = DOMAIN;
  194.     xcpt.name = funcname;
  195.     xcpt.arg1 = x;
  196.     if (!matherr (&xcpt)) {
  197.         fprintf (stderr, "%s: DOMAIN error\n", funcname);
  198.         errno = EDOM;
  199.         xcpt.retval = 0.0;
  200.     }
  201.     } else {
  202.     m = frexp (x, &k);
  203. #ifdef OLD
  204.     u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
  205.     for (count = 0, y = u; count < ITERATIONS; count++) {
  206.         y = 0.5 * (y + (m / y));
  207.     }
  208.     if ((kmod2 = (k % 2)) < 0) {
  209.         y /= SQRT2;
  210.     } else if (kmod2 > 0) {
  211.         y *= SQRT2;
  212.     }
  213.     bugfix = 2;
  214.     xcpt.retval = ldexp (y, k/bugfix);
  215. #else
  216.     u = (((P4 * m + P3) * m + P2) * m + P1) * m + P0;
  217.     u = ldexp((u + (m / u)), -1);    /* Heron's iteration */
  218.     u += m / u;            /* and a part of the second one */
  219.     if (k & 1) {
  220.         u *= SQRT2;
  221.     }
  222.     /* 
  223.      * here we rely on the fact that -3/2 and (-3 >> 1)
  224.      * do give different results
  225.      */
  226.     xcpt.retval = ldexp (u, (k >> 1) - 1);
  227. #endif
  228.     }
  229.     DBUG_3 ("sqrtout", "result %le", xcpt.retval);
  230.     DBUG_RETURN (xcpt.retval);
  231. }
  232.