Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/m9x/remquo.c
+++ new/usr/src/lib/libm/common/m9x/remquo.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
|
↓ open down ↓ |
14 lines elided |
↑ open up ↑ |
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
20 20 */
21 21
22 22 /*
23 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 +
25 26 /*
26 27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 28 * Use is subject to license terms.
28 29 */
29 30
30 31 #pragma weak __remquo = remquo
31 32
32 -/* INDENT OFF */
33 +
33 34 /*
34 35 * double remquo(double x, double y, int *quo) return remainder(x,y) and an
35 36 * integer pointer quo such that *quo = N mod {2**31}, where N is the
36 37 * exact integral part of x/y rounded to nearest even.
37 38 *
38 39 * remquo call internal fmodquo
39 40 */
40 -/* INDENT ON */
41 41
42 42 #include "libm.h"
43 43 #include "libm_protos.h"
44 -#include <math.h> /* fabs() */
44 +#include <math.h> /* fabs() */
45 45 #include <sys/isa_defs.h>
46 46
47 47 #if defined(_BIG_ENDIAN)
48 -#define HIWORD 0
49 -#define LOWORD 1
48 +#define HIWORD 0
49 +#define LOWORD 1
50 50 #else
51 -#define HIWORD 1
52 -#define LOWORD 0
51 +#define HIWORD 1
52 +#define LOWORD 0
53 53 #endif
54 -#define __HI(x) ((int *) &x)[HIWORD]
55 -#define __LO(x) ((int *) &x)[LOWORD]
56 -
57 -static const double one = 1.0, Zero[] = {0.0, -0.0};
54 +#define __HI(x) ((int *)&x)[HIWORD]
55 +#define __LO(x) ((int *)&x)[LOWORD]
58 56
57 +static const double one = 1.0, Zero[] = { 0.0, -0.0 };
59 58 static double
60 -fmodquo(double x, double y, int *quo) {
59 +fmodquo(double x, double y, int *quo)
60 +{
61 61 int n, hx, hy, hz, ix, iy, sx, sq, i, m;
62 62 unsigned lx, ly, lz;
63 63
64 - hx = __HI(x); /* high word of x */
65 - lx = __LO(x); /* low word of x */
66 - hy = __HI(y); /* high word of y */
67 - ly = __LO(y); /* low word of y */
68 - sx = hx & 0x80000000; /* sign of x */
64 + hx = __HI(x); /* high word of x */
65 + lx = __LO(x); /* low word of x */
66 + hy = __HI(y); /* high word of y */
67 + ly = __LO(y); /* low word of y */
68 + sx = hx & 0x80000000; /* sign of x */
69 69 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */
70 - hx ^= sx; /* |x| */
71 - hy &= 0x7fffffff; /* |y| */
70 + hx ^= sx; /* |x| */
71 + hy &= 0x7fffffff; /* |y| */
72 72
73 73 /* purge off exception values */
74 74 *quo = 0;
75 +
75 76 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */
76 77 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */
77 78 return ((x * y) / (x * y));
79 +
78 80 if (hx <= hy) {
79 81 if (hx < hy || lx < ly)
80 82 return (x); /* |x|<|y| return x */
83 +
81 84 if (lx == ly) {
82 85 *quo = 1 + (sq >> 30);
83 86 /* |x|=|y| return x*0 */
84 - return (Zero[(unsigned) sx >> 31]);
87 + return (Zero[(unsigned)sx >> 31]);
85 88 }
86 89 }
87 90
88 91 /* determine ix = ilogb(x) */
89 - if (hx < 0x00100000) { /* subnormal x */
92 + if (hx < 0x00100000) { /* subnormal x */
90 93 if (hx == 0) {
91 94 for (ix = -1043, i = lx; i > 0; i <<= 1)
92 95 ix -= 1;
93 96 } else {
94 97 for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
95 98 ix -= 1;
96 99 }
97 - } else
100 + } else {
98 101 ix = (hx >> 20) - 1023;
102 + }
99 103
100 104 /* determine iy = ilogb(y) */
101 - if (hy < 0x00100000) { /* subnormal y */
105 + if (hy < 0x00100000) { /* subnormal y */
102 106 if (hy == 0) {
103 107 for (iy = -1043, i = ly; i > 0; i <<= 1)
104 108 iy -= 1;
105 109 } else {
106 110 for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
107 111 iy -= 1;
108 112 }
109 - } else
113 + } else {
110 114 iy = (hy >> 20) - 1023;
115 + }
111 116
112 117 /* set up {hx,lx}, {hy,ly} and align y to x */
113 - if (ix >= -1022)
118 + if (ix >= -1022) {
114 119 hx = 0x00100000 | (0x000fffff & hx);
115 - else { /* subnormal x, shift x to normal */
120 + } else { /* subnormal x, shift x to normal */
116 121 n = -1022 - ix;
122 +
117 123 if (n <= 31) {
118 124 hx = (hx << n) | (lx >> (32 - n));
119 125 lx <<= n;
120 126 } else {
121 127 hx = lx << (n - 32);
122 128 lx = 0;
123 129 }
124 130 }
125 - if (iy >= -1022)
131 +
132 + if (iy >= -1022) {
126 133 hy = 0x00100000 | (0x000fffff & hy);
127 - else { /* subnormal y, shift y to normal */
134 + } else { /* subnormal y, shift y to normal */
128 135 n = -1022 - iy;
136 +
129 137 if (n <= 31) {
130 138 hy = (hy << n) | (ly >> (32 - n));
131 139 ly <<= n;
132 140 } else {
133 141 hy = ly << (n - 32);
134 142 ly = 0;
135 143 }
136 144 }
137 145
138 146 /* fix point fmod */
139 147 n = ix - iy;
140 148 m = 0;
149 +
141 150 while (n--) {
142 151 hz = hx - hy;
143 152 lz = lx - ly;
153 +
144 154 if (lx < ly)
145 155 hz -= 1;
156 +
146 157 if (hz < 0) {
147 158 hx = hx + hx + (lx >> 31);
148 159 lx = lx + lx;
149 160 } else {
150 161 m += 1;
162 +
151 163 if ((hz | lz) == 0) { /* return sign(x)*0 */
152 164 if (n < 31)
153 165 m <<= 1 + n;
154 166 else
155 167 m = 0;
168 +
156 169 m &= 0x7fffffff;
157 170 *quo = sq >= 0 ? m : -m;
158 - return (Zero[(unsigned) sx >> 31]);
171 + return (Zero[(unsigned)sx >> 31]);
159 172 }
173 +
160 174 hx = hz + hz + (lz >> 31);
161 175 lx = lz + lz;
162 176 }
177 +
163 178 m += m;
164 179 }
180 +
165 181 hz = hx - hy;
166 182 lz = lx - ly;
183 +
167 184 if (lx < ly)
168 185 hz -= 1;
186 +
169 187 if (hz >= 0) {
170 188 hx = hz;
171 189 lx = lz;
172 190 m += 1;
173 191 }
192 +
174 193 m &= 0x7fffffff;
175 194 *quo = sq >= 0 ? m : -m;
176 195
177 196 /* convert back to floating value and restore the sign */
178 - if ((hx | lx) == 0) { /* return sign(x)*0 */
179 - return (Zero[(unsigned) sx >> 31]);
180 - }
197 + if ((hx | lx) == 0) /* return sign(x)*0 */
198 + return (Zero[(unsigned)sx >> 31]);
199 +
181 200 while (hx < 0x00100000) { /* normalize x */
182 201 hx = hx + hx + (lx >> 31);
183 202 lx = lx + lx;
184 203 iy -= 1;
185 204 }
186 - if (iy >= -1022) { /* normalize output */
205 +
206 + if (iy >= -1022) { /* normalize output */
187 207 hx = (hx - 0x00100000) | ((iy + 1023) << 20);
188 208 __HI(x) = hx | sx;
189 209 __LO(x) = lx;
190 210 } else { /* subnormal output */
191 211 n = -1022 - iy;
212 +
192 213 if (n <= 20) {
193 - lx = (lx >> n) | ((unsigned) hx << (32 - n));
214 + lx = (lx >> n) | ((unsigned)hx << (32 - n));
194 215 hx >>= n;
195 216 } else if (n <= 31) {
196 217 lx = (hx << (32 - n)) | (lx >> n);
197 218 hx = sx;
198 219 } else {
199 220 lx = hx >> (n - 32);
200 221 hx = sx;
201 222 }
223 +
202 224 __HI(x) = hx | sx;
203 225 __LO(x) = lx;
204 - x *= one; /* create necessary signal */
226 + x *= one; /* create necessary signal */
205 227 }
206 - return (x); /* exact output */
228 +
229 + return (x); /* exact output */
207 230 }
208 231
209 232 #define zero Zero[0]
210 233
211 234 double
212 -remquo(double x, double y, int *quo) {
235 +remquo(double x, double y, int *quo)
236 +{
213 237 int hx, hy, sx, sq;
214 238 double v;
215 239 unsigned ly;
216 240
217 - hx = __HI(x); /* high word of x */
218 - hy = __HI(y); /* high word of y */
219 - ly = __LO(y); /* low word of y */
220 - sx = hx & 0x80000000; /* sign of x */
241 + hx = __HI(x); /* high word of x */
242 + hy = __HI(y); /* high word of y */
243 + ly = __LO(y); /* low word of y */
244 + sx = hx & 0x80000000; /* sign of x */
221 245 sq = (hx ^ hy) & 0x80000000; /* sign of x/y */
222 - hx ^= sx; /* |x| */
223 - hy &= 0x7fffffff; /* |y| */
246 + hx ^= sx; /* |x| */
247 + hy &= 0x7fffffff; /* |y| */
224 248
225 249 /* purge off exception values */
226 250 *quo = 0;
251 +
227 252 if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */
228 253 (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */
229 254 return ((x * y) / (x * y));
230 255
231 256 y = fabs(y);
232 257 x = fabs(x);
258 +
233 259 if (hy <= 0x7fdfffff) {
234 260 x = fmodquo(x, y + y, quo);
235 261 *quo = ((*quo) & 0x3fffffff) << 1;
236 262 }
263 +
237 264 if (hy < 0x00200000) {
238 265 if (x + x > y) {
239 266 *quo += 1;
267 +
240 268 if (x == y)
241 269 x = zero;
242 270 else
243 271 x -= y;
272 +
244 273 if (x + x >= y) {
245 274 x -= y;
246 275 *quo += 1;
247 276 }
248 277 }
249 278 } else {
250 279 v = 0.5 * y;
280 +
251 281 if (x > v) {
252 282 *quo += 1;
283 +
253 284 if (x == y)
254 285 x = zero;
255 286 else
256 287 x -= y;
288 +
257 289 if (x >= v) {
258 290 x -= y;
259 291 *quo += 1;
260 292 }
261 293 }
262 294 }
295 +
263 296 if (sq != 0)
264 297 *quo = -(*quo);
298 +
265 299 return (sx == 0 ? x : -x);
266 300 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX