Module: BigMath

# File 'lib/bigdecimal/math.rb', line 54

private_class_method def _sin_periodic_reduction(x, prec, add_half_pi: false) # :nodoc:
  return [1, x] if -Math::PI/2 <= x && x <= Math::PI/2 && !add_half_pi

  mod_prec = prec + BigDecimal.double_fig
  pi_extra_prec = [x.exponent, 0].max + BigDecimal.double_fig
  while true
    pi = PI(mod_prec + pi_extra_prec)
    half_pi = pi / 2
    div, mod = (add_half_pi ? x + pi : x + half_pi).divmod(pi)
    mod -= half_pi
    if mod.zero? || mod_prec + mod.exponent <= 0
      # mod is too small to estimate required pi precision
      mod_prec = mod_prec * 3 / 2 + BigDecimal.double_fig
    elsif mod_prec + mod.exponent < prec
      # Estimate required precision of pi
      mod_prec = prec - mod.exponent + BigDecimal.double_fig
    else
      return [div % 2 == 0 ? 1 : -1, mod.mult(1, prec)]
    end
  end
end

# File 'lib/bigdecimal/math.rb', line 162

def atan(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :atan)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  n = prec + BigDecimal.double_fig
  pi = PI(n)
  x = -x if neg = x < 0
  return pi.div(neg ? -2 : 2, prec) if x.infinite?
  return pi.div(neg ? -4 : 4, prec) if x.round(prec) == 1
  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.mult(1, prec)
end

# File 'lib/bigdecimal/math.rb', line 124

def cos(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :cos)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cos)
  return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan?
  sign, x = _sin_periodic_reduction(x, prec + BigDecimal.double_fig, add_half_pi: true)
  sign * sin(x, prec)
end

# File 'lib/bigdecimal/math.rb', line 245

def E(prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :E)
  BigMath.exp(1, prec)
end

# File 'lib/bigdecimal.rb', line 328

def self.exp(x, prec)
  prec = BigDecimal::Internal.coerce_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?

  # exp(x * 10**cnt) = exp(x)**(10**cnt)
  cnt = x < -1 || x > 1 ? x.exponent : 0
  prec2 = prec + BigDecimal.double_fig + cnt
  x = x._decimal_shift(-cnt)

  # Calculation of exp(small_prec) is fast because calculation of x**n is fast
  # Calculation of exp(small_abs) converges fast.
  # exp(x) = exp(small_prec_part + small_abs_part) = exp(small_prec_part) * exp(small_abs_part)
  x_small_prec = x.round(Integer.sqrt(prec2))
  y = _exp_taylor(x_small_prec, prec2).mult(_exp_taylor(x.sub(x_small_prec, prec2), prec2), prec2)

  # 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

# File 'lib/bigdecimal.rb', line 251

def self.log(x, prec)
  prec = BigDecimal::Internal.coerce_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, 'Negative argument for log' if x < 0
  return -BigDecimal::Internal.infinity_computation_result if x.zero?
  return BigDecimal::Internal.infinity_computation_result if x.infinite?
  return BigDecimal(0) if x == 1

  prec2 = prec + BigDecimal.double_fig
  BigDecimal.save_limit do
    BigDecimal.limit(0)
    if x > 10 || x < 0.1
      log10 = log(BigDecimal(10), prec2)
      exponent = x.exponent
      x = x._decimal_shift(-exponent)
      if x < 0.3
        x *= 10
        exponent -= 1
      end
      return (log10 * exponent).add(log(x, prec2), prec)
    end

    x_minus_one_exponent = (x - 1).exponent

    # log(x) = log(sqrt(sqrt(sqrt(sqrt(x))))) * 2**sqrt_steps
    sqrt_steps = [Integer.sqrt(prec2) + 3 * x_minus_one_exponent, 0].max

    lg2 = 0.3010299956639812
    sqrt_prec = prec2 + [-x_minus_one_exponent, 0].max + (sqrt_steps * lg2).ceil

    sqrt_steps.times do
      x = x.sqrt(sqrt_prec)
    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, sqrt_prec)
    y = x
    x2 = x.mult(x, prec2)
    1.step do |i|
      n = prec2 + 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), prec2)
    end

    y.mult(2 ** (sqrt_steps + 1), prec)
  end
end

# File 'lib/bigdecimal/math.rb', line 200

def PI(prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :PI)
  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.mult(1, prec)
end

# File 'lib/bigdecimal/math.rb', line 87

def sin(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :sin)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sin)
  return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan?
  n    = prec + BigDecimal.double_fig
  one  = BigDecimal("1")
  two  = BigDecimal("2")
  sign, x = _sin_periodic_reduction(x, n)
  x1   = x
  x2   = x.mult(x,n)
  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
    x1  = -x2.mult(x1,n)
    i  += two
    z  *= (i-one) * i
    d   = x1.div(z,m)
    y  += d
  end
  y = BigDecimal("1") if y > 1
  y.mult(sign, prec)
end

# File 'lib/bigdecimal/math.rb', line 43

def sqrt(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :sqrt)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sqrt)
  x.sqrt(prec)
end

# File 'lib/bigdecimal/math.rb', line 146

def tan(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :tan)
  sin(x, prec + BigDecimal.double_fig).div(cos(x, prec + BigDecimal.double_fig), prec)
end