module BigMath

static VALUE
BigMath_s_exp(VALUE klass, VALUE x, VALUE vprec)
{
    ssize_t prec, n, i;
    Real* vx = NULL;
    VALUE one, d, y;
    int negative = 0;
    int infinite = 0;
    int nan = 0;
    double flo;

    prec = NUM2SSIZET(vprec);
    if (prec <= 0) {
        rb_raise(rb_eArgError, "Zero or negative precision for exp");
    }

    /* TODO: the following switch statement is almost same as one in the
     *       BigDecimalCmp function. */
    switch (TYPE(x)) {
      case T_DATA:
        if (!is_kind_of_BigDecimal(x)) break;
        vx = DATA_PTR(x);
        negative = BIGDECIMAL_NEGATIVE_P(vx);
        infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
        nan = VpIsNaN(vx);
        break;

      case T_FIXNUM:
        /* fall through */
      case T_BIGNUM:
        vx = GetVpValue(x, 0);
        break;

      case T_FLOAT:
        flo = RFLOAT_VALUE(x);
        negative = flo < 0;
        infinite = isinf(flo);
        nan = isnan(flo);
        if (!infinite && !nan) {
            vx = GetVpValueWithPrec(x, DBLE_FIG, 0);
        }
        break;

      case T_RATIONAL:
        vx = GetVpValueWithPrec(x, prec, 0);
        break;

      default:
        break;
    }
    if (infinite) {
        if (negative) {
            return ToValue(GetVpValueWithPrec(INT2FIX(0), prec, 1));
        }
        else {
            Real* vy;
            vy = VpCreateRbObject(prec, "#0");
            VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
            RB_GC_GUARD(vy->obj);
            return ToValue(vy);
        }
    }
    else if (nan) {
        Real* vy;
        vy = VpCreateRbObject(prec, "#0");
        VpSetNaN(vy);
        RB_GC_GUARD(vy->obj);
        return ToValue(vy);
    }
    else if (vx == NULL) {
        cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
    }
    x = vx->obj;

    n = prec + rmpd_double_figures();
    negative = BIGDECIMAL_NEGATIVE_P(vx);
    if (negative) {
        VALUE x_zero = INT2NUM(1);
        VALUE x_copy = f_BigDecimal(1, &x_zero, klass);
        x = BigDecimal_initialize_copy(x_copy, x);
        vx = DATA_PTR(x);
        VpSetSign(vx, 1);
    }

    one = ToValue(VpCreateRbObject(1, "1"));
    y   = one;
    d   = y;
    i   = 1;

    while (!VpIsZero((Real*)DATA_PTR(d))) {
        SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
        SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
        ssize_t m = n - vabs(ey - ed);

        rb_thread_check_ints();

        if (m <= 0) {
            break;
        }
        else if ((size_t)m < rmpd_double_figures()) {
            m = rmpd_double_figures();
        }

        d = BigDecimal_mult(d, x);                             /* d <- d * x */
        d = BigDecimal_div2(d, SSIZET2NUM(i), SSIZET2NUM(m));  /* d <- d / i */
        y = BigDecimal_add(y, d);                              /* y <- y + d  */
        ++i;                                                   /* i  <- i + 1 */
    }

    if (negative) {
        return BigDecimal_div2(one, y, vprec);
    }
    else {
        vprec = SSIZET2NUM(prec - VpExponent10(DATA_PTR(y)));
        return BigDecimal_round(1, &vprec, y);
    }

    RB_GC_GUARD(one);
    RB_GC_GUARD(x);
    RB_GC_GUARD(y);
    RB_GC_GUARD(d);
}

