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