Feat: Added exp, log and pow functions as requested in #1632 (#1781)

* Feat: Added exp, log and pow functions as requested in #1632
---------

Co-authored-by: Christoffer Lerno <christoffer@aegik.com>

lib/std/math/math_nolibc/exp.c3

@@ -0,0 +1,60 @@
+module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
+
+const double EXP_LN2_HI = 6.93147180369123816490e-01;
+const double EXP_LN2_LO = 1.90821492927058770002e-10;
+const double EXP_INV_LN2 = 1.44269504088896338700e+00;
+const double EXP_P1 = 1.66666666666666019037e-01;
+const double EXP_P2 = -2.77777777770155933842e-03;
+const double EXP_P3 = 6.61375632143793436117e-05;
+const double EXP_P4 = -1.65339022054652515390e-06;
+const double EXP_P5 = 4.13813679705723846039e-08;
+/*--------------------------------------------*/
+const float EXPF_LN2_HI = 6.9314575195e-01f;
+const float EXPF_LN2_LO = 1.4286067653e-06f;
+const float EXPF_INV_LN2 = 1.4426950216e+00f;
+const float EXPF_P1 = 1.6666667163e-01f;
+const float EXPF_P2 = -2.7777778450e-03f;
+const float EXPF_P3 = 6.6137559770e-05f;
+const float EXPF_P4 = -1.6533901999e-06f;
+
+fn double exp(double x) @extern("exp")
+{
+    if (x != x) return x;
+    if (x == double.inf) return double.inf;
+    if (x == -double.inf) return 0.0; // IEEE 754 spec
+    // Overflow threshold for exp (approx +709.78 for double)
+    if (x > 709.782712893384) return double.inf;
+    // Underflow threshold for exp (approx -745.13 for double)
+    if (x < -745.133219101941) return 0.0;
+
+    double px = x * EXP_INV_LN2;
+    double k = _floor(px + 0.5);
+    double r = x - k * EXP_LN2_HI - k * EXP_LN2_LO;
+
+    double r2 = r * r;
+    double p = r2 * (EXP_P1 + r2 * (EXP_P2 + r2 * (EXP_P3 + r2 * (EXP_P4 + r2 * EXP_P5))));
+    double exp_r = 1.0 + r + r * p;
+
+    return ldexp(exp_r, (int)k);
+}
+
+fn float expf(float x) @extern("expf")
+{
+    if (x != x) return x;
+    if (x == float.inf) return float.inf;
+    if (x == -float.inf) return 0.0f; // IEEE 754 spec
+    // Overflow threshold (approx +88.72 for float)
+    if (x > 88.7228f) return float.inf;
+    // Underflow threshold (approx -103.97 for float)
+    if (x < -103.972084f) return 0.0f;
+
+    float px = x * EXPF_INV_LN2;
+    float k = _floorf(px + 0.5f);
+    float r = x - k * EXPF_LN2_HI - k * EXPF_LN2_LO;
+
+    float r2 = r * r;
+    float p = r2 * (EXPF_P1 + r2 * (EXPF_P2 + r2 * (EXPF_P3 + r2 * EXPF_P4)));
+    float exp_r = 1.0f + r + r * p;
+
+    return ldexpf(exp_r, (int)k);
+}

lib/std/math/math_nolibc/fabs.c3

@@ -0,0 +1,15 @@
+module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
+
+fn double _fabs(double x) @weak @extern("fabs") @nostrip
+{
+    ulong ix = bitcast(x, ulong);
+    ix &= ~(1ul << 63);
+    return bitcast(ix, double);
+}
+
+fn float _fabsf(float x) @weak @extern("fabsf") @nostrip
+{
+    uint ix = bitcast(x, uint);
+    ix &= 0x7fffffff;
+    return bitcast(ix, float);
+}

lib/std/math/math_nolibc/frexp.c3

@@ -0,0 +1,59 @@
+module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
+
+fn double frexp(double x, int* exp) @extern("frexp")
+{
+    uint hx = x.high_word();
+    uint ix = hx & 0x7fffffff;
+    uint lx = x.low_word();
+
+    if (ix >= 0x7ff00000 || (ix | lx) == 0)
+    {
+        *exp = 0;
+        return x;
+    }
+
+    // exponent extraction and normalization
+    int e = (int)((ix >> 20) & 0x7ff);
+    if (e == 0)
+    {
+        // subnormal number
+        x *= 0x1p64;
+        hx = x.high_word();
+        e = (int)((hx >> 20) & 0x7ff) - 64;
+    }
+    *exp = e - 1022;
+
+    // set exponent to -1 (fraction in [0.5, 1))
+    hx = (hx & 0x800fffff) | 0x3fe00000;
+    {
+        ulong rep = ((ulong)hx << 32) | lx;
+        return bitcast(rep, double);
+    }
+}
+
+fn float frexpf(float x, int* exp) @extern("frexpf")
+{
+    uint ix = x.word();
+    uint hx = ix & 0x7fffffff;
+
+    if (hx >= 0x7f800000 || hx == 0)
+    {
+        *exp = 0;
+        return x;
+    }
+
+    // exponent extraction and normalization
+    int e = (int)((hx >> 23) & 0xff);
+    if (e == 0)
+    {
+        // subnormal number
+        x *= 0x1p64f;
+        ix = x.word();
+        e = (int)((ix >> 23) & 0xff) - 64;
+    }
+    *exp = e - 126;
+
+    // set exponent to -1 (fraction in [0.5, 1))
+    ix = (ix & 0x807fffff) | 0x3f000000;
+    return bitcast(ix, float);
+}

