Added tests for math functions and fixed bugs

Author: ZERICO2005

src/libc/erfcf.c

@@ -6,7 +6,7 @@
  *
  * @remarks Minimum ulp:
  * ulp of +16     at +0x1.a13834p-1 with ideal expf (x < 1.0f)
- * ulp of +594482 at +0x1.251634p+3 with ideal expf (x >= 1.0f)
+ * ulp of +594482 at +0x1.251634p+3 with current expf (x >= 1.0f)
  */
 float erfcf(float x)
 {

src/libc/erfcl.c

@@ -7,7 +7,7 @@
  * @remarks Minimum ulp:
  * ulp of -516 at +2.373300261e+01 with ideal expl
  */
-long double _erfcl_c(long double arg) {
+long double erfcl(long double arg) {
     long double x, x_sqr;
     long double t0, t1, t2, t3, t4, t5;
     long double ret;

src/libc/expf.c

@@ -49,7 +49,7 @@ float _expf_c(float arg) {
     if ( arg == 0.0f){
         return 1.0f;
     }
-    if ( arg < -exp_min_arg ){
+    if ( arg < exp_min_arg ){
         return 0.0f;
     }
     if ( arg > exp_max_arg ){

src/libc/include/math.h

@@ -27,7 +27,6 @@ extern "C" {
 #define M_2_SQRTPI    1.12837916709551257390     /* 2/sqrt(pi)     */
 #define M_SQRT2       1.41421356237309504880     /* sqrt(2)        */
 #define M_SQRT1_2     0.707106781186547524401    /* 1/sqrt(2)      */
-// /* deprecated */ #define M_LOG_2M_PI   1.83787706640934548        /* log(2*M_PI)    */
 
 #define FP_ILOGB0     (~__INT_MAX__)
 #define FP_ILOGBNAN     __INT_MAX__

src/libc/include/tgmath.h

@@ -9,7 +9,6 @@
 
 #include <math.h>
 
-
 #define __tgmath_promote(x) _Generic((x), \
     float: ((float)0.f), \
     default: ((double)0.), \
@@ -378,8 +377,6 @@
     float: __hypot3f \
 )((x), (y), (z))
 
-#define __tgmath_hypot2(x,y)2
-#define __tgmath_hypot3(x,y,z)3
 #define __tgmath_hypot_count(_1,_2,_3,count,...)count
 #define __tgmath_JOIN2(X,Y)X##Y
 #define __tgmath_CONCAT(X,Y)__tgmath_JOIN2(X,Y)

src/libc/isnormall.src

@@ -20,7 +20,7 @@ __isnormall:
 	ld	d, $40
 	; test for zero/subnormal
 	sbc	hl, de
-	jr	c, .exp_all_zero
+;	jr	c, .exp_all_zero
 
 .exp_all_ones:
 .exp_all_zero:

src/libc/lgammaf.c

@@ -42,7 +42,7 @@ float lgammaf(float x) { /* the natural logarithm of the Gamma function. */
      */
     const int maximum_iter = 36 + N;
     for (int iter = 0; iter < maximum_iter; iter++) {
-        if (x < (float)N) {
+        if (!(x < (float)N)) {
             break;
         }
         v *= x;
@@ -53,7 +53,7 @@ float lgammaf(float x) { /* the natural logarithm of the Gamma function. */
                 + (B12 / (12.0f * 11.0f))) * w + (B10 / (10.0f *  9.0f))) * w
                 + (B8  / ( 8.0f *  7.0f))) * w + (B6  / ( 6.0f *  5.0f))) * w
                 + (B4  / ( 4.0f *  3.0f))) * w + (B2  / ( 2.0f *  1.0f))) / x
-                + ln_pi_div_2 - logf(v) - x + (x - 0.5f) * logf(x);
+                + ln_pi_div_2 - logf(fabsf(v)) - x + (x - 0.5f) * logf(fabsf(x));
 }
 
 double lgamma(double) __attribute__((alias("lgammaf")));

src/libc/sinf.c

@@ -11,16 +11,16 @@
 */
 #include <math.h>
 
-#define two_over_pi  0.636619772367581
-#define p0           0.135788409787738e8
-#define p1           -0.494290810090284e7
-#define p2           0.440103053537527e6
-#define p3           -0.138472724998245e5
-#define p4           0.145968840666577e3
-#define q0           0.864455865292253e7
-#define q1           0.408179225234330e6
-#define q2           0.946309610153821e4
-#define q3           0.132653490878614e3
+#define two_over_pi  0.63661977236758134f
+#define p0           0.135788409787738e8f
+#define p1           -0.494290810090284e7f
+#define p2           0.440103053537527e6f
+#define p3           -0.138472724998245e5f
+#define p4           0.145968840666577e3f
+#define q0           0.864455865292253e7f
+#define q1           0.408179225234330e6f
+#define q2           0.946309610153821e4f
+#define q3           0.132653490878614e3f
 
 float sinus(float arg, int quad)
 {
@@ -31,34 +31,44 @@ float sinus(float arg, int quad)
     float temp1, temp2;
 
     x = arg;
-    if(x<0) {
+    if (x<0.0f) {
         x = -x;
         quad = quad + 2;
     }
-    x = x*two_over_pi; /*underflow?*/
-    if(x>32764){
+    x = x * two_over_pi; /* underflow? */
+    if (x > 32764.0f) {
         y = modff(x,&e);
         e = e + quad;
-        modff(0.25*e,&f);
-        quad = e - 4*f;
-    }else{
+        modff(0.25f * e,&f);
+        quad = e - 4.0f * f;
+    } else {
         k = x;
         y = x - k;
-        quad = (quad + k) & 03;
+        quad = (quad + k) & 0x3;
     }
-    if (quad & 01)
-        y = 1-y;
-    if(quad > 1)
+    if (quad & 0x1) {
+        y = 1.0f - y;
+    }
+    if (quad > 1) {
         y = -y;
+    }
 
     ysq = y*y;
     temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
     temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
     return(temp1/temp2);
 }
 
+/**
+ * @remarks Minimum ulp:
+ * ulp of -5 at +0x1.fe2dd0p-9 (|x| < pi/2)
+ */
 float _sinf_c(float arg) {
+    if (fabsf(arg) < 0x1.0p-11f) {
+        return arg;
+    }
     return sinus(arg, 0);
 }
 
+
 double _sin_c(double) __attribute__((alias("_sinf_c")));

src/libc/tgammaf.c

@@ -29,7 +29,7 @@ float tgammaf(float x) { /* Gamma function */
             return zero/zero; /* probably better to return NAN here */
         }
         ret = (float)M_PI / (sinf((float)M_PI * f) * expf(lgammaf(1.0f - x)));
-        if (((unsigned int)i & 0x1) == 0) {
+        if (((unsigned char)i & 0x1) == 0) {
             ret = -ret;
         }
         return ret;

src/libcxx/include/cmath

@@ -239,7 +239,7 @@ template<typename _Tp, typename _Up> inline constexpr enable_if_t<is_integral_v<
 hypot(_Tp __x, _Up __y) { return hypot(__x, __y); }
 
 inline float hypot(float __x, float __y, float __z) { return __hypot3f(__x, __y, __z); }
-inline long double hypot(long double __x, long double __y, long double __z) { __hypot3l(__x, __y, __z); }
+inline long double hypot(long double __x, long double __y, long double __z) { return __hypot3l(__x, __y, __z); }
 template<typename _Tp, typename _Up, typename _Vp>
 inline enable_if_t<is_integral_v<_Tp> && is_integral_v<_Up> && is_integral_v<_Vp>, double>
 hypot(_Tp __x, _Up __y, _Vp __z) { return __hypot3(__x, __y, __z); }