Implemented math.h

libc/math/abs.c

@@ -0,0 +1,10 @@
+#include <math.h>
+
+int abs(int n) {
+    return (n < 0) ? -n : n;
+}
+
+double fabs(double x) {
+    return (x < 0.0) ? -x : x;
+}
+

libc/math/const.c

@@ -0,0 +1,43 @@
+#include <math.h>
+
+double nan(const char* tagp) {
+    (void)tagp;
+
+    // Create a NaN bit pattern in IEEE 754
+    uint64_t nanBits = 0x7ff8000000000000ULL;
+    double result;
+    memcpy(&result, &nanBits, sizeof(result));
+    return result;
+}
+
+double inf() {
+    // Create an infinity bit pattern in IEEE 754
+    uint64_t infBits = 0x7ff0000000000000ULL;
+    double result;
+    memcpy(&result, &infBits, sizeof(result));
+    return result;
+}
+
+int isfinite(double x) {
+    uint64_t bits;
+    memcpy(&bits, &x, sizeof(bits));
+    uint64_t exponent = (bits >> 52) & 0x7ff;
+    return exponent != 0x7ff;
+}
+
+int isinf(double x) {
+    uint64_t bits;
+    memcpy(&bits, &x, sizeof(bits));
+    uint64_t exponent = (bits >> 52) & 0x7ff;
+    uint64_t mantissa = bits & 0xfffffffffffff;
+    return (exponent == 0x7ff && mantissa == 0);
+}
+
+int isnan(double x) {
+    uint64_t bits;
+    memcpy(&bits, &x, sizeof(bits));
+    uint64_t exponent = (bits >> 52) & 0x7ff;
+    uint64_t mantissa = bits & 0xfffffffffffff;
+    return (exponent == 0x7ff && mantissa != 0);
+}
+

libc/math/exp.c

@@ -0,0 +1,46 @@
+#include <math.h>
+
+#define _EXP_TAYLOR_TERMS   10
+#define _EXP2_TAYLOR_TERMS  5
+
+double exp(double x) {
+    if (x == 0.0) return 1.0;
+    if (x < 0.0) return 1.0 / exp(-x);
+
+    double sum = 1.0;   // Result
+    double term = 1.0;  // Curent term in the series
+    
+    /*
+     Calculate the Taylor series:
+     e^x = 1 + x + x^2/2! + x^3/3! + ...
+     */
+    for (int i = 1; i < _EXP_TAYLOR_TERMS; i++) {
+        term *= x / (double)i;
+        sum += term;
+    }
+
+    return sum;
+}
+
+double exp2(double x) {
+    if (x == 0.0) return 1.0;
+    if (x < 0.0) return 1.0 / exp2(-x);
+
+    int n = (int)x;     // Integer part
+    double frac = x - n;// Fractional part
+
+    // Integer power of 2
+    double result = (double)(1 << n);   // 2^n is just 1 << n
+
+    // Approximate the fractional part
+    double sum = 1.0;
+    double term = frac;
+
+    for (int i = 1; i < _EXP2_TAYLOR_TERMS; i++) {
+        sum += term;
+        term *= frac / (i+1);
+    }
+
+    return result * sum;
+}
+

libc/math/htrig.c

@@ -0,0 +1,16 @@
+#include <math.h>
+
+double sinh(double x) {
+    return (exp(x) - exp(-x)) / 2.0;
+}
+
+double cosh(double x) {
+    return (exp(x) + exp(-x)) / 2.0;
+}
+
+double tanh(double x) {
+    double ex = exp(x);
+    double emx = exp(-x);
+    return (ex - exp) / (ex + exp);
+}
+

libc/math/log.c

@@ -0,0 +1,35 @@
+#include <math.h>
+
+#define _LOG_TERMS  20
+
+double log(double x) {
+    if (x <= 0.0) return -INFINITY;
+
+    double result = 0.0;
+    double term = (x - 1.0) / (x + 1.0);
+    double termSquared = term*term;
+    double numerator = term;
+
+    for (int i = 1; i < _LOG_TERMS; i++) {
+        result += numerator / (double)i;
+        numerator *= termSquared;
+    }
+
+    return 2.0 * result;
+}
+
+double log2(double x) {
+    if (x <= 0.0) return -INFINITY;
+    return log(x) / M_LN2;
+}
+
+double log10(double x) {
+    if (x <= 0.0) return -INFINITY;
+    return log(x) / M_LN10;
+}
+
+double logn(double x, double n) {
+    if (x <= 0.0) return -INFINITY;
+    return log(x) / log(n);
+}
+

