Go to most recent revision | Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
4349 | Serge | 1 | |
2 | /* |
||
3 | * ==================================================== |
||
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
||
5 | * |
||
6 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
||
7 | * Permission to use, copy, modify, and distribute this |
||
8 | * software is freely granted, provided that this notice |
||
9 | * is preserved. |
||
10 | * ==================================================== |
||
11 | */ |
||
12 | |||
13 | |||
14 | #include |
||
15 | #include |
||
16 | #include |
||
17 | |||
18 | |||
19 | #define _XOPEN_MODE |
||
20 | |||
21 | |||
22 | number, and many need to know whether the result of an operation will |
||
23 | overflow. These conditions depend on whether the largest exponent is |
||
24 | used for NaNs & infinities, or whether it's used for finite numbers. The |
||
25 | macros below wrap up that kind of information: |
||
26 | |||
27 | |||
28 | True if a positive float with bitmask X is finite. |
||
29 | |||
30 | |||
31 | True if a positive float with bitmask X is not a number. |
||
32 | |||
33 | |||
34 | True if a positive float with bitmask X is +infinity. |
||
35 | |||
36 | |||
37 | The bitmask of FLT_MAX. |
||
38 | |||
39 | |||
40 | The bitmask of FLT_MAX/2. |
||
41 | |||
42 | |||
43 | The bitmask of the largest finite exponent (129 if the largest |
||
44 | exponent is used for finite numbers, 128 otherwise). |
||
45 | |||
46 | |||
47 | The bitmask of log(FLT_MAX), rounded down. This value is the largest |
||
48 | input that can be passed to exp() without producing overflow. |
||
49 | |||
50 | |||
51 | The bitmask of log(2*FLT_MAX), rounded down. This value is the |
||
52 | largest input than can be passed to cosh() without producing |
||
53 | overflow. |
||
54 | |||
55 | |||
56 | The largest biased exponent that can be used for finite numbers |
||
57 | (255 if the largest exponent is used for finite numbers, 254 |
||
58 | otherwise) */ |
||
59 | |||
60 | |||
61 | #define FLT_UWORD_IS_FINITE(x) 1 |
||
62 | #define FLT_UWORD_IS_NAN(x) 0 |
||
63 | #define FLT_UWORD_IS_INFINITE(x) 0 |
||
64 | #define FLT_UWORD_MAX 0x7fffffff |
||
65 | #define FLT_UWORD_EXP_MAX 0x43010000 |
||
66 | #define FLT_UWORD_LOG_MAX 0x42b2d4fc |
||
67 | #define FLT_UWORD_LOG_2MAX 0x42b437e0 |
||
68 | #define HUGE ((float)0X1.FFFFFEP128) |
||
69 | #else |
||
70 | #define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L) |
||
71 | #define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L) |
||
72 | #define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L) |
||
73 | #define FLT_UWORD_MAX 0x7f7fffffL |
||
74 | #define FLT_UWORD_EXP_MAX 0x43000000 |
||
75 | #define FLT_UWORD_LOG_MAX 0x42b17217 |
||
76 | #define FLT_UWORD_LOG_2MAX 0x42b2d4fc |
||
77 | #define HUGE ((float)3.40282346638528860e+38) |
||
78 | #endif |
||
79 | #define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23)) |
||
80 | #define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23) |
||
81 | |||
82 | |||
83 | on whether the target supports denormals or not: |
||
84 | |||
85 | |||
86 | True if a positive float with bitmask X is +0. Without denormals, |
||
87 | any float with a zero exponent is a +0 representation. With |
||
88 | denormals, the only +0 representation is a 0 bitmask. |
||
89 | |||
90 | |||
91 | True if a non-zero positive float with bitmask X is subnormal. |
||
92 | (Routines should check for zeros first.) |
||
93 | |||
94 | |||
95 | The bitmask of the smallest float above +0. Call this number |
||
96 | REAL_FLT_MIN... |
||
97 | |||
98 | |||
99 | The bitmask of the float representation of REAL_FLT_MIN's exponent. |
||
100 | |||
101 | |||
102 | The bitmask of |log(REAL_FLT_MIN)|, rounding down. |
||
103 | |||
104 | |||
105 | REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported, |
||
106 | -22 if they are). |
||
107 | */ |
||
108 | |||
109 | |||
110 | #define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L) |
||
111 | #define FLT_UWORD_IS_SUBNORMAL(x) 0 |
||
112 | #define FLT_UWORD_MIN 0x00800000 |
||
113 | #define FLT_UWORD_EXP_MIN 0x42fc0000 |
||
114 | #define FLT_UWORD_LOG_MIN 0x42aeac50 |
||
115 | #define FLT_SMALLEST_EXP 1 |
||
116 | #else |
||
117 | #define FLT_UWORD_IS_ZERO(x) ((x)==0) |
||
118 | #define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L) |
||
119 | #define FLT_UWORD_MIN 0x00000001 |
||
120 | #define FLT_UWORD_EXP_MIN 0x43160000 |
||
121 | #define FLT_UWORD_LOG_MIN 0x42cff1b5 |
||
122 | #define FLT_SMALLEST_EXP -22 |
||
123 | #endif |
||
124 | |||
125 | |||
126 | #undef __P |
||
127 | #define __P(p) p |
||
128 | #else |
||
129 | #define __P(p) () |
||
130 | #endif |
||
131 | |||
132 | |||
133 | * set X_TLOSS = pi*2**52, which is possibly defined in |
||
134 | * (one may replace the following line by "#include |
||
135 | */ |
||
136 | |||
137 | |||
138 | |||
139 | |||
140 | |||
141 | |||
142 | extern double scalb __P((double, int)); |
||
143 | #else |
||
144 | extern double scalb __P((double, double)); |
||
145 | #endif |
||
146 | extern double significand __P((double)); |
||
147 | |||
148 | |||
149 | extern double __ieee754_sqrt __P((double)); |
||
150 | extern double __ieee754_acos __P((double)); |
||
151 | extern double __ieee754_acosh __P((double)); |
||
152 | extern double __ieee754_log __P((double)); |
||
153 | extern double __ieee754_atanh __P((double)); |
||
154 | extern double __ieee754_asin __P((double)); |
||
155 | extern double __ieee754_atan2 __P((double,double)); |
||
156 | extern double __ieee754_exp __P((double)); |
||
157 | extern double __ieee754_cosh __P((double)); |
||
158 | extern double __ieee754_fmod __P((double,double)); |
||
159 | extern double __ieee754_pow __P((double,double)); |
||
160 | extern double __ieee754_lgamma_r __P((double,int *)); |
||
161 | extern double __ieee754_gamma_r __P((double,int *)); |
||
162 | extern double __ieee754_log10 __P((double)); |
||
163 | extern double __ieee754_sinh __P((double)); |
||
164 | extern double __ieee754_hypot __P((double,double)); |
||
165 | extern double __ieee754_j0 __P((double)); |
||
166 | extern double __ieee754_j1 __P((double)); |
||
167 | extern double __ieee754_y0 __P((double)); |
||
168 | extern double __ieee754_y1 __P((double)); |
||
169 | extern double __ieee754_jn __P((int,double)); |
||
170 | extern double __ieee754_yn __P((int,double)); |
||
171 | extern double __ieee754_remainder __P((double,double)); |
||
172 | extern __int32_t __ieee754_rem_pio2 __P((double,double*)); |
||
173 | #ifdef _SCALB_INT |
||
174 | extern double __ieee754_scalb __P((double,int)); |
||
175 | #else |
||
176 | extern double __ieee754_scalb __P((double,double)); |
||
177 | #endif |
||
178 | |||
179 | |||
180 | extern double __kernel_standard __P((double,double,int)); |
||
181 | extern double __kernel_sin __P((double,double,int)); |
||
182 | extern double __kernel_cos __P((double,double)); |
||
183 | extern double __kernel_tan __P((double,double,int)); |
||
184 | extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*)); |
||
185 | |||
186 | |||
187 | #ifdef _SCALB_INT |
||
188 | extern float scalbf __P((float, int)); |
||
189 | #else |
||
190 | extern float scalbf __P((float, float)); |
||
191 | #endif |
||
192 | extern float significandf __P((float)); |
||
193 | |||
194 | |||
195 | extern float __ieee754_sqrtf __P((float)); |
||
196 | extern float __ieee754_acosf __P((float)); |
||
197 | extern float __ieee754_acoshf __P((float)); |
||
198 | extern float __ieee754_logf __P((float)); |
||
199 | extern float __ieee754_atanhf __P((float)); |
||
200 | extern float __ieee754_asinf __P((float)); |
||
201 | extern float __ieee754_atan2f __P((float,float)); |
||
202 | extern float __ieee754_expf __P((float)); |
||
203 | extern float __ieee754_coshf __P((float)); |
||
204 | extern float __ieee754_fmodf __P((float,float)); |
||
205 | extern float __ieee754_powf __P((float,float)); |
||
206 | extern float __ieee754_lgammaf_r __P((float,int *)); |
||
207 | extern float __ieee754_gammaf_r __P((float,int *)); |
||
208 | extern float __ieee754_log10f __P((float)); |
||
209 | extern float __ieee754_sinhf __P((float)); |
||
210 | extern float __ieee754_hypotf __P((float,float)); |
||
211 | extern float __ieee754_j0f __P((float)); |
||
212 | extern float __ieee754_j1f __P((float)); |
||
213 | extern float __ieee754_y0f __P((float)); |
||
214 | extern float __ieee754_y1f __P((float)); |
||
215 | extern float __ieee754_jnf __P((int,float)); |
||
216 | extern float __ieee754_ynf __P((int,float)); |
||
217 | extern float __ieee754_remainderf __P((float,float)); |
||
218 | extern __int32_t __ieee754_rem_pio2f __P((float,float*)); |
||
219 | #ifdef _SCALB_INT |
||
220 | extern float __ieee754_scalbf __P((float,int)); |
||
221 | #else |
||
222 | extern float __ieee754_scalbf __P((float,float)); |
||
223 | #endif |
||
224 | |||
225 | |||
226 | extern float __kernel_sinf __P((float,float,int)); |
||
227 | extern float __kernel_cosf __P((float,float)); |
||
228 | extern float __kernel_tanf __P((float,float,int)); |
||
229 | extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*)); |
||
230 | |||
231 | |||
232 | n0 = ((*(int*)&one)>>29)^1; * index of high word * |
||
233 | ix0 = *(n0+(int*)&x); * high word of x * |
||
234 | ix1 = *((1-n0)+(int*)&x); * low word of x * |
||
235 | to dig two 32 bit words out of the 64 bit IEEE floating point |
||
236 | value. That is non-ANSI, and, moreover, the gcc instruction |
||
237 | scheduler gets it wrong. We instead use the following macros. |
||
238 | Unlike the original code, we determine the endianness at compile |
||
239 | time, not at run time; I don't see much benefit to selecting |
||
240 | endianness at run time. */ |
||
241 | |||
242 | |||
243 | #ifndef __IEEE_LITTLE_ENDIAN |
||
244 | #error Must define endianness |
||
245 | #endif |
||
246 | #endif |
||
247 | |||
248 | |||
249 | ints. */ |
||
250 | |||
251 | |||
252 | |||
253 | |||
254 | { |
||
255 | double value; |
||
256 | struct |
||
257 | { |
||
258 | __uint32_t msw; |
||
259 | __uint32_t lsw; |
||
260 | } parts; |
||
261 | } ieee_double_shape_type; |
||
262 | |||
263 | |||
264 | |||
265 | |||
266 | |||
267 | |||
268 | { |
||
269 | double value; |
||
270 | struct |
||
271 | { |
||
272 | __uint32_t lsw; |
||
273 | __uint32_t msw; |
||
274 | } parts; |
||
275 | } ieee_double_shape_type; |
||
276 | |||
277 | |||
278 | |||
279 | |||
280 | |||
281 | |||
282 | do { \ |
||
283 | ieee_double_shape_type ew_u; \ |
||
284 | ew_u.value = (d); \ |
||
285 | (ix0) = ew_u.parts.msw; \ |
||
286 | (ix1) = ew_u.parts.lsw; \ |
||
287 | } while (0) |
||
288 | |||
289 | |||
290 | |||
291 | |||
292 | do { \ |
||
293 | ieee_double_shape_type gh_u; \ |
||
294 | gh_u.value = (d); \ |
||
295 | (i) = gh_u.parts.msw; \ |
||
296 | } while (0) |
||
297 | |||
298 | |||
299 | |||
300 | |||
301 | do { \ |
||
302 | ieee_double_shape_type gl_u; \ |
||
303 | gl_u.value = (d); \ |
||
304 | (i) = gl_u.parts.lsw; \ |
||
305 | } while (0) |
||
306 | |||
307 | |||
308 | |||
309 | |||
310 | do { \ |
||
311 | ieee_double_shape_type iw_u; \ |
||
312 | iw_u.parts.msw = (ix0); \ |
||
313 | iw_u.parts.lsw = (ix1); \ |
||
314 | (d) = iw_u.value; \ |
||
315 | } while (0) |
||
316 | |||
317 | |||
318 | |||
319 | |||
320 | do { \ |
||
321 | ieee_double_shape_type sh_u; \ |
||
322 | sh_u.value = (d); \ |
||
323 | sh_u.parts.msw = (v); \ |
||
324 | (d) = sh_u.value; \ |
||
325 | } while (0) |
||
326 | |||
327 | |||
328 | |||
329 | |||
330 | do { \ |
||
331 | ieee_double_shape_type sl_u; \ |
||
332 | sl_u.value = (d); \ |
||
333 | sl_u.parts.lsw = (v); \ |
||
334 | (d) = sl_u.value; \ |
||
335 | } while (0) |
||
336 | |||
337 | |||
338 | int. */ |
||
339 | |||
340 | |||
341 | { |
||
342 | float value; |
||
343 | __uint32_t word; |
||
344 | } ieee_float_shape_type; |
||
345 | |||
346 | |||
347 | |||
348 | |||
349 | do { \ |
||
350 | ieee_float_shape_type gf_u; \ |
||
351 | gf_u.value = (d); \ |
||
352 | (i) = gf_u.word; \ |
||
353 | } while (0) |
||
354 | |||
355 | |||
356 | |||
357 | |||
358 | do { \ |
||
359 | ieee_float_shape_type sf_u; \ |
||
360 | sf_u.word = (i); \ |
||
361 | (d) = sf_u.value; \ |
||
362 | } while (0) |
||
363 | |||
364 | |||
365 | of a shift is exactly equal to the size of the shifted operand. */ |
||
366 | |||
367 | |||
368 | (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0) |
||
369 | |||
370 | |||
371 | (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0) |
||
372 | |||
373 | |||
374 | |||
375 | |||
376 | * Quoting from ISO/IEC 9899:TC2: |
||
377 | * |
||
378 | * 6.2.5.13 Types |
||
379 | * Each complex type has the same representation and alignment requirements as |
||
380 | * an array type containing exactly two elements of the corresponding real type; |
||
381 | * the first element is equal to the real part, and the second element to the |
||
382 | * imaginary part, of the complex number. |
||
383 | */ |
||
384 | typedef union { |
||
385 | float complex z; |
||
386 | float parts[2]; |
||
387 | } float_complex; |
||
388 | |||
389 | |||
390 | double complex z; |
||
391 | double parts[2]; |
||
392 | } double_complex; |
||
393 | |||
394 | |||
395 | long double complex z; |
||
396 | long double parts[2]; |
||
397 | } long_double_complex; |
||
398 | |||
399 | |||
400 | #define IMAG_PART(z) ((z).parts[1]) |
||
401 | |||
402 | |||
403 | #define>0x00800000L) |
||
404 |