lib/std/math/math_nolibc/ldexp.c3

@@ -0,0 +1,67 @@
+module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
+
+fn double ldexp(double x, int exp) @extern("ldexp") 
+{
+    uint hx = x.high_word();
+    int hexp = (int)((hx & 0x7ff00000) >> 20);
+    
+ 
+    if (hexp == 0x7ff) return x;
+    if (hexp == 0)
+    {
+        // subnormal number handling
+        x *= 0x1p64;
+        hx = x.high_word();
+        hexp = (int)((hx & 0x7ff00000) >> 20) - 64;
+    }
+
+    // new exponent calculation
+    hexp += exp;
+
+    if (hexp > 0x7fe) return x * double.inf;
+    if (hexp < 1)
+    {
+        x *= 0x1p-1022;
+        hexp += 1022;
+        if (hexp < 1) x *= 0x1p-1022;
+        return x;
+    }
+
+    // set new exponent
+    hx = ((ulong)hx & 0x800fffff) | ((ulong)hexp << 20);
+    {
+        ulong rep = ((ulong)hx << 32) | x.low_word();
+        return bitcast(rep, double);
+    }
+}
+
+fn float ldexpf(float x, int exp) @extern("ldexpf")
+{
+    uint ix = x.word();
+    int hexp = (int)((ix & 0x7f800000) >> 23);
+
+    if (hexp == 0xff) return x;
+    if (hexp == 0)
+    {
+        // subnormal number handling
+        x *= 0x1p64f;
+        ix = x.word();
+        hexp = (int)((ix & 0x7f800000) >> 23) - 64;
+    }
+
+    // new exponent calculation
+    hexp += exp;
+
+    if (hexp > 0xfe) return x * float.inf;
+    if (hexp < 1)
+    {
+        x *= 0x1p-126f;
+        hexp += 126;
+        if (hexp < 1) x *= 0x1p-126f;
+        return x;
+    }
+
+    // set new exponent
+    ix = (ix & 0x807fffff) | (hexp << 23);
+    return bitcast(ix, float);
+}

lib/std/math/math_nolibc/log.c3

@@ -0,0 +1,86 @@
+module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
+
+const double LOG_LN2_HI = 6.93147180369123816490e-01;
+const double LOG_LN2_LO = 1.90821492927058770002e-10;
+const double LOG_L1 = 6.666666666666735130e-01;
+const double LOG_L2 = 3.999999999940941908e-01;
+const double LOG_L3 = 2.857142874366239149e-01;
+const double LOG_L4 = 2.222219843214978396e-01;
+const double LOG_L5 = 1.818357216161805012e-01;
+const double LOG_L6 = 1.531383769920937332e-01;
+/*--------------------------------------------*/
+const float LOGF_LN2_HI = 6.9313812256e-01f;
+const float LOGF_LN2_LO = 9.0580006145e-06f;
+const float LOGF_L1 = 6.6666662693e-01f;
+const float LOGF_L2 = 4.0000972152e-01f;
+const float LOGF_L3 = 2.8498786688e-01f;
+const float LOGF_L4 = 2.4279078841e-01f;
+
+const double SQRT2 = 1.41421356237309504880;
+const float SQRT2F = 1.41421356237309504880f;
+
+fn double log(double x) @extern("log")
+{
+    if (x != x) return x;
+    if (x < 0.0) return double.nan;
+    if (x == 0.0) return -double.inf;
+    if (x == double.inf) return double.inf;
+
+    int k;
+    double f = frexp(x, &k);
+    if (f < SQRT2 * 0.5)
+    {
+        f *= 2.0;
+        k--;
+    }
+
+    // polynomial approximation of log(1 + f), with f in [0, sqrt(2) - 1]
+    f -= 1.0;
+    double s = f / (2.0 + f);
+    double z = s * s;
+    double w = z * z;
+
+    // even-part polynomial terms (t1) and odd-part polynomial terms (t2)
+    double t1 = w * (LOG_L1 + w * (LOG_L3 + w * LOG_L5));
+    double t2 = z * (LOG_L2 + w * (LOG_L4 + w * LOG_L6));
+    double r = t1 + t2;
+
+    double hfsq = 0.5 * f * f;
+
+    return k * LOG_LN2_HI - ((hfsq - (s * (hfsq + r) + k * LOG_LN2_LO)) - f);
+}
+
+fn float logf(float x) @extern("logf")
+{
+    if (x != x) return x;
+    if (x < 0.0f) return float.nan;
+    if (x == 0.0f) return -float.inf;
+    if (x == float.inf) return float.inf;
+
+    int k;
+    float f = frexpf(x, &k);
+    if (f < SQRT2F * 0.5f)
+    {
+        f *= 2.0f;
+        k--;
+    }
+
+    // polynomial approximation for log(1 + f)
+    f -= 1.0f;
+    float s = f / (2.0f + f);
+    float z = s * s;
+    float w = z * z;
+
+    /*
+        logf uses fewer terms in its polynomial approximation
+        compared to the double precision version because
+        single-precision floating point doesn't benefit from the additional terms.
+        LOGF_L1, ... ,LOGF_L4 provide sufficient accuracy for 32-bit float calculations.
+    */
+    float t1 = w * (LOGF_L1 + w * LOGF_L3);
+    float t2 = z * (LOGF_L2 + w * LOGF_L4);
+    float r = t1 + t2;
+
+    float hfsq = 0.5f * f * f;
+    return k * LOGF_LN2_HI - ((hfsq - (s * (hfsq + r) + k * LOGF_LN2_LO)) - f);
+}

