Implementation of some math.h functions

math.c

// uses only elementary methods
// guarantees none of speed, precision or accuracy
// I didn't read the specs
// does not handle any errors
// might be educational; definitely not practical

#include "math.h"

// double has 53 bits mantissa
// precision is log_10(2^53) = 15.9546
#define PI 3.141592653589793
#define E 2.718281828459045
#define LN_2 0.6931471805599453
#define SQRT_2 1.414213562373095

#define NAN (0.0 / 0.0)
#define POS_INF (1.0 / 0.0)
#define NEG_INF (-1.0 / 0.0)

double pow_int(double x, int n);
double log_any(double x, double y);

double acos(double x) {
  if (x < 0) {
    return PI - acos(-x);
  }
  if (fabs(x) > 1) {
    return NAN;
  }
  if (x == 0) {
    return PI / 2;
  }
  if (x == -1) {
    return PI;
  }
  if (x == 1) {
    return 0;
  }
  return atan(sqrt(1 - x * x) / x);
}

double asin(double x) {
  if (fabs(x) > 1) {
    return NAN;
  }
  if (x == -1) {
    return -PI / 2;
  }
  if (x == 1) {
    return PI / 2;
  }
  return atan(x / sqrt(1 - x * x));
}

// let tan(a) = x
// arctan((sqrt(1 + x^2) - 1) / x)
// = arctan((sec(a) - 1) / tan(a) * (cos(a) / cos(a)))
// = arctan((1 - cos(a)) / sin(a))
// = arctan((2 * (sin(a/2))^2) / (2 * sin(a/2) * cos(a / 2)))
// = arctan(sin(a/2) / cos(a/2)) = arctan(tan(a/2))
// = a/2 = 1/2 * arctan(x)
#define EPSILON 0.1
double atan(double x) {
  double s = 1;
  while (fabs(x) > EPSILON) {
    x = (sqrt(1 + x * x) - 1) / x;
    s *= 2;
  }
  double approx = x - pow_int(x, 3) / 3 + pow_int(x, 5) / 5 -
                  pow_int(x, 7) / 7 + pow_int(x, 9) / 9 - pow_int(x, 11) / 11 +
                  pow_int(x, 13) / 13 - pow_int(x, 15) / 15 +
                  pow_int(x, 17) / 17 - pow_int(x, 19) / 19;
  return s * approx;
}
#undef EPSILON

double atan2(double y, double x) {
  if (x == 0) {
    if (y == 0) {
      return 0;
    }
    if (y > 0) {
      return PI / 2;
    } else {
      return -PI / 2;
    }
  }
  if (x > 0) {
    return atan(y / x);
  } else {
    if (y >= 0) {
      return atan(y / x) + PI;
    } else {
      return atan(y / x) - PI;
    }
  }
  return 0;
}

double cos(double x) {
  int sign = 1;
  // cos(-a) = cos(a)
  // a in [0, \infinity)
  if (x < 0) {
    x = -x;
  }
  // cos(a + 2k * PI) = cos(a)
  // a in [0, 2 * PI]
  if (x > 2 * PI) {
    x -= (int)(x / 2 / PI) * 2 * PI;
  }
  // cos(PI + a) = -cos(a)
  // a in [0, PI]
  if (x > PI) {
    x -= PI;
    sign *= -1;
  }
  // cos(PI - a) = -cos(a)
  // a in [0, PI / 2]
  if (x > PI / 2) {
    x = PI - x;
    sign *= -1;
  }
  // sin(PI / 2 - a) = cos(a)
  // a in [0, PI / 4]
  if (x > PI / 4) {
    return sign * sin(PI / 2 - x);
  }
  double approx = 1 - pow_int(x, 2) / 2 + pow_int(x, 4) / 24 -
                  pow_int(x, 6) / 720 + pow_int(x, 8) / 40320 -
                  pow_int(x, 10) / 3628800 + pow_int(x, 12) / 479001600 -
                  pow_int(x, 14) / 87178291200;
  return sign * approx;
}

double cosh(double x) { return (exp(x) + exp(-x)) / 2; }

