Subversion Repositories Kolibri OS

/*

 * ====================================================

 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.

 *

 * Developed at SunPro, a Sun Microsystems, Inc. business.

 * Permission to use, copy, modify, and distribute this

 * software is freely granted, provided that this notice

 * is preserved.

 * ====================================================

 */

#include 

#include 

#include 

#define _XOPEN_MODE

   number, and many need to know whether the result of an operation will

   overflow.  These conditions depend on whether the largest exponent is

   used for NaNs & infinities, or whether it's used for finite numbers.  The

   macros below wrap up that kind of information:

	True if a positive float with bitmask X is finite.

	True if a positive float with bitmask X is not a number.

	True if a positive float with bitmask X is +infinity.

	The bitmask of FLT_MAX.

	The bitmask of FLT_MAX/2.

	The bitmask of the largest finite exponent (129 if the largest

	exponent is used for finite numbers, 128 otherwise).

	The bitmask of log(FLT_MAX), rounded down.  This value is the largest

	input that can be passed to exp() without producing overflow.

	The bitmask of log(2*FLT_MAX), rounded down.  This value is the

	largest input than can be passed to cosh() without producing

	overflow.

	The largest biased exponent that can be used for finite numbers

	(255 if the largest exponent is used for finite numbers, 254

	otherwise) */

#define FLT_UWORD_IS_FINITE(x) 1

#define FLT_UWORD_IS_NAN(x) 0

#define FLT_UWORD_IS_INFINITE(x) 0

#define FLT_UWORD_MAX 0x7fffffff

#define FLT_UWORD_EXP_MAX 0x43010000

#define FLT_UWORD_LOG_MAX 0x42b2d4fc

#define FLT_UWORD_LOG_2MAX 0x42b437e0

#define HUGE ((float)0X1.FFFFFEP128)

#else

#define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)

#define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)

#define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)

#define FLT_UWORD_MAX 0x7f7fffffL

#define FLT_UWORD_EXP_MAX 0x43000000

#define FLT_UWORD_LOG_MAX 0x42b17217

#define FLT_UWORD_LOG_2MAX 0x42b2d4fc

#define HUGE ((float)3.40282346638528860e+38)

#endif

#define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23))

#define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)

   on whether the target supports denormals or not:

	True if a positive float with bitmask X is +0.	Without denormals,

	any float with a zero exponent is a +0 representation.	With

	denormals, the only +0 representation is a 0 bitmask.

	True if a non-zero positive float with bitmask X is subnormal.

	(Routines should check for zeros first.)

	The bitmask of the smallest float above +0.  Call this number

	REAL_FLT_MIN...

	The bitmask of the float representation of REAL_FLT_MIN's exponent.

	The bitmask of |log(REAL_FLT_MIN)|, rounding down.

	REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,

	-22 if they are).

*/

#define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)

#define FLT_UWORD_IS_SUBNORMAL(x) 0

#define FLT_UWORD_MIN 0x00800000

#define FLT_UWORD_EXP_MIN 0x42fc0000

#define FLT_UWORD_LOG_MIN 0x42aeac50

#define FLT_SMALLEST_EXP 1

#else

#define FLT_UWORD_IS_ZERO(x) ((x)==0)

#define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)

#define FLT_UWORD_MIN 0x00000001

#define FLT_UWORD_EXP_MIN 0x43160000

#define FLT_UWORD_LOG_MIN 0x42cff1b5

#define FLT_SMALLEST_EXP -22

#endif

#undef __P

#define	__P(p)	p

#else

#define	__P(p)	()

#endif

 * set X_TLOSS = pi*2**52, which is possibly defined in 

 * (one may replace the following line by "#include ")

 */

extern double scalb __P((double, int));

#else

extern double scalb __P((double, double));

#endif

extern double significand __P((double));

extern double __ieee754_sqrt __P((double));

extern double __ieee754_acos __P((double));

extern double __ieee754_acosh __P((double));

extern double __ieee754_log __P((double));

extern double __ieee754_atanh __P((double));

extern double __ieee754_asin __P((double));

extern double __ieee754_atan2 __P((double,double));

extern double __ieee754_exp __P((double));

extern double __ieee754_cosh __P((double));

extern double __ieee754_fmod __P((double,double));

extern double __ieee754_pow __P((double,double));

extern double __ieee754_lgamma_r __P((double,int *));

extern double __ieee754_gamma_r __P((double,int *));

extern double __ieee754_log10 __P((double));

extern double __ieee754_sinh __P((double));

extern double __ieee754_hypot __P((double,double));

extern double __ieee754_j0 __P((double));

extern double __ieee754_j1 __P((double));

extern double __ieee754_y0 __P((double));

extern double __ieee754_y1 __P((double));

extern double __ieee754_jn __P((int,double));

extern double __ieee754_yn __P((int,double));

extern double __ieee754_remainder __P((double,double));

extern __int32_t __ieee754_rem_pio2 __P((double,double*));

#ifdef _SCALB_INT

extern double __ieee754_scalb __P((double,int));

#else

