sysdeps/ieee754/sqrt.c

   1 /*
   2  * Copyright (c) 1985 Regents of the University of California.
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms are permitted provided
   6  * that: (1) source distributions retain this entire copyright notice and
   7  * comment, and (2) distributions including binaries display the following
   8  * acknowledgement:  ``This product includes software developed by the
   9  * University of California, Berkeley and its contributors'' in the
  10  * documentation or other materials provided with the distribution and in
  11  * all advertising materials mentioning features or use of this software.
  12  * Neither the name of the University nor the names of its contributors may
  13  * be used to endorse or promote products derived from this software without
  14  * specific prior written permission.
  15  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
  16  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
  17  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  18  */
  19
  20 #include <ansidecl.h>
  21 #include <errno.h>
  22 #include <math.h>
  23
  24 /* Return the square root of X.  */
  25 double
  26 DEFUN (sqrt, (x), double x)
  27 {
  28   double q, s, b, r, t;
  29   CONST double zero = 0.0;
  30   int m, n, i;
  31
  32   /* sqrt (NaN) is NaN; sqrt (+-0) is +-0.  */
  33   if (__isnan (x) || x == zero)
  34     return x;
  35
  36   if (x < zero)
  37     return zero / zero;
  38
  39   /* sqrt (Inf) is Inf.  */
  40   if (__isinf (x))
  41     return x;
  42
  43   /* Scale X to [1,4).  */
  44   n = __logb (x);
  45   x = __scalb (x, -n);
  46   m = __logb (x);
  47   if (m != 0)
  48     /* Subnormal number.  */
  49     x = __scalb (x, -m);
  50
  51   m += n;
  52   n = m / 2;
  53
  54   if ((n + n) != m)
  55     {
  56       x *= 2;
  57       --m;
  58       n = m / 2;
  59     }
  60
  61   /* Generate sqrt (X) bit by bit (accumulating in Q).  */
  62   q = 1.0;
  63   s = 4.0;
  64   x -= 1.0;
  65   r = 1;
  66   for (i = 1; i <= 51; i++)
  67     {
  68       t = s + 1;
  69       x *= 4;
  70       r /= 2;
  71       if (t <= x)
  72         {
  73           s = t + t + 2, x -= t;
  74           q += r;
  75         }
  76       else
  77         s *= 2;
  78     }
  79
  80   /* Generate the last bit and determine the final rounding.  */
  81   r /= 2;
  82   x *= 4;
  83   if (x == zero)
  84     goto end;
  85   (void) (100 + r);             /* Trigger inexact flag.  */
  86   if (s < x)
  87     {
  88       q += r;
  89       x -= s;
  90       s += 2;
  91       s *= 2;
  92       x *= 4;
  93       t = (x - s) - 5;
  94       b = 1.0 + 3 * r / 4;
  95       if (b == 1.0)
  96         goto end;               /* B == 1: Round to zero.  */
  97       b = 1.0 + r / 4;
  98       if (b > 1.0)
  99         t = 1;                  /* B > 1: Round to +Inf.  */
 100       if (t >= 0)
 101         q += r;
 102     }                           /* Else round to nearest.  */
 103   else
 104     {
 105       s *= 2;
 106       x *= 4;
 107       t = (x - s) - 1;
 108       b = 1.0 + 3 * r / 4;
 109       if (b == 1.0)
 110         goto end;
 111       b = 1.0 + r / 4;
 112       if (b > 1.0)
 113         t = 1;
 114       if (t >= 0)
 115         q += r;
 116     }
 117
 118 end:
 119   return __scalb (q, n);
 120 }