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 __remquof = remquof
31
32 /* INDENT OFF */
33 /*
34 * float remquof(float x, float y, int *quo) return remainderf(x,y) and an
35 * integer pointer quo such that *quo = N mod (2**31), where N is the
36 * exact integeral part of x/y rounded to nearest even.
37 *
38 * remquof call internal fmodquof
39 */
40
41 #include "libm.h"
42 #include "libm_protos.h"
43 #include <math.h>
44 extern float fabsf(float);
45
46 static const int
47 is = (int) 0x80000000,
48 im = 0x007fffff,
49 ii = 0x7f800000,
50 iu = 0x00800000;
51
52 static const float zero = 0.0F, half = 0.5F;
53 /* INDENT ON */
54
55 static float
56 fmodquof(float x, float y, int *quo) {
57 float w;
58 int hx, ix, iy, iz, k, ny, nd, m, sq;
59
60 hx = *(int *) &x;
61 ix = hx & 0x7fffffff;
62 iy = *(int *) &y;
63 sq = (iy ^ hx) & is; /* sign of x/y */
64 iy &= 0x7fffffff;
65
66 /* purge off exception values */
67 *quo = 0;
68 if (ix >= ii || iy > ii || iy == 0) {
69 w = x * y;
70 w = w / w;
71 } else if (ix <= iy) {
72 if (ix < iy)
73 w = x; /* return x if |x|<|y| */
74 else {
75 *quo = 1 + (sq >> 30);
76 w = zero * x; /* return sign(x)*0.0 */
77 }
78 } else {
79 /* INDENT OFF */
80 /*
81 * scale x,y to "normal" with
82 * ny = exponent of y
83 * nd = exponent of x minus exponent of y
84 */
85 /* INDENT ON */
86 ny = iy >> 23;
87 k = ix >> 23;
88
89 /* special case for subnormal y or x */
90 if (ny == 0) {
91 ny = 1;
92 while (iy < iu) {
93 ny -= 1;
94 iy += iy;
95 }
96 nd = k - ny;
97 if (k == 0) {
98 nd += 1;
99 while (ix < iu) {
100 nd -= 1;
101 ix += ix;
102 }
103 } else
104 ix = iu | (ix & im);
105 } else {
106 nd = k - ny;
107 ix = iu | (ix & im);
108 iy = iu | (iy & im);
109 }
110 /* INDENT OFF */
111 /* fix point fmod for normalized ix and iy */
112 /*
113 * while (nd--) {
114 * iz = ix - iy;
115 * if (iz < 0)
116 * ix = ix + ix;
117 * else if (iz == 0) {
118 * *(int *) &w = is & hx;
119 * return w;
120 * } else
121 * ix = iz + iz;
122 * }
123 */
124 /* INDENT ON */
125 /* unroll the above loop 4 times to gain performance */
126 m = 0;
127 k = nd >> 2;
128 nd -= (k << 2);
129 while (k--) {
130 iz = ix - iy;
131 if (iz >= 0) {
132 m += 1;
133 ix = iz + iz;
134 } else
135 ix += ix;
136 m += m;
137 iz = ix - iy;
138 if (iz >= 0) {
139 m += 1;
140 ix = iz + iz;
141 } else
142 ix += ix;
143 m += m;
144 iz = ix - iy;
145 if (iz >= 0) {
146 m += 1;
147 ix = iz + iz;
148 } else
149 ix += ix;
150 m += m;
151 iz = ix - iy;
152 if (iz >= 0) {
153 m += 1;
154 ix = iz + iz;
155 } else
156 ix += ix;
157 m += m;
158 if (iz == 0) {
159 iz = (k << 2) + nd;
160 if (iz < 32)
161 m <<= iz;
162 else
163 m = 0;
164 m &= 0x7fffffff;
165 *quo = sq >= 0 ? m : -m;
166 *(int *) &w = is & hx;
167 return (w);
168 }
169 }
170 while (nd--) {
171 iz = ix - iy;
172 if (iz >= 0) {
173 m += 1;
174 ix = iz + iz;
175 } else
176 ix += ix;
177 m += m;
178 }
179 /* end of unrolling */
180
181 iz = ix - iy;
182 if (iz >= 0) {
183 m += 1;
184 ix = iz;
185 }
186 m &= 0x7fffffff;
187 *quo = sq >= 0 ? m : -m;
188
189 /* convert back to floating value and restore the sign */
190 if (ix == 0) {
191 *(int *) &w = is & hx;
192 return (w);
193 }
194 while (ix < iu) {
195 ix += ix;
196 ny -= 1;
197 }
198 while (ix > (iu + iu)) {
199 ny += 1;
200 ix >>= 1;
201 }
202 if (ny > 0)
203 *(int *) &w = (is & hx) | (ix & im) | (ny << 23);
204 else { /* subnormal output */
205 k = -ny + 1;
206 ix >>= k;
207 *(int *) &w = (is & hx) | ix;
208 }
209 }
210 return (w);
211 }
212
213 float
214 remquof(float x, float y, int *quo) {
215 int hx, hy, sx, sq;
216 float v;
217
218 hx = *(int *) &x; /* high word of x */
219 hy = *(int *) &y; /* high word of y */
220 sx = hx & is; /* sign of x */
221 sq = (hx ^ hy) & is; /* sign of x/y */
222 hx ^= sx; /* |x| */
223 hy &= 0x7fffffff; /* |y| */
224
225 /* purge off exception values: y is 0 or NaN, x is Inf or NaN */
226 *quo = 0;
227 if (hx >= ii || hy > ii || hy == 0) {
228 v = x * y;
229 return (v / v);
230 }
231
232 y = fabsf(y);
233 x = fabsf(x);
234 if (hy <= 0x7f7fffff) {
235 x = fmodquof(x, y + y, quo);
236 *quo = ((*quo) & 0x3fffffff) << 1;
237 }
238 if (hy < 0x01000000) {
239 if (x + x > y) {
240 *quo += 1;
241 if (x == y)
242 x = zero;
243 else
244 x -= y;
245 if (x + x >= y) {
246 x -= y;
247 *quo += 1;
248 }
249 }
250 } else {
251 v = half * y;
252 if (x > v) {
253 *quo += 1;
254 if (x == y)
255 x = zero;
256 else
257 x -= y;
258 if (x >= v) {
259 x -= y;
260 *quo += 1;
261 }
262 }
263 }
264 if (sq != 0)
265 *quo = -(*quo);
266 return (sx == 0 ? x : -x);
267 }
|
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 __remquof = remquof
32
33
34 /*
35 * float remquof(float x, float y, int *quo) return remainderf(x,y) and an
36 * integer pointer quo such that *quo = N mod (2**31), where N is the
37 * exact integeral part of x/y rounded to nearest even.
38 *
39 * remquof call internal fmodquof
40 */
41
42 #include "libm.h"
43 #include "libm_protos.h"
44 #include <math.h>
45
46 extern float fabsf(float);
47 static const int is = (int)0x80000000,
48 im = 0x007fffff,
49 ii = 0x7f800000,
50 iu = 0x00800000;
51 static const float zero = 0.0F, half = 0.5F;
52
53 static float
54 fmodquof(float x, float y, int *quo)
55 {
56 float w;
57 int hx, ix, iy, iz, k, ny, nd, m, sq;
58
59 hx = *(int *)&x;
60 ix = hx & 0x7fffffff;
61 iy = *(int *)&y;
62 sq = (iy ^ hx) & is; /* sign of x/y */
63 iy &= 0x7fffffff;
64
65 /* purge off exception values */
66 *quo = 0;
67
68 if (ix >= ii || iy > ii || iy == 0) {
69 w = x * y;
70 w = w / w;
71 } else if (ix <= iy) {
72 if (ix < iy) {
73 w = x; /* return x if |x|<|y| */
74 } else {
75 *quo = 1 + (sq >> 30);
76 w = zero * x; /* return sign(x)*0.0 */
77 }
78 } else {
79
80 /*
81 * scale x,y to "normal" with
82 * ny = exponent of y
83 * nd = exponent of x minus exponent of y
84 */
85 ny = iy >> 23;
86 k = ix >> 23;
87
88 /* special case for subnormal y or x */
89 if (ny == 0) {
90 ny = 1;
91
92 while (iy < iu) {
93 ny -= 1;
94 iy += iy;
95 }
96
97 nd = k - ny;
98
99 if (k == 0) {
100 nd += 1;
101
102 while (ix < iu) {
103 nd -= 1;
104 ix += ix;
105 }
106 } else {
107 ix = iu | (ix & im);
108 }
109 } else {
110 nd = k - ny;
111 ix = iu | (ix & im);
112 iy = iu | (iy & im);
113 }
114
115 /*
116 * fix point fmod for normalized ix and iy
117 */
118
119 /*
120 * while (nd--) {
121 * iz = ix - iy;
122 * if (iz < 0)
123 * ix = ix + ix;
124 * else if (iz == 0) {
125 * *(int *) &w = is & hx;
126 * return w;
127 * } else
128 * ix = iz + iz;
129 * }
130 */
131
132 /*
133 * unroll the above loop 4 times to gain performance
134 */
135 m = 0;
136 k = nd >> 2;
137 nd -= (k << 2);
138
139 while (k--) {
140 iz = ix - iy;
141
142 if (iz >= 0) {
143 m += 1;
144 ix = iz + iz;
145 } else {
146 ix += ix;
147 }
148
149 m += m;
150 iz = ix - iy;
151
152 if (iz >= 0) {
153 m += 1;
154 ix = iz + iz;
155 } else {
156 ix += ix;
157 }
158
159 m += m;
160 iz = ix - iy;
161
162 if (iz >= 0) {
163 m += 1;
164 ix = iz + iz;
165 } else {
166 ix += ix;
167 }
168
169 m += m;
170 iz = ix - iy;
171
172 if (iz >= 0) {
173 m += 1;
174 ix = iz + iz;
175 } else {
176 ix += ix;
177 }
178
179 m += m;
180
181 if (iz == 0) {
182 iz = (k << 2) + nd;
183
184 if (iz < 32)
185 m <<= iz;
186 else
187 m = 0;
188
189 m &= 0x7fffffff;
190 *quo = sq >= 0 ? m : -m;
191 *(int *)&w = is & hx;
192 return (w);
193 }
194 }
195
196 while (nd--) {
197 iz = ix - iy;
198
199 if (iz >= 0) {
200 m += 1;
201 ix = iz + iz;
202 } else {
203 ix += ix;
204 }
205
206 m += m;
207 }
208
209 /* end of unrolling */
210
211 iz = ix - iy;
212
213 if (iz >= 0) {
214 m += 1;
215 ix = iz;
216 }
217
218 m &= 0x7fffffff;
219 *quo = sq >= 0 ? m : -m;
220
221 /* convert back to floating value and restore the sign */
222 if (ix == 0) {
223 *(int *)&w = is & hx;
224 return (w);
225 }
226
227 while (ix < iu) {
228 ix += ix;
229 ny -= 1;
230 }
231
232 while (ix > (iu + iu)) {
233 ny += 1;
234 ix >>= 1;
235 }
236
237 if (ny > 0) {
238 *(int *)&w = (is & hx) | (ix & im) | (ny << 23);
239 } else { /* subnormal output */
240 k = -ny + 1;
241 ix >>= k;
242 *(int *)&w = (is & hx) | ix;
243 }
244 }
245
246 return (w);
247 }
248
249 float
250 remquof(float x, float y, int *quo)
251 {
252 int hx, hy, sx, sq;
253 float v;
254
255 hx = *(int *)&x; /* high word of x */
256 hy = *(int *)&y; /* high word of y */
257 sx = hx & is; /* sign of x */
258 sq = (hx ^ hy) & is; /* sign of x/y */
259 hx ^= sx; /* |x| */
260 hy &= 0x7fffffff; /* |y| */
261
262 /* purge off exception values: y is 0 or NaN, x is Inf or NaN */
263 *quo = 0;
264
265 if (hx >= ii || hy > ii || hy == 0) {
266 v = x * y;
267 return (v / v);
268 }
269
270 y = fabsf(y);
271 x = fabsf(x);
272
273 if (hy <= 0x7f7fffff) {
274 x = fmodquof(x, y + y, quo);
275 *quo = ((*quo) & 0x3fffffff) << 1;
276 }
277
278 if (hy < 0x01000000) {
279 if (x + x > y) {
280 *quo += 1;
281
282 if (x == y)
283 x = zero;
284 else
285 x -= y;
286
287 if (x + x >= y) {
288 x -= y;
289 *quo += 1;
290 }
291 }
292 } else {
293 v = half * y;
294
295 if (x > v) {
296 *quo += 1;
297
298 if (x == y)
299 x = zero;
300 else
301 x -= y;
302
303 if (x >= v) {
304 x -= y;
305 *quo += 1;
306 }
307 }
308 }
309
310 if (sq != 0)
311 *quo = -(*quo);
312
313 return (sx == 0 ? x : -x);
314 }
|