Print this page
5261 libm should stop using synonyms.h
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/Q/fmodl.c
+++ new/usr/src/lib/libm/common/Q/fmodl.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
|
↓ open down ↓ |
19 lines elided |
↑ open up ↑ |
20 20 */
21 21
22 22 /*
23 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 25 /*
26 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 27 * Use is subject to license terms.
28 28 */
29 29
30 -#pragma weak fmodl = __fmodl
30 +#pragma weak __fmodl = fmodl
31 31
32 32 #include "libm.h"
33 33
34 34 static const int
35 35 is = -0x7fffffff - 1,
36 36 im = 0x0000ffff,
37 37 iu = 0x00010000;
38 38
39 39 static const long double
40 40 zero = 0.0L,
41 41 one = 1.0L;
42 42
43 43 #ifdef __LITTLE_ENDIAN
44 44 #define __H0(x) *(3 + (int *) &x)
45 45 #define __H1(x) *(2 + (int *) &x)
46 46 #define __H2(x) *(1 + (int *) &x)
47 47 #define __H3(x) *(0 + (int *) &x)
48 48 #else
49 49 #define __H0(x) *(0 + (int *) &x)
50 50 #define __H1(x) *(1 + (int *) &x)
51 51 #define __H2(x) *(2 + (int *) &x)
52 52 #define __H3(x) *(3 + (int *) &x)
53 53 #endif
54 54
55 55 long double
56 56 fmodl(long double x, long double y) {
57 57 long double a, b;
58 58 int n, ix, iy, k, sx;
59 59 int hx;
60 60 int x0, y0, z0, carry;
61 61 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
62 62
63 63 hx = __H0(x);
64 64 x1 = __H1(x);
65 65 x2 = __H2(x);
66 66 x3 = __H3(x);
67 67 y0 = __H0(y);
68 68 y1 = __H1(y);
69 69 y2 = __H2(y);
70 70 y3 = __H3(y);
71 71
72 72 sx = hx & 0x80000000;
73 73 x0 = hx ^ sx;
74 74 y0 &= 0x7fffffff;
75 75
76 76 /* purge off exception values */
77 77 if (x0 >= 0x7fff0000 || /* !finitel(x) */
78 78 (y0 > 0x7fff0000) || (y0 == 0x7fff0000 && ((y1 | y2 | y3) != 0)) ||
79 79 (y0 | y1 | y2 | y3) == 0) /* isnanl(y) || y = 0 */
80 80 return ((x * y) / (x * y));
81 81 a = fabsl(x);
82 82 b = fabsl(y);
83 83 if (a <= b) {
84 84 if (a < b)
85 85 return (x);
86 86 else
87 87 return (zero * x);
88 88 }
89 89 /* determine ix = ilogbl(x) */
90 90 if (x0 < iu) { /* subnormal x */
91 91 ix = -16382;
92 92 while (x0 == 0) {
93 93 ix -= 16;
94 94 x0 = x1 >> 16;
95 95 x1 = (x1 << 16) | (x2 >> 16);
96 96 x2 = (x2 << 16) | (x3 >> 16);
97 97 x3 = (x3 << 16);
98 98 }
99 99 while (x0 < iu) {
100 100 ix -= 1;
101 101 x0 = (x0 << 1) | (x1 >> 31);
102 102 x1 = (x1 << 1) | (x2 >> 31);
103 103 x2 = (x2 << 1) | (x3 >> 31);
104 104 x3 <<= 1;
105 105 }
106 106 } else {
107 107 ix = (x0 >> 16) - 16383;
108 108 x0 = iu | (x0 & im);
109 109 }
110 110
111 111 /* determine iy = ilogbl(y) */
112 112 if (y0 < iu) { /* subnormal y */
113 113 iy = -16382;
114 114 while (y0 == 0) {
115 115 iy -= 16;
116 116 y0 = y1 >> 16;
117 117 y1 = (y1 << 16) | (y2 >> 16);
118 118 y2 = (y2 << 16) | (y3 >> 16);
119 119 y3 = (y3 << 16);
120 120 }
121 121 while (y0 < iu) {
122 122 iy -= 1;
123 123 y0 = (y0 << 1) | (y1 >> 31);
124 124 y1 = (y1 << 1) | (y2 >> 31);
125 125 y2 = (y2 << 1) | (y3 >> 31);
126 126 y3 <<= 1;
127 127 }
128 128 } else {
129 129 iy = (y0 >> 16) - 16383;
130 130 y0 = iu | (y0 & im);
131 131 }
132 132
133 133 /* fix point fmod */
134 134 n = ix - iy;
135 135 while (n--) {
136 136 while (x0 == 0 && n >= 16) {
137 137 n -= 16;
138 138 x0 = x1 >> 16;
139 139 x1 = (x1 << 16) | (x2 >> 16);
140 140 x2 = (x2 << 16) | (x3 >> 16);
141 141 x3 = (x3 << 16);
142 142 }
143 143 while (x0 < iu && n >= 1) {
144 144 n -= 1;
145 145 x0 = (x0 << 1) | (x1 >> 31);
146 146 x1 = (x1 << 1) | (x2 >> 31);
147 147 x2 = (x2 << 1) | (x3 >> 31);
148 148 x3 = (x3 << 1);
149 149 }
150 150 carry = 0;
151 151 z3 = x3 - y3;
152 152 carry = (z3 > x3);
153 153 if (carry == 0) {
154 154 z2 = x2 - y2;
155 155 carry = (z2 > x2);
156 156 } else {
157 157 z2 = x2 - y2 - 1;
158 158 carry = (z2 >= x2);
159 159 }
160 160 if (carry == 0) {
161 161 z1 = x1 - y1;
162 162 carry = (z1 > x1);
163 163 } else {
164 164 z1 = x1 - y1 - 1;
165 165 carry = (z1 >= x1);
166 166 }
167 167 z0 = x0 - y0 - carry;
168 168 if (z0 < 0) { /* double x */
169 169 x0 = x0 + x0 + ((x1 & is) != 0);
170 170 x1 = x1 + x1 + ((x2 & is) != 0);
171 171 x2 = x2 + x2 + ((x3 & is) != 0);
172 172 x3 = x3 + x3;
173 173 } else {
174 174 if (z0 == 0) {
175 175 if ((z1 | z2 | z3) == 0) { /* 0: done */
176 176 __H0(a) = hx & is;
177 177 __H1(a) = __H2(a) = __H3(a) = 0;
178 178 return (a);
179 179 }
180 180 }
181 181 /* x = z << 1 */
182 182 z0 = z0 + z0 + ((z1 & is) != 0);
183 183 z1 = z1 + z1 + ((z2 & is) != 0);
184 184 z2 = z2 + z2 + ((z3 & is) != 0);
185 185 z3 = z3 + z3;
186 186 x0 = z0;
187 187 x1 = z1;
188 188 x2 = z2;
189 189 x3 = z3;
190 190 }
191 191 }
192 192
193 193 carry = 0;
194 194 z3 = x3 - y3;
195 195 carry = (z3 > x3);
196 196 if (carry == 0) {
197 197 z2 = x2 - y2;
198 198 carry = (z2 > x2);
199 199 } else {
200 200 z2 = x2 - y2 - 1;
201 201 carry = (z2 >= x2);
202 202 }
203 203 if (carry == 0) {
204 204 z1 = x1 - y1;
205 205 carry = (z1 > x1);
206 206 } else {
207 207 z1 = x1 - y1 - 1;
208 208 carry = (z1 >= x1);
209 209 }
210 210 z0 = x0 - y0 - carry;
211 211 if (z0 >= 0) {
212 212 x0 = z0;
213 213 x1 = z1;
214 214 x2 = z2;
215 215 x3 = z3;
216 216 }
217 217 /* convert back to floating value and restore the sign */
218 218 if ((x0 | x1 | x2 | x3) == 0) {
219 219 __H0(a) = hx & is;
220 220 __H1(a) = __H2(a) = __H3(a) = 0;
221 221 return (a);
222 222 }
223 223 while (x0 < iu) {
224 224 if (x0 == 0) {
225 225 iy -= 16;
226 226 x0 = x1 >> 16;
227 227 x1 = (x1 << 16) | (x2 >> 16);
228 228 x2 = (x2 << 16) | (x3 >> 16);
229 229 x3 = (x3 << 16);
230 230 } else {
231 231 x0 = x0 + x0 + ((x1 & is) != 0);
232 232 x1 = x1 + x1 + ((x2 & is) != 0);
233 233 x2 = x2 + x2 + ((x3 & is) != 0);
234 234 x3 = x3 + x3;
235 235 iy -= 1;
236 236 }
237 237 }
238 238
239 239 /* normalize output */
240 240 if (iy >= -16382) {
241 241 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
242 242 __H1(a) = x1;
243 243 __H2(a) = x2;
244 244 __H3(a) = x3;
245 245 } else { /* subnormal output */
246 246 n = -16382 - iy;
247 247 k = n & 31;
248 248 if (k != 0) {
249 249 if (k <= 16) {
250 250 x3 = (x2 << (32 - k)) | (x3 >> k);
251 251 x2 = (x1 << (32 - k)) | (x2 >> k);
252 252 x1 = (x0 << (32 - k)) | (x1 >> k);
253 253 x0 >>= k;
254 254 } else {
255 255 x3 = (x2 << (32 - k)) | (x3 >> k);
256 256 x2 = (x1 << (32 - k)) | (x2 >> k);
257 257 x1 = (x0 << (32 - k)) | (x1 >> k);
258 258 x0 = 0;
259 259 }
260 260 }
261 261 while (n >= 32) {
262 262 n -= 32;
263 263 x3 = x2;
264 264 x2 = x1;
265 265 x1 = x0;
266 266 x0 = 0;
267 267 }
268 268 __H0(a) = x0 | sx;
269 269 __H1(a) = x1;
270 270 __H2(a) = x2;
271 271 __H3(a) = x3;
272 272 a *= one;
273 273 }
274 274 return (a);
275 275 }
|
↓ open down ↓ |
235 lines elided |
↑ open up ↑ |
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX