newlib/libm/mathfp/s_sqrt.c

   1
   2 /* @(#)z_sqrt.c 1.0 98/08/13 */
   3 /*****************************************************************
   4  * The following routines are coded directly from the algorithms
   5  * and coefficients given in "Software Manual for the Elementary
   6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
   7  * Hall, 1980.
   8  *****************************************************************/
   9
  10 /*
  11 FUNCTION
  12         <<sqrt>>, <<sqrtf>>---positive square root
  13
  14 INDEX
  15         sqrt
  16 INDEX
  17         sqrtf
  18
  19 SYNOPSIS
  20         #include <math.h>
  21         double sqrt(double <[x]>);
  22         float  sqrtf(float <[x]>);
  23
  24 DESCRIPTION
  25         <<sqrt>> computes the positive square root of the argument.
  26
  27 RETURNS
  28         On success, the square root is returned. If <[x]> is real and
  29         positive, then the result is positive.  If <[x]> is real and
  30         negative, the global value <<errno>> is set to <<EDOM>> (domain error).
  31
  32
  33 PORTABILITY
  34         <<sqrt>> is ANSI C.  <<sqrtf>> is an extension.
  35 */
  36
  37 /******************************************************************
  38  * Square Root
  39  *
  40  * Input:
  41  *   x - floating point value
  42  *
  43  * Output:
  44  *   square-root of x
  45  *
  46  * Description:
  47  *   This routine performs floating point square root.
  48  *
  49  *   The initial approximation is computed as
  50  *     y0 = 0.41731 + 0.59016 * f
  51  *   where f is a fraction such that x = f * 2^exp.
  52  *
  53  *   Three Newton iterations in the form of Heron's formula
  54  *   are then performed to obtain the final value:
  55  *     y[i] = (y[i-1] + f / y[i-1]) / 2, i = 1, 2, 3.
  56  *
  57  *****************************************************************/
  58
  59 #include "fdlibm.h"
  60 #include "zmath.h"
  61
  62 #ifndef _DOUBLE_IS_32BITS
  63
  64 double
  65 sqrt (double x)
  66 {
  67   double f, y;
  68   int exp, i, odd;
  69
  70   /* Check for special values. */
  71   switch (numtest (x))
  72     {
  73       case NAN:
  74         errno = EDOM;
  75         return (x);
  76       case INF:
  77         if (ispos (x))
  78           {
  79             errno = EDOM;
  80             return (z_notanum.d);
  81           }
  82         else
  83           {
  84             errno = ERANGE;
  85             return (z_infinity.d);
  86           }
  87     }
  88
  89   /* Initial checks are performed here. */
  90   if (x == 0.0)
  91     return (0.0);
  92   if (x < 0)
  93     {
  94       errno = EDOM;
  95       return (z_notanum.d);
  96     }
  97
  98   /* Find the exponent and mantissa for the form x = f * 2^exp. */
  99   f = frexp (x, &exp);
 100
 101   odd = exp & 1;
 102
 103   /* Get the initial approximation. */
 104   y = 0.41731 + 0.59016 * f;
 105
 106   f /= 2.0;
 107   /* Calculate the remaining iterations. */
 108   for (i = 0; i < 3; ++i)
 109     y = y / 2.0 + f / y;
 110
 111   /* Calculate the final value. */
 112   if (odd)
 113     {
 114       y *= __SQRT_HALF;
 115       exp++;
 116     }
 117   exp >>= 1;
 118   y = ldexp (y, exp);
 119
 120   return (y);
 121 }
 122
 123 #endif /* _DOUBLE_IS_32BITS */