2 /* @(#)fdlibm.h 5.1 93/09/24 */
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
11 * ====================================================
17 #define __int32_t int32_t
18 #define __uint32_t uint32_t
20 #define __ieee754_acosf acosf
21 #define __ieee754_sqrtf sqrtf
23 /* REDHAT LOCAL: Include files. */
25 /* #include <sys/types.h> */
26 #include <machine/ieeefp.h>
28 /* REDHAT LOCAL: Default to XOPEN_MODE. */
31 /* Most routines need to check whether a float is finite, infinite, or not a
32 number, and many need to know whether the result of an operation will
33 overflow. These conditions depend on whether the largest exponent is
34 used for NaNs & infinities, or whether it's used for finite numbers. The
35 macros below wrap up that kind of information:
37 FLT_UWORD_IS_FINITE(X)
38 True if a positive float with bitmask X is finite.
41 True if a positive float with bitmask X is not a number.
43 FLT_UWORD_IS_INFINITE(X)
44 True if a positive float with bitmask X is +infinity.
47 The bitmask of FLT_MAX.
50 The bitmask of FLT_MAX/2.
53 The bitmask of the largest finite exponent (129 if the largest
54 exponent is used for finite numbers, 128 otherwise).
57 The bitmask of log(FLT_MAX), rounded down. This value is the largest
58 input that can be passed to exp() without producing overflow.
61 The bitmask of log(2*FLT_MAX), rounded down. This value is the
62 largest input than can be passed to cosh() without producing
66 The largest biased exponent that can be used for finite numbers
67 (255 if the largest exponent is used for finite numbers, 254
70 #ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
71 #define FLT_UWORD_IS_FINITE(x) 1
72 #define FLT_UWORD_IS_NAN(x) 0
73 #define FLT_UWORD_IS_INFINITE(x) 0
74 #define FLT_UWORD_MAX 0x7fffffff
75 #define FLT_UWORD_EXP_MAX 0x43010000
76 #define FLT_UWORD_LOG_MAX 0x42b2d4fc
77 #define FLT_UWORD_LOG_2MAX 0x42b437e0
78 #define HUGE ((float)0X1.FFFFFEP128)
80 #define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
81 #define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
82 #define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
83 #define FLT_UWORD_MAX 0x7f7fffffL
84 #define FLT_UWORD_EXP_MAX 0x43000000
85 #define FLT_UWORD_LOG_MAX 0x42b17217
86 #define FLT_UWORD_LOG_2MAX 0x42b2d4fc
87 #define HUGE ((float)3.40282346638528860e+38)
89 #define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23))
90 #define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
92 /* Many routines check for zero and subnormal numbers. Such things depend
93 on whether the target supports denormals or not:
96 True if a positive float with bitmask X is +0. Without denormals,
97 any float with a zero exponent is a +0 representation. With
98 denormals, the only +0 representation is a 0 bitmask.
100 FLT_UWORD_IS_SUBNORMAL(X)
101 True if a non-zero positive float with bitmask X is subnormal.
102 (Routines should check for zeros first.)
105 The bitmask of the smallest float above +0. Call this number
109 The bitmask of the float representation of REAL_FLT_MIN's exponent.
112 The bitmask of |log(REAL_FLT_MIN)|, rounding down.
115 REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
119 #ifdef _FLT_NO_DENORMALS
120 #define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
121 #define FLT_UWORD_IS_SUBNORMAL(x) 0
122 #define FLT_UWORD_MIN 0x00800000
123 #define FLT_UWORD_EXP_MIN 0x42fc0000
124 #define FLT_UWORD_LOG_MIN 0x42aeac50
125 #define FLT_SMALLEST_EXP 1
127 #define FLT_UWORD_IS_ZERO(x) ((x)==0)
128 #define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
129 #define FLT_UWORD_MIN 0x00000001
130 #define FLT_UWORD_EXP_MIN 0x43160000
131 #define FLT_UWORD_LOG_MIN 0x42cff1b5
132 #define FLT_SMALLEST_EXP -22
143 * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
144 * (one may replace the following line by "#include <values.h>")
147 #define X_TLOSS 1.41484755040568800000e+16
149 /* Functions that are not documented, and are not in <math.h>. */
152 extern double scalb __P((double, int));
154 extern double scalb __P((double, double));
156 extern double significand __P((double));
158 /* ieee style elementary functions */
159 extern double __ieee754_sqrt __P((double));
160 extern double __ieee754_acos __P((double));
161 extern double __ieee754_acosh __P((double));
162 extern double __ieee754_log __P((double));
163 extern double __ieee754_atanh __P((double));
164 extern double __ieee754_asin __P((double));
165 extern double __ieee754_atan2 __P((double,double));
166 extern double __ieee754_exp __P((double));
167 extern double __ieee754_cosh __P((double));
168 extern double __ieee754_fmod __P((double,double));
169 extern double __ieee754_pow __P((double,double));
170 extern double __ieee754_lgamma_r __P((double,int *));
171 extern double __ieee754_gamma_r __P((double,int *));
172 extern double __ieee754_log10 __P((double));
173 extern double __ieee754_sinh __P((double));
174 extern double __ieee754_hypot __P((double,double));
175 extern double __ieee754_j0 __P((double));
176 extern double __ieee754_j1 __P((double));
177 extern double __ieee754_y0 __P((double));
178 extern double __ieee754_y1 __P((double));
179 extern double __ieee754_jn __P((int,double));
180 extern double __ieee754_yn __P((int,double));
181 extern double __ieee754_remainder __P((double,double));
182 extern __int32_t __ieee754_rem_pio2 __P((double,double*));
184 extern double __ieee754_scalb __P((double,int));
186 extern double __ieee754_scalb __P((double,double));
189 /* fdlibm kernel function */
190 extern double __kernel_standard __P((double,double,int));
191 extern double __kernel_sin __P((double,double,int));
192 extern double __kernel_cos __P((double,double));
193 extern double __kernel_tan __P((double,double,int));
194 extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
196 /* Undocumented float functions. */
198 extern float scalbf __P((float, int));
200 extern float scalbf __P((float, float));
202 extern float significandf __P((float));
204 /* ieee style elementary float functions */
205 extern float __ieee754_sqrtf __P((float));
206 extern float __ieee754_acosf __P((float));
207 extern float __ieee754_acoshf __P((float));
208 extern float __ieee754_logf __P((float));
209 extern float __ieee754_atanhf __P((float));
210 extern float __ieee754_asinf __P((float));
211 extern float __ieee754_atan2f __P((float,float));
212 extern float __ieee754_expf __P((float));
213 extern float __ieee754_coshf __P((float));
214 extern float __ieee754_fmodf __P((float,float));
215 extern float __ieee754_powf __P((float,float));
216 extern float __ieee754_lgammaf_r __P((float,int *));
217 extern float __ieee754_gammaf_r __P((float,int *));
218 extern float __ieee754_log10f __P((float));
219 extern float __ieee754_sinhf __P((float));
220 extern float __ieee754_hypotf __P((float,float));
221 extern float __ieee754_j0f __P((float));
222 extern float __ieee754_j1f __P((float));
223 extern float __ieee754_y0f __P((float));
224 extern float __ieee754_y1f __P((float));
225 extern float __ieee754_jnf __P((int,float));
226 extern float __ieee754_ynf __P((int,float));
227 extern float __ieee754_remainderf __P((float,float));
228 extern __int32_t __ieee754_rem_pio2f __P((float,float*));
230 extern float __ieee754_scalbf __P((float,int));
232 extern float __ieee754_scalbf __P((float,float));
235 /* float versions of fdlibm kernel functions */
236 extern float __kernel_sinf __P((float,float,int));
237 extern float __kernel_cosf __P((float,float));
238 extern float __kernel_tanf __P((float,float,int));
239 extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
241 /* The original code used statements like
242 n0 = ((*(int*)&one)>>29)^1; * index of high word *
243 ix0 = *(n0+(int*)&x); * high word of x *
244 ix1 = *((1-n0)+(int*)&x); * low word of x *
245 to dig two 32 bit words out of the 64 bit IEEE floating point
246 value. That is non-ANSI, and, moreover, the gcc instruction
247 scheduler gets it wrong. We instead use the following macros.
248 Unlike the original code, we determine the endianness at compile
249 time, not at run time; I don't see much benefit to selecting
250 endianness at run time. */
252 #ifndef __IEEE_BIG_ENDIAN
253 #ifndef __IEEE_LITTLE_ENDIAN
254 #error Must define endianness
258 /* A union which permits us to convert between a double and two 32 bit
261 #ifdef __IEEE_BIG_ENDIAN
271 } ieee_double_shape_type;
275 #ifdef __IEEE_LITTLE_ENDIAN
285 } ieee_double_shape_type;
289 /* Get two 32 bit ints from a double. */
291 #define EXTRACT_WORDS(ix0,ix1,d) \
293 ieee_double_shape_type ew_u; \
295 (ix0) = ew_u.parts.msw; \
296 (ix1) = ew_u.parts.lsw; \
299 /* Get the more significant 32 bit int from a double. */
301 #define GET_HIGH_WORD(i,d) \
303 ieee_double_shape_type gh_u; \
305 (i) = gh_u.parts.msw; \
308 /* Get the less significant 32 bit int from a double. */
310 #define GET_LOW_WORD(i,d) \
312 ieee_double_shape_type gl_u; \
314 (i) = gl_u.parts.lsw; \
317 /* Set a double from two 32 bit ints. */
319 #define INSERT_WORDS(d,ix0,ix1) \
321 ieee_double_shape_type iw_u; \
322 iw_u.parts.msw = (ix0); \
323 iw_u.parts.lsw = (ix1); \
327 /* Set the more significant 32 bits of a double from an int. */
329 #define SET_HIGH_WORD(d,v) \
331 ieee_double_shape_type sh_u; \
333 sh_u.parts.msw = (v); \
337 /* Set the less significant 32 bits of a double from an int. */
339 #define SET_LOW_WORD(d,v) \
341 ieee_double_shape_type sl_u; \
343 sl_u.parts.lsw = (v); \
347 /* A union which permits us to convert between a float and a 32 bit
354 } ieee_float_shape_type;
356 /* Get a 32 bit int from a float. */
358 #define GET_FLOAT_WORD(i,d) \
360 ieee_float_shape_type gf_u; \
365 /* Set a float from a 32 bit int. */
367 #define SET_FLOAT_WORD(d,i) \
369 ieee_float_shape_type sf_u; \
374 /* Macros to avoid undefined behaviour that can arise if the amount
375 of a shift is exactly equal to the size of the shifted operand. */
377 #define SAFE_LEFT_SHIFT(op,amt) \
378 (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0)
380 #define SAFE_RIGHT_SHIFT(op,amt) \
381 (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0)
386 * Quoting from ISO/IEC 9899:TC2:
389 * Each complex type has the same representation and alignment requirements as
390 * an array type containing exactly two elements of the corresponding real type;
391 * the first element is equal to the real part, and the second element to the
392 * imaginary part, of the complex number.
405 long double complex z;
406 long double parts[2];
407 } long_double_complex;
409 #define REAL_PART(z) ((z).parts[0])
410 #define IMAG_PART(z) ((z).parts[1])
412 #endif /* _COMPLEX_H */