mmath

#include "erl_nif.h"
#include "mmath.h"


ffloat
float_mul(decimal v, double m) {
  // TODO fix in situation when m is so big, that it will overflow coefficient
  return (ffloat) {
    .value = v.value * m,
      .confidence = v.confidence
      };
}

decimal
dec_div(decimal v, int64_t m) {
  // TODO: Improve accuracy by adjusting coefficient before dividing
  dec_inflate(&v);
  decimal r = {
    .coefficient = v.coefficient / m,
    .exponent = v.exponent,
    .confidence = v.confidence
  };
  dec_reduce(&r);
  return r;
}

static decimal
dec_add_aligned(int64_t a_coef, int64_t b_coef, int8_t e, uint32_t confidence) {
  decimal r = {
    .coefficient = a_coef + b_coef,
    .exponent = e,
    .confidence = confidence
  };
  dec_reduce(&r);
  return r;
}

/* Add decimal b to a. Decimal a has always bigger coefficient */
static decimal
dec_add_not_aligned(decimal big, decimal small, uint32_t confidence) {
  int8_t over_digits = 0;
  int8_t e, de;
  int64_t a_coef, b_coef;

  if (big.coefficient == 0) {
    small.confidence = confidence;
    return small;
  }
  if (small.coefficient == 0) {
    big.confidence = confidence;
    return big;
  }

  e = small.exponent;
  de = big.exponent - e;
  b_coef = small.coefficient;
  over_digits = qlog10(llabs(big.coefficient)) + de - MAX_DIGITS - 1;

  if (over_digits > 0) {
    e += over_digits;
    de -= over_digits;
    b_coef /= (int64_t)qipow10(over_digits);
  }

  a_coef = big.coefficient * qipow10(de);
  return dec_add_aligned(a_coef, b_coef, e, confidence);
}

inline decimal
dec_add(decimal a, decimal b) {
  uint32_t confidence = (a.confidence + b.confidence) / 2;
  if (a.exponent == b.exponent) {
    return dec_add_aligned(a.coefficient, b.coefficient, a.exponent,
                           confidence);
  } else {
    if (a.exponent >= b.exponent) {
      return dec_add_not_aligned(a, b, confidence);
    } else {
      return dec_add_not_aligned(b, a, confidence);
    }
  }
}

inline decimal
dec_add3(decimal a, decimal b, decimal c) {
  return dec_add(dec_add(a, b), c);
}

inline decimal /* a - b */
dec_sub(decimal a, decimal b) {
  return dec_add(a, dec_neg(b));
}

inline decimal
dec_neg(decimal a) {
  return (decimal) {
    .coefficient = -a.coefficient,
      .exponent = a.exponent,
      .confidence = a.confidence
  };
}

/* return 1 when 1 > b, 0 when equal and -1 otherwise */
int
dec_cmp(decimal a, decimal b) {
  //TODO: optimise, not going into actual computation
  decimal r = dec_sub(a, b);
  if (r.coefficient == 0 && r.exponent == 0) {
    return 0;
  } else if (r.coefficient < 0) {
    return -1;
  } else {
    return 1;
  }
}