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