extern double __ieee754_scalb __P((double,double));

#endif

extern double __kernel_standard __P((double,double,int));

extern double __kernel_sin __P((double,double,int));

extern double __kernel_cos __P((double,double));

extern double __kernel_tan __P((double,double,int));

extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));

#ifdef _SCALB_INT

extern float scalbf __P((float, int));

#else

extern float scalbf __P((float, float));

#endif

extern float significandf __P((float));

extern float __ieee754_sqrtf __P((float));

extern float __ieee754_acosf __P((float));

extern float __ieee754_acoshf __P((float));

extern float __ieee754_logf __P((float));

extern float __ieee754_atanhf __P((float));

extern float __ieee754_asinf __P((float));

extern float __ieee754_atan2f __P((float,float));

extern float __ieee754_expf __P((float));

extern float __ieee754_coshf __P((float));

extern float __ieee754_fmodf __P((float,float));

extern float __ieee754_powf __P((float,float));

extern float __ieee754_lgammaf_r __P((float,int *));

extern float __ieee754_gammaf_r __P((float,int *));

extern float __ieee754_log10f __P((float));

extern float __ieee754_sinhf __P((float));

extern float __ieee754_hypotf __P((float,float));

extern float __ieee754_j0f __P((float));

extern float __ieee754_j1f __P((float));

extern float __ieee754_y0f __P((float));

extern float __ieee754_y1f __P((float));

extern float __ieee754_jnf __P((int,float));

extern float __ieee754_ynf __P((int,float));

extern float __ieee754_remainderf __P((float,float));

extern __int32_t __ieee754_rem_pio2f __P((float,float*));

#ifdef _SCALB_INT

extern float __ieee754_scalbf __P((float,int));

#else

extern float __ieee754_scalbf __P((float,float));

#endif

extern float __kernel_sinf __P((float,float,int));

extern float __kernel_cosf __P((float,float));

extern float __kernel_tanf __P((float,float,int));

extern int   __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));

	n0 = ((*(int*)&one)>>29)^1;		* index of high word *

	ix0 = *(n0+(int*)&x);			* high word of x *

	ix1 = *((1-n0)+(int*)&x);		* low word of x *

   to dig two 32 bit words out of the 64 bit IEEE floating point

   value.  That is non-ANSI, and, moreover, the gcc instruction

   scheduler gets it wrong.  We instead use the following macros.

   Unlike the original code, we determine the endianness at compile

   time, not at run time; I don't see much benefit to selecting

   endianness at run time.  */

#ifndef __IEEE_LITTLE_ENDIAN

 #error Must define endianness

#endif

#endif

   ints.  */

{

  double value;

  struct

  {

    __uint32_t msw;

    __uint32_t lsw;

  } parts;

} ieee_double_shape_type;

{

  double value;

  struct

  {

    __uint32_t lsw;

    __uint32_t msw;

  } parts;

} ieee_double_shape_type;

do {								\

  ieee_double_shape_type ew_u;					\

  ew_u.value = (d);						\

  (ix0) = ew_u.parts.msw;					\

  (ix1) = ew_u.parts.lsw;					\

} while (0)

do {								\

  ieee_double_shape_type gh_u;					\

  gh_u.value = (d);						\

  (i) = gh_u.parts.msw;						\

} while (0)

do {								\

  ieee_double_shape_type gl_u;					\

  gl_u.value = (d);						\

  (i) = gl_u.parts.lsw;						\

} while (0)

do {								\

  ieee_double_shape_type iw_u;					\

  iw_u.parts.msw = (ix0);					\

  iw_u.parts.lsw = (ix1);					\

  (d) = iw_u.value;						\

} while (0)

do {								\

  ieee_double_shape_type sh_u;					\

  sh_u.value = (d);						\

  sh_u.parts.msw = (v);						\

  (d) = sh_u.value;						\

} while (0)

do {								\

  ieee_double_shape_type sl_u;					\

  sl_u.value = (d);						\

  sl_u.parts.lsw = (v);						\

  (d) = sl_u.value;						\

} while (0)

   int.  */

{

  float value;

  __uint32_t word;

} ieee_float_shape_type;

do {								\

  ieee_float_shape_type gf_u;					\

  gf_u.value = (d);						\

  (i) = gf_u.word;						\

} while (0)

do {								\

  ieee_float_shape_type sf_u;					\

  sf_u.word = (i);						\

  (d) = sf_u.value;						\

} while (0)

   of a shift is exactly equal to the size of the shifted operand.  */

  (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0)

  (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0)

 * Quoting from ISO/IEC 9899:TC2:

 *

 * 6.2.5.13 Types

 * Each complex type has the same representation and alignment requirements as

 * an array type containing exactly two elements of the corresponding real type;

 * the first element is equal to the real part, and the second element to the

 * imaginary part, of the complex number.

 */

typedef union {

        float complex z;

        float parts[2];

} float_complex;

        double complex z;

        double parts[2];

} double_complex;

        long double complex z;

        long double parts[2];

} long_double_complex;

#define IMAG_PART(z)    ((z).parts[1])

#define>