double sin(double x) {
  int sign = 1;
  // sin(-a) = -sin(a)
  // a in [0, \infinity)
  if (x < 0) {
    x = -x;
    sign *= -1;
  }
  // sin(a + 2k * PI) = sin(a)
  // a in [0, 2 * PI]
  if (x > 2 * PI) {
    x -= (int)(x / 2 / PI) * 2 * PI;
  }
  // sin(PI + a) = -sin(a)
  // a in [0, PI]
  if (x > PI) {
    x -= PI;
    sign *= -1;
  }
  // sin(PI - a) = sin(a)
  // a in [0, PI / 2]
  if (x > PI / 2) {
    x = PI - x;
  }
  // cos(PI / 2 - a) = sin(a)
  // a in [0, PI / 4]
  if (x > PI / 4) {
    return sign * cos(PI / 2 - x);
  }
  double approx = x - pow_int(x, 3) / 6 + pow_int(x, 5) / 120 -
                  pow_int(x, 7) / 5040 + pow_int(x, 9) / 362880 -
                  pow_int(x, 11) / 39916800 + pow_int(x, 13) / 6227020800;
  return sign * approx;
}

double sinh(double x) { return (exp(x) - exp(-x)) / 2; }

double tan(double x) { return sin(x) / cos(x); }

double tanh(double x) {
  double base = exp(x);
  return 1 - 2 / (base * base + 1);
}

// e^x = 2^(x/ln(2)) = 2^(int exponent [n]) + 2^(fractional exponent [x])
double exp(double x) {
  if (x == 0) {
    return 1;
  }
  if (x < 0) {
    return 1 / exp(-x);
  }
  x = x / LN_2;
  int n = (int)x;
  x -= (n + 0.5);
  // powers of ln(2)
  double ln2[13] = {1,
                    LN_2,
                    0.480453013918201,
                    0.333024651988929,
                    0.230835098583083,
                    0.160002697757141,
                    0.110905418832348,
                    0.0768737783724616,
                    0.0532848427378619,
                    0.0369342385103290,
                    0.0256008632895631,
                    0.0177451662090613,
                    0.0123000119263784};
  // Taylor approximate 2^x near 0.5
  // at least power of 2 is fast
  double frac_result = SQRT_2 + (SQRT_2 * x * LN_2) +
                       (pow_int(x, 2) * ln2[2] / SQRT_2) +
                       (pow_int(x, 3) * ln2[3] / (3 * SQRT_2)) +
                       (pow_int(x, 4) * ln2[4] / (12 * SQRT_2)) +
                       (pow_int(x, 5) * ln2[5] / (60 * SQRT_2)) +
                       (pow_int(x, 6) * ln2[6] / (360 * SQRT_2)) +
                       (pow_int(x, 7) * ln2[7] / (2520 * SQRT_2)) +
                       (pow_int(x, 8) * ln2[8] / (20160 * SQRT_2)) +
                       (pow_int(x, 9) * ln2[9] / (181440 * SQRT_2)) +
                       (pow_int(x, 10) * ln2[10] / (1814400 * SQRT_2)) +
                       (pow_int(x, 11) * ln2[11] / (19958400 * SQRT_2)) +
                       (pow_int(x, 12) * ln2[12] / (239500800 * SQRT_2));
  return (1 << n) * frac_result;
}

double frexp(double x, int *exponent) {
  if (x < 0) {
    return -frexp(-x, exponent);
  }
  int n = 1;
  double a = x;
  while (a >= 1.0) {
    a = x / (1 << n++);
  }
  *exponent = n - 1;
  return a;
}

double ldexp(double x, int exponent) { return x * (1 << exponent); }

double ln(double x) {
  if (x < 0) {
    return NAN;
  }
  if (x == 0) {
    return NEG_INF;
  }
  if (x < 1) {
    return -ln(1 / x);
  }
  // factor out a power of 2
  // such that ln(a * 2^n) = ln(a) + n * ln(2)
  // https://math.stackexchange.com/a/3383716
  int n = 1;
  double a = x;
  while (a > 3.0 / 4.0) {
    a = x / (1 << n++);
  }
  a *= 2;
  n -= 2;
  // expand at (1 + a) / 2
  // 1/x = \sum_{k=0}^{\infty} (-1)^k * (2/(a+1))^{k+1} * (x-(1+a)/(2))^k
  // \ln(a) = \int_{1}^{a} 1/x dx = \sum_{k=0}^{\infty}
  // (-1)^k * (2/(a+1))^{k+1} * \int_{1}^{a}(x-(1+a)/2)^k dx
  double base = (a - 1) / (a + 1);
  double ln_a = (2 * base) + (2 * pow_int(base, 3) / 3) +
                (2 * pow_int(base, 5) / 5) + (2 * pow_int(base, 7) / 7) +
                (2 * pow_int(base, 9) / 9) + (2 * pow_int(base, 11) / 11) +
                (2 * pow_int(base, 13) / 13) + (2 * pow_int(base, 15) / 15) +
                (2 * pow_int(base, 17) / 17) + (2 * pow_int(base, 19) / 19);
  return ln_a + n * LN_2;
}

