math.asa

// ------------ math.asa ------------
// Contains functions for doing math
import "logical.asa";

// precision used in math.asa
def precision: push 0.000000000000001f; end

// round tolerance if that close to integer => round
//   => not used by default! 
def round_tolerance: push 0.0000000001f; end

// mathematical constant pi
def pi: push 3.141592653589793f; end

// mathematical constant e (euler's number)
def e:  push 2.718281828459045f; end

// modulo: computes the smallest positive representation of a mod b
def modulo: pop b; pop a;
  push a; push 0; cmp; pop a_to_0;

  push a_to_0; ifgoto 0 rem_0;        // if a == 0, jump to rem_0
  push a_to_0; ifgoto 1 a_greater_0;  // if a  > 0, jump to a_greater_0
  // else: here means a < 0 (negative)
  // => add b to a until a >= 0
  label a_lesser_0;
    push a; push b; add; pop a;
    push a; push 0; cmp; pop a_to_0;
  push a_to_0; ifgoto -1 a_lesser_0;

  // If a > 0, jump to a_greater_0
  push a; push 0; cmp; pop a_to_0;
  push a_to_0; ifgoto 1 a_greater_0;

  // If a == 0
  // or a == b
  // push 0 to rem and return
  label rem_0;
    push 0; pop rem; goto return;

  // If a > 0, compare a with b
  label a_greater_0;
    push a; push b; cmp; pop a_to_b;
    push a_to_b; ifgoto -1 a_lesser_b; // if a  < b, jump to a_lesser_b
    push a_to_b; ifgoto 0 rem_0;       // if a == b, jump to rem_0

    // If a > b, subtract b from a until a < b
    label a_greater_b;
      push a; push b; sub; pop a;
      push a; push b; cmp; pop a_to_b;
    push a_to_b; ifgoto 1 a_greater_b;
    push a_to_b; ifgoto -1 a_lesser_b;
    goto rem_0;                        // because at this point a == b

  // If a < b, set remainder to a and return
  label a_lesser_b;
    push a; pop rem; goto return;

  label return;
    push rem;
end

// max: returns the maximum of two values
def max: pop a; pop b;
  push a; push b; call logical/lesser_or_eq?; pop a_leq_b;
  push a_leq_b; call logical/not; pop b_gt_a;
  push a; push b_gt_a; mul;  // a * (b > a)
  push b; push a_leq_b; mul; // b * (a <= b)
  add;                       // (a * (b > a)) + (b * (a <= b))
end

// min: returns the minimum of two values
def min: pop a; pop b;
  push a; push b; call logical/lesser_or_eq?; pop a_leq_b;
  push a_leq_b; call logical/not; pop b_gt_a;
  push a; push a_leq_b; mul;  // a * (a <= b)
  push b; push b_gt_a; mul;   // b * (b > a)
  add;                        // (a * (a <= b)) + (b * (b > a))
end

// abs: pushes the absolute value of the element on top of the stack
def abs: pop n;
  push n; push 0; cmp; ifgoto -1 negative;        // if n < 0 return -1 * n
  push n; pop res; goto return;                   // else     return n
  label negative; push n; push -1; mul; pop res;
  label return; push res;
end

// exp: computes e^x, where x => 0
def exp: pop x; call math/e; push x; call math/ipow; end

// ln: computes natural log of x for x > 0 using the taylor sequence
//   => more info on the taylor series here: https://en.wikipedia.org/wiki/Taylor_series
def ln: pop x;
  push x; push 0; call logical/greater?; ifgoto 1 compute;
  push "math/ln: input x must be greater than 0!"; raise;

  label compute;
  push x; push 1; call logical/equals?;      ifgoto 1 ret0; // if x == 1 return 0
  push x; call math/e; call logical/equals?; ifgoto 1 ret1; // if x == e return 1
  push 0.693147180559945bd; pop ln2;                        // yes.. I actually saved it in a variable :)
  push x; push 2; call logical/equals?; ifgoto 1 retln2;    // if x == 2 return ln(2)
  // else if x > 2:
  push x; push 0.0bd; add; pop x;                          // convert x to BigDouble

  // simplify x so that 1 <= x <= 2:
  //   so that we can later use logarithmic identity:
  //   ln(x) = ln(x/2) + ln(2)
  //   example:
  //    ln(10)	= ln(5)    + ln(2)
  //            = ln(2.5)  + ln(2) + ln(2)
  //            = ln(1.25) + ln(2) + ln(2) + ln(2) = 2.30259...
  //    and then we compute ln(1.25) and add the ln(2)'s
  //    to get the final result
  // push 0; pop n;
  var n 0;
  label simplify;
    incr n;
    push x; push 2; div; pop x;            // x = x/2
    push x; push 2; call logical/lesser?;
    ifgoto 1 compute_taylor;               // if x < 2 goto compute_taylor
    goto simplify;                         // else repeat the loop

  // taylor series for ln(x):
  // (x-1)^1 - 1/2 * (x-1)^2 + 1/3 * (x-1)^3 - ...
  label compute_taylor;
    var k 1.0bd;                               // k is used for alternating sign
    var i 1.0bd;                               // current iteration counter (BigDouble for coefficient computation)
    var res 0.0bd;                             // result init to 0.0bd (BigDouble)
    push n; push ln2; mul; pop res;            // res = ln(2) * n
    var max_iter 100;                           // max amount of iterations to do
    var term 1.0bd;                            // term = 1.0bd
    var subterm 1.0bd;                         // subterm = 1.0bd
    label loop;
      push k; push i; div; pop coeff;                   // coeff = k/i
      push x; push 1; sub; push i; call math/ipow;
      pop subterm;                                      // subterm sb = (x-1)^i
      push coeff; push subterm; mul; pop term;          // term = coeff * sb

      // if term < precision return current res
      push term; call math/abs; call math/precision;
      call logical/lesser?;
      ifgoto 1 return;

      push res; push term; add; pop res;           // res = res + term
      incr i;
      push k; push -1; mul; pop k;                 // k = -k
      push i; push max_iter; call logical/lesser?;
      ifgoto 1 loop;                               // if i < max_iter repeat loop
  goto return;

  label ret0;    push 0; pop res; goto return;
  label ret1;    push 1; pop res; goto return;
  label retln2;  push ln2;  pop res; goto return;
  label return;  push res;
