]>
Commit | Line | Data |
---|---|---|
f964490f RM |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
3 | * written by International Business Machines Corp. | |
568035b7 | 4 | * Copyright (C) 2001-2013 Free Software Foundation, Inc. |
f964490f RM |
5 | * |
6 | * This program is free software; you can redistribute it and/or modify | |
7 | * it under the terms of the GNU Lesser General Public License as published by | |
8 | * the Free Software Foundation; either version 2.1 of the License, or | |
9 | * (at your option) any later version. | |
10 | * | |
11 | * This program is distributed in the hope that it will be useful, | |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
14 | * GNU Lesser General Public License for more details. | |
15 | * | |
16 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
f964490f RM |
18 | */ |
19 | /*********************************************************************/ | |
20 | /* MODULE_NAME: uroot.c */ | |
21 | /* */ | |
22 | /* FUNCTION: usqrt */ | |
23 | /* */ | |
24 | /* FILES NEEDED: dla.h endian.h mydefs.h uroot.h */ | |
25 | /* uroot.tbl */ | |
26 | /* */ | |
27 | /* An ultimate sqrt routine. Given an IEEE double machine number x */ | |
28 | /* it computes the correctly rounded (to nearest) value of square */ | |
29 | /* root of x. */ | |
30 | /* Assumption: Machine arithmetic operations are performed in */ | |
31 | /* round to nearest mode of IEEE 754 standard. */ | |
32 | /* */ | |
33 | /*********************************************************************/ | |
34 | ||
35 | #include <math_private.h> | |
36 | ||
37 | typedef unsigned int int4; | |
38 | typedef union {int4 i[4]; long double x; double d[2]; } mynumber; | |
39 | ||
40 | static const mynumber | |
41 | t512 = {{0x5ff00000, 0x00000000, 0x00000000, 0x00000000 }}, /* 2^512 */ | |
42 | tm256 = {{0x2ff00000, 0x00000000, 0x00000000, 0x00000000 }}; /* 2^-256 */ | |
43 | static const double | |
44 | two54 = 1.80143985094819840000e+16, /* 0x4350000000000000 */ | |
45 | twom54 = 5.55111512312578270212e-17; /* 0x3C90000000000000 */ | |
46 | ||
47 | /*********************************************************************/ | |
48 | /* An ultimate sqrt routine. Given an IEEE double machine number x */ | |
49 | /* it computes the correctly rounded (to nearest) value of square */ | |
50 | /* root of x. */ | |
51 | /*********************************************************************/ | |
0ac5ae23 | 52 | long double __ieee754_sqrtl(long double x) |
f964490f RM |
53 | { |
54 | static const long double big = 134217728.0, big1 = 134217729.0; | |
55 | long double t,s,i; | |
56 | mynumber a,c; | |
57 | int4 k, l, m; | |
58 | int n; | |
59 | double d; | |
60 | ||
61 | a.x=x; | |
62 | k=a.i[0] & 0x7fffffff; | |
63 | /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ | |
64 | if (k>0x000fffff && k<0x7ff00000) { | |
65 | if (x < 0) return (big1-big1)/(big-big); | |
66 | l = (k&0x001fffff)|0x3fe00000; | |
67 | if (((a.i[2] & 0x7fffffff) | a.i[3]) != 0) { | |
68 | n = (int) ((l - k) * 2) >> 21; | |
69 | m = (a.i[2] >> 20) & 0x7ff; | |
70 | if (m == 0) { | |
71 | a.d[1] *= two54; | |
72 | m = ((a.i[2] >> 20) & 0x7ff) - 54; | |
73 | } | |
74 | m += n; | |
da93d214 | 75 | if ((int) m > 0) |
f964490f | 76 | a.i[2] = (a.i[2] & 0x800fffff) | (m << 20); |
da93d214 | 77 | else if ((int) m <= -54) { |
f964490f RM |
78 | a.i[2] &= 0x80000000; |
79 | a.i[3] = 0; | |
80 | } else { | |
81 | m += 54; | |
82 | a.i[2] = (a.i[2] & 0x800fffff) | (m << 20); | |
83 | a.d[1] *= twom54; | |
84 | } | |
85 | } | |
86 | a.i[0] = l; | |
87 | s = a.x; | |
88 | d = __ieee754_sqrt (a.d[0]); | |
89 | c.i[0] = 0x20000000+((k&0x7fe00000)>>1); | |
90 | c.i[1] = 0; | |
91 | c.i[2] = 0; | |
92 | c.i[3] = 0; | |
93 | i = d; | |
94 | t = 0.5L * (i + s / i); | |
95 | i = 0.5L * (t + s / t); | |
96 | return c.x * i; | |
97 | } | |
98 | else { | |
99 | if (k>=0x7ff00000) { | |
100 | if (a.i[0] == 0xfff00000 && a.i[1] == 0) | |
101 | return (big1-big1)/(big-big); /* sqrt (-Inf) = NaN. */ | |
102 | return x; /* sqrt (NaN) = NaN, sqrt (+Inf) = +Inf. */ | |
103 | } | |
104 | if (x == 0) return x; | |
105 | if (x < 0) return (big1-big1)/(big-big); | |
106 | return tm256.x*__ieee754_sqrtl(x*t512.x); | |
107 | } | |
108 | } | |
0ac5ae23 | 109 | strong_alias (__ieee754_sqrtl, __sqrtl_finite) |