diff --git a/lib/std/math/math_nolibc/exp.c3 b/lib/std/math/math_nolibc/exp.c3 new file mode 100644 index 000000000..ce19a4d4e --- /dev/null +++ b/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); +} diff --git a/lib/std/math/math_nolibc/fabs.c3 b/lib/std/math/math_nolibc/fabs.c3 new file mode 100644 index 000000000..d434ca8af --- /dev/null +++ b/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); +} diff --git a/lib/std/math/math_nolibc/frexp.c3 b/lib/std/math/math_nolibc/frexp.c3 new file mode 100644 index 000000000..e0565bba8 --- /dev/null +++ b/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); +} diff --git a/lib/std/math/math_nolibc/ldexp.c3 b/lib/std/math/math_nolibc/ldexp.c3 new file mode 100644 index 000000000..4b0388718 --- /dev/null +++ b/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); +} diff --git a/lib/std/math/math_nolibc/log.c3 b/lib/std/math/math_nolibc/log.c3 new file mode 100644 index 000000000..77301348f --- /dev/null +++ b/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); +} \ No newline at end of file diff --git a/lib/std/math/math_nolibc/pow.c3 b/lib/std/math/math_nolibc/pow.c3 index 3637834b9..a039e5998 100644 --- a/lib/std/math/math_nolibc/pow.c3 +++ b/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; +} \ No newline at end of file