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_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
fmodquof(float x,float y,int * quo)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
remquof(float x,float y,int * quo)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 }
268