Overhauled WIN32 math and added missing functions:
1) Cleanup: moved function implementations from win32/include/math.h to win32/include/tcc/tcc_libm.h 2) Added missing math functions: MUSL: asinh(), acosh(), atanh(), scalbn(). My impl: log1p(), expm1(), ilogb(), scalbln(), nexttoward() - now only a few are missing: remquo(), remainder(), fma(), nan(), erf(), erfc(), nearbyint(). 3) Added/defined all missing *f() and *l() math functions. 4) Added functions have short but accurate/fast implementations. 5) Added <tgmath.h> for all platforms. (not too useful, IMO, but is C99 standard).
This commit is contained in:
		
							parent
							
								
									82611f5e6d
								
							
						
					
					
						commit
						b11144d69c
					
				
					 3 changed files with 484 additions and 321 deletions
				
			
		
							
								
								
									
										89
									
								
								include/tgmath.h
									
										
									
									
									
										Normal file
									
								
							
							
						
						
									
										89
									
								
								include/tgmath.h
									
										
									
									
									
										Normal file
									
								
							|  | @ -0,0 +1,89 @@ | |||
| /*
 | ||||
|  * ISO C Standard:  7.22  Type-generic math <tgmath.h> | ||||
|  */ | ||||
| 
 | ||||
| #ifndef _TGMATH_H | ||||
| #define _TGMATH_H | ||||
| 
 | ||||
| #include <math.h> | ||||
| 
 | ||||
| #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 <math.h> and <complex.h> (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 <math.h> 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 <complex.h> 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 */ | ||||
|  | @ -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<double>(_Number)); } | ||||
|   __CRT_INLINE long double _copysignl(long double _Number,long double _Sign) { return _copysign(static_cast<double>(_Number),static_cast<double>(_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<class _Ty> 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 */ | ||||
|  |  | |||
|  | @ -123,16 +123,8 @@ __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;                       \ | ||||
|   volatile typeof(x) __x;         \ | ||||
|   __x = (x);                      \ | ||||
| } else if (sizeof(x) == sizeof(double)) {   \ | ||||
|   volatile double __x;                      \ | ||||
|   __x = (x);                                \ | ||||
| } else {                                    \ | ||||
|   volatile long double __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 <stdlib.h>: 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_ */ | ||||
|  |  | |||
		Loading…
	
	Add table
		
		Reference in a new issue