1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
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 |
17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
18 | */ |
19 | /************************************************************************/ |
20 | /* */ |
21 | /* MODULE_NAME:halfulp.c */ |
22 | /* */ |
23 | /* FUNCTIONS:halfulp */ |
24 | /* FILES NEEDED: mydefs.h dla.h endian.h */ |
25 | /* uroot.c */ |
26 | /* */ |
27 | /*Routine halfulp(double x, double y) computes x^y where result does */ |
28 | /*not need rounding. If the result is closer to 0 than can be */ |
29 | /*represented it returns 0. */ |
30 | /* In the following cases the function does not compute anything */ |
31 | /*and returns a negative number: */ |
32 | /*1. if the result needs rounding, */ |
33 | /*2. if y is outside the interval [0, 2^20-1], */ |
34 | /*3. if x can be represented by x=2**n for some integer n. */ |
35 | /************************************************************************/ |
36 | |
37 | #include "endian.h" |
38 | #include "mydefs.h" |
39 | #include <dla.h> |
40 | #include <math_private.h> |
41 | |
42 | #ifndef SECTION |
43 | # define SECTION |
44 | #endif |
45 | |
46 | static const int4 tab54[32] = { |
47 | 262143, 11585, 1782, 511, 210, 107, 63, 42, |
48 | 30, 22, 17, 14, 12, 10, 9, 7, |
49 | 7, 6, 5, 5, 5, 4, 4, 4, |
50 | 3, 3, 3, 3, 3, 3, 3, 3 |
51 | }; |
52 | |
53 | |
54 | double |
55 | SECTION |
56 | __halfulp (double x, double y) |
57 | { |
58 | mynumber v; |
59 | double z, u, uu; |
60 | #ifndef DLA_FMS |
61 | double j1, j2, j3, j4, j5; |
62 | #endif |
63 | int4 k, l, m, n; |
64 | if (y <= 0) /*if power is negative or zero */ |
65 | { |
66 | v.x = y; |
67 | if (v.i[LOW_HALF] != 0) |
68 | return -10.0; |
69 | v.x = x; |
70 | if (v.i[LOW_HALF] != 0) |
71 | return -10.0; |
72 | if ((v.i[HIGH_HALF] & 0x000fffff) != 0) |
73 | return -10; /* if x =2 ^ n */ |
74 | k = ((v.i[HIGH_HALF] & 0x7fffffff) >> 20) - 1023; /* find this n */ |
75 | z = (double) k; |
76 | return (z * y == -1075.0) ? 0 : -10.0; |
77 | } |
78 | /* if y > 0 */ |
79 | v.x = y; |
80 | if (v.i[LOW_HALF] != 0) |
81 | return -10.0; |
82 | |
83 | v.x = x; |
84 | /* case where x = 2**n for some integer n */ |
85 | if (((v.i[HIGH_HALF] & 0x000fffff) | v.i[LOW_HALF]) == 0) |
86 | { |
87 | k = (v.i[HIGH_HALF] >> 20) - 1023; |
88 | return (((double) k) * y == -1075.0) ? 0 : -10.0; |
89 | } |
90 | |
91 | v.x = y; |
92 | k = v.i[HIGH_HALF]; |
93 | m = k << 12; |
94 | l = 0; |
95 | while (m) |
96 | { |
97 | m = m << 1; l++; |
98 | } |
99 | n = (k & 0x000fffff) | 0x00100000; |
100 | n = n >> (20 - l); /* n is the odd integer of y */ |
101 | k = ((k >> 20) - 1023) - l; /* y = n*2**k */ |
102 | if (k > 5) |
103 | return -10.0; |
104 | if (k > 0) |
105 | for (; k > 0; k--) |
106 | n *= 2; |
107 | if (n > 34) |
108 | return -10.0; |
109 | k = -k; |
110 | if (k > 5) |
111 | return -10.0; |
112 | |
113 | /* now treat x */ |
114 | while (k > 0) |
115 | { |
116 | z = __ieee754_sqrt (x); |
117 | EMULV (z, z, u, uu, j1, j2, j3, j4, j5); |
118 | if (((u - x) + uu) != 0) |
119 | break; |
120 | x = z; |
121 | k--; |
122 | } |
123 | if (k) |
124 | return -10.0; |
125 | |
126 | /* it is impossible that n == 2, so the mantissa of x must be short */ |
127 | |
128 | v.x = x; |
129 | if (v.i[LOW_HALF]) |
130 | return -10.0; |
131 | k = v.i[HIGH_HALF]; |
132 | m = k << 12; |
133 | l = 0; |
134 | while (m) |
135 | { |
136 | m = m << 1; l++; |
137 | } |
138 | m = (k & 0x000fffff) | 0x00100000; |
139 | m = m >> (20 - l); /* m is the odd integer of x */ |
140 | |
141 | /* now check whether the length of m**n is at most 54 bits */ |
142 | |
143 | if (m > tab54[n - 3]) |
144 | return -10.0; |
145 | |
146 | /* yes, it is - now compute x**n by simple multiplications */ |
147 | |
148 | u = x; |
149 | for (k = 1; k < n; k++) |
150 | u = u * x; |
151 | return u; |
152 | } |
153 | |