Initial revision
[kopensolaris-gnu/glibc.git] / 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 <= k; 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 }