double log(double x) { return ln(x); }

double log10(double x) { return log_any(10, x); }

double log_any(double x, double y) { return ln(y) / ln(x); }

double modf(double x, double *integer) {
  *integer = (int)x;
  return x - (int)x;
}

double pow_int_rec(double a, double x, int n) {
  if (n < 0) {
    return pow_int_rec(a, 1 / x, -n);
  }
  if (n == 0) {
    return a;
  }
  if (n & 1) {
    return pow_int_rec(x * a, x * x, (n - 1) / 2);
  } else {
    return pow_int_rec(a, x * x, n / 2);
  }
}

// used by Taylor series
double pow_int(double x, int n) {
  return pow_int_rec(1, x, n);
}

// a^x = e^ln(a)^x = e^(xln(a))
double pow(double x, double n) {
  if (n == 0) {
    return 1;
  }
  return exp(n * ln(x));
}

// Newton's (Babylonian) method
// x_n+1 = x_n - f(x_n) / f'(x_n)
// let f(x) = x^2 - S = 0
// x_n+1 = x_n - (x_n^2 - S) / (2 * x_n)
//       = 1/2 * (x_n + S/x_n)
// intuition is that if estimate e is low,
// S / e is high & correct is in between
double sqrt(double x) {
  // keep floating point precision
  double s = 1;
  while (x > 100) {
    x /= 100;
    s *= 10;
  }
  int n = 10;
  double res = x / 2;
  while (n--) {
    res = (res + x / res) / 2;
  }
  return s * res;
}

double ceil(double x) {
  if (x == (int)x) {
    return x;
  }
  return x > 0 ? (int)x + 1 : (int)x;
}

double fabs(double x) { return (x > 0 ? x : -x); }

double floor(double x) {
  if (x == (int)x) {
    return x;
  }
  return x > 0 ? (int)x : (int)x - 1;
}

double fmod(double x, double y) { return x - (int)(x / y) * y; }

math.h

// https://www.tutorialspoint.com/c_standard_library/math_h.htm

double acos(double x);
// Returns the arc cosine of x in radians.

double asin(double x);
// Returns the arc sine of x in radians.

double atan(double x);
// Returns the arc tangent of x in radians.

double atan2(double y, double x);
// Returns the arc tangent in radians of y/x based on the signs of both values
// to determine the correct quadrant.

double cos(double x);
// Returns the cosine of a radian angle x.

double cosh(double x);
// Returns the hyperbolic cosine of x.

double sin(double x);
// Returns the sine of a radian angle x.

double sinh(double x);
// Returns the hyperbolic sine of x.

double tan(double x);
// Returns the tangent of a radian angle x.

double tanh(double x);
// Returns the hyperbolic tangent of x.

double exp(double x);
// Returns the value of e raised to the xth power.

double frexp(double x, int *exponent);
// The returned value is the mantissa and the integer pointed to by exponent is
// the exponent. The resultant value is x = mantissa * 2 ^ exponent.

double ldexp(double x, int exponent);
// Returns x multiplied by 2 raised to the power of exponent.

double log(double x);
// Returns the natural logarithm (base-e logarithm) of x.

double log10(double x);
// Returns the common logarithm (base-10 logarithm) of x.

double modf(double x, double *integer);
// The returned value is the fraction component (part after the decimal), and
// sets integer to the integer component.

double pow(double x, double n);
// Returns x raised to the power of y.

double sqrt(double x);
// Returns the square root of x.

double ceil(double x);
// Returns the smallest integer value greater than or equal to x.

double fabs(double x);
// Returns the absolute value of x.

double floor(double x);
// Returns the largest integer value less than or equal to x.

double fmod(double x, double y);
// Returns the remainder of x divided by y.