5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak __remquol = remquol
31
32 #include "libm.h"
33 #if defined(__SUNPRO_C)
34 #include <sunmath.h> /* fabsl */
35 #endif
36 /* INDENT OFF */
37 static const int
38 is = -0x7fffffff - 1,
39 im = 0x0000ffff,
40 iu = 0x00010000;
41
42 static const long double zero = 0.0L, one = 1.0L;
43 /* INDENT ON */
44
45 #if defined(__sparc)
46 #define __H0(x) ((int *) &x)[0]
47 #define __H1(x) ((int *) &x)[1]
48 #define __H2(x) ((int *) &x)[2]
49 #define __H3(x) ((int *) &x)[3]
50 #else
51 #error Unsupported architecture
52 #endif
53
54 /*
55 * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered
56 */
57 static long double
58 fmodquol(long double x, long double y, int *quo) {
59 long double a, b;
60 int n, ix, iy, k, sx, sq, m;
61 int hx;
62 int x0, y0, z0, carry;
63 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
64
65 hx = __H0(x);
66 x1 = __H1(x);
67 x2 = __H2(x);
68 x3 = __H3(x);
69 y0 = __H0(y);
70 y1 = __H1(y);
71 y2 = __H2(y);
72 y3 = __H3(y);
73
74 sx = hx & is;
75 sq = (hx ^ y0) & is;
76 x0 = hx ^ sx;
77 y0 &= ~0x80000000;
78
79 a = fabsl(x);
80 b = fabsl(y);
81 if (a <= b) {
82 if (a < b)
83 return (x);
84 else {
85 *quo = 1 + (sq >> 30);
86 return (zero * x);
87 }
88 }
89 /* determine ix = ilogbl(x) */
90 if (x0 < iu) { /* subnormal x */
91 ix = 0;
92 ix = -16382;
93 while (x0 == 0) {
94 ix -= 16;
95 x0 = x1 >> 16;
96 x1 = (x1 << 16) | (x2 >> 16);
97 x2 = (x2 << 16) | (x3 >> 16);
98 x3 = (x3 << 16);
99 }
100 while (x0 < iu) {
101 ix -= 1;
102 x0 = (x0 << 1) | (x1 >> 31);
103 x1 = (x1 << 1) | (x2 >> 31);
104 x2 = (x2 << 1) | (x3 >> 31);
105 x3 <<= 1;
106 }
107 } else {
108 ix = (x0 >> 16) - 16383;
109 x0 = iu | (x0 & im);
110 }
111
112 /* determine iy = ilogbl(y) */
113 if (y0 < iu) { /* subnormal y */
114 iy = -16382;
115 while (y0 == 0) {
116 iy -= 16;
117 y0 = y1 >> 16;
118 y1 = (y1 << 16) | (y2 >> 16);
119 y2 = (y2 << 16) | (y3 >> 16);
120 y3 = (y3 << 16);
121 }
122 while (y0 < iu) {
123 iy -= 1;
124 y0 = (y0 << 1) | (y1 >> 31);
125 y1 = (y1 << 1) | (y2 >> 31);
126 y2 = (y2 << 1) | (y3 >> 31);
127 y3 <<= 1;
128 }
129 } else {
130 iy = (y0 >> 16) - 16383;
131 y0 = iu | (y0 & im);
132 }
133
134
135 /* fix point fmod */
136 n = ix - iy;
137 m = 0;
138 while (n--) {
139 while (x0 == 0 && n >= 16) {
140 m <<= 16;
141 n -= 16;
142 x0 = x1 >> 16;
143 x1 = (x1 << 16) | (x2 >> 16);
144 x2 = (x2 << 16) | (x3 >> 16);
145 x3 = (x3 << 16);
146 }
147 while (x0 < iu && n >= 1) {
148 m += m;
149 n -= 1;
150 x0 = (x0 << 1) | (x1 >> 31);
151 x1 = (x1 << 1) | (x2 >> 31);
152 x2 = (x2 << 1) | (x3 >> 31);
153 x3 = (x3 << 1);
154 }
155 carry = 0;
156 z3 = x3 - y3;
157 carry = z3 > x3;
158 if (carry == 0) {
159 z2 = x2 - y2;
160 carry = z2 > x2;
161 } else {
162 z2 = x2 - y2 - 1;
163 carry = z2 >= x2;
164 }
165 if (carry == 0) {
166 z1 = x1 - y1;
167 carry = z1 > x1;
168 } else {
169 z1 = x1 - y1 - 1;
170 carry = z1 >= x1;
171 }
172 z0 = x0 - y0 - carry;
173 if (z0 < 0) { /* double x */
174 x0 = x0 + x0 + ((x1 & is) != 0);
175 x1 = x1 + x1 + ((x2 & is) != 0);
176 x2 = x2 + x2 + ((x3 & is) != 0);
177 x3 = x3 + x3;
178 m += m;
179 } else {
180 m += 1;
181 if (z0 == 0) {
182 if ((z1 | z2 | z3) == 0) {
183 /* 0: we are done */
184 if (n < 31)
185 m <<= (1 + n);
186 else
187 m = 0;
188 m &= ~0x80000000;
189 *quo = sq >= 0 ? m : -m;
190 __H0(a) = hx & is;
191 __H1(a) = __H2(a) = __H3(a) = 0;
192 return (a);
193 }
194 }
195 /* x = z << 1 */
196 z0 = z0 + z0 + ((z1 & is) != 0);
197 z1 = z1 + z1 + ((z2 & is) != 0);
198 z2 = z2 + z2 + ((z3 & is) != 0);
199 z3 = z3 + z3;
200 x0 = z0;
201 x1 = z1;
202 x2 = z2;
203 x3 = z3;
204 m += m;
205 }
206 }
207 carry = 0;
208 z3 = x3 - y3;
209 carry = z3 > x3;
210 if (carry == 0) {
211 z2 = x2 - y2;
212 carry = z2 > x2;
213 } else {
214 z2 = x2 - y2 - 1;
215 carry = z2 >= x2;
216 }
217 if (carry == 0) {
218 z1 = x1 - y1;
219 carry = z1 > x1;
220 } else {
221 z1 = x1 - y1 - 1;
222 carry = z1 >= x1;
223 }
224 z0 = x0 - y0 - carry;
225 if (z0 >= 0) {
226 x0 = z0;
227 x1 = z1;
228 x2 = z2;
229 x3 = z3;
230 m += 1;
231 }
232 m &= ~0x80000000;
233 *quo = sq >= 0 ? m : -m;
234
235 /* convert back to floating value and restore the sign */
236 if ((x0 | x1 | x2 | x3) == 0) {
237 __H0(a) = hx & is;
238 __H1(a) = __H2(a) = __H3(a) = 0;
239 return (a);
240 }
241 while (x0 < iu) {
242 if (x0 == 0) {
243 iy -= 16;
244 x0 = x1 >> 16;
245 x1 = (x1 << 16) | (x2 >> 16);
246 x2 = (x2 << 16) | (x3 >> 16);
247 x3 = (x3 << 16);
248 } else {
249 x0 = x0 + x0 + ((x1 & is) != 0);
250 x1 = x1 + x1 + ((x2 & is) != 0);
251 x2 = x2 + x2 + ((x3 & is) != 0);
252 x3 = x3 + x3;
253 iy -= 1;
254 }
255 }
256
257 /* normalize output */
258 if (iy >= -16382) {
259 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
260 __H1(a) = x1;
261 __H2(a) = x2;
262 __H3(a) = x3;
263 } else { /* subnormal output */
264 n = -16382 - iy;
265 k = n & 31;
266 if (k <= 16) {
267 x3 = (x2 << (32 - k)) | (x3 >> k);
268 x2 = (x1 << (32 - k)) | (x2 >> k);
269 x1 = (x0 << (32 - k)) | (x1 >> k);
270 x0 >>= k;
271 } else {
272 x3 = (x2 << (32 - k)) | (x3 >> k);
273 x2 = (x1 << (32 - k)) | (x2 >> k);
274 x1 = (x0 << (32 - k)) | (x1 >> k);
275 x0 = 0;
276 }
277 while (n >= 32) {
278 n -= 32;
279 x3 = x2;
280 x2 = x1;
281 x1 = x0;
282 x0 = 0;
283 }
284 __H0(a) = x0 | sx;
285 __H1(a) = x1;
286 __H2(a) = x2;
287 __H3(a) = x3;
288 a *= one;
289 }
290 return (a);
291 }
292
293 long double
294 remquol(long double x, long double y, int *quo) {
295 int hx, hy, sx, sq;
296 long double v;
297
298 hx = __H0(x); /* high word of x */
299 hy = __H0(y); /* high word of y */
300 sx = hx & is; /* sign of x */
301 sq = (hx ^ hy) & is; /* sign of x/y */
302 hx ^= sx; /* |x| */
303 hy &= ~0x80000000;
304
305 /* purge off exception values */
306 *quo = 0;
307 /* y=0, y is NaN, x is NaN or inf */
308 if (y == 0.0L || y != y || hx >= 0x7fff0000)
309 return ((x * y) / (x * y));
310
311 y = fabsl(y);
312 x = fabsl(x);
313 if (hy <= 0x7ffdffff) {
314 x = fmodquol(x, y + y, quo);
315 *quo = ((*quo) & 0x3fffffff) << 1;
316 }
317 if (hy < 0x00020000) {
318 if (x + x > y) {
319 *quo += 1;
320 if (x == y)
321 x = zero;
322 else
323 x -= y;
324 if (x + x >= y) {
325 x -= y;
326 *quo += 1;
327 }
328 }
329 } else {
330 v = 0.5L * y;
331 if (x > v) {
332 *quo += 1;
333 if (x == y)
334 x = zero;
335 else
336 x -= y;
337 if (x >= v) {
338 x -= y;
339 *quo += 1;
340 }
341 }
342 }
343 if (sq != 0)
344 *quo = -(*quo);
345 return (sx == 0 ? x : -x);
346 }
|
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __remquol = remquol
32
33 #include "libm.h"
34 #if defined(__SUNPRO_C)
35 #include <sunmath.h> /* fabsl */
36 #endif
37
38 static const int is = -0x7fffffff - 1, im = 0x0000ffff, iu = 0x00010000;
39 static const long double zero = 0.0L, one = 1.0L;
40
41
42 #if defined(__sparc)
43 #define __H0(x) ((int *)&x)[0]
44 #define __H1(x) ((int *)&x)[1]
45 #define __H2(x) ((int *)&x)[2]
46 #define __H3(x) ((int *)&x)[3]
47 #else
48 #error Unsupported architecture
49 #endif
50
51 /*
52 * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered
53 */
54 static long double
55 fmodquol(long double x, long double y, int *quo)
56 {
57 long double a, b;
58 int n, ix, iy, k, sx, sq, m;
59 int hx;
60 int x0, y0, z0, carry;
61 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
62
63 hx = __H0(x);
64 x1 = __H1(x);
65 x2 = __H2(x);
66 x3 = __H3(x);
67 y0 = __H0(y);
68 y1 = __H1(y);
69 y2 = __H2(y);
70 y3 = __H3(y);
71
72 sx = hx & is;
73 sq = (hx ^ y0) & is;
74 x0 = hx ^ sx;
75 y0 &= ~0x80000000;
76
77 a = fabsl(x);
78 b = fabsl(y);
79
80 if (a <= b) {
81 if (a < b) {
82 return (x);
83 } else {
84 *quo = 1 + (sq >> 30);
85 return (zero * x);
86 }
87 }
88
89 /* determine ix = ilogbl(x) */
90 if (x0 < iu) { /* subnormal x */
91 ix = 0;
92 ix = -16382;
93
94 while (x0 == 0) {
95 ix -= 16;
96 x0 = x1 >> 16;
97 x1 = (x1 << 16) | (x2 >> 16);
98 x2 = (x2 << 16) | (x3 >> 16);
99 x3 = (x3 << 16);
100 }
101
102 while (x0 < iu) {
103 ix -= 1;
104 x0 = (x0 << 1) | (x1 >> 31);
105 x1 = (x1 << 1) | (x2 >> 31);
106 x2 = (x2 << 1) | (x3 >> 31);
107 x3 <<= 1;
108 }
109 } else {
110 ix = (x0 >> 16) - 16383;
111 x0 = iu | (x0 & im);
112 }
113
114 /* determine iy = ilogbl(y) */
115 if (y0 < iu) { /* subnormal y */
116 iy = -16382;
117
118 while (y0 == 0) {
119 iy -= 16;
120 y0 = y1 >> 16;
121 y1 = (y1 << 16) | (y2 >> 16);
122 y2 = (y2 << 16) | (y3 >> 16);
123 y3 = (y3 << 16);
124 }
125
126 while (y0 < iu) {
127 iy -= 1;
128 y0 = (y0 << 1) | (y1 >> 31);
129 y1 = (y1 << 1) | (y2 >> 31);
130 y2 = (y2 << 1) | (y3 >> 31);
131 y3 <<= 1;
132 }
133 } else {
134 iy = (y0 >> 16) - 16383;
135 y0 = iu | (y0 & im);
136 }
137
138 /* fix point fmod */
139 n = ix - iy;
140 m = 0;
141
142 while (n--) {
143 while (x0 == 0 && n >= 16) {
144 m <<= 16;
145 n -= 16;
146 x0 = x1 >> 16;
147 x1 = (x1 << 16) | (x2 >> 16);
148 x2 = (x2 << 16) | (x3 >> 16);
149 x3 = (x3 << 16);
150 }
151
152 while (x0 < iu && n >= 1) {
153 m += m;
154 n -= 1;
155 x0 = (x0 << 1) | (x1 >> 31);
156 x1 = (x1 << 1) | (x2 >> 31);
157 x2 = (x2 << 1) | (x3 >> 31);
158 x3 = (x3 << 1);
159 }
160
161 carry = 0;
162 z3 = x3 - y3;
163 carry = z3 > x3;
164
165 if (carry == 0) {
166 z2 = x2 - y2;
167 carry = z2 > x2;
168 } else {
169 z2 = x2 - y2 - 1;
170 carry = z2 >= x2;
171 }
172
173 if (carry == 0) {
174 z1 = x1 - y1;
175 carry = z1 > x1;
176 } else {
177 z1 = x1 - y1 - 1;
178 carry = z1 >= x1;
179 }
180
181 z0 = x0 - y0 - carry;
182
183 if (z0 < 0) { /* double x */
184 x0 = x0 + x0 + ((x1 & is) != 0);
185 x1 = x1 + x1 + ((x2 & is) != 0);
186 x2 = x2 + x2 + ((x3 & is) != 0);
187 x3 = x3 + x3;
188 m += m;
189 } else {
190 m += 1;
191
192 if (z0 == 0) {
193 if ((z1 | z2 | z3) == 0) {
194 /* 0: we are done */
195 if (n < 31)
196 m <<= (1 + n);
197 else
198 m = 0;
199
200 m &= ~0x80000000;
201 *quo = sq >= 0 ? m : -m;
202 __H0(a) = hx & is;
203 __H1(a) = __H2(a) = __H3(a) = 0;
204 return (a);
205 }
206 }
207
208 /* x = z << 1 */
209 z0 = z0 + z0 + ((z1 & is) != 0);
210 z1 = z1 + z1 + ((z2 & is) != 0);
211 z2 = z2 + z2 + ((z3 & is) != 0);
212 z3 = z3 + z3;
213 x0 = z0;
214 x1 = z1;
215 x2 = z2;
216 x3 = z3;
217 m += m;
218 }
219 }
220
221 carry = 0;
222 z3 = x3 - y3;
223 carry = z3 > x3;
224
225 if (carry == 0) {
226 z2 = x2 - y2;
227 carry = z2 > x2;
228 } else {
229 z2 = x2 - y2 - 1;
230 carry = z2 >= x2;
231 }
232
233 if (carry == 0) {
234 z1 = x1 - y1;
235 carry = z1 > x1;
236 } else {
237 z1 = x1 - y1 - 1;
238 carry = z1 >= x1;
239 }
240
241 z0 = x0 - y0 - carry;
242
243 if (z0 >= 0) {
244 x0 = z0;
245 x1 = z1;
246 x2 = z2;
247 x3 = z3;
248 m += 1;
249 }
250
251 m &= ~0x80000000;
252 *quo = sq >= 0 ? m : -m;
253
254 /* convert back to floating value and restore the sign */
255 if ((x0 | x1 | x2 | x3) == 0) {
256 __H0(a) = hx & is;
257 __H1(a) = __H2(a) = __H3(a) = 0;
258 return (a);
259 }
260
261 while (x0 < iu) {
262 if (x0 == 0) {
263 iy -= 16;
264 x0 = x1 >> 16;
265 x1 = (x1 << 16) | (x2 >> 16);
266 x2 = (x2 << 16) | (x3 >> 16);
267 x3 = (x3 << 16);
268 } else {
269 x0 = x0 + x0 + ((x1 & is) != 0);
270 x1 = x1 + x1 + ((x2 & is) != 0);
271 x2 = x2 + x2 + ((x3 & is) != 0);
272 x3 = x3 + x3;
273 iy -= 1;
274 }
275 }
276
277 /* normalize output */
278 if (iy >= -16382) {
279 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
280 __H1(a) = x1;
281 __H2(a) = x2;
282 __H3(a) = x3;
283 } else { /* subnormal output */
284 n = -16382 - iy;
285 k = n & 31;
286
287 if (k <= 16) {
288 x3 = (x2 << (32 - k)) | (x3 >> k);
289 x2 = (x1 << (32 - k)) | (x2 >> k);
290 x1 = (x0 << (32 - k)) | (x1 >> k);
291 x0 >>= k;
292 } else {
293 x3 = (x2 << (32 - k)) | (x3 >> k);
294 x2 = (x1 << (32 - k)) | (x2 >> k);
295 x1 = (x0 << (32 - k)) | (x1 >> k);
296 x0 = 0;
297 }
298
299 while (n >= 32) {
300 n -= 32;
301 x3 = x2;
302 x2 = x1;
303 x1 = x0;
304 x0 = 0;
305 }
306
307 __H0(a) = x0 | sx;
308 __H1(a) = x1;
309 __H2(a) = x2;
310 __H3(a) = x3;
311 a *= one;
312 }
313
314 return (a);
315 }
316
317 long double
318 remquol(long double x, long double y, int *quo)
319 {
320 int hx, hy, sx, sq;
321 long double v;
322
323 hx = __H0(x); /* high word of x */
324 hy = __H0(y); /* high word of y */
325 sx = hx & is; /* sign of x */
326 sq = (hx ^ hy) & is; /* sign of x/y */
327 hx ^= sx; /* |x| */
328 hy &= ~0x80000000;
329
330 /* purge off exception values */
331 *quo = 0;
332
333 /* y=0, y is NaN, x is NaN or inf */
334 if (y == 0.0L || y != y || hx >= 0x7fff0000)
335 return ((x * y) / (x * y));
336
337 y = fabsl(y);
338 x = fabsl(x);
339
340 if (hy <= 0x7ffdffff) {
341 x = fmodquol(x, y + y, quo);
342 *quo = ((*quo) & 0x3fffffff) << 1;
343 }
344
345 if (hy < 0x00020000) {
346 if (x + x > y) {
347 *quo += 1;
348
349 if (x == y)
350 x = zero;
351 else
352 x -= y;
353
354 if (x + x >= y) {
355 x -= y;
356 *quo += 1;
357 }
358 }
359 } else {
360 v = 0.5L * y;
361
362 if (x > v) {
363 *quo += 1;
364
365 if (x == y)
366 x = zero;
367 else
368 x -= y;
369
370 if (x >= v) {
371 x -= y;
372 *quo += 1;
373 }
374 }
375 }
376
377 if (sq != 0)
378 *quo = -(*quo);
379
380 return (sx == 0 ? x : -x);
381 }
|