/*****************************************************************************/
/* @file:    math.h                                                          */
/* @version: V1.01                                                           */
/*****************************************************************************/

#ifndef _MATH
#define _MATH

#ifdef __cplusplus
extern "C" {
#endif

#ifndef _ERRNO
#define _ERRNO

#define EDOM   1
#define ERANGE 2
#define ENOENT 3
#define EFPOS  5
#define EILSEQ 6

#define __set_errno(val) (errno = (val))

#endif  /* _ERRNO */


/*
 * Macros for our inline functions down below.
 * unsigned& __FLT(float x) - returns the bit pattern of x
 * unsigned& __HI(double x) - returns the bit pattern of the high part of x
 *                            (high part has exponent & sign bit in it)
 * unsigned& __LO(double x) - returns the bit pattern of the low part of x
 *
 * We can assign to __FLT, __HI, and __LO and the appropriate bits get set in
 * the floating point variable used.
 *
 * __HI & __LO are affected by the little endian.
 */
#define __FLT(x) (*(unsigned int *)&(x))
#define __HI(x)  (*(1 + (unsigned int *)&(x)))
#define __LO(x)  (*(unsigned int *)&(x))

#define	signbit(x) \
  (sizeof(x) == sizeof(float) ? (__FLT(x) >> 31) : (__HI(x) >> 31))

/***************************************************************/
/* CONSTANT.                                                   */
/***************************************************************/
#define HUGE_VALF \
  (((union { unsigned int __l ; float __d; })    \
    { __l: 0x7f800000UL }).__d)
#define HUGE_VAL \
  (((union { unsigned long long __l; double __d; })  \
    { __l: 0x7ff0000000000000ULL }).__d)
#define HUGE_VALL  HUGE_VAL
#define NANF \
  (((union { unsigned int __l; float __d; })  \
    { __l: 0x7fc00000UL }).__d)
#define NAN \
  (((union { unsigned long long __l; double __d; })  \
    { __l: 0x7ff8000000000000ULL }).__d)
#define INFINITY   HUGE_VALF

#define M_E			2.71828182845904523540e+00 /* 0x4005BF0A, 0x8B145769 */	/* e */
#define M_LOG2E		1.44269504088896340740e+00 /* 0x3FF71547, 0x652B82FE */	/* log_2 e */
#define M_LOG10E	0.43429448190325182765e+00 /* 0x3FDBCB7B, 0x1526E50E */	/* log_10 e */
#define M_LN2		0.69314718055994530942e+00 /* 0x3FE62E42, 0xFEFA39EF */	/* log_e 2 */
#define M_LN10		2.30258509299404568402e+00 /* 0x40026BB1, 0xBBB55516 */	/* log_e 10 */
#define M_PI		3.14159265358979323846e+00 /* 0x400921FB, 0x54442D18 */	/* pi */
#define M_PI_2		1.57079632679489655800e+00 /* 0x3FF921FB, 0x54442D18 */	/* pi/2 */
#define M_PI_4		0.78539816339744830962e+00 /* 0x3FE921FB, 0x54442D18 */	/* pi/4 */
#define M_1_PI		0.31830988618379067154e+00 /* 0x3FD45F30, 0x6DC9C883 */	/* 1/pi */
#define M_2_PI		0.63661977236758134308e+00 /* 0x3FE45F30, 0x6DC9C883 */	/* 2/pi */
#define M_2_SQRTPI	1.12837916709551257390e+00 /* 0x3FF20DD7, 0x50429B6D */	/* 2/sqrt(pi) */
#define M_SQRT2		1.41421356237309504880e+00 /* 0x3FF6A09E, 0x667F3BCD */	/* sqrt(2) */
#define M_SQRT1_2	0.70710678118654752440e+00 /* 0x3FE6A09E, 0x667F3BCD */	/* 1/sqrt(2) */

enum
{
	FP_INFINITE,
	FP_NAN,
	FP_NORMAL,
	FP_SUBNORMAL,
	FP_ZERO
};

/* simplistic isinf and isnan */
#define isinf(x)      ((x) == INFINITY || (x) == -INFINITY)
#define isnan(x)      ((x) != (x))
#define isfinite(x)   (sizeof (x) == sizeof (float) ? finitef (x) : finite (x))
#define isnormal(x)   (fpclassify(x) == FP_NORMAL)
#define fpclassify(x) (sizeof (x) == sizeof (float) ? __fpclassifyf (x) : __fpclassify (x))
//#define signbit(x)    (sizeof (x) == sizeof (float) ? __signbitf (x) : __signbit (x))

/* floating-point comparisons without exception for NaN */
#define isgreater(x, y)         __builtin_isgreater(x, y)
#define isgreaterequal(x, y)    __builtin_isgreaterequal(x, y)
#define isless(x, y)            __builtin_isless(x, y)
#define islessequal(x, y)       __builtin_islessequal(x, y)
#define islessgreater(x, y)     __builtin_islessgreater(x, y)
#define isunordered(x, y)       __builtin_isunordered(x, y)

#ifndef _IEEE_LIBM
#define _POSIX_ 3
#define _XOPEN_ 2
#define _SVID_ 1
#define _IEEE_ 0
typedef int _LIB_VERSION_TYPE;
extern _LIB_VERSION_TYPE _LIB_VERSION;
#endif
extern int signgam;
extern int errno;
/***************************************************************/
/* FUNCTION DEFINITIONS.                                       */
/***************************************************************/
/* trigonometric function */
extern  float       sinf            ( float );
extern  double      sin             ( double );
extern  long double sinl            ( long double );
extern  float       cosf            ( float );
extern  double      cos             ( double );
extern  long double cosl            ( long double );
extern  float       tanf            ( float );
extern  double      tan             ( double );
extern  long double tanl            ( long double );
extern  void        sincosf         ( float, float *, float * );

/* hyperbolic trigonometric function */
extern  float       sinhf           ( float );
extern  double      sinh            ( double );
extern  long double sinhl           ( long double );
extern  float       coshf           ( float );
extern  double      cosh            ( double );
extern  long double coshl           ( long double );
extern  float       tanhf           ( float );
extern  double      tanh            ( double );
extern  long double tanhl           ( long double );

/* arc trigonometric function */
extern  float       asinf           ( float );
extern  double      asin            ( double );
extern  long double asinl           ( long double );
extern  float       acosf           ( float );
extern  double      acos            ( double );
extern  long double acosl           ( long double );
extern  float       atanf           ( float );
extern  double      atan            ( double );
extern  long double atanl           ( long double );

/* arc tangent function of two variables */
extern  float       atan2f          ( float, float );
extern  double      atan2           ( double, double );
extern  long double atan2l          ( long double, long double );

/* inverse hyperbolic trigonometric function */
extern  float       asinhf          ( float );
extern  double      asinh           ( double );
extern  long double asinhl          ( long double );
extern  float       acoshf          ( float );
extern  double      acosh           ( double );
extern  long double acoshl          ( long double );
extern  float       atanhf          ( float );
extern  double      atanh           ( double );
extern  long double atanhl          ( long double );

/* base-e exponential function */
extern  float       expf            ( float );
extern  double      exp             ( double );
extern  long double expl            ( long double );
/* base-2 exponential function */
extern  float       exp2f           ( float );
extern  double      exp2            ( double );
extern  long double exp2l           ( long double );
/* exponential minus 1 */
extern  float       expm1f          ( float );
extern  double      expm1           ( double );
extern  long double expm1l          ( long double );

/* natural logarithmic function */
extern  float       logf            ( float );
extern  double      log             ( double );
extern  long double logl            ( long double );
/* base-2 logarithmic function */
extern  float       log2f           ( float );
extern  double      log2            ( double );
extern  long double log2l           ( long double );
/* base-10 logarithmic function */
extern  float       log10f          ( float );
extern  double      log10           ( double );
extern  long double log10l          ( long double );
/* get exponent of a floating-point value */
extern  float       logbf           ( float );
extern  double      logb            ( double );
extern  long double logbl           ( long double );
/* logarithm of 1 plus argument: log (1 + x) */
extern  float       log1pf          ( float );
extern  double      log1p           ( double );
extern  long double log1pl          ( long double );

/* power functions */
extern  float       powf            ( float, float );
extern  double      pow             ( double, double );
extern  long double powl            ( long double, long double );

/* square root function */
extern  float       sqrtf           ( float );
extern  double      sqrt            ( double );
extern  long double sqrtl           ( long double );

/* cube root function */
extern  float       cbrtf           ( float );
extern  double      cbrt            ( double );
extern  long double cbrtl           ( long double );

/*  ceiling function: smallest integral value not less than argument */
extern  float       ceilf           ( float );
extern  double      ceil            ( double );
extern  long double ceill           ( long double );

/* largest integral value not greater than argument */
extern  float       floorf          ( float );
extern  double      floor           ( double );
extern  long double floorl          ( long double );

/* absolute value of floating-point number */
extern  float       fabsf           ( float );
extern  double      fabs            ( double );
extern  long double fabsl           ( long double );

/* multiply floating-point number by integral power of 2 */
extern  float       ldexpf          ( float, int );
extern  double      ldexp           ( double, int );
extern  long double ldexpl          ( long double, int );

/* convert floating-point number to fractional and integral components */
extern  float       frexpf          ( float, int * );
extern  double      frexp           ( double, int * );
extern  long double frexpl          ( long double, int * );

/* extract signed integral and fractional values from floating-point number */
extern  float       modff           ( float, float * );
extern  double      modf            ( double, double * );
extern  long double modfl           ( long double, long double * );

/* floating-point remainder function */
extern  float       fmodf           ( float, float );
extern  double      fmod            ( double, double );
extern  long double fmodl           ( long double, long double );

/* Euclidean distance function */
extern  float       hypotf          ( float, float );
extern  double      hypot           ( double, double );
extern  long double hypotl          ( long double, long double );

/* Bessel functions of the first kind */
extern  float       j0f             ( float );
extern  double      j0              ( double );
extern  long double j0l             ( long double );
extern  float       j1f             ( float );
extern  double      j1              ( double );
extern  long double j1l             ( long double );
extern  float       jnf             ( int, float );
extern  double      jn              ( int, double );
extern  long double jnl             ( int, long double );

/* Bessel functions of the second kind */
extern  float       y0f             ( float );
extern  double      y0              ( double );
extern  long double y0l             ( long double );
extern  float       y1f             ( float );
extern  double      y1              ( double );
extern  long double y1l             ( long double );
extern  float       ynf             ( int, float );
extern  double      yn              ( int, double );
extern  long double ynl             ( int, long double );

/* copy sign of a number */
extern  float       copysignf       ( float, float );
extern  double      copysign        ( double, double );
extern  long double copysignl       ( long double, long double );

/* complementary error function */
extern  float       erfcf           ( float );
extern  double      erfc            ( double );
extern  long double erfcl           ( long double );

/* error function */
extern  float       erff            ( float );
extern  double      erf             ( double );
extern  long double erfl            ( long double );

/* positive difference */
extern  float       fdimf           ( float, float );
extern  double      fdim            ( double, double );
extern  long double fdiml           ( long double, long double );

/* floating-point multiply and add */
extern  float       fmaf            ( float, float, float );
extern  double      fma             ( double, double, double );
extern  long double fmal            ( long double, long double, long double );

/* determine maximum of two floating-point numbers */
extern  float       fmaxf           ( float, float );
extern  double      fmax            ( double, double );
extern  long double fmaxl           ( long double, long double );

/* determine minimum of two floating-point numbers */
extern  float       fminf           ( float, float );
extern  double      fmin            ( double, double );
extern  long double fminl           ( long double, long double );

/* get integer exponent of a floating-point value */
#define FP_ILOGB0                   (-2147483647)
#define FP_ILOGBNAN                 (2147483647)
extern  int         ilogbf          ( float );
extern  int         ilogb           ( double );
extern  int         ilogbl          ( long double );

/* log gamma function */
extern  float       lgammaf         ( float );
extern  double      lgamma          ( double );
extern  long double lgammal         ( long double );
extern  float       lgammaf_r       ( float, int * );
extern  double      lgamma_r        ( double, int * );
extern  long double lgammal_r       ( long double, int * );

/* round to nearest integer */
extern  float       rintf           ( float );
extern  double      rint            ( double );
extern  long double rintl           ( long double );
extern  long        lrintf          ( float );
extern  long        lrint           ( double );
extern  long        lrintl          ( long double );
extern  long long   llrintf         ( float );
extern  long long   llrint          ( double );
extern  long long   llrintl         ( long double );

/* round to nearest integer, away from zero */
extern  float       roundf          ( float );
extern  double      round           ( double );
extern  long double roundl          ( long double );
extern  long        lroundf         ( float );
extern  long        lround          ( double );
extern  long        lroundl         ( long double );
extern  long long   llroundf        ( float );
extern  long long   llround         ( double );
extern  long long   llroundl        ( long double );

/* round to nearest integer */
extern  float       nearbyintf      ( float );
extern  double      nearbyint       ( double );
extern  long double nearbyintl      ( long double );

/* floating-point number manipulation */
extern  float       nextafterf      ( float, float );
extern  double      nextafter       ( double, double );
extern  long double nextafterl      ( long double, long double );
extern  float       nexttowardf     ( float, double );
extern  double      nexttoward      ( double, double );
extern  long double nexttowardl     ( long double, long double );
/* Return the least floating-point number greater than ARG */
extern  float       nextupf         ( float );

/* floating-point remainder function */
extern  float       remainderf      ( float, float );
extern  double      remainder       ( double, double );
extern  long double remainderl      ( long double, long double );
extern  double      drem            ( double, double );
extern  long double dreml           ( long double, long double );

/* remainder and part of quotient */
extern  float       remquof         ( float, float, int * );
extern  double      remquo          ( double, double, int * );
extern  long double remquol         ( long double, long double, int * );

/* multiply floating-point number by integral power of radix */
extern  float       scalbnf         ( float, int );
extern  double      scalbn          ( double, long );
extern  long double scalbnl         ( long double, int );
extern  float       scalblnf        ( float, int );
extern  double      scalbln         ( double, long );
extern  long double scalblnl        ( long double, long );

/* true gamma function */
extern  float       tgammaf         ( float );
extern  double      tgamma          ( double );
extern  long double tgammal         ( long double );

/* round to integer, toward zero */
extern  float       truncf          ( float );
extern  double      trunc           ( double );
extern  long double truncl          ( long double );

/* get mantissa of floating-point number */
extern  double      significand     ( double );
extern  long double significandl    ( long double );

/* BSD floating-point classification functions */
extern  int         finitef         ( float );
extern  int         finite          ( double );
extern  int         finitel         ( double );
extern  int         __isinff        ( float );
extern  int         __isinf         ( double );
extern  int         __isinfl        ( long double );
extern  int         __isnanf        ( float );
extern  int         __isnan         ( double );
extern  int         __isnanl        ( long double );

/* test sign of a real floating-point number */
extern  int         __signbitf      ( float );
extern  int         __signbit       ( double );
extern  int         __signbitl      ( long double );

extern  int         __fpclassifyf   ( float );
extern  int         __fpclassify    ( double );
extern  int         __fpclassifyl   ( long double );

/* logarithm of the gamma function */
extern  double      gamma           ( double );
extern  long double gammal          ( long double );
extern  float       gammaf_r        ( float, int * );
extern  double      gamma_r         ( double, int * );
extern  long double gammal_r        ( long double, int * );

/* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N
   - 1, in the form R * (1 + *EPS) where the return value R is an
   approximation to the product and *EPS is set to indicate the
   approximate error in the return value.  X is such that all the
   values X + 1, ..., X + N - 1 are exactly representable, and X_EPS /
   X is small enough that factors quadratic in it can be
   neglected.  */
extern  float       __gamma_productf (float x, float x_eps, int n, float *eps);

/* lgammaf expanding around zeros.
   Compute lgamma of a negative argument -15 < X < -2, setting
   *SIGNGAMP accordingly.  */
extern  float       __lgamma_negf   ( float x, int *signgamp );

#ifdef __cplusplus
}
#endif
#endif /* _MATH */
