diff --git a/include/tgmath.h b/include/tgmath.h new file mode 100644 index 00000000..5d3e3578 --- /dev/null +++ b/include/tgmath.h @@ -0,0 +1,89 @@ +/* + * ISO C Standard: 7.22 Type-generic math + */ + +#ifndef _TGMATH_H +#define _TGMATH_H + +#include + +#ifndef __cplusplus +#define __tgmath_real(x, F) \ + _Generic ((x), float: F##f, long double: F##l, default: F)(x) +#define __tgmath_real_2_1(x, y, F) \ + _Generic ((x), float: F##f, long double: F##l, default: F)(x, y) +#define __tgmath_real_2(x, y, F) \ + _Generic ((x)+(y), float: F##f, long double: F##l, default: F)(x, y) +#define __tgmath_real_3_2(x, y, z, F) \ + _Generic ((x)+(y), float: F##f, long double: F##l, default: F)(x, y, z) +#define __tgmath_real_3(x, y, z, F) \ + _Generic ((x)+(y)+(z), float: F##f, long double: F##l, default: F)(x, y, z) + +/* Functions defined in both and (7.22p4) */ +#define acos(z) __tgmath_real(z, acos) +#define asin(z) __tgmath_real(z, asin) +#define atan(z) __tgmath_real(z, atan) +#define acosh(z) __tgmath_real(z, acosh) +#define asinh(z) __tgmath_real(z, asinh) +#define atanh(z) __tgmath_real(z, atanh) +#define cos(z) __tgmath_real(z, cos) +#define sin(z) __tgmath_real(z, sin) +#define tan(z) __tgmath_real(z, tan) +#define cosh(z) __tgmath_real(z, cosh) +#define sinh(z) __tgmath_real(z, sinh) +#define tanh(z) __tgmath_real(z, tanh) +#define exp(z) __tgmath_real(z, exp) +#define log(z) __tgmath_real(z, log) +#define pow(z1,z2) __tgmath_real_2(z1, z2, pow) +#define sqrt(z) __tgmath_real(z, sqrt) +#define fabs(z) __tgmath_real(z, fabs) + +/* Functions defined in only (7.22p5) */ +#define atan2(x,y) __tgmath_real_2(x, y, atan2) +#define cbrt(x) __tgmath_real(x, cbrt) +#define ceil(x) __tgmath_real(x, ceil) +#define copysign(x,y) __tgmath_real_2(x, y, copysign) +#define erf(x) __tgmath_real(x, erf) +#define erfc(x) __tgmath_real(x, erfc) +#define exp2(x) __tgmath_real(x, exp2) +#define expm1(x) __tgmath_real(x, expm1) +#define fdim(x,y) __tgmath_real_2(x, y, fdim) +#define floor(x) __tgmath_real(x, floor) +#define fma(x,y,z) __tgmath_real_3(x, y, z, fma) +#define fmax(x,y) __tgmath_real_2(x, y, fmax) +#define fmin(x,y) __tgmath_real_2(x, y, fmin) +#define fmod(x,y) __tgmath_real_2(x, y, fmod) +#define frexp(x,y) __tgmath_real_2_1(x, y, frexp) +#define hypot(x,y) __tgmath_real_2(x, y, hypot) +#define ilogb(x) __tgmath_real(x, ilogb) +#define ldexp(x,y) __tgmath_real_2_1(x, y, ldexp) +#define lgamma(x) __tgmath_real(x, lgamma) +#define llrint(x) __tgmath_real(x, llrint) +#define llround(x) __tgmath_real(x, llround) +#define log10(x) __tgmath_real(x, log10) +#define log1p(x) __tgmath_real(x, log1p) +#define log2(x) __tgmath_real(x, log2) +#define logb(x) __tgmath_real(x, logb) +#define lrint(x) __tgmath_real(x, lrint) +#define lround(x) __tgmath_real(x, lround) +#define nearbyint(x) __tgmath_real(x, nearbyint) +#define nextafter(x,y) __tgmath_real_2(x, y, nextafter) +#define nexttoward(x,y) __tgmath_real_2(x, y, nexttoward) +#define remainder(x,y) __tgmath_real_2(x, y, remainder) +#define remquo(x,y,z) __tgmath_real_3_2(x, y, z, remquo) +#define rint(x) __tgmath_real(x, rint) +#define round(x) __tgmath_real(x, round) +#define scalbln(x,y) __tgmath_real_2_1(x, y, scalbln) +#define scalbn(x,y) __tgmath_real_2_1(x, y, scalbn) +#define tgamma(x) __tgmath_real(x, tgamma) +#define trunc(x) __tgmath_real(x, trunc) + +/* Functions defined in only (7.22p6) +#define carg(z) __tgmath_cplx_only(z, carg) +#define cimag(z) __tgmath_cplx_only(z, cimag) +#define conj(z) __tgmath_cplx_only(z, conj) +#define cproj(z) __tgmath_cplx_only(z, cproj) +#define creal(z) __tgmath_cplx_only(z, creal) +*/ +#endif /* __cplusplus */ +#endif /* _TGMATH_H */ diff --git a/win32/include/math.h b/win32/include/math.h index c270c57c..745ed376 100644 --- a/win32/include/math.h +++ b/win32/include/math.h @@ -197,105 +197,6 @@ extern "C" { int __cdecl _fpclassf(float _X); #endif -#ifndef __cplusplus -#define _hypotl(x,y) ((long double)_hypot((double)(x),(double)(y))) -#define _matherrl _matherr - __CRT_INLINE long double _chgsignl(long double _Number) { return _chgsign((double)(_Number)); } - __CRT_INLINE long double _copysignl(long double _Number,long double _Sign) { return _copysign((double)(_Number),(double)(_Sign)); } - __CRT_INLINE float frexpf(float _X,int *_Y) { return ((float)frexp((double)_X,_Y)); } - -#if !defined (__ia64__) - __CRT_INLINE float __cdecl fabsf (float x) - { -#ifdef _WIN64 - *((int *) &x) &= 0x7fffffff; - return x; -#else - float res; - __asm__ ("fabs;" : "=t" (res) : "0" (x)); - return res; -#endif - } - - __CRT_INLINE float __cdecl ldexpf (float x, int expn) { return (float) ldexp (x, expn); } -#endif -#if defined (_WIN32) && !defined(_WIN64) - __CRT_INLINE float acosf(float x) { return (float) acos(x); } - __CRT_INLINE float asinf(float x) { return (float) asin(x); } - __CRT_INLINE float atanf(float x) { return (float) atan(x); } - __CRT_INLINE float atan2f(float x, float y) { return (float) atan2(x, y); } - __CRT_INLINE float ceilf(float x) { return (float) ceil(x); } - __CRT_INLINE float cosf(float x) { return (float) cos(x); } - __CRT_INLINE float coshf(float x) { return (float) cosh(x); } - __CRT_INLINE float expf(float x) { return (float) exp(x); } - __CRT_INLINE float floorf(float x) { return (float) floor(x); } - __CRT_INLINE float fmodf(float x, float y) { return (float) fmod(x, y); } - __CRT_INLINE float logf(float x) { return (float) log(x); } - __CRT_INLINE float logbf(float x) { return (float) logb(x); } - __CRT_INLINE float log10f(float x) { return (float) log10(x); } - __CRT_INLINE float modff(float x, float *y) { - double di, df = modf(x, &di); - *y = (float) di; return (float) df; - } - __CRT_INLINE float powf(float x, float y) { return (float) pow(x, y); } - __CRT_INLINE float sinf(float x) { return (float) sin(x); } - __CRT_INLINE float sinhf(float x) { return (float) sinh(x); } - __CRT_INLINE float sqrtf(float x) { return (float) sqrt(x); } - __CRT_INLINE float tanf(float x) { return (float) tan(x); } - __CRT_INLINE float tanhf(float x) { return (float) tanh(x); } -#endif -#else - // cplusplus - __CRT_INLINE long double __cdecl fabsl (long double x) - { - long double res; - __asm__ ("fabs;" : "=t" (res) : "0" (x)); - return res; - } - __CRT_INLINE long double modfl(long double _X,long double *_Y) { - double _Di,_Df = modf((double)_X,&_Di); - *_Y = (long double)_Di; - return (_Df); - } - __CRT_INLINE long double _chgsignl(long double _Number) { return _chgsign(static_cast(_Number)); } - __CRT_INLINE long double _copysignl(long double _Number,long double _Sign) { return _copysign(static_cast(_Number),static_cast(_Sign)); } - __CRT_INLINE float frexpf(float _X,int *_Y) { return ((float)frexp((double)_X,_Y)); } -#ifndef __ia64__ - __CRT_INLINE float __cdecl fabsf (float x) - { - float res; - __asm__ ("fabs;" : "=t" (res) : "0" (x)); - return res; - } - __CRT_INLINE float __cdecl ldexpf (float x, int expn) { return (float) ldexp (x, expn); } -#ifndef __x86_64 - __CRT_INLINE float acosf(float _X) { return ((float)acos((double)_X)); } - __CRT_INLINE float asinf(float _X) { return ((float)asin((double)_X)); } - __CRT_INLINE float atanf(float _X) { return ((float)atan((double)_X)); } - __CRT_INLINE float atan2f(float _X,float _Y) { return ((float)atan2((double)_X,(double)_Y)); } - __CRT_INLINE float ceilf(float _X) { return ((float)ceil((double)_X)); } - __CRT_INLINE float cosf(float _X) { return ((float)cos((double)_X)); } - __CRT_INLINE float coshf(float _X) { return ((float)cosh((double)_X)); } - __CRT_INLINE float expf(float _X) { return ((float)exp((double)_X)); } - __CRT_INLINE float floorf(float _X) { return ((float)floor((double)_X)); } - __CRT_INLINE float fmodf(float _X,float _Y) { return ((float)fmod((double)_X,(double)_Y)); } - __CRT_INLINE float logf(float _X) { return ((float)log((double)_X)); } - __CRT_INLINE float log10f(float _X) { return ((float)log10((double)_X)); } - __CRT_INLINE float modff(float _X,float *_Y) { - double _Di,_Df = modf((double)_X,&_Di); - *_Y = (float)_Di; - return ((float)_Df); - } - __CRT_INLINE float powf(float _X,float _Y) { return ((float)pow((double)_X,(double)_Y)); } - __CRT_INLINE float sinf(float _X) { return ((float)sin((double)_X)); } - __CRT_INLINE float sinhf(float _X) { return ((float)sinh((double)_X)); } - __CRT_INLINE float sqrtf(float _X) { return ((float)sqrt((double)_X)); } - __CRT_INLINE float tanf(float _X) { return ((float)tan((double)_X)); } - __CRT_INLINE float tanhf(float _X) { return ((float)tanh((double)_X)); } -#endif -#endif -#endif - #ifndef NO_OLDNAMES #define matherr _matherr @@ -339,10 +240,13 @@ extern "C" { extern int __cdecl __fpclassify (double); extern int __cdecl __fpclassifyl (long double); -/* Implemented at tcc/tcc_libm.h */ +/* Implemented at tcc/tcc_libm.h #define fpclassify(x) (sizeof (x) == sizeof (float) ? __fpclassifyf (x) \ : sizeof (x) == sizeof (double) ? __fpclassify (x) \ : __fpclassifyl (x)) +*/ +#define fpclassify(x) \ + _Generic(x, float: __fpclassifyf, double: __fpclassify, long double: __fpclassifyl)(x) /* 7.12.3.2 */ #define isfinite(x) ((fpclassify(x) & FP_NAN) == 0) @@ -364,10 +268,13 @@ extern "C" { extern int __cdecl __signbit (double); extern int __cdecl __signbitl (long double); -/* Implemented at tcc/tcc_libm.h */ +/* Implemented at tcc/tcc_libm.h #define signbit(x) (sizeof (x) == sizeof (float) ? __signbitf (x) \ : sizeof (x) == sizeof (double) ? __signbit (x) \ : __signbitl (x)) +*/ +#define signbit(x) \ + _Generic(x, float: __signbitf, double: __signbit, long double: signbitl)(x) extern double __cdecl exp2(double); extern float __cdecl exp2f(float); @@ -390,31 +297,6 @@ extern "C" { extern double __cdecl logb (double); extern float __cdecl logbf (float); extern long double __cdecl logbl (long double); -#ifndef _WIN32 - __CRT_INLINE double __cdecl logb (double x) - { - double res; - __asm__ ("fxtract\n\t" - "fstp %%st" : "=t" (res) : "0" (x)); - return res; - } - - __CRT_INLINE float __cdecl logbf (float x) - { - float res; - __asm__ ("fxtract\n\t" - "fstp %%st" : "=t" (res) : "0" (x)); - return res; - } - - __CRT_INLINE long double __cdecl logbl (long double x) - { - long double res; - __asm__ ("fxtract\n\t" - "fstp %%st" : "=t" (res) : "0" (x)); - return res; - } -#endif extern long double __cdecl modfl (long double, long double*); @@ -433,8 +315,8 @@ extern "C" { extern float __cdecl cbrtf (float); extern long double __cdecl cbrtl (long double); - __CRT_INLINE float __cdecl hypotf (float x, float y) - { return (float) hypot (x, y);} + extern double __cdecl hypot (double, double); + extern float __cdecl hypotf (float, float); extern long double __cdecl hypotl (long double, long double); extern long double __cdecl powl (long double, long double); @@ -490,117 +372,23 @@ extern "C" { /* 7.12.9.4 */ /* round, using fpu control word settings */ - __CRT_INLINE double __cdecl rint (double x) - { - double retval; - __asm__ ( - "fldl %1\n" - "frndint \n" - "fstpl %0\n" : "=m" (retval) : "m" (x)); - return retval; - } + extern double __cdecl rint (double); + extern float __cdecl rintf (float); + extern long double __cdecl rintl (long double); - __CRT_INLINE float __cdecl rintf (float x) - { - float retval; - __asm__ ( - "flds %1\n" - "frndint \n" - "fstps %0\n" : "=m" (retval) : "m" (x)); - return retval; - } + extern long __cdecl lrint (double); + extern long __cdecl lrintf (float); + extern long __cdecl lrintl (long double); - __CRT_INLINE long double __cdecl rintl (long double x) - { -#ifdef _WIN32 - // on win32 'long double' is double internally - return rint(x); -#else - long double retval; - __asm__ ( - "fldt %1\n" - "frndint \n" - "fstpt %0\n" : "=m" (retval) : "m" (x)); - return retval; -#endif - } - - /* 7.12.9.5 */ - __CRT_INLINE long __cdecl lrint (double x) - { - long retval; - __asm__ __volatile__ \ - ("fldl %1\n" \ - "fistpl %0" : "=m" (retval) : "m" (x)); \ - return retval; - } - - __CRT_INLINE long __cdecl lrintf (float x) - { - long retval; - __asm__ __volatile__ \ - ("flds %1\n" \ - "fistpl %0" : "=m" (retval) : "m" (x)); \ - return retval; - } - - __CRT_INLINE long __cdecl lrintl (long double x) - { - long retval; - __asm__ __volatile__ \ - ("fldt %1\n" \ - "fistpl %0" : "=m" (retval) : "m" (x)); \ - return retval; - } - - __CRT_INLINE long long __cdecl llrint (double x) - { - long long retval; - __asm__ __volatile__ \ - ("fldl %1\n" \ - "fistpll %0" : "=m" (retval) : "m" (x)); \ - return retval; - } - - __CRT_INLINE long long __cdecl llrintf (float x) - { - long long retval; - __asm__ __volatile__ \ - ("flds %1\n" \ - "fistpll %0" : "=m" (retval) : "m" (x)); \ - return retval; - } - - __CRT_INLINE long long __cdecl llrintl (long double x) - { - long long retval; - __asm__ __volatile__ \ - ("fldt %1\n" \ - "fistpll %0" : "=m" (retval) : "m" (x)); \ - return retval; - } + extern long long __cdecl llrint (double); + extern long long __cdecl llrintf (float); + extern long long __cdecl llrintl (long double); #define FE_TONEAREST 0x0000 #define FE_DOWNWARD 0x0400 #define FE_UPWARD 0x0800 #define FE_TOWARDZERO 0x0c00 - __CRT_INLINE double trunc (double _x) - { - double retval; - unsigned short saved_cw; - unsigned short tmp_cw; - __asm__ ("fnstcw %0;" : "=m" (saved_cw)); /* save FPU control word */ - tmp_cw = (saved_cw & ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO)) - | FE_TOWARDZERO; - __asm__ ("fldcw %0;" : : "m" (tmp_cw)); - __asm__ ("fldl %1;" - "frndint;" - "fstpl %0;" : "=m" (retval) : "m" (_x)); /* round towards zero */ - __asm__ ("fldcw %0;" : : "m" (saved_cw) ); /* restore saved control word */ - return retval; - } - /* 7.12.9.6 */ /* round away from zero, regardless of fpu control word settings */ extern double __cdecl round (double); @@ -685,70 +473,12 @@ extern "C" { extern long double __cdecl fmal (long double, long double, long double); -#if 0 // gr: duplicate, see below - /* 7.12.14 */ - /* - * With these functions, comparisons involving quiet NaNs set the FP - * condition code to "unordered". The IEEE floating-point spec - * dictates that the result of floating-point comparisons should be - * false whenever a NaN is involved, with the exception of the != op, - * which always returns true: yes, (NaN != NaN) is true). - */ - -#if __GNUC__ >= 3 - -#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) - -#else - /* helper */ - __CRT_INLINE int __cdecl - __fp_unordered_compare (long double x, long double y){ - unsigned short retval; - __asm__ ("fucom %%st(1);" - "fnstsw;": "=a" (retval) : "t" (x), "u" (y)); - return retval; - } - -#define isgreater(x, y) ((__fp_unordered_compare(x, y) \ - & 0x4500) == 0) -#define isless(x, y) ((__fp_unordered_compare (y, x) \ - & 0x4500) == 0) -#define isgreaterequal(x, y) ((__fp_unordered_compare (x, y) \ - & FP_INFINITE) == 0) -#define islessequal(x, y) ((__fp_unordered_compare(y, x) \ - & FP_INFINITE) == 0) -#define islessgreater(x, y) ((__fp_unordered_compare(x, y) \ - & FP_SUBNORMAL) == 0) -#define isunordered(x, y) ((__fp_unordered_compare(x, y) \ - & 0x4500) == 0x4500) - -#endif -#endif //0 - - #endif /* __STDC_VERSION__ >= 199901L */ #endif /* __NO_ISOCEXT */ #ifdef __cplusplus } -extern "C++" { - template inline _Ty _Pow_int(_Ty _X,int _Y) { - unsigned int _N; - if(_Y >= 0) _N = (unsigned int)_Y; - else _N = (unsigned int)(-_Y); - for(_Ty _Z = _Ty(1);;_X *= _X) { - if((_N & 1)!=0) _Z *= _X; - if((_N >>= 1)==0) return (_Y < 0 ? _Ty(1) / _Z : _Z); - } - } -} #endif - #pragma pack(pop) /* 7.12.14 */ diff --git a/win32/include/tcc/tcc_libm.h b/win32/include/tcc/tcc_libm.h index 90e26c47..2c6065ef 100644 --- a/win32/include/tcc/tcc_libm.h +++ b/win32/include/tcc/tcc_libm.h @@ -122,17 +122,9 @@ __CRT_INLINE long double __cdecl fmaxl (long double x, long double y) { /* *round* */ -#define TCCFP_FORCE_EVAL(x) do { \ -if (sizeof(x) == sizeof(float)) { \ - volatile float __x; \ - __x = (x); \ -} else if (sizeof(x) == sizeof(double)) { \ - volatile double __x; \ - __x = (x); \ -} else { \ - volatile long double __x; \ - __x = (x); \ -} \ +#define TCCFP_FORCE_EVAL(x) do { \ + volatile typeof(x) __x; \ + __x = (x); \ } while(0) __CRT_INLINE double __cdecl round (double x) { @@ -150,15 +142,8 @@ __CRT_INLINE double __cdecl round (double x) { return 0*u.f; } y = (double)(x + 0x1p52) - 0x1p52 - x; - if (y > 0.5) - y = y + x - 1; - else if (y <= -0.5) - y = y + x + 1; - else - y = y + x; - if (u.i >> 63) - y = -y; - return y; + y = y + x - (y > 0.5) + (y <= -0.5); /* branchless */ + return (u.i >> 63) ? -y : y; } __CRT_INLINE long __cdecl lround (double x) { @@ -194,17 +179,334 @@ __CRT_INLINE long long __cdecl llroundl (long double x) { } +/* MUSL asinh, acosh, atanh */ + +__CRT_INLINE double __cdecl asinh(double x) { + union {double f; uint64_t i;} u = {x}; + unsigned e = u.i >> 52 & 0x7ff, s = u.i >> 63; + u.i &= -1ull / 2, x = u.f; + if (e >= 0x3ff + 26) x = log(x) + 0.69314718055994530941723212145818; + else if (e >= 0x3ff + 1) x = log(2*x + 1 / (sqrt(x*x + 1) + x)); + else if (e >= 0x3ff - 26) x = log1p(x + x*x / (sqrt(x*x + 1) + 1)); + else TCCFP_FORCE_EVAL(x + 0x1p120f); + return s ? -x : x; +} + +__CRT_INLINE double __cdecl acosh(double x) { + union {double f; uint64_t i;} u = {x}; + unsigned e = u.i >> 52 & 0x7ff; + if (e < 0x3ff + 1) return log1p(x - 1 + sqrt((x - 1)*(x - 1) + 2*(x - 1))); + if (e < 0x3ff + 26) return log(2*x - 1 / (x + sqrt(x*x - 1))); + return log(x) + 0.69314718055994530941723212145818; +} + +__CRT_INLINE double __cdecl atanh(double x) { + union {double f; uint64_t i;} u = {x}; + unsigned e = u.i >> 52 & 0x7ff, s = u.i >> 63; + u.i &= -1ull / 2, x = u.f; + if (e < 0x3ff - 1) { + if (e < 0x3ff - 32) { if (e == 0) TCCFP_FORCE_EVAL((float)x); } + else x = 0.5 * log1p(2*x + 2*x*x / (1 - x)); /* |x| < 0.5 */ + } else x = 0.5 * log1p(2*(x / (1 - x))); /* avoid overflow */ + return s ? -x : x; +} + +/* MUSL scalbn */ + +__CRT_INLINE double __cdecl scalbn(double x, int n) { + union {double f; uint64_t i;} u; + if (n > 1023) { + x *= 0x1p1023, n -= 1023; + if (n > 1023) { + x *= 0x1p1023, n -= 1023; + if (n > 1023) n = 1023; + } + } else if (n < -1022) { + x *= 0x1p-1022 * 0x1p53, n += 1022 - 53; + if (n < -1022) { + x *= 0x1p-1022 * 0x1p53, n += 1022 - 53; + if (n < -1022) n = -1022; + } + } + u.i = (0x3ffull + n) << 52; + return x * u.f; +} + + /******************************************************************************* End of code based on MUSL *******************************************************************************/ +/* Following are math functions missing from msvcrt.dll, and not defined + * in math.h or above. Functions still remaining: + * remquo(), remainder(), fma(), nan(), erf(), erfc(), nearbyint(). + * In : lldiv(). + */ + +__CRT_INLINE float __cdecl scalbnf(float x, int n) { + return scalbn(x, n); +} +__CRT_INLINE long double __cdecl scalbnl(long double x, int n) { + return scalbn(x, n); +} + +__CRT_INLINE double __cdecl scalbln(double x, long n) { + int m = n > INT_MAX ? INT_MAX : (n < INT_MIN ? INT_MIN : n); + return scalbn(x, m); +} +__CRT_INLINE float __cdecl scalblnf(float x, long n) { + int m = n > INT_MAX ? INT_MAX : (n < INT_MIN ? INT_MIN : n); + return scalbn(x, m); +} + + +__CRT_INLINE float __cdecl ldexpf(float x, int expn) { + return scalbn(x, expn); +} +__CRT_INLINE long double __cdecl ldexpl(long double x, int expn) { + return scalbn(x, expn); +} + +__CRT_INLINE float __cdecl frexpf(float x, int *y) { + return frexp(x, y); +} +__CRT_INLINE long double __cdecl frexpl (long double x, int* y) { + return frexp(x, y); +} + + +__CRT_INLINE double __cdecl rint(double x) { +double retval; + __asm__ ( + "fldl %1\n" + "frndint \n" + "fstpl %0\n" : "=m" (retval) : "m" (x)); + return retval; +} + +__CRT_INLINE float __cdecl rintf(float x) { + float retval; + __asm__ ( + "flds %1\n" + "frndint \n" + "fstps %0\n" : "=m" (retval) : "m" (x)); + return retval; +} + + +/* 7.12.9.5 */ +__CRT_INLINE long __cdecl lrint(double x) { + long retval; + __asm__ __volatile__ + ("fldl %1\n" + "fistpl %0" : "=m" (retval) : "m" (x)); + return retval; +} + +__CRT_INLINE long __cdecl lrintf(float x) { + long retval; + __asm__ __volatile__ + ("flds %1\n" + "fistpl %0" : "=m" (retval) : "m" (x)); + return retval; +} + +__CRT_INLINE long long __cdecl llrint(double x) { +long long retval; + __asm__ __volatile__ + ("fldl %1\n" + "fistpll %0" : "=m" (retval) : "m" (x)); + return retval; +} + +__CRT_INLINE long long __cdecl llrintf(float x) { + long long retval; + __asm__ __volatile__ + ("flds %1\n" + "fistpll %0" : "=m" (retval) : "m" (x)); + return retval; +} + +/* + __CRT_INLINE long double __cdecl rintl (long double x) { + long double retval; + __asm__ ( + "fldt %1\n" + "frndint \n" + "fstpt %0\n" : "=m" (retval) : "m" (x)); + return retval; + } + __CRT_INLINE long __cdecl lrintl (long double x) { + long retval; + __asm__ __volatile__ + ("fldt %1\n" + "fistpl %0" : "=m" (retval) : "m" (x)); + return retval; + } + __CRT_INLINE long long __cdecl llrintl (long double x) { + long long retval; + __asm__ __volatile__ + ("fldt %1\n" + "fistpll %0" : "=m" (retval) : "m" (x)); + return retval; + } + __CRT_INLINE float __cdecl fabsf (float x) { + float res; + __asm__ ("fabs;" : "=t" (res) : "0" (x)); + return res; + } +*/ + +__CRT_INLINE long double __cdecl rintl(long double x) { + return rint(x); +} +__CRT_INLINE long __cdecl lrintl(long double x) { + return lrint(x); +} +__CRT_INLINE long long __cdecl llrintl(long double x) { + return llrint(x); +} + + +__CRT_INLINE double trunc(double _x) { + double retval; + unsigned short saved_cw; + unsigned short tmp_cw; + __asm__ ("fnstcw %0;" : "=m" (saved_cw)); /* save FPU control word */ + tmp_cw = (saved_cw & ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO)) + | FE_TOWARDZERO; + __asm__ ("fldcw %0;" : : "m" (tmp_cw)); + __asm__ ("fldl %1;" + "frndint;" + "fstpl %0;" : "=m" (retval) : "m" (_x)); /* round towards zero */ + __asm__ ("fldcw %0;" : : "m" (saved_cw) ); /* restore saved control word */ + return retval; +} + +__CRT_INLINE float __cdecl truncf(float x) { + return (float) ((int) x); +} +__CRT_INLINE long double __cdecl truncl(long double x) { + return trunc(x); +} + + +__CRT_INLINE long double __cdecl nextafterl(long double x, long double to) { + return nextafter(x, to); +} + +__CRT_INLINE double __cdecl nexttoward(double x, long double to) { + return nextafter(x, to); +} +__CRT_INLINE float __cdecl nexttowardf(float x, long double to) { + return nextafterf(x, to); +} +__CRT_INLINE long double __cdecl nexttowardl(long double x, long double to) { + return nextafter(x, to); +} + + +__CRT_INLINE float __cdecl fabsf (float x) { + union {float f; uint32_t i;} u = {x}; + u.i &= 0x7fffffff; + return u.f; +} +__CRT_INLINE long double __cdecl fabsl (long double x) { + return fabs(x); +} + + +#if defined(_WIN32) && !defined(_WIN64) && !defined(__ia64__) + __CRT_INLINE float acosf(float x) { return acos(x); } + __CRT_INLINE float asinf(float x) { return asin(x); } + __CRT_INLINE float atanf(float x) { return atan(x); } + __CRT_INLINE float atan2f(float x, float y) { return atan2(x, y); } + __CRT_INLINE float ceilf(float x) { return ceil(x); } + __CRT_INLINE float cosf(float x) { return cos(x); } + __CRT_INLINE float coshf(float x) { return cosh(x); } + __CRT_INLINE float expf(float x) { return exp(x); } + __CRT_INLINE float floorf(float x) { return floor(x); } + __CRT_INLINE float fmodf(float x, float y) { return fmod(x, y); } + __CRT_INLINE float logf(float x) { return log(x); } + __CRT_INLINE float logbf(float x) { return logb(x); } + __CRT_INLINE float log10f(float x) { return log10(x); } + __CRT_INLINE float modff(float x, float *y) { double di, df = modf(x, &di); *y = di; return df; } + __CRT_INLINE float powf(float x, float y) { return pow(x, y); } + __CRT_INLINE float sinf(float x) { return sin(x); } + __CRT_INLINE float sinhf(float x) { return sinh(x); } + __CRT_INLINE float sqrtf(float x) { return sqrt(x); } + __CRT_INLINE float tanf(float x) { return tan(x); } + __CRT_INLINE float tanhf(float x) { return tanh(x); } +#endif +__CRT_INLINE float __cdecl asinhf(float x) { return asinh(x); } +__CRT_INLINE float __cdecl acoshf(float x) { return acosh(x); } +__CRT_INLINE float __cdecl atanhf(float x) { return atanh(x); } + +__CRT_INLINE long double __cdecl asinhl(long double x) { return asinh(x); } +__CRT_INLINE long double __cdecl acoshl(long double x) { return acosh(x); } +__CRT_INLINE long double __cdecl atanhl(long double x) { return atanh(x); } +__CRT_INLINE long double __cdecl asinl(long double x) { return asin(x); } +__CRT_INLINE long double __cdecl acosl(long double x) { return acos(x); } +__CRT_INLINE long double __cdecl atanl(long double x) { return atan(x); } +__CRT_INLINE long double __cdecl ceill(long double x) { return ceil(x); } +__CRT_INLINE long double __cdecl coshl(long double x) { return cosh(x); } +__CRT_INLINE long double __cdecl cosl(long double x) { return cos(x); } +__CRT_INLINE long double __cdecl expl(long double x) { return exp(x); } +__CRT_INLINE long double __cdecl floorl(long double x) { return floor(x); } +__CRT_INLINE long double __cdecl fmodl(long double x, long double y) { return fmod(x, y); } +__CRT_INLINE long double __cdecl hypotl(long double x, long double y) { return hypot(x, y); } +__CRT_INLINE long double __cdecl logl(long double x) { return log(x); } +__CRT_INLINE long double __cdecl logbl(long double x) { return logb(x); } +__CRT_INLINE long double __cdecl log10l(long double x) { return log10(x); } +__CRT_INLINE long double __cdecl modfl(long double x, long double* y) { return modf(x, y); } +__CRT_INLINE long double __cdecl powl(long double x, long double y) { return pow(x, y); } +__CRT_INLINE long double __cdecl sinhl(long double x) { return sinh(x); } +__CRT_INLINE long double __cdecl sinl(long double x) { return sin(x); } +__CRT_INLINE long double __cdecl sqrtl(long double x) { return sqrt(x); } +__CRT_INLINE long double __cdecl tanhl(long double x) { return tanh(x); } +__CRT_INLINE long double __cdecl tanl(long double x) { return tan(x); } + +/* Following are accurate, but much shorter implementations than MUSL lib. */ + +__CRT_INLINE double __cdecl log1p(double x) { + double u = 1.0 + x; + return u == 1.0 ? x : log(u)*(x / (u - 1.0)); +} +__CRT_INLINE float __cdecl log1pf(float x) { + float u = 1.0f + x; + return u == 1.0f ? x : logf(u)*(x / (u - 1.0f)); +} +__CRT_INLINE long double __cdecl log1pl(long double x) { + return log1p(x); +} + + +__CRT_INLINE double __cdecl expm1(double x) { + if (fabs(x) > 0.0024) return exp(x) - 1.0; + double u, v = x*x, t = x + 0.5*v + 0.1666666666666666667*(u = v*x); + return t + 0.04166666666666666667*u*x; +} +__CRT_INLINE float __cdecl expm1f(float x) { + if (fabsf(x) > 0.085f) return expf(x) - 1.0f; + float u, v = x*x, t = x + 0.5f*v + 0.1666666666666666667f*(u = v*x); + return t + 0.04166666666666666667f*u*x; +} +__CRT_INLINE long double __cdecl expm1l(long double x) { + return expm1(x); +} + + __CRT_INLINE double __cdecl cbrt(double x) { - return (1.0 - ((x < 0.0) << 1)) * pow(fabs(x), 1.0 / 3.0); + return (1 - ((x < 0.0) << 1)) * pow(fabs(x), 1/3.0); } __CRT_INLINE float __cdecl cbrtf(float x) { - return (1.0f - ((x < 0.0f) << 1)) * (float) pow(fabs(x), 1.0 / 3.0); + return (1 - ((x < 0.0f) << 1)) * powf(fabsf((float) x), 1/3.0f); } +__CRT_INLINE long double __cdecl cbrtl(long double x) { + return (1 - ((x < 0.0) << 1)) * pow(fabs(x), 1/3.0); +} + __CRT_INLINE double __cdecl log2(double x) { return log(x) * 1.4426950408889634073599246810019; @@ -212,6 +514,10 @@ __CRT_INLINE double __cdecl log2(double x) { __CRT_INLINE float __cdecl log2f(float x) { return logf(x) * 1.4426950408889634073599246810019f; } +__CRT_INLINE long double __cdecl log2l(long double x) { + return log(x) * 1.4426950408889634073599246810019; +} + __CRT_INLINE double __cdecl exp2(double x) { return exp(x * 0.69314718055994530941723212145818); @@ -219,48 +525,86 @@ __CRT_INLINE double __cdecl exp2(double x) { __CRT_INLINE float __cdecl exp2f(float x) { return expf(x * 0.69314718055994530941723212145818f); } +__CRT_INLINE long double __cdecl exp2l(long double x) { + return exp(x * 0.69314718055994530941723212145818); +} + + +__CRT_INLINE int __cdecl ilogb(double x) { + return (int) logb(x); +} +__CRT_INLINE int __cdecl ilogbf(float x) { + return (int) logbf(x); +} +__CRT_INLINE int __cdecl ilogbl(long double x) { + return (int) logb(x); +} + + +__CRT_INLINE double __cdecl fdim(double x, double y) { + if (isnan(x) || isnan(y)) return NAN; + return x > y ? x - y : 0; +} +__CRT_INLINE float __cdecl fdimf(float x, float y) { + if (isnan(x) || isnan(y)) return NAN; + return x > y ? x - y : 0; +} +__CRT_INLINE long double __cdecl fdiml(long double x, long double y) { + if (isnan(x) || isnan(y)) return NAN; + return x > y ? x - y : 0; +} + /* tgamma and lgamma: Lanczos approximation * https://rosettacode.org/wiki/Gamma_function * https://www.johndcook.com/blog/cpp_gamma */ -__CRT_INLINE double __cdecl lgamma(double x); __CRT_INLINE double __cdecl tgamma(double x) { double m = 1.0, t = 3.14159265358979323; - if (x == (int)x) { + if (x > 172.0) + return INFINITY; + if (x == floor(x)) { for (int k = 2; k < x; ++k) m *= k; return x > 0.0 ? m : (x == 0.0 ? INFINITY : NAN); } if (x < 0.5) - return t / (sin(t * x) * tgamma(1.0 - x)); + return t / (sin(t*x)*tgamma(1.0 - x)); if (x > 12.0) return exp(lgamma(x)); + static const double c[8] = {676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7}; - m = 0.99999999999980993, t = x + (8 - 1.5); + m = 0.99999999999980993, t = x + 6.5; /* x-1+8-.5 */ for (int k = 0; k < 8; ++k) m += c[k] / (x + k); - return 2.50662827463100050 * pow(t, x - 0.5) * exp(-t) * m; /* sqrt(2pi) */ + return 2.50662827463100050 * pow(t, x - 0.5)*exp(-t)*m; /* C=sqrt(2pi) */ } + __CRT_INLINE double __cdecl lgamma(double x) { if (x < 12.0) - return x <= 0.0 && x == (int)x ? INFINITY : log(fabs(tgamma(x))); + return x > 0.0 || x != floor(x) ? log(fabs(tgamma(x))) : INFINITY; + static const double c[7] = {1.0/12.0, -1.0/360.0, 1.0/1260.0, -1.0/1680.0, 1.0/1188.0, -691.0/360360.0, 1.0/156.0}; - double m = -3617.0/122400.0, z = 1.0 / (x * x); - for (int k = 6; k >= 0; --k) m = m * z + c[k]; - return (x - 0.5) * log(x) - x + 0.918938533204672742 + m / x; /* log(2pi)/2 */ + double m = -3617.0/122400.0, t = 1.0 / (x*x); + for (int k = 6; k >= 0; --k) m = m*t + c[k]; + return (x - 0.5)*log(x) - x + 0.918938533204672742 + m / x; /* C=log(2pi)/2 */ } __CRT_INLINE float __cdecl tgammaf(float x) { return tgamma(x); } - __CRT_INLINE float __cdecl lgammaf(float x) { return lgamma(x); } +__CRT_INLINE long double __cdecl tgammal(long double x) { + return tgamma(x); +} +__CRT_INLINE long double __cdecl lgammal(long double x) { + return lgamma(x); +} #endif /* _TCC_LIBM_H_ */