static VALUE
BigMath_s_log(VALUE klass, VALUE x, VALUE vprec)
{
    ssize_t prec, n, i;
    SIGNED_VALUE expo;
    Real* vx = NULL;
    VALUE vn, one, two, w, x2, y, d;
    int zero = 0;
    int negative = 0;
    int infinite = 0;
    int nan = 0;
    double flo;
    long fix;

    if (!is_integer(vprec)) {
        rb_raise(rb_eArgError, "precision must be an Integer");
    }

    prec = NUM2SSIZET(vprec);
    if (prec <= 0) {
        rb_raise(rb_eArgError, "Zero or negative precision for exp");
    }

    /* TODO: the following switch statement is almost same as one in the
     *       BigDecimalCmp function. */
    switch (TYPE(x)) {
      case T_DATA:
          if (!is_kind_of_BigDecimal(x)) break;
          vx = DATA_PTR(x);
          zero = VpIsZero(vx);
          negative = BIGDECIMAL_NEGATIVE_P(vx);
          infinite = VpIsPosInf(vx) || VpIsNegInf(vx);
          nan = VpIsNaN(vx);
          break;

      case T_FIXNUM:
        fix = FIX2LONG(x);
        zero = fix == 0;
        negative = fix < 0;
        goto get_vp_value;

      case T_BIGNUM:
        i = FIX2INT(rb_big_cmp(x, INT2FIX(0)));
        zero = i == 0;
        negative = i < 0;
get_vp_value:
        if (zero || negative) break;
        vx = GetVpValue(x, 0);
        break;

      case T_FLOAT:
        flo = RFLOAT_VALUE(x);
        zero = flo == 0;
        negative = flo < 0;
        infinite = isinf(flo);
        nan = isnan(flo);
        if (!zero && !negative && !infinite && !nan) {
            vx = GetVpValueWithPrec(x, DBLE_FIG, 1);
        }
        break;

      case T_RATIONAL:
        zero = RRATIONAL_ZERO_P(x);
        negative = RRATIONAL_NEGATIVE_P(x);
        if (zero || negative) break;
        vx = GetVpValueWithPrec(x, prec, 1);
        break;

      case T_COMPLEX:
        rb_raise(rb_eMathDomainError,
                 "Complex argument for BigMath.log");

      default:
        break;
    }
    if (infinite && !negative) {
        Real* vy;
        vy = VpCreateRbObject(prec, "#0");
        RB_GC_GUARD(vy->obj);
        VpSetInf(vy, VP_SIGN_POSITIVE_INFINITE);
        return ToValue(vy);
    }
    else if (nan) {
        Real* vy;
        vy = VpCreateRbObject(prec, "#0");
        RB_GC_GUARD(vy->obj);
        VpSetNaN(vy);
        return ToValue(vy);
    }
    else if (zero || negative) {
        rb_raise(rb_eMathDomainError,
                 "Zero or negative argument for log");
    }
    else if (vx == NULL) {
        cannot_be_coerced_into_BigDecimal(rb_eArgError, x);
    }
    x = ToValue(vx);

    RB_GC_GUARD(one) = ToValue(VpCreateRbObject(1, "1"));
    RB_GC_GUARD(two) = ToValue(VpCreateRbObject(1, "2"));

    n = prec + rmpd_double_figures();
    RB_GC_GUARD(vn) = SSIZET2NUM(n);
    expo = VpExponent10(vx);
    if (expo < 0 || expo >= 3) {
        char buf[DECIMAL_SIZE_OF_BITS(SIZEOF_VALUE * CHAR_BIT) + 4];
        snprintf(buf, sizeof(buf), "1E%"PRIdVALUE, -expo);
        x = BigDecimal_mult2(x, ToValue(VpCreateRbObject(1, buf)), vn);
    }
    else {
        expo = 0;
    }
    w = BigDecimal_sub(x, one);
    x = BigDecimal_div2(w, BigDecimal_add(x, one), vn);
    RB_GC_GUARD(x2) = BigDecimal_mult2(x, x, vn);
    RB_GC_GUARD(y)  = x;
    RB_GC_GUARD(d)  = y;
    i = 1;
    while (!VpIsZero((Real*)DATA_PTR(d))) {
        SIGNED_VALUE const ey = VpExponent10(DATA_PTR(y));
        SIGNED_VALUE const ed = VpExponent10(DATA_PTR(d));
        ssize_t m = n - vabs(ey - ed);
        if (m <= 0) {
            break;
        }
        else if ((size_t)m < rmpd_double_figures()) {
            m = rmpd_double_figures();
        }

        x = BigDecimal_mult2(x2, x, vn);
        i += 2;
        d = BigDecimal_div2(x, SSIZET2NUM(i), SSIZET2NUM(m));
        y = BigDecimal_add(y, d);
    }

    y = BigDecimal_mult(y, two);
    if (expo != 0) {
        VALUE log10, vexpo, dy;
        log10 = BigMath_s_log(klass, INT2FIX(10), vprec);
        vexpo = ToValue(GetVpValue(SSIZET2NUM(expo), 1));
        dy = BigDecimal_mult(log10, vexpo);
        y = BigDecimal_add(y, dy);
    }

    return y;
}

