module BigMath

def self.exp(x, prec)
  BigDecimal::Internal.validate_prec(prec, :exp)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :exp)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return x.positive? ? BigDecimal::Internal.infinity_computation_result : BigDecimal(0) if x.infinite?
  return BigDecimal(1) if x.zero?
  return BigDecimal(1).div(exp(-x, prec), prec) if x < 0
  # exp(x * 10**cnt) = exp(x)**(10**cnt)
  cnt = x > 1 ? x.exponent : 0
  prec2 = prec + BigDecimal.double_fig + cnt
  x = x._decimal_shift(-cnt)
  xn = BigDecimal(1)
  y = BigDecimal(1)
  # Taylor series for exp(x) around 0
  1.step do |i|
    n = prec2 + xn.exponent
    break if n <= 0 || xn.zero?
    x = x.mult(1, n)
    xn = xn.mult(x, n).div(i, n)
    y = y.add(xn, prec2)
  end
  # calculate exp(x * 10**cnt) from exp(x)
  # exp(x * 10**k) = exp(x * 10**(k - 1)) ** 10
  cnt.times do
    y2 = y.mult(y, prec2)
    y5 = y2.mult(y2, prec2).mult(y, prec2)
    y = y5.mult(y5, prec2)
  end
  y.mult(1, prec)
end

def self.log(x, prec)
  BigDecimal::Internal.validate_prec(prec, :log)
  raise Math::DomainError, 'Complex argument for BigMath.log' if Complex === x
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  raise Math::DomainError, 'Zero or negative argument for log' if x <= 0
  return BigDecimal::Internal.infinity_computation_result if x.infinite?
  return BigDecimal(0) if x == 1
  BigDecimal.save_limit do
    BigDecimal.limit(0)
    if x > 10 || x < 0.1
      log10 = log(BigDecimal(10), prec)
      exponent = x.exponent
      x = x._decimal_shift(-exponent)
      if x < 0.3
        x *= 10
        exponent -= 1
      end
      return log10 * exponent + log(x, prec)
    end
    x_minus_one_exponent = (x - 1).exponent
    prec += BigDecimal.double_fig
    # log(x) = log(sqrt(sqrt(sqrt(sqrt(x))))) * 2**sqrt_steps
    sqrt_steps = [Integer.sqrt(prec) + 3 * x_minus_one_exponent, 0].max
    lg2 = 0.3010299956639812
    prec2 = prec + [-x_minus_one_exponent, 0].max + (sqrt_steps * lg2).ceil
    sqrt_steps.times do
      x = x.sqrt(prec2)
      # Workaround for https://github.com/ruby/bigdecimal/issues/354
      x = x.mult(1, prec2 + BigDecimal.double_fig)
    end
    # Taylor series for log(x) around 1
    # log(x) = -log((1 + X) / (1 - X)) where X = (x - 1) / (x + 1)
    # log(x) = 2 * (X + X**3 / 3 + X**5 / 5 + X**7 / 7 + ...)
    x = (x - 1).div(x + 1, prec2)
    y = x
    x2 = x.mult(x, prec)
    1.step do |i|
      n = prec + x.exponent - y.exponent + x2.exponent
      break if n <= 0 || x.zero?
      x = x.mult(x2.round(n - x2.exponent), n)
      y = y.add(x.div(2 * i + 1, n), prec)
    end
    y.mult(2 ** (sqrt_steps + 1), prec)
  end
end

def E(prec)
  raise ArgumentError, "Zero or negative precision for E" if prec <= 0
  BigMath.exp(1, prec)
end

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

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
  n = prec + BigDecimal.double_fig
  x = BigDecimal("1").div(x, n) if inv = x > 1
  x = (-1 + sqrt(1 + x.mult(x, n), n)).div(x, n) if dbl = x > 0.5
  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

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 > 6
    twopi = two * BigMath.PI(prec + x.exponent)
    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 < -1 ? BigDecimal("-1") : y > 1 ? BigDecimal("1") : y
end

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 > 6
    twopi = two * BigMath.PI(prec + x.exponent)
    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
  y = BigDecimal("1") if y > 1
  neg ? -y : y
end

def sqrt(x, prec)
  x.sqrt(prec)
end