end

// even?: pushes 1 to the stack if the number is even
def even?: pop n; 
  push n; push 2; call math/modulo; push 0; call logical/equals?;
end

// factorial: calculates the nth factorial
def factorial: pop n;
    var ans 1bi; 
    var i 1;
    label fact;
      push i; push n; call logical/greater?;
      ifgoto 1 return;                        // if i > n return ans
      push ans; push i; mul; pop ans;         // ans = ans * i
      incr i;                                 // increment i
      goto fact;
    label return; push ans;
end

// binomial_coeff: calculates a choose b <=> a! / ((a-b)! * b!)
def binomial_coeff: pop b; pop a;
  push a; push b; call logical/equals?;
  ifgoto 1 ret1;                             // if a == b return 1
  // else...
  push a; call math/factorial;               // a!
  push a; push b; sub; call math/factorial;  // (a-b)!
  push b; call math/factorial;               // b!
  mul;                                       // (a-b)! * b!
  div;                                       // a! / ((a-b)! * b!)
  pop res; goto return;

  label ret1;   var res 1;
  label return; push res;
end

// ipow: computes base^n where n is an integer
//   => If you have fractional exponents, compute result in
//      combination with iroot_nth function, since
//      b^(n/m) <=> iroot(ipow(b, n), m)
def ipow: pop n; pop base;
  push n; push 1; call logical/equals?; ifgoto 1 retbase; // if n == 1 return base
  push n; push 0; call logical/lesser?; pop n_lt_0;       // if n < 0 then n_lt_0 = 1 else 0

  push n; call math/abs; pop n_abs;
  var res 1.0f;                                                   // initialize result at 1 for x == 0
  push n_abs; push 0; call logical/equals?; ifgoto 1 ret;         // if x == 0 return 1
  var i 1;
  label loop;
    push i; push n_abs; call logical/lesser_or_eq?; ifgoto 0 done; // if i > n goto return
    push res; push base; mul; pop res;                             // res = res * base
    incr i;                                                        // increment i
    goto loop;                                                     // repeat...

  label done;
    push n_lt_0; ifgoto 0 ret;                 // if n > 0 goto return
    push 1; push res; div; pop res; goto ret;  // else: res = 1/res

  label retbase; push base; pop res; goto ret;
  label ret;     push res;
end

// iroot: computes the nth root of base
//   => using newton raphson method to approximate integer root
// def iroot:
//   pop n; pop base;
//   var corr 0.0bd; // correction
//   var res 1.0bd;  // res
  
//   push base; push 0.0bd; call logical/equals?; 
//   ifgoto 1 ret0;
  
//   push n; push 1; call logical/lesser?; 
//   ifgoto 1 return_nan;
  
//   push base; push 0.0bd; call logical/lesser?; ifgoto 1 check_even;
//   goto compute;

//   label check_even;
//     push n; call math/even?; ifgoto 1 return_nan;

//   label compute;
//     label loop;
//       push base; push res; push n; push 1; sub; call math/ipow; div; push res; sub; push n; div; pop corr;
//       push res; push corr; add; pop res;
//       push corr; call math/abs; call math/precision; call logical/greater_or_eq?; ifgoto 1 loop;
  
//   label ret0;
//     push 0.0bd; goto ret;

//   label return_nan;
//     push "math/iroot: exponent n must be integer and n >= 1"; raise;

//   label ret;
//     push res;
// end


// iroot: computes the nth root of base using the Newton-Raphson method
def iroot: pop n; pop base;
  var corr 0.0bd; // Correction term for Newton-Raphson update
  var res 1.0bd;  // Initial approximation of the nth root
  
  // If base is 0, return 0
  push base; push 0.0bd; call logical/equals?; 
  ifgoto 1 ret_0;
  
  // If n is less than 1, raise an error
  push n; push 1; call logical/lesser?; 
  ifgoto 1 ret_nan;
  
  // If base is less than 0 and n is even, raise an error (no real root)
  push base; push 0.0bd; call logical/lesser?; ifgoto 1 check_even;
  goto compute;

  label check_even;
    push n; call math/even?; ifgoto 1 ret_no_real_root;

  // Compute the nth root using the Newton-Raphson method
  label compute;
    label loop;
      // corr = ((base / (res^(n-1))) - res) / n
      push base; 
      push res; push n; push 1; sub; call math/ipow; div;
      push res; sub; push n; div; pop corr;

      // res = res + corr
      push res; push corr; add; pop res;

      // continue until corr < precision, else return res
      push corr; call math/abs; call math/precision; 
      call logical/greater_or_eq?; ifgoto 1 loop; goto ret;
  
  label ret_0; push 0.0bd; pop res; goto ret;

  label ret_nan; 
    push "math/iroot: cannot compute nth root because n must be integer and n >= 1"; raise;
    
  label ret_no_real_root; 
    push "math/iroot: cannot compute nth root because b < 0 and n is even => no real root"; raise;

  label ret; push res;
end