lib/std/math/math_nolibc/pow.c3

@@ -1,11 +1,109 @@
 module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
 
-fn float powf_broken(float x, float f) @extern("powf") @weak @nostrip
+fn double pow(double x, double y) @extern("pow")
 {
-	unreachable("'powf' not supported");
+    if (x != x || y != y) return double.nan;
+
+    if (y == double.inf)
+    {
+        if (x == 1.0 || x == -1.0) return 1.0;
+        return (_fabs(x) < 1.0) ? 0.0 : double.inf;
+    }
+    if (y == -double.inf)
+    {
+        if (x == 1.0 || x == -1.0) return 1.0;
+        return (_fabs(x) < 1.0) ? double.inf : 0.0;
+    }
+    if (x == double.inf)
+    {
+        return (y < 0.0) ? 0.0 : double.inf;
+    }
+    if (x == -double.inf)
+    {
+        if (y != _floor(y)) return -double.nan;
+        if (y < 0.0) return 0.0;
+        return ((int)y & 1) ? -double.inf : double.inf;
+    }
+
+    if (y == 0.0) return 1.0;
+
+    if (x == 0.0)
+    {
+        if (y < 0.0) return double.inf;
+        if (y > 0.0) return 0.0;
+        return 1.0;
+    }
+
+    if (y == 1.0) return x;
+    if (x == 1.0) return 1.0;
+
+    if (x < 0.0)
+    {
+        if (y != _floor(y)) return double.nan;
+        return ((int)y & 1) ? -pow(-x, y) : pow(-x, y);
+    }
+
+    double result = exp(y * log(x));
+
+    if (result == double.inf || result == -double.inf)
+    {
+        return (y < 0.0) ? 0.0 : double.inf;
+    }
+    if (result == 0.0) return 0.0;
+
+    return result;
 }
 
-fn double pow_broken(double x, double y) @extern("pow") @weak @nostrip
+fn float powf(float x, float y) @extern("powf")
 {
-	unreachable("'pow' not supported");
+    if (x != x || y != y) return float.nan;
-}
+
+    if (y == float.inf)
+    {
+        if (x == 1.0f || x == -1.0f) return 1.0f;
+        return (_fabsf(x) < 1.0f) ? 0.0f : float.inf;
+    }
+    if (y == -float.inf)
+    {
+        if (x == 1.0f || x == -1.0f) return 1.0f;
+        return (_fabsf(x) < 1.0f) ? float.inf : 0.0f;
+    }
+    if (x == float.inf)
+    {
+        return (y < 0.0f) ? 0.0f : float.inf;
+    }
+    if (x == -float.inf)
+    {
+        if (y != _floorf(y)) return float.nan;
+        if (y < 0.0f) return 0.0f;
+        return ((int)y & 1) ? -float.inf : float.inf;
+    }
+
+    if (y == 0.0f) return 1.0f;
+
+    if (x == 0.0f)
+    {
+        if (y < 0.0f) return float.inf;
+        if (y > 0.0f) return 0.0f;
+        return 1.0f;
+    }
+
+    if (y == 1.0f) return x;
+    if (x == 1.0f) return 1.0f;
+
+    if (x < 0.0f)
+    {
+        if (y != _floorf(y)) return float.nan;
+        return ((int)y & 1) ? -powf(-x, y) : powf(-x, y);
+    }
+
+    float result = expf(y * logf(x));
+
+    if (result == float.inf || result == -float.inf)
+    {
+        return (y < 0.0f) ? 0.0f : float.inf;
+    }
+    if (result == 0.0f) return 0.0f;
+
+    return result;
+}