# File ext/bigdecimal/lib/bigdecimal/math.rb, line 228
def E(prec)
  raise ArgumentError, "Zero or negative precision for E" if prec <= 0
  BigMath.exp(1, prec)
end

# File ext/bigdecimal/lib/bigdecimal/math.rb, line 183
def PI(prec)
  raise ArgumentError, "Zero or negative precision for PI" if prec <= 0
  n      = prec + BigDecimal.double_fig
  zero   = BigDecimal("0")
  one    = BigDecimal("1")
  two    = BigDecimal("2")

  m25    = BigDecimal("-0.04")
  m57121 = BigDecimal("-57121")

  pi     = zero

  d = one
  k = one
  t = BigDecimal("-80")
  while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    t   = t*m25
    d   = t.div(k,m)
    k   = k+two
    pi  = pi + d
  end

  d = one
  k = one
  t = BigDecimal("956")
  while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    t   = t.div(m57121,n)
    d   = t.div(k,m)
    pi  = pi + d
    k   = k+two
  end
  pi
end

# File ext/bigdecimal/lib/bigdecimal/math.rb, line 146
def atan(x, prec)
  raise ArgumentError, "Zero or negative precision for atan" if prec <= 0
  return BigDecimal("NaN") if x.nan?
  pi = PI(prec)
  x = -x if neg = x < 0
  return pi.div(neg ? -2 : 2, prec) if x.infinite?
  return pi / (neg ? -4 : 4) if x.round(prec) == 1
  x = BigDecimal("1").div(x, prec) if inv = x > 1
  x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5
  n    = prec + BigDecimal.double_fig
  y = x
  d = y
  t = x
  r = BigDecimal("3")
  x2 = x.mult(x,n)
  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    t = -t.mult(x2,n)
    d = t.div(r,m)
    y += d
    r += 2
  end
  y *= 2 if dbl
  y = pi / 2 - y if inv
  y = -y if neg
  y
end

# File ext/bigdecimal/lib/bigdecimal/math.rb, line 102
def cos(x, prec)
  raise ArgumentError, "Zero or negative precision for cos" if prec <= 0
  return BigDecimal("NaN") if x.infinite? || x.nan?
  n    = prec + BigDecimal.double_fig
  one  = BigDecimal("1")
  two  = BigDecimal("2")
  x = -x if x < 0
  if x > (twopi = two * BigMath.PI(prec))
    if x > 30
      x %= twopi
    else
      x -= twopi while x > twopi
    end
  end
  x1 = one
  x2 = x.mult(x,n)
  sign = 1
  y = one
  d = y
  i = BigDecimal("0")
  z = one
  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    sign = -sign
    x1  = x2.mult(x1,n)
    i  += two
    z  *= (i-one) * i
    d   = sign * x1.div(z,m)
    y  += d
  end
  y
end

# File ext/bigdecimal/lib/bigdecimal/math.rb, line 58
def sin(x, prec)
  raise ArgumentError, "Zero or negative precision for sin" if prec <= 0
  return BigDecimal("NaN") if x.infinite? || x.nan?
  n    = prec + BigDecimal.double_fig
  one  = BigDecimal("1")
  two  = BigDecimal("2")
  x = -x if neg = x < 0
  if x > (twopi = two * BigMath.PI(prec))
    if x > 30
      x %= twopi
    else
      x -= twopi while x > twopi
    end
  end
  x1   = x
  x2   = x.mult(x,n)
  sign = 1
  y    = x
  d    = y
  i    = one
  z    = one
  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    sign = -sign
    x1  = x2.mult(x1,n)
    i  += two
    z  *= (i-one) * i
    d   = sign * x1.div(z,m)
    y  += d
  end
  neg ? -y : y
end

# File ext/bigdecimal/lib/bigdecimal/math.rb, line 43
def sqrt(x, prec)
  x.sqrt(prec)
end