Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/m9x/remquol.c
+++ new/usr/src/lib/libm/common/m9x/remquol.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 __remquol = remquol
31 32
32 33 #include "libm.h"
33 34 #if defined(__SUNPRO_C)
34 35 #include <sunmath.h> /* fabsl */
35 36 #endif
36 -/* INDENT OFF */
37 -static const int
38 - is = -0x7fffffff - 1,
39 - im = 0x0000ffff,
40 - iu = 0x00010000;
41 37
38 +static const int is = -0x7fffffff - 1, im = 0x0000ffff, iu = 0x00010000;
42 39 static const long double zero = 0.0L, one = 1.0L;
43 -/* INDENT ON */
40 +
44 41
45 42 #if defined(__sparc)
46 -#define __H0(x) ((int *) &x)[0]
47 -#define __H1(x) ((int *) &x)[1]
48 -#define __H2(x) ((int *) &x)[2]
49 -#define __H3(x) ((int *) &x)[3]
43 +#define __H0(x) ((int *)&x)[0]
44 +#define __H1(x) ((int *)&x)[1]
45 +#define __H2(x) ((int *)&x)[2]
46 +#define __H3(x) ((int *)&x)[3]
50 47 #else
51 48 #error Unsupported architecture
52 49 #endif
53 50
54 51 /*
55 52 * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered
56 53 */
57 54 static long double
58 -fmodquol(long double x, long double y, int *quo) {
55 +fmodquol(long double x, long double y, int *quo)
56 +{
59 57 long double a, b;
60 58 int n, ix, iy, k, sx, sq, m;
61 59 int hx;
62 60 int x0, y0, z0, carry;
63 61 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
64 62
65 63 hx = __H0(x);
66 64 x1 = __H1(x);
67 65 x2 = __H2(x);
68 66 x3 = __H3(x);
69 67 y0 = __H0(y);
70 68 y1 = __H1(y);
|
↓ open down ↓ |
2 lines elided |
↑ open up ↑ |
71 69 y2 = __H2(y);
72 70 y3 = __H3(y);
73 71
74 72 sx = hx & is;
75 73 sq = (hx ^ y0) & is;
76 74 x0 = hx ^ sx;
77 75 y0 &= ~0x80000000;
78 76
79 77 a = fabsl(x);
80 78 b = fabsl(y);
79 +
81 80 if (a <= b) {
82 - if (a < b)
81 + if (a < b) {
83 82 return (x);
84 - else {
83 + } else {
85 84 *quo = 1 + (sq >> 30);
86 85 return (zero * x);
87 86 }
88 87 }
88 +
89 89 /* determine ix = ilogbl(x) */
90 - if (x0 < iu) { /* subnormal x */
90 + if (x0 < iu) { /* subnormal x */
91 91 ix = 0;
92 92 ix = -16382;
93 +
93 94 while (x0 == 0) {
94 95 ix -= 16;
95 96 x0 = x1 >> 16;
96 97 x1 = (x1 << 16) | (x2 >> 16);
97 98 x2 = (x2 << 16) | (x3 >> 16);
98 99 x3 = (x3 << 16);
99 100 }
101 +
100 102 while (x0 < iu) {
101 103 ix -= 1;
102 104 x0 = (x0 << 1) | (x1 >> 31);
103 105 x1 = (x1 << 1) | (x2 >> 31);
104 106 x2 = (x2 << 1) | (x3 >> 31);
105 107 x3 <<= 1;
106 108 }
107 109 } else {
108 110 ix = (x0 >> 16) - 16383;
109 111 x0 = iu | (x0 & im);
110 112 }
111 113
112 114 /* determine iy = ilogbl(y) */
113 - if (y0 < iu) { /* subnormal y */
115 + if (y0 < iu) { /* subnormal y */
114 116 iy = -16382;
117 +
115 118 while (y0 == 0) {
116 119 iy -= 16;
117 120 y0 = y1 >> 16;
118 121 y1 = (y1 << 16) | (y2 >> 16);
119 122 y2 = (y2 << 16) | (y3 >> 16);
120 123 y3 = (y3 << 16);
121 124 }
125 +
122 126 while (y0 < iu) {
123 127 iy -= 1;
124 128 y0 = (y0 << 1) | (y1 >> 31);
125 129 y1 = (y1 << 1) | (y2 >> 31);
126 130 y2 = (y2 << 1) | (y3 >> 31);
127 131 y3 <<= 1;
128 132 }
129 133 } else {
130 134 iy = (y0 >> 16) - 16383;
131 135 y0 = iu | (y0 & im);
132 136 }
133 137
134 -
135 138 /* fix point fmod */
136 139 n = ix - iy;
137 140 m = 0;
141 +
138 142 while (n--) {
139 143 while (x0 == 0 && n >= 16) {
140 144 m <<= 16;
141 145 n -= 16;
142 146 x0 = x1 >> 16;
143 147 x1 = (x1 << 16) | (x2 >> 16);
144 148 x2 = (x2 << 16) | (x3 >> 16);
145 149 x3 = (x3 << 16);
146 150 }
151 +
147 152 while (x0 < iu && n >= 1) {
148 153 m += m;
149 154 n -= 1;
150 155 x0 = (x0 << 1) | (x1 >> 31);
151 156 x1 = (x1 << 1) | (x2 >> 31);
152 157 x2 = (x2 << 1) | (x3 >> 31);
153 158 x3 = (x3 << 1);
154 159 }
160 +
155 161 carry = 0;
156 162 z3 = x3 - y3;
157 163 carry = z3 > x3;
164 +
158 165 if (carry == 0) {
159 166 z2 = x2 - y2;
160 167 carry = z2 > x2;
161 168 } else {
162 169 z2 = x2 - y2 - 1;
163 170 carry = z2 >= x2;
164 171 }
172 +
165 173 if (carry == 0) {
166 174 z1 = x1 - y1;
167 175 carry = z1 > x1;
168 176 } else {
169 177 z1 = x1 - y1 - 1;
170 178 carry = z1 >= x1;
171 179 }
180 +
172 181 z0 = x0 - y0 - carry;
173 - if (z0 < 0) { /* double x */
182 +
183 + if (z0 < 0) { /* double x */
174 184 x0 = x0 + x0 + ((x1 & is) != 0);
175 185 x1 = x1 + x1 + ((x2 & is) != 0);
176 186 x2 = x2 + x2 + ((x3 & is) != 0);
177 187 x3 = x3 + x3;
178 188 m += m;
179 189 } else {
180 190 m += 1;
191 +
181 192 if (z0 == 0) {
182 193 if ((z1 | z2 | z3) == 0) {
183 194 /* 0: we are done */
184 195 if (n < 31)
185 196 m <<= (1 + n);
186 197 else
187 198 m = 0;
199 +
188 200 m &= ~0x80000000;
189 201 *quo = sq >= 0 ? m : -m;
190 202 __H0(a) = hx & is;
191 203 __H1(a) = __H2(a) = __H3(a) = 0;
192 204 return (a);
193 205 }
194 206 }
207 +
195 208 /* x = z << 1 */
196 209 z0 = z0 + z0 + ((z1 & is) != 0);
197 210 z1 = z1 + z1 + ((z2 & is) != 0);
198 211 z2 = z2 + z2 + ((z3 & is) != 0);
199 212 z3 = z3 + z3;
200 213 x0 = z0;
201 214 x1 = z1;
202 215 x2 = z2;
203 216 x3 = z3;
204 217 m += m;
205 218 }
206 219 }
220 +
207 221 carry = 0;
208 222 z3 = x3 - y3;
209 223 carry = z3 > x3;
224 +
210 225 if (carry == 0) {
211 226 z2 = x2 - y2;
212 227 carry = z2 > x2;
213 228 } else {
214 229 z2 = x2 - y2 - 1;
215 230 carry = z2 >= x2;
216 231 }
232 +
217 233 if (carry == 0) {
218 234 z1 = x1 - y1;
219 235 carry = z1 > x1;
220 236 } else {
221 237 z1 = x1 - y1 - 1;
222 238 carry = z1 >= x1;
223 239 }
240 +
224 241 z0 = x0 - y0 - carry;
242 +
225 243 if (z0 >= 0) {
226 244 x0 = z0;
227 245 x1 = z1;
228 246 x2 = z2;
229 247 x3 = z3;
230 248 m += 1;
231 249 }
250 +
232 251 m &= ~0x80000000;
233 252 *quo = sq >= 0 ? m : -m;
234 253
235 254 /* convert back to floating value and restore the sign */
236 255 if ((x0 | x1 | x2 | x3) == 0) {
237 256 __H0(a) = hx & is;
238 257 __H1(a) = __H2(a) = __H3(a) = 0;
239 258 return (a);
240 259 }
260 +
241 261 while (x0 < iu) {
242 262 if (x0 == 0) {
243 263 iy -= 16;
244 264 x0 = x1 >> 16;
245 265 x1 = (x1 << 16) | (x2 >> 16);
246 266 x2 = (x2 << 16) | (x3 >> 16);
247 267 x3 = (x3 << 16);
248 268 } else {
249 269 x0 = x0 + x0 + ((x1 & is) != 0);
250 270 x1 = x1 + x1 + ((x2 & is) != 0);
251 271 x2 = x2 + x2 + ((x3 & is) != 0);
252 272 x3 = x3 + x3;
|
↓ open down ↓ |
2 lines elided |
↑ open up ↑ |
253 273 iy -= 1;
254 274 }
255 275 }
256 276
257 277 /* normalize output */
258 278 if (iy >= -16382) {
259 279 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
260 280 __H1(a) = x1;
261 281 __H2(a) = x2;
262 282 __H3(a) = x3;
263 - } else { /* subnormal output */
283 + } else { /* subnormal output */
264 284 n = -16382 - iy;
265 285 k = n & 31;
286 +
266 287 if (k <= 16) {
267 288 x3 = (x2 << (32 - k)) | (x3 >> k);
268 289 x2 = (x1 << (32 - k)) | (x2 >> k);
269 290 x1 = (x0 << (32 - k)) | (x1 >> k);
270 291 x0 >>= k;
271 292 } else {
272 293 x3 = (x2 << (32 - k)) | (x3 >> k);
273 294 x2 = (x1 << (32 - k)) | (x2 >> k);
274 295 x1 = (x0 << (32 - k)) | (x1 >> k);
275 296 x0 = 0;
276 297 }
298 +
277 299 while (n >= 32) {
278 300 n -= 32;
279 301 x3 = x2;
280 302 x2 = x1;
281 303 x1 = x0;
282 304 x0 = 0;
283 305 }
306 +
284 307 __H0(a) = x0 | sx;
285 308 __H1(a) = x1;
286 309 __H2(a) = x2;
287 310 __H3(a) = x3;
288 311 a *= one;
289 312 }
313 +
290 314 return (a);
291 315 }
292 316
293 317 long double
294 -remquol(long double x, long double y, int *quo) {
318 +remquol(long double x, long double y, int *quo)
319 +{
295 320 int hx, hy, sx, sq;
296 321 long double v;
297 322
298 - hx = __H0(x); /* high word of x */
299 - hy = __H0(y); /* high word of y */
300 - sx = hx & is; /* sign of x */
301 - sq = (hx ^ hy) & is; /* sign of x/y */
302 - hx ^= sx; /* |x| */
323 + hx = __H0(x); /* high word of x */
324 + hy = __H0(y); /* high word of y */
325 + sx = hx & is; /* sign of x */
326 + sq = (hx ^ hy) & is; /* sign of x/y */
327 + hx ^= sx; /* |x| */
303 328 hy &= ~0x80000000;
304 329
305 330 /* purge off exception values */
306 331 *quo = 0;
332 +
307 333 /* y=0, y is NaN, x is NaN or inf */
308 334 if (y == 0.0L || y != y || hx >= 0x7fff0000)
309 335 return ((x * y) / (x * y));
310 336
311 337 y = fabsl(y);
312 338 x = fabsl(x);
339 +
313 340 if (hy <= 0x7ffdffff) {
314 341 x = fmodquol(x, y + y, quo);
315 342 *quo = ((*quo) & 0x3fffffff) << 1;
316 343 }
344 +
317 345 if (hy < 0x00020000) {
318 346 if (x + x > y) {
319 347 *quo += 1;
348 +
320 349 if (x == y)
321 350 x = zero;
322 351 else
323 352 x -= y;
353 +
324 354 if (x + x >= y) {
325 355 x -= y;
326 356 *quo += 1;
327 357 }
328 358 }
329 359 } else {
330 360 v = 0.5L * y;
361 +
331 362 if (x > v) {
332 363 *quo += 1;
364 +
333 365 if (x == y)
334 366 x = zero;
335 367 else
336 368 x -= y;
369 +
337 370 if (x >= v) {
338 371 x -= y;
339 372 *quo += 1;
340 373 }
341 374 }
342 375 }
376 +
343 377 if (sq != 0)
344 378 *quo = -(*quo);
379 +
345 380 return (sx == 0 ? x : -x);
346 381 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX