newlib/libm/mathfp/sf_sqrt.c

   1
   2 /* @(#)z_sqrtf.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  * Square Root
  11  *
  12  * Input:
  13  *   x - floating point value
  14  *
  15  * Output:
  16  *   square-root of x
  17  *
  18  * Description:
  19  *   This routine performs floating point square root.
  20  *
  21  *   The initial approximation is computed as
  22  *     y0 = 0.41731 + 0.59016 * f
  23  *   where f is a fraction such that x = f * 2^exp.
  24  *
  25  *   Three Newton iterations in the form of Heron's formula
  26  *   are then performed to obtain the final value:
  27  *     y[i] = (y[i-1] + f / y[i-1]) / 2, i = 1, 2, 3.
  28  *
  29  *****************************************************************/
  30
  31 #include "fdlibm.h"
  32 #include "zmath.h"
  33
  34 float
  35 sqrtf (float x)
  36 {
  37   float f, y;
  38   int exp, i, odd;
  39
  40   /* Check for special values. */
  41   switch (numtestf (x))
  42     {
  43       case NAN:
  44         errno = EDOM;
  45         return (x);
  46       case INF:
  47         if (isposf (x))
  48           {
  49             errno = EDOM;
  50             return (z_notanum_f.f);
  51           }
  52         else
  53           {
  54             errno = ERANGE;
  55             return (z_infinity_f.f);
  56           }
  57     }
  58
  59   /* Initial checks are performed here. */
  60   if (x == 0.0)
  61     return (0.0);
  62   if (x < 0)
  63     {
  64       errno = EDOM;
  65       return (z_notanum_f.f);
  66     }
  67
  68   /* Find the exponent and mantissa for the form x = f * 2^exp. */
  69   f = frexpf (x, &exp);
  70   odd = exp & 1;
  71
  72   /* Get the initial approximation. */
  73   y = 0.41731 + 0.59016 * f;
  74
  75   f *= 0.5;
  76   /* Calculate the remaining iterations. */
  77   for (i = 0; i < 2; ++i)
  78     y = y * 0.5 + f / y;
  79
  80   /* Calculate the final value. */
  81   if (odd)
  82     {
  83       y *= __SQRT_HALF;
  84       exp++;
  85     }
  86   exp >>= 1;
  87   y = ldexpf (y, exp);
  88
  89   return (y);
  90 }
  91
  92 #ifdef _DOUBLE_IS_32BITS
  93
  94 double sqrt (double x)
  95 {
  96   return (double) sqrtf ((float) x);
  97 }
  98
  99 #endif /* _DOUBLE_IS_32BITS */