libc/math/pow.c

@@ -0,0 +1,42 @@
+#include <math.h>
+
+#define _POW_LOOP_MAX   60
+
+double sqrt(double x) {
+    if (x < 0.0) return -NAN;
+
+    double result = x;
+    double oldResult = 0.0;
+
+    while (result != oldResult) {
+        oldResult = result;
+        result = 0.5 * (result + x/result);
+    }
+
+    return result;
+}
+
+double cbrt(double x) {
+    double result = x;
+    double oldResult = 0.0;
+
+    while (result != oldResult) {
+        oldResult = result;
+        result = (2.0 * result + x/(result * result)) / 3.0;
+    }
+
+    return result;
+}
+
+double hypot(double x, double y) {
+    return sqrt(x*x + y*y);
+}
+
+double pow(double base, double expn) { // exp is a function so we have to use expn
+    if (base == 0.0 && expn == 0.0) return NAN;
+    if (base == 0.0) return 0.0;
+    if (expn == 0.0) return 1.0;
+
+    return exp(expn * log(base));
+}
+

libc/math/round.c

@@ -0,0 +1,46 @@
+#include <math.h>
+
+double round(double x) {
+    double res;
+    __asm__ volatile (
+        "fldl %1\n"       // Load x to the FPU stack
+        "frndint\n"       // Round to nearest integer
+        "fstpl %0"        // Store result
+        : "=m"(res)
+        : "m"(x)
+    );
+    if ((x - res) == 0.5 || (x - res) == -0.5) {
+        // Halfway cases: round away from zero
+        res += (x > 0) ? 1.0 : -1.0;
+    }
+    return res;
+}
+
+double floor(double x) {
+    double res = round(x);
+    if (res > x) { res -= 1.0; }
+    return res;
+}
+
+double ceil(double x){
+    double res = round(x);
+    if (res < x) { res += 1.0; }
+    return res;
+}
+
+double trunc(double x) {
+    double res = round(x);
+    if (x < 0 && x != res) { res -= 1.0; }
+    return res;
+}
+
+double fmod(double x, double y) {
+    return x - y * floor(x / y);
+}
+
+double modf(double x, double* iptr) {
+    double whole = trunc(x);
+    *iptr = whole;
+    return x - whole;
+}
+

libc/math/trig.c

@@ -0,0 +1,61 @@
+#include <math.h>
+
+double sin(double x) {
+    double res;
+    __asm__ volatile (
+        "fldl %1\n"     // Load x to the FPU stack
+        "fsin\n"        // Compute sin(x)
+        "fstpl %0"      // Store result
+        : "=m"(res) : "m"(x)
+    );
+    return res;
+}
+
+double cos(double x) {
+    double res;
+    __asm__ volatile (
+        "fldl %1\n"     // Load x to the FPU stack
+        "fcos\n"        // Compute cos(x)
+        "fstpl %0"      // Store result
+        : "=m"(res) : "m"(x)
+    );
+    return res;
+}
+
+
+double tan(double x) {
+    return sin(x) / cos(x);
+}
+
+double asin(double x) {
+    return atan2(x, sqrt(1 - x*x));
+}
+
+double acos(double x) {
+    return M_PI_2 - asin(x);
+}
+
+double atan(double x) {
+    double res;
+    __asm__ volatile (
+        "fld1\n"        // Load 1.0 to the FPU stack
+        "fldl %1\n"     // Load x to the FPU stack
+        "fpatan\n"      // Compute atan(x)
+        "fstpl %0"      // Store result
+        : "=m"(res) : "m"(x)
+    );
+    return res;
+}
+
+double atan2(double y, double x) {
+    double res;
+    __asm__ volatile (
+        "fldl %2\n"       // Load x to the FPU stack
+        "fldl %1\n"       // Load y to the FPU stack
+        "fpatan\n"        // Compute atan2(y, x)
+        "fstpl %0"        // Store the result
+        : "=m"(res) : "m"(y), "m"(x)
+    );
+    return res;
+}
+