1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
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 }