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 remquo = __remquo
31
32 /* INDENT OFF */
33 /*
34 * double remquo(double x, double y, int *quo) return remainder(x,y) and an
35 * integer pointer quo such that *quo = N mod {2**31}, where N is the
36 * exact integral part of x/y rounded to nearest even.
37 *
38 * remquo call internal fmodquo
39 */
40 /* INDENT ON */
41
42 #include "libm.h"
43 #include "libm_synonyms.h"
44 #include "libm_protos.h"
45 #include <math.h> /* fabs() */
46 #include <sys/isa_defs.h>
47
48 #if defined(_BIG_ENDIAN)
49 #define HIWORD 0
50 #define LOWORD 1
51 #else
52 #define HIWORD 1
53 #define LOWORD 0
54 #endif
55 #define __HI(x) ((int *) &x)[HIWORD]
56 #define __LO(x) ((int *) &x)[LOWORD]
57
58 static const double one = 1.0, Zero[] = {0.0, -0.0};
59
60 static double
61 fmodquo(double x, double y, int *quo) {
62 int n, hx, hy, hz, ix, iy, sx, sq, i, m;
63 unsigned lx, ly, lz;
64
65 hx = __HI(x); /* high word of x */
66 lx = __LO(x); /* low word of x */
67 hy = __HI(y); /* high word of y */
68 ly = __LO(y); /* low word of y */
69 sx = hx & 0x80000000; /* sign of x */
70 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */
71 hx ^= sx; /* |x| */
72 hy &= 0x7fffffff; /* |y| */
73
74 /* purge off exception values */
75 *quo = 0;
76 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */
77 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */
78 return ((x * y) / (x * y));
79 if (hx <= hy) {
80 if (hx < hy || lx < ly)
81 return (x); /* |x|<|y| return x */
82 if (lx == ly) {
83 *quo = 1 + (sq >> 30);
84 /* |x|=|y| return x*0 */
85 return (Zero[(unsigned) sx >> 31]);
86 }
87 }
88
89 /* determine ix = ilogb(x) */
90 if (hx < 0x00100000) { /* subnormal x */
91 if (hx == 0) {
92 for (ix = -1043, i = lx; i > 0; i <<= 1)
93 ix -= 1;
94 } else {
95 for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
96 ix -= 1;
97 }
98 } else
99 ix = (hx >> 20) - 1023;
100
101 /* determine iy = ilogb(y) */
102 if (hy < 0x00100000) { /* subnormal y */
103 if (hy == 0) {
104 for (iy = -1043, i = ly; i > 0; i <<= 1)
105 iy -= 1;
106 } else {
107 for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
108 iy -= 1;
109 }
110 } else
111 iy = (hy >> 20) - 1023;
112
113 /* set up {hx,lx}, {hy,ly} and align y to x */
114 if (ix >= -1022)
115 hx = 0x00100000 | (0x000fffff & hx);
116 else { /* subnormal x, shift x to normal */
117 n = -1022 - ix;
118 if (n <= 31) {
119 hx = (hx << n) | (lx >> (32 - n));
120 lx <<= n;
121 } else {
122 hx = lx << (n - 32);
123 lx = 0;
124 }
125 }
126 if (iy >= -1022)
127 hy = 0x00100000 | (0x000fffff & hy);
128 else { /* subnormal y, shift y to normal */
129 n = -1022 - iy;
130 if (n <= 31) {
131 hy = (hy << n) | (ly >> (32 - n));
132 ly <<= n;
133 } else {
134 hy = ly << (n - 32);
135 ly = 0;
136 }
137 }
138
139 /* fix point fmod */
140 n = ix - iy;
141 m = 0;
142 while (n--) {
143 hz = hx - hy;
144 lz = lx - ly;
145 if (lx < ly)
146 hz -= 1;
147 if (hz < 0) {
148 hx = hx + hx + (lx >> 31);
149 lx = lx + lx;
150 } else {
151 m += 1;
152 if ((hz | lz) == 0) { /* return sign(x)*0 */
153 if (n < 31)
154 m <<= 1 + n;
155 else
156 m = 0;
157 m &= 0x7fffffff;
158 *quo = sq >= 0 ? m : -m;
159 return (Zero[(unsigned) sx >> 31]);
160 }
161 hx = hz + hz + (lz >> 31);
162 lx = lz + lz;
163 }
164 m += m;
165 }
166 hz = hx - hy;
167 lz = lx - ly;
168 if (lx < ly)
169 hz -= 1;
170 if (hz >= 0) {
171 hx = hz;
172 lx = lz;
173 m += 1;
174 }
175 m &= 0x7fffffff;
176 *quo = sq >= 0 ? m : -m;
177
178 /* convert back to floating value and restore the sign */
179 if ((hx | lx) == 0) { /* return sign(x)*0 */
180 return (Zero[(unsigned) sx >> 31]);
181 }
182 while (hx < 0x00100000) { /* normalize x */
183 hx = hx + hx + (lx >> 31);
184 lx = lx + lx;
185 iy -= 1;
186 }
187 if (iy >= -1022) { /* normalize output */
188 hx = (hx - 0x00100000) | ((iy + 1023) << 20);
189 __HI(x) = hx | sx;
190 __LO(x) = lx;
191 } else { /* subnormal output */
192 n = -1022 - iy;
193 if (n <= 20) {
194 lx = (lx >> n) | ((unsigned) hx << (32 - n));
195 hx >>= n;
196 } else if (n <= 31) {
197 lx = (hx << (32 - n)) | (lx >> n);
198 hx = sx;
199 } else {
200 lx = hx >> (n - 32);
201 hx = sx;
202 }
203 __HI(x) = hx | sx;
204 __LO(x) = lx;
205 x *= one; /* create necessary signal */
206 }
207 return (x); /* exact output */
208 }
209
210 #define zero Zero[0]
211
212 double
213 remquo(double x, double y, int *quo) {
214 int hx, hy, sx, sq;
215 double v;
216 unsigned ly;
217
218 hx = __HI(x); /* high word of x */
219 hy = __HI(y); /* high word of y */
220 ly = __LO(y); /* low word of y */
221 sx = hx & 0x80000000; /* sign of x */
222 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */
223 hx ^= sx; /* |x| */
224 hy &= 0x7fffffff; /* |y| */
225
226 /* purge off exception values */
227 *quo = 0;
228 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */
229 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */
230 return ((x * y) / (x * y));
231
232 y = fabs(y);
233 x = fabs(x);
234 if (hy <= 0x7fdfffff) {
235 x = fmodquo(x, y + y, quo);
236 *quo = ((*quo) & 0x3fffffff) << 1;
237 }
238 if (hy < 0x00200000) {
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 = 0.5 * 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 }