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 __fmodl = fmodl
31
32 #include "libm.h"
33
34 static const int
35 is = -0x7fffffff - 1,
36 im = 0x0000ffff,
37 iu = 0x00010000;
38
39 static const long double
40 zero = 0.0L,
41 one = 1.0L;
42
43 #ifdef __LITTLE_ENDIAN
44 #define __H0(x) *(3 + (int *) &x)
45 #define __H1(x) *(2 + (int *) &x)
46 #define __H2(x) *(1 + (int *) &x)
47 #define __H3(x) *(0 + (int *) &x)
48 #else
49 #define __H0(x) *(0 + (int *) &x)
50 #define __H1(x) *(1 + (int *) &x)
51 #define __H2(x) *(2 + (int *) &x)
52 #define __H3(x) *(3 + (int *) &x)
53 #endif
54
55 long double
56 fmodl(long double x, long double y) {
57 long double a, b;
58 int n, ix, iy, k, sx;
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 & 0x80000000;
73 x0 = hx ^ sx;
74 y0 &= 0x7fffffff;
75
76 /* purge off exception values */
77 if (x0 >= 0x7fff0000 || /* !finitel(x) */
78 (y0 > 0x7fff0000) || (y0 == 0x7fff0000 && ((y1 | y2 | y3) != 0)) ||
79 (y0 | y1 | y2 | y3) == 0) /* isnanl(y) || y = 0 */
80 return ((x * y) / (x * y));
81 a = fabsl(x);
82 b = fabsl(y);
83 if (a <= b) {
84 if (a < b)
85 return (x);
86 else
87 return (zero * x);
88 }
89 /* determine ix = ilogbl(x) */
90 if (x0 < iu) { /* subnormal x */
91 ix = -16382;
92 while (x0 == 0) {
93 ix -= 16;
94 x0 = x1 >> 16;
95 x1 = (x1 << 16) | (x2 >> 16);
96 x2 = (x2 << 16) | (x3 >> 16);
97 x3 = (x3 << 16);
98 }
99 while (x0 < iu) {
100 ix -= 1;
101 x0 = (x0 << 1) | (x1 >> 31);
102 x1 = (x1 << 1) | (x2 >> 31);
103 x2 = (x2 << 1) | (x3 >> 31);
104 x3 <<= 1;
105 }
106 } else {
107 ix = (x0 >> 16) - 16383;
108 x0 = iu | (x0 & im);
109 }
110
111 /* determine iy = ilogbl(y) */
112 if (y0 < iu) { /* subnormal y */
113 iy = -16382;
114 while (y0 == 0) {
115 iy -= 16;
116 y0 = y1 >> 16;
117 y1 = (y1 << 16) | (y2 >> 16);
118 y2 = (y2 << 16) | (y3 >> 16);
119 y3 = (y3 << 16);
120 }
121 while (y0 < iu) {
122 iy -= 1;
123 y0 = (y0 << 1) | (y1 >> 31);
124 y1 = (y1 << 1) | (y2 >> 31);
125 y2 = (y2 << 1) | (y3 >> 31);
126 y3 <<= 1;
127 }
128 } else {
129 iy = (y0 >> 16) - 16383;
130 y0 = iu | (y0 & im);
131 }
132
133 /* fix point fmod */
134 n = ix - iy;
135 while (n--) {
136 while (x0 == 0 && n >= 16) {
137 n -= 16;
138 x0 = x1 >> 16;
139 x1 = (x1 << 16) | (x2 >> 16);
140 x2 = (x2 << 16) | (x3 >> 16);
141 x3 = (x3 << 16);
142 }
143 while (x0 < iu && n >= 1) {
144 n -= 1;
145 x0 = (x0 << 1) | (x1 >> 31);
146 x1 = (x1 << 1) | (x2 >> 31);
147 x2 = (x2 << 1) | (x3 >> 31);
148 x3 = (x3 << 1);
149 }
150 carry = 0;
151 z3 = x3 - y3;
152 carry = (z3 > x3);
153 if (carry == 0) {
154 z2 = x2 - y2;
155 carry = (z2 > x2);
156 } else {
157 z2 = x2 - y2 - 1;
158 carry = (z2 >= x2);
159 }
160 if (carry == 0) {
161 z1 = x1 - y1;
162 carry = (z1 > x1);
163 } else {
164 z1 = x1 - y1 - 1;
165 carry = (z1 >= x1);
166 }
167 z0 = x0 - y0 - carry;
168 if (z0 < 0) { /* double x */
169 x0 = x0 + x0 + ((x1 & is) != 0);
170 x1 = x1 + x1 + ((x2 & is) != 0);
171 x2 = x2 + x2 + ((x3 & is) != 0);
172 x3 = x3 + x3;
173 } else {
174 if (z0 == 0) {
175 if ((z1 | z2 | z3) == 0) { /* 0: done */
176 __H0(a) = hx & is;
177 __H1(a) = __H2(a) = __H3(a) = 0;
178 return (a);
179 }
180 }
181 /* x = z << 1 */
182 z0 = z0 + z0 + ((z1 & is) != 0);
183 z1 = z1 + z1 + ((z2 & is) != 0);
184 z2 = z2 + z2 + ((z3 & is) != 0);
185 z3 = z3 + z3;
186 x0 = z0;
187 x1 = z1;
188 x2 = z2;
189 x3 = z3;
190 }
191 }
192
193 carry = 0;
194 z3 = x3 - y3;
195 carry = (z3 > x3);
196 if (carry == 0) {
197 z2 = x2 - y2;
198 carry = (z2 > x2);
199 } else {
200 z2 = x2 - y2 - 1;
201 carry = (z2 >= x2);
202 }
203 if (carry == 0) {
204 z1 = x1 - y1;
205 carry = (z1 > x1);
206 } else {
207 z1 = x1 - y1 - 1;
208 carry = (z1 >= x1);
209 }
210 z0 = x0 - y0 - carry;
211 if (z0 >= 0) {
212 x0 = z0;
213 x1 = z1;
214 x2 = z2;
215 x3 = z3;
216 }
217 /* convert back to floating value and restore the sign */
218 if ((x0 | x1 | x2 | x3) == 0) {
219 __H0(a) = hx & is;
220 __H1(a) = __H2(a) = __H3(a) = 0;
221 return (a);
222 }
223 while (x0 < iu) {
224 if (x0 == 0) {
225 iy -= 16;
226 x0 = x1 >> 16;
227 x1 = (x1 << 16) | (x2 >> 16);
228 x2 = (x2 << 16) | (x3 >> 16);
229 x3 = (x3 << 16);
230 } else {
231 x0 = x0 + x0 + ((x1 & is) != 0);
232 x1 = x1 + x1 + ((x2 & is) != 0);
233 x2 = x2 + x2 + ((x3 & is) != 0);
234 x3 = x3 + x3;
235 iy -= 1;
236 }
237 }
238
239 /* normalize output */
240 if (iy >= -16382) {
241 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
242 __H1(a) = x1;
243 __H2(a) = x2;
244 __H3(a) = x3;
245 } else { /* subnormal output */
246 n = -16382 - iy;
247 k = n & 31;
248 if (k != 0) {
249 if (k <= 16) {
250 x3 = (x2 << (32 - k)) | (x3 >> k);
251 x2 = (x1 << (32 - k)) | (x2 >> k);
252 x1 = (x0 << (32 - k)) | (x1 >> k);
253 x0 >>= k;
254 } else {
255 x3 = (x2 << (32 - k)) | (x3 >> k);
256 x2 = (x1 << (32 - k)) | (x2 >> k);
257 x1 = (x0 << (32 - k)) | (x1 >> k);
258 x0 = 0;
259 }
260 }
261 while (n >= 32) {
262 n -= 32;
263 x3 = x2;
264 x2 = x1;
265 x1 = x0;
266 x0 = 0;
267 }
268 __H0(a) = x0 | sx;
269 __H1(a) = x1;
270 __H2(a) = x2;
271 __H3(a) = x3;
272 a *= one;
273 }
274 return (a);
275 }
|
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 __fmodl = fmodl
32
33 #include "libm.h"
34
35 static const int is = -0x7fffffff - 1, im = 0x0000ffff, iu = 0x00010000;
36 static const long double zero = 0.0L, one = 1.0L;
37
38 #ifdef __LITTLE_ENDIAN
39 #define __H0(x) *(3 + (int *)&x)
40 #define __H1(x) *(2 + (int *)&x)
41 #define __H2(x) *(1 + (int *)&x)
42 #define __H3(x) *(0 + (int *)&x)
43 #else
44 #define __H0(x) *(0 + (int *)&x)
45 #define __H1(x) *(1 + (int *)&x)
46 #define __H2(x) *(2 + (int *)&x)
47 #define __H3(x) *(3 + (int *)&x)
48 #endif
49
50 long double
51 fmodl(long double x, long double y)
52 {
53 long double a, b;
54 int n, ix, iy, k, sx;
55 int hx;
56 int x0, y0, z0, carry;
57 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
58
59 hx = __H0(x);
60 x1 = __H1(x);
61 x2 = __H2(x);
62 x3 = __H3(x);
63 y0 = __H0(y);
64 y1 = __H1(y);
65 y2 = __H2(y);
66 y3 = __H3(y);
67
68 sx = hx & 0x80000000;
69 x0 = hx ^ sx;
70 y0 &= 0x7fffffff;
71
72 /* purge off exception values */
73 if (x0 >= 0x7fff0000 || /* !finitel(x) */
74 (y0 > 0x7fff0000) || (y0 == 0x7fff0000 && ((y1 | y2 | y3) != 0)) ||
75 (y0 | y1 | y2 | y3) == 0) /* isnanl(y) || y = 0 */
76 return ((x * y) / (x * y));
77
78 a = fabsl(x);
79 b = fabsl(y);
80
81 if (a <= b) {
82 if (a < b)
83 return (x);
84 else
85 return (zero * x);
86 }
87
88 /* determine ix = ilogbl(x) */
89 if (x0 < iu) { /* subnormal x */
90 ix = -16382;
91
92 while (x0 == 0) {
93 ix -= 16;
94 x0 = x1 >> 16;
95 x1 = (x1 << 16) | (x2 >> 16);
96 x2 = (x2 << 16) | (x3 >> 16);
97 x3 = (x3 << 16);
98 }
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
116 while (y0 == 0) {
117 iy -= 16;
118 y0 = y1 >> 16;
119 y1 = (y1 << 16) | (y2 >> 16);
120 y2 = (y2 << 16) | (y3 >> 16);
121 y3 = (y3 << 16);
122 }
123
124 while (y0 < iu) {
125 iy -= 1;
126 y0 = (y0 << 1) | (y1 >> 31);
127 y1 = (y1 << 1) | (y2 >> 31);
128 y2 = (y2 << 1) | (y3 >> 31);
129 y3 <<= 1;
130 }
131 } else {
132 iy = (y0 >> 16) - 16383;
133 y0 = iu | (y0 & im);
134 }
135
136 /* fix point fmod */
137 n = ix - iy;
138
139 while (n--) {
140 while (x0 == 0 && n >= 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
148 while (x0 < iu && n >= 1) {
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
156 carry = 0;
157 z3 = x3 - y3;
158 carry = (z3 > x3);
159
160 if (carry == 0) {
161 z2 = x2 - y2;
162 carry = (z2 > x2);
163 } else {
164 z2 = x2 - y2 - 1;
165 carry = (z2 >= x2);
166 }
167
168 if (carry == 0) {
169 z1 = x1 - y1;
170 carry = (z1 > x1);
171 } else {
172 z1 = x1 - y1 - 1;
173 carry = (z1 >= x1);
174 }
175
176 z0 = x0 - y0 - carry;
177
178 if (z0 < 0) { /* double x */
179 x0 = x0 + x0 + ((x1 & is) != 0);
180 x1 = x1 + x1 + ((x2 & is) != 0);
181 x2 = x2 + x2 + ((x3 & is) != 0);
182 x3 = x3 + x3;
183 } else {
184 if (z0 == 0) {
185 if ((z1 | z2 | z3) == 0) { /* 0: done */
186 __H0(a) = hx & is;
187 __H1(a) = __H2(a) = __H3(a) = 0;
188 return (a);
189 }
190 }
191
192 /* x = z << 1 */
193 z0 = z0 + z0 + ((z1 & is) != 0);
194 z1 = z1 + z1 + ((z2 & is) != 0);
195 z2 = z2 + z2 + ((z3 & is) != 0);
196 z3 = z3 + z3;
197 x0 = z0;
198 x1 = z1;
199 x2 = z2;
200 x3 = z3;
201 }
202 }
203
204 carry = 0;
205 z3 = x3 - y3;
206 carry = (z3 > x3);
207
208 if (carry == 0) {
209 z2 = x2 - y2;
210 carry = (z2 > x2);
211 } else {
212 z2 = x2 - y2 - 1;
213 carry = (z2 >= x2);
214 }
215
216 if (carry == 0) {
217 z1 = x1 - y1;
218 carry = (z1 > x1);
219 } else {
220 z1 = x1 - y1 - 1;
221 carry = (z1 >= x1);
222 }
223
224 z0 = x0 - y0 - carry;
225
226 if (z0 >= 0) {
227 x0 = z0;
228 x1 = z1;
229 x2 = z2;
230 x3 = z3;
231 }
232
233 /* convert back to floating value and restore the sign */
234 if ((x0 | x1 | x2 | x3) == 0) {
235 __H0(a) = hx & is;
236 __H1(a) = __H2(a) = __H3(a) = 0;
237 return (a);
238 }
239
240 while (x0 < iu) {
241 if (x0 == 0) {
242 iy -= 16;
243 x0 = x1 >> 16;
244 x1 = (x1 << 16) | (x2 >> 16);
245 x2 = (x2 << 16) | (x3 >> 16);
246 x3 = (x3 << 16);
247 } else {
248 x0 = x0 + x0 + ((x1 & is) != 0);
249 x1 = x1 + x1 + ((x2 & is) != 0);
250 x2 = x2 + x2 + ((x3 & is) != 0);
251 x3 = x3 + x3;
252 iy -= 1;
253 }
254 }
255
256 /* normalize output */
257 if (iy >= -16382) {
258 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
259 __H1(a) = x1;
260 __H2(a) = x2;
261 __H3(a) = x3;
262 } else { /* subnormal output */
263 n = -16382 - iy;
264 k = n & 31;
265
266 if (k != 0) {
267 if (k <= 16) {
268 x3 = (x2 << (32 - k)) | (x3 >> k);
269 x2 = (x1 << (32 - k)) | (x2 >> k);
270 x1 = (x0 << (32 - k)) | (x1 >> k);
271 x0 >>= k;
272 } else {
273 x3 = (x2 << (32 - k)) | (x3 >> k);
274 x2 = (x1 << (32 - k)) | (x2 >> k);
275 x1 = (x0 << (32 - k)) | (x1 >> k);
276 x0 = 0;
277 }
278 }
279
280 while (n >= 32) {
281 n -= 32;
282 x3 = x2;
283 x2 = x1;
284 x1 = x0;
285 x0 = 0;
286 }
287
288 __H0(a) = x0 | sx;
289 __H1(a) = x1;
290 __H2(a) = x2;
291 __H3(a) = x3;
292 a *= one;
293 }
294
295 return (a);
296 }
|