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
143extern double scalb __P((double, int));
144#else
145extern double scalb __P((double, double));
146#endif
147extern double significand __P((double));
148
149/* ieee style elementary functions */
150extern double __ieee754_sqrt __P((double));
151extern double __ieee754_acos __P((double));
152extern double __ieee754_acosh __P((double));
153extern double __ieee754_log __P((double));
154extern double __ieee754_atanh __P((double));
155extern double __ieee754_asin __P((double));
156extern double __ieee754_atan2 __P((double,double));
157extern double __ieee754_exp __P((double));
158extern double __ieee754_cosh __P((double));
159extern double __ieee754_fmod __P((double,double));
160extern double __ieee754_pow __P((double,double));
161extern double __ieee754_lgamma_r __P((double,int *));
162extern double __ieee754_gamma_r __P((double,int *));
163extern double __ieee754_log10 __P((double));
164extern double __ieee754_sinh __P((double));
165extern double __ieee754_hypot __P((double,double));
166extern double __ieee754_j0 __P((double));
167extern double __ieee754_j1 __P((double));
168extern double __ieee754_y0 __P((double));
169extern double __ieee754_y1 __P((double));
170extern double __ieee754_jn __P((int,double));
171extern double __ieee754_yn __P((int,double));
172extern double __ieee754_remainder __P((double,double));
173extern __int32_t __ieee754_rem_pio2 __P((double,double*));
174#ifdef _SCALB_INT
175extern double __ieee754_scalb __P((double,int));
176#else
177extern double __ieee754_scalb __P((double,double));
178#endif
179
180/* fdlibm kernel function */
181extern double __kernel_standard __P((double,double,int));
182extern double __kernel_sin __P((double,double,int));
183extern double __kernel_cos __P((double,double));
184extern double __kernel_tan __P((double,double,int));
185extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
186
187/* Undocumented float functions.  */
188#ifdef _SCALB_INT
189extern float scalbf __P((float, int));
190#else
191extern float scalbf __P((float, float));
192#endif
193extern float significandf __P((float));
194
195/* ieee style elementary float functions */
196extern float __ieee754_sqrtf __P((float));
197extern float __ieee754_acosf __P((float));
198extern float __ieee754_acoshf __P((float));
199extern float __ieee754_logf __P((float));
200extern float __ieee754_atanhf __P((float));
201extern float __ieee754_asinf __P((float));
202extern float __ieee754_atan2f __P((float,float));
203extern float __ieee754_expf __P((float));
204extern float __ieee754_coshf __P((float));
205extern float __ieee754_fmodf __P((float,float));
206extern float __ieee754_powf __P((float,float));
207extern float __ieee754_lgammaf_r __P((float,int *));
208extern float __ieee754_gammaf_r __P((float,int *));
209extern float __ieee754_log10f __P((float));
210extern float __ieee754_sinhf __P((float));
211extern float __ieee754_hypotf __P((float,float));
212extern float __ieee754_j0f __P((float));
213extern float __ieee754_j1f __P((float));
214extern float __ieee754_y0f __P((float));
215extern float __ieee754_y1f __P((float));
216extern float __ieee754_jnf __P((int,float));
217extern float __ieee754_ynf __P((int,float));
218extern float __ieee754_remainderf __P((float,float));
219extern __int32_t __ieee754_rem_pio2f __P((float,float*));
220#ifdef _SCALB_INT
221extern float __ieee754_scalbf __P((float,int));
222#else
223extern float __ieee754_scalbf __P((float,float));
224#endif
225
226/* float versions of fdlibm kernel functions */
227extern float __kernel_sinf __P((float,float,int));
228extern float __kernel_cosf __P((float,float));
229extern float __kernel_tanf __P((float,float,int));
230extern int   __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
231
232/* The original code used statements like
233	n0 = ((*(int*)&one)>>29)^1;		* index of high word *
234	ix0 = *(n0+(int*)&x);			* high word of x *
235	ix1 = *((1-n0)+(int*)&x);		* low word of x *
236   to dig two 32 bit words out of the 64 bit IEEE floating point
237   value.  That is non-ANSI, and, moreover, the gcc instruction
238   scheduler gets it wrong.  We instead use the following macros.
239   Unlike the original code, we determine the endianness at compile
240   time, not at run time; I don't see much benefit to selecting
241   endianness at run time.  */
242
243#ifndef __IEEE_BIG_ENDIAN
244#ifndef __IEEE_LITTLE_ENDIAN
245 #error Must define endianness
246#endif
247#endif
248
249/* A union which permits us to convert between a double and two 32 bit
250   ints.  */
251
252#ifdef __IEEE_BIG_ENDIAN
253
254typedef union
255{
256  double value;
257  struct
258  {
259    __uint32_t msw;
260    __uint32_t lsw;
261  } parts;
262} ieee_double_shape_type;
263
264#endif
265
266#ifdef __IEEE_LITTLE_ENDIAN
267
268typedef union
269{
270  double value;
271  struct
272  {
273    __uint32_t lsw;
274    __uint32_t msw;
275  } parts;
276} ieee_double_shape_type;
277
278#endif
279
280/* Get two 32 bit ints from a double.  */
281
282#define EXTRACT_WORDS(ix0,ix1,d)				\
283do {								\
284  ieee_double_shape_type ew_u;					\
285  ew_u.value = (d);						\
286  (ix0) = ew_u.parts.msw;					\
287  (ix1) = ew_u.parts.lsw;					\
288} while (0)
289
290/* Get the more significant 32 bit int from a double.  */
291
292#define GET_HIGH_WORD(i,d)					\
293do {								\
294  ieee_double_shape_type gh_u;					\
295  gh_u.value = (d);						\
296  (i) = gh_u.parts.msw;						\
297} while (0)
298
299/* Get the less significant 32 bit int from a double.  */
300
301#define GET_LOW_WORD(i,d)					\
302do {								\
303  ieee_double_shape_type gl_u;					\
304  gl_u.value = (d);						\
305  (i) = gl_u.parts.lsw;						\
306} while (0)
307
308/* Set a double from two 32 bit ints.  */
309
310#define INSERT_WORDS(d,ix0,ix1)					\
311do {								\
312  ieee_double_shape_type iw_u;					\
313  iw_u.parts.msw = (ix0);					\
314  iw_u.parts.lsw = (ix1);					\
315  (d) = iw_u.value;						\
316} while (0)
317
318/* Set the more significant 32 bits of a double from an int.  */
319
320#define SET_HIGH_WORD(d,v)					\
321do {								\
322  ieee_double_shape_type sh_u;					\
323  sh_u.value = (d);						\
324  sh_u.parts.msw = (v);						\
325  (d) = sh_u.value;						\
326} while (0)
327
328/* Set the less significant 32 bits of a double from an int.  */
329
330#define SET_LOW_WORD(d,v)					\
331do {								\
332  ieee_double_shape_type sl_u;					\
333  sl_u.value = (d);						\
334  sl_u.parts.lsw = (v);						\
335  (d) = sl_u.value;						\
336} while (0)
337
338/* A union which permits us to convert between a float and a 32 bit
339   int.  */
340
341typedef union
342{
343  float value;
344  __uint32_t word;
345} ieee_float_shape_type;
346
347/* Get a 32 bit int from a float.  */
348
349#define GET_FLOAT_WORD(i,d)					\
350do {								\
351  ieee_float_shape_type gf_u;					\
352  gf_u.value = (d);						\
353  (i) = gf_u.word;						\
354} while (0)
355
356/* Set a float from a 32 bit int.  */
357
358#define SET_FLOAT_WORD(d,i)					\
359do {								\
360  ieee_float_shape_type sf_u;					\
361  sf_u.word = (i);						\
362  (d) = sf_u.value;						\
363} while (0)
364
365/* Macros to avoid undefined behaviour that can arise if the amount
366   of a shift is exactly equal to the size of the shifted operand.  */
367
368#define SAFE_LEFT_SHIFT(op,amt)					\
369  (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0)
370
371#define SAFE_RIGHT_SHIFT(op,amt)				\
372  (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0)
373
374#ifdef  _COMPLEX_H
375
376/*
377 * Quoting from ISO/IEC 9899:TC2:
378 *
379 * 6.2.5.13 Types
380 * Each complex type has the same representation and alignment requirements as
381 * an array type containing exactly two elements of the corresponding real type;
382 * the first element is equal to the real part, and the second element to the
383 * imaginary part, of the complex number.
384 */
385typedef union {
386        float complex z;
387        float parts[2];
388} float_complex;
389
390typedef union {
391        double complex z;
392        double parts[2];
393} double_complex;
394
395typedef union {
396        long double complex z;
397        long double parts[2];
398} long_double_complex;
399
400#define REAL_PART(z)    ((z).parts[0])
401#define IMAG_PART(z)    ((z).parts[1])
402
403#endif  /* _COMPLEX_H */
404
405