Print this page
5261 libm should stop using synonyms.h
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/m9x/remquol.c
+++ new/usr/src/lib/libm/common/m9x/remquol.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 remquol = __remquol
30 +#pragma weak __remquol = remquol
31 31
32 32 #include "libm.h"
33 -#include "libm_synonyms.h"
34 33 #if defined(__SUNPRO_C)
35 34 #include <sunmath.h> /* fabsl */
36 35 #endif
37 36 /* INDENT OFF */
38 37 static const int
39 38 is = -0x7fffffff - 1,
40 39 im = 0x0000ffff,
41 40 iu = 0x00010000;
42 41
43 42 static const long double zero = 0.0L, one = 1.0L;
44 43 /* INDENT ON */
45 44
46 45 #if defined(__sparc)
47 46 #define __H0(x) ((int *) &x)[0]
48 47 #define __H1(x) ((int *) &x)[1]
49 48 #define __H2(x) ((int *) &x)[2]
50 49 #define __H3(x) ((int *) &x)[3]
51 50 #else
52 51 #error Unsupported architecture
53 52 #endif
54 53
55 54 /*
56 55 * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered
57 56 */
58 57 static long double
59 58 fmodquol(long double x, long double y, int *quo) {
60 59 long double a, b;
61 60 int n, ix, iy, k, sx, sq, m;
62 61 int hx;
63 62 int x0, y0, z0, carry;
64 63 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
65 64
66 65 hx = __H0(x);
67 66 x1 = __H1(x);
68 67 x2 = __H2(x);
69 68 x3 = __H3(x);
70 69 y0 = __H0(y);
71 70 y1 = __H1(y);
72 71 y2 = __H2(y);
73 72 y3 = __H3(y);
74 73
75 74 sx = hx & is;
76 75 sq = (hx ^ y0) & is;
77 76 x0 = hx ^ sx;
78 77 y0 &= ~0x80000000;
79 78
80 79 a = fabsl(x);
81 80 b = fabsl(y);
82 81 if (a <= b) {
83 82 if (a < b)
84 83 return (x);
85 84 else {
86 85 *quo = 1 + (sq >> 30);
87 86 return (zero * x);
88 87 }
89 88 }
90 89 /* determine ix = ilogbl(x) */
91 90 if (x0 < iu) { /* subnormal x */
92 91 ix = 0;
93 92 ix = -16382;
94 93 while (x0 == 0) {
95 94 ix -= 16;
96 95 x0 = x1 >> 16;
97 96 x1 = (x1 << 16) | (x2 >> 16);
98 97 x2 = (x2 << 16) | (x3 >> 16);
99 98 x3 = (x3 << 16);
100 99 }
101 100 while (x0 < iu) {
102 101 ix -= 1;
103 102 x0 = (x0 << 1) | (x1 >> 31);
104 103 x1 = (x1 << 1) | (x2 >> 31);
105 104 x2 = (x2 << 1) | (x3 >> 31);
106 105 x3 <<= 1;
107 106 }
108 107 } else {
109 108 ix = (x0 >> 16) - 16383;
110 109 x0 = iu | (x0 & im);
111 110 }
112 111
113 112 /* determine iy = ilogbl(y) */
114 113 if (y0 < iu) { /* subnormal y */
115 114 iy = -16382;
116 115 while (y0 == 0) {
117 116 iy -= 16;
118 117 y0 = y1 >> 16;
119 118 y1 = (y1 << 16) | (y2 >> 16);
120 119 y2 = (y2 << 16) | (y3 >> 16);
121 120 y3 = (y3 << 16);
122 121 }
123 122 while (y0 < iu) {
124 123 iy -= 1;
125 124 y0 = (y0 << 1) | (y1 >> 31);
126 125 y1 = (y1 << 1) | (y2 >> 31);
127 126 y2 = (y2 << 1) | (y3 >> 31);
128 127 y3 <<= 1;
129 128 }
130 129 } else {
131 130 iy = (y0 >> 16) - 16383;
132 131 y0 = iu | (y0 & im);
133 132 }
134 133
135 134
136 135 /* fix point fmod */
137 136 n = ix - iy;
138 137 m = 0;
139 138 while (n--) {
140 139 while (x0 == 0 && n >= 16) {
141 140 m <<= 16;
142 141 n -= 16;
143 142 x0 = x1 >> 16;
144 143 x1 = (x1 << 16) | (x2 >> 16);
145 144 x2 = (x2 << 16) | (x3 >> 16);
146 145 x3 = (x3 << 16);
147 146 }
148 147 while (x0 < iu && n >= 1) {
149 148 m += m;
150 149 n -= 1;
151 150 x0 = (x0 << 1) | (x1 >> 31);
152 151 x1 = (x1 << 1) | (x2 >> 31);
153 152 x2 = (x2 << 1) | (x3 >> 31);
154 153 x3 = (x3 << 1);
155 154 }
156 155 carry = 0;
157 156 z3 = x3 - y3;
158 157 carry = z3 > x3;
159 158 if (carry == 0) {
160 159 z2 = x2 - y2;
161 160 carry = z2 > x2;
162 161 } else {
163 162 z2 = x2 - y2 - 1;
164 163 carry = z2 >= x2;
165 164 }
166 165 if (carry == 0) {
167 166 z1 = x1 - y1;
168 167 carry = z1 > x1;
169 168 } else {
170 169 z1 = x1 - y1 - 1;
171 170 carry = z1 >= x1;
172 171 }
173 172 z0 = x0 - y0 - carry;
174 173 if (z0 < 0) { /* double x */
175 174 x0 = x0 + x0 + ((x1 & is) != 0);
176 175 x1 = x1 + x1 + ((x2 & is) != 0);
177 176 x2 = x2 + x2 + ((x3 & is) != 0);
178 177 x3 = x3 + x3;
179 178 m += m;
180 179 } else {
181 180 m += 1;
182 181 if (z0 == 0) {
183 182 if ((z1 | z2 | z3) == 0) {
184 183 /* 0: we are done */
185 184 if (n < 31)
186 185 m <<= (1 + n);
187 186 else
188 187 m = 0;
189 188 m &= ~0x80000000;
190 189 *quo = sq >= 0 ? m : -m;
191 190 __H0(a) = hx & is;
192 191 __H1(a) = __H2(a) = __H3(a) = 0;
193 192 return (a);
194 193 }
195 194 }
196 195 /* x = z << 1 */
197 196 z0 = z0 + z0 + ((z1 & is) != 0);
198 197 z1 = z1 + z1 + ((z2 & is) != 0);
199 198 z2 = z2 + z2 + ((z3 & is) != 0);
200 199 z3 = z3 + z3;
201 200 x0 = z0;
202 201 x1 = z1;
203 202 x2 = z2;
204 203 x3 = z3;
205 204 m += m;
206 205 }
207 206 }
208 207 carry = 0;
209 208 z3 = x3 - y3;
210 209 carry = z3 > x3;
211 210 if (carry == 0) {
212 211 z2 = x2 - y2;
213 212 carry = z2 > x2;
214 213 } else {
215 214 z2 = x2 - y2 - 1;
216 215 carry = z2 >= x2;
217 216 }
218 217 if (carry == 0) {
219 218 z1 = x1 - y1;
220 219 carry = z1 > x1;
221 220 } else {
222 221 z1 = x1 - y1 - 1;
223 222 carry = z1 >= x1;
224 223 }
225 224 z0 = x0 - y0 - carry;
226 225 if (z0 >= 0) {
227 226 x0 = z0;
228 227 x1 = z1;
229 228 x2 = z2;
230 229 x3 = z3;
231 230 m += 1;
232 231 }
233 232 m &= ~0x80000000;
234 233 *quo = sq >= 0 ? m : -m;
235 234
236 235 /* convert back to floating value and restore the sign */
237 236 if ((x0 | x1 | x2 | x3) == 0) {
238 237 __H0(a) = hx & is;
239 238 __H1(a) = __H2(a) = __H3(a) = 0;
240 239 return (a);
241 240 }
242 241 while (x0 < iu) {
243 242 if (x0 == 0) {
244 243 iy -= 16;
245 244 x0 = x1 >> 16;
246 245 x1 = (x1 << 16) | (x2 >> 16);
247 246 x2 = (x2 << 16) | (x3 >> 16);
248 247 x3 = (x3 << 16);
249 248 } else {
250 249 x0 = x0 + x0 + ((x1 & is) != 0);
251 250 x1 = x1 + x1 + ((x2 & is) != 0);
252 251 x2 = x2 + x2 + ((x3 & is) != 0);
253 252 x3 = x3 + x3;
254 253 iy -= 1;
255 254 }
256 255 }
257 256
258 257 /* normalize output */
259 258 if (iy >= -16382) {
260 259 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
261 260 __H1(a) = x1;
262 261 __H2(a) = x2;
263 262 __H3(a) = x3;
264 263 } else { /* subnormal output */
265 264 n = -16382 - iy;
266 265 k = n & 31;
267 266 if (k <= 16) {
268 267 x3 = (x2 << (32 - k)) | (x3 >> k);
269 268 x2 = (x1 << (32 - k)) | (x2 >> k);
270 269 x1 = (x0 << (32 - k)) | (x1 >> k);
271 270 x0 >>= k;
272 271 } else {
273 272 x3 = (x2 << (32 - k)) | (x3 >> k);
274 273 x2 = (x1 << (32 - k)) | (x2 >> k);
275 274 x1 = (x0 << (32 - k)) | (x1 >> k);
276 275 x0 = 0;
277 276 }
278 277 while (n >= 32) {
279 278 n -= 32;
280 279 x3 = x2;
281 280 x2 = x1;
282 281 x1 = x0;
283 282 x0 = 0;
284 283 }
285 284 __H0(a) = x0 | sx;
286 285 __H1(a) = x1;
287 286 __H2(a) = x2;
288 287 __H3(a) = x3;
289 288 a *= one;
290 289 }
291 290 return (a);
292 291 }
293 292
294 293 long double
295 294 remquol(long double x, long double y, int *quo) {
296 295 int hx, hy, sx, sq;
297 296 long double v;
298 297
299 298 hx = __H0(x); /* high word of x */
300 299 hy = __H0(y); /* high word of y */
301 300 sx = hx & is; /* sign of x */
302 301 sq = (hx ^ hy) & is; /* sign of x/y */
303 302 hx ^= sx; /* |x| */
304 303 hy &= ~0x80000000;
305 304
306 305 /* purge off exception values */
307 306 *quo = 0;
308 307 /* y=0, y is NaN, x is NaN or inf */
309 308 if (y == 0.0L || y != y || hx >= 0x7fff0000)
310 309 return ((x * y) / (x * y));
311 310
312 311 y = fabsl(y);
313 312 x = fabsl(x);
314 313 if (hy <= 0x7ffdffff) {
315 314 x = fmodquol(x, y + y, quo);
316 315 *quo = ((*quo) & 0x3fffffff) << 1;
317 316 }
318 317 if (hy < 0x00020000) {
319 318 if (x + x > y) {
320 319 *quo += 1;
321 320 if (x == y)
322 321 x = zero;
323 322 else
324 323 x -= y;
325 324 if (x + x >= y) {
326 325 x -= y;
327 326 *quo += 1;
328 327 }
329 328 }
330 329 } else {
331 330 v = 0.5L * y;
332 331 if (x > v) {
333 332 *quo += 1;
334 333 if (x == y)
335 334 x = zero;
336 335 else
337 336 x -= y;
338 337 if (x >= v) {
339 338 x -= y;
340 339 *quo += 1;
341 340 }
342 341 }
343 342 }
344 343 if (sq != 0)
345 344 *quo = -(*quo);
346 345 return (sx == 0 ? x : -x);
347 346 }
|
↓ open down ↓ |
304 lines elided |
↑ open up ↑ |
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX