sllibc/math.c

#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include "internal.h"

unsigned int _libc_rand_seed = 12345;

int rand() {
  _libc_rand_seed ^= _libc_rand_seed << 13;
  _libc_rand_seed ^= _libc_rand_seed >> 17;
  _libc_rand_seed ^= _libc_rand_seed << 5;
  return _libc_rand_seed % (RAND_MAX+1);
}

void srand(unsigned int r) {
  if (r == 0) {
    _libc_rand_seed = 1;
  } else {
    _libc_rand_seed = r;
  }
}

double fabs(double n) {
  if (n < 0) {
    return 0.0 - n;
  } else {
    return n;
  }
}

double ldexp(double x, int exp) {
    UNIMPLEMENTED();
}

double frexp(double x, int* expptr) {
  // TODO: Check special cases?
  uint64_t bits = *((uint64_t*)((void*)&x));
  uint64_t exp = (bits>>52)&0x7FF;
  *expptr = (int) exp - 1023;
  bits &= 0x800FFFFFFFFFFFFFULL;
  bits |= 1023ULL << 52;
  return *((double*)((void*)&bits));
}

/*
double _libc_factorial(int n) {
  double x = 1.0;
  for (int i = 1; i <= n; i++) {
    x *= i;
  }
  return x;
}

// This implementation is very very dumb.
double pow(double x, double p) {
  double result = x;
  for (int i = 2; i < p; i++) {
    result *= x;
  }
  return result;
}

#define EXP_EPSILON 0.0000000001

// NOTE: This function was written with the help of Qwen AI
double exp(double x) {
  double r = 1.0;
  double t = 1.0;
  int n = 1;
  
  while (fabs(t) < EXP_EPSILON) {
    t = pow(x, n) / _libc_factorial(n);
    r += t;
    n++;
  }
  
  return r;
}
*/

////////////////////////////////////////////////////////////////////
// NOTE: OLD CODE taken from PDPC400's math.c ("public domain")
////////////////////////////////////////////////////////////////////

#define HUGE_VAL 9.999999999999999999999E72

/*

  Some constants to make life easier elsewhere
  (These should I guess be in math.h)

*/
double _libc_pi   = 3.1415926535897932384626433832795;
double _libc_ln10 = 2.3025850929940456840179914546844;
double _libc_ln2 = 0.69314718055994530941723212145818;

/*

exp(x) = 1 + x + x2/2 + x3/6 + x4/24 + x5/120 + ... + xn/n! + ...

*/
double exp (double x)
{
    int i;
    double term,answer,work;

    i=2;
    term=x;
    answer=x;

    while (1)
    {
        work = i;
        term =  (term * x)/work;
        if ( answer == (answer + term) )break;
        answer = answer + (term);
        i++;
    }

    answer=answer+1.0;
    return(answer);
}

/*

   Calculate LOG using Taylor series.

   log(1+ x) = x - x**2  +  x**3  - x**4 + x**5
                   ====     ====    ====   ====    .........
                    2         3       4     8

   Note this only works for small x so we scale....

*/
double log (double x)
{
    int i,scale;
    double term,answer,work,xs;

    if (x <= 0 )
    {
        /* need to set signal */
        // TODO: Proper errno errno=EDOM;
        return (HUGE_VAL);
    }
    if( x == 1.0)return(0.0);

/*
  Scale arguments to be in range 1 < x <= 10
*/

/*
    scale = 0;
    xs = x;
    while ( xs > 10.0 ) { scale ++; xs=xs/10.0;}
    while ( xs < 1.0 ) { scale --; xs=xs*10.0;}
*/
    xs = frexp(x,&scale);
    xs = (1.0 * xs) - 1.0;
    scale = scale - 0;

    i=2;
    term=answer=xs;

    while (1)
    {
        work = i;
        term = - (term * xs);
        if ( answer == (answer + (term/work)) )break;
        answer = answer + (term/work);
        i++;
    }

    answer = answer + (double)scale * _libc_ln2;
    return(answer);
}


double log10(double x)
{
    return ( log(x) / _libc_ln10 );
}


/*

   This code uses log and exp to calculate x to the power y.
   If

*/

double pow(double x,double y)
{
    int j,neg;
    double yy,xx;
    neg=0;
    j=y;
    yy=j;
    if( yy == y) {
        xx = x;
        if ( y < 0 ){neg = 1; j = -j;}
        if ( y == 0) return (1.0);
        --j;
        while(j>0){ xx=xx * x; j--;}
        if(neg)xx=1.0/xx;
        return (xx);
    }
    if (x < 0.0)
    {
         // TODO: Proper errno errno=EDOM;
         return(0.0);
    }
    if (y == 0.0) return (1.0);

    return (exp(y*log(x)));
}