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