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