module BigMath

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

def PI(prec)
  # Gauss–Legendre algorithm
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :PI)
  n = prec + BigDecimal::Internal::EXTRA_PREC
  a = BigDecimal(1)
  b = BigDecimal(0.5, 0).sqrt(n)
  s = BigDecimal(0.25, 0)
  t = 1
  while a != b && (a - b).exponent > 1 - n
    c = (a - b).div(2, n)
    a, b = (a + b).div(2, n), (a * b).sqrt(n)
    s = s.sub(c * c * t, n)
    t *= 2
  end
  (a * b).div(s, prec)
end

def _erf_taylor(x, a, erf_a, prec) # :nodoc:
ero?
x+a)*exp((x+a)**2)*sqrt(pi)/2
 c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...
a)*f(x+a)
c2*x + 3*c3*x**2 + 4*c4*x**3 + 5*c5*x**4 + ...
a)*(c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...)
 1
* a * c1) / 2
* a * c2) / 3
* a * c3) / 4
are positive when a >= 0
2).div(sqrt(PI(prec), prec), prec)
scale.mult(exp(-a*a, prec), prec), prec)
_prev).add(1, prec).mult(x, prec)
next, prec)
(x, prec) + a * c_next).mult(2, prec).mult(x, prec).div(k, prec)
 prec)
c_next, cn
 c_next].all? { |c| c.zero?  || (c.exponent < sum.exponent - prec) }
ale.mult(exp(-(x + a).mult(x + a, prec), prec), prec), prec)
mal(1) : value

def _erfc_asymptotic(x, prec) # :nodoc:
)*sqrt(pi)*exp(x**2)/2
e following differential equation:
+ 1
uation, we can derive the following asymptotic expansion:
sum { (-1)**k * (2*k)! / 4**k / k! / x**(2*k)) } / x
xpansion does not converge.
 k that satisfies (2*k)! / 4**k / k! / x**(2*k) < 10**(-prec),
alculate erfc within the given precision.
approximation, we can simplify this condition to:
(k) - k - 2*k*log(x) < -prec*log(10)
 is minimized when k = x**2.
:Internal::EXTRA_PREC
ternal.fast_to_f(x)
).floor).bsearch do |k|
 k * Math.log(k) - k - 2 * k * Math.log(xf) < -prec * Math.log(10)

(x2, prec)`, x2 needs extra log10(x**2) digits of precision
 + (2 * Math.log10(xf)).ceil)
k|
c).mult(1 - 2 * k, prec).div(2, prec)
prec)
c).mult(PI(prec).sqrt(prec), prec), prec).div(x, prec)

def _exp_binary_splitting(x, prec) # :nodoc:
) if x.zero?
fies x**k / k! < 10**(-prec)
)
Internal.float_log(x.abs)
h { |k| Math.lgamma(k + 1)[0] - k * logx > prec * log10 }
/2*(1+x/3*(1+x/4*(1+x/5*(1+...)))))
ernal.taylor_sum_binary_splitting(x, [*1..step], prec)

def _gamma_positive_integer(x, prec) # :nodoc:
).map {|i| BigDecimal(i) }
> 1
.each_slice(2).map {|a, b| b ? a.mult(b, prec) : a }

def _gamma_spouge_sum_part(x, prec) # :nodoc:
ation
+ 0.5) * exp(-x - a) * (sqrt(2 * pi)  + (1..a - 1).sum{|k| c[k] / (x + k) } + epsilon)
**k * (a - k)**(k - 0.5) * exp(a - k) / (k - 1)!
unded by a**(-0.5) * (2 * pi) ** (-a - 0.5)
 a for given precision
g10(2 * Math::PI)).ceil
t of c[k] in low precision to estimate required precision
0)
(1, low_prec)
a-1).map do |k|
h.log10(k - 1) if k > 1
 log10f + (k - 0.5) * Math.log10(a - k) - BigDecimal::Internal.float_log(x_low_prec.add(k, low_prec)) / log10f
 of sum by Stirling's approximation
 = x < 1 ? -Math.log10(a) / 2 : Math.log10(2 * Math::PI) / 2 + x_low_prec.add(0.5, low_prec) * Math.log10(x_low_prec / x_low_prec.add(a, low_prec))
d precision of c[k]
ts.max.ceil - approx_sum_exponent.floor, 0].max + prec
-1, prec2)
).sqrt(prec).mult(BigMath.exp(-a, prec), prec)
|k|
* (a - k)**(k - 0.5) * exp(-k) / (k-1)! / (x + k)
prec2) if k > 1
prec2)
imal((a - k) ** k), prec2).div(BigDecimal(a - k).sqrt(prec2).mult(x.add(k, prec2), prec2), prec2)
x + k)
prec2)

def _sin_around_zero(x, prec) # :nodoc:
eral parts
) = sin(x.xx + 0.00xx + 0.0000xxxx + ...)
each part and restore sin(0.xxxxxxxx...) using addition theorem.


ncate(n)
plitting(partial_x, prec)
qrt(prec)
 c).add(cos * s, prec), (cos * c).sub(sin * s, prec)
l(-1), BigDecimal(1))

def _sin_binary_splitting(x, prec) # :nodoc:

)
fies x2**k / (2k+1)! < 10**(-prec)
)
Internal.float_log(x.abs)
h { |k| Math.lgamma(2 * k + 1)[0] - 2 * k * logx > prec * log10 }
ator sequence for binary splitting
(2*3)*(1-x2/(4*5)*(1-x2/(6*7)*(1-x2/(8*9)*(1-...)))))
{|i| -(2 * i) * (2 * i + 1) }
al::Internal.taylor_sum_binary_splitting(x2, ds, prec), prec)

def _sin_periodic_reduction(x, prec, add_half_pi: false) # :nodoc:
ath::PI/2 <= x && x <= Math::PI/2 && !add_half_pi
igDecimal::Internal::EXTRA_PREC
exponent, 0].max + BigDecimal::Internal::EXTRA_PREC
+ pi_extra_prec)
alf_pi ? x + pi : x + half_pi).divmod(pi)
od_prec + mod.exponent <= 0
all to estimate required pi precision
prec * 3 / 2 + BigDecimal::Internal::EXTRA_PREC
mod.exponent < prec
ired precision of pi
 - mod.exponent + BigDecimal::Internal::EXTRA_PREC
 == 0 ? 1 : -1, mod.mult(1, prec)]

def _sinpix(x, pi, prec) # :nodoc:
 1
5 # to avoid sin(pi*x) loss of precision for x close to 1
i, prec), prec)

def acos(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :acos)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :acos)
  raise Math::DomainError, "Out of domain argument for acos" if x < -1 || x > 1
  return BigDecimal::Internal.nan_computation_result if x.nan?
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  return (PI(prec2) / 2).sub(asin(x, prec2), prec) if x < 0
  return PI(prec2).div(2, prec) if x.zero?
  sin = (1 - x**2).sqrt(prec2)
  atan(sin.div(x, prec2), prec)
end

def acosh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :acosh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :acosh)
  raise Math::DomainError, "Out of domain argument for acosh" if x < 1
  return BigDecimal::Internal.infinity_computation_result if x.infinite?
  return BigDecimal::Internal.nan_computation_result if x.nan?
  log(x + sqrt(x**2 - 1, prec + BigDecimal::Internal::EXTRA_PREC), prec)
end

def asin(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :asin)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :asin)
  raise Math::DomainError, "Out of domain argument for asin" if x < -1 || x > 1
  return BigDecimal::Internal.nan_computation_result if x.nan?
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  cos = (1 - x**2).sqrt(prec2)
  if cos.zero?
    PI(prec2).div(x > 0 ? 2 : -2, prec)
  else
    atan(x.div(cos, prec2), prec)
  end
end

def asinh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :asinh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :asinh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?
  return -asinh(-x, prec) if x < 0
  sqrt_prec = prec + [-x.exponent, 0].max + BigDecimal::Internal::EXTRA_PREC
  log(x + sqrt(x**2 + 1, sqrt_prec), prec)
end

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::Internal::EXTRA_PREC
  return PI(n).div(2 * x.infinite?, prec) if x.infinite?
  x = -x if neg = x < 0
  x = BigDecimal(1).div(x, n) if inv = x < -1 || x > 1
  # Solve tan(y) - x = 0 with Newton's method
  # Repeat: y -= (tan(y) - x) * cos(y)**2
  y = BigDecimal(Math.atan(BigDecimal::Internal.fast_to_f(x)), 0)
  BigDecimal::Internal.newton_loop(n) do |p|
    s = sin(y, p)
    c = (1 - s * s).sqrt(p)
    y = y.sub(c * (s.sub(c * x.mult(1, p), p)), p)
  end
  y = PI(n) / 2 - y if inv
  y.mult(neg ? -1 : 1, prec)
end

def atan2(y, x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :atan2)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan2)
  y = BigDecimal::Internal.coerce_to_bigdecimal(y, prec, :atan2)
  return BigDecimal::Internal.nan_computation_result if x.nan? || y.nan?
  if x.infinite? || y.infinite?
    one = BigDecimal(1)
    zero = BigDecimal(0)
    x = x.infinite? ? (x > 0 ? one : -one) : zero
    y = y.infinite? ? (y > 0 ? one : -one) : y.sign * zero
  end
  return x.sign >= 0 ? BigDecimal(0) : y.sign * PI(prec) if y.zero?
  y = -y if neg = y < 0
  xlarge = y.abs < x.abs
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  if x > 0
    v = xlarge ? atan(y.div(x, prec2), prec) : PI(prec2) / 2 - atan(x.div(y, prec2), prec2)
  else
    v = xlarge ? PI(prec2) - atan(-y.div(x, prec2), prec2) : PI(prec2) / 2 + atan(x.div(-y, prec2), prec2)
  end
  v.mult(neg ? -1 : 1, prec)
end

def atanh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :atanh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atanh)
  raise Math::DomainError, "Out of domain argument for atanh" if x < -1 || x > 1
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x == 1
  return -BigDecimal::Internal.infinity_computation_result if x == -1
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  (log(x + 1, prec2) - log(1 - x, prec2)).div(2, prec)
end

def cbrt(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :cbrt)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cbrt)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?
  return BigDecimal(0) if x.zero?
  x = -x if neg = x < 0
  ex = x.exponent / 3
  x = x._decimal_shift(-3 * ex)
  y = BigDecimal(Math.cbrt(BigDecimal::Internal.fast_to_f(x)), 0)
  BigDecimal::Internal.newton_loop(prec + BigDecimal::Internal::EXTRA_PREC) do |p|
    y = (2 * y + x.div(y, p).div(y, p)).div(3, p)
  end
  y._decimal_shift(ex).mult(neg ? -1 : 1, prec)
end

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?
  n = prec + BigDecimal::Internal::EXTRA_PREC
  sign, x = _sin_periodic_reduction(x, n, add_half_pi: true)
  _sin_around_zero(x, n).mult(sign, prec)
end

def cosh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :cosh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cosh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite?
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  e = exp(x, prec2)
  (e + BigDecimal(1).div(e, prec2)).div(2, prec)
end

def erf(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :erf)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erf)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal(x.infinite?) if x.infinite?
  return BigDecimal(0) if x == 0
  return -erf(-x, prec) if x < 0
  return BigDecimal(1) if x > 5000000000 # erf(5000000000) > 1 - 1e-10000000000000000000
  if x > 8
    xf = BigDecimal::Internal.fast_to_f(x)
    log10_erfc = -xf ** 2 / Math.log(10) - Math.log10(xf * Math::PI ** 0.5)
    erfc_prec = [prec + log10_erfc.ceil, 1].max
    erfc = _erfc_asymptotic(x, erfc_prec)
    return BigDecimal(1).sub(erfc, prec) if erfc
  end
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  x_smallprec = x.mult(1, Integer.sqrt(prec2) / 2)
  # Taylor series of x with small precision is fast
  erf1 = _erf_taylor(x_smallprec, BigDecimal(0), BigDecimal(0), prec2)
  # Taylor series converges quickly for small x
  _erf_taylor(x - x_smallprec, x_smallprec, erf1, prec2).mult(1, prec)
end

def erfc(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :erfc)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erfc)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal(1 - x.infinite?) if x.infinite?
  return BigDecimal(1).sub(erf(x, prec + BigDecimal::Internal::EXTRA_PREC), prec) if x < 0.5
  return BigDecimal(0) if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow)
  if x >= 8
    y = _erfc_asymptotic(x, prec)
    return y.mult(1, prec) if y
  end
  # erfc(x) = 1 - erf(x) < exp(-x**2)/x/sqrt(pi)
  # Precision of erf(x) needs about log10(exp(-x**2)/x/sqrt(pi)) extra digits
  log10 = 2.302585092994046
  xf = BigDecimal::Internal.fast_to_f(x)
  high_prec = prec + BigDecimal::Internal::EXTRA_PREC + ((xf**2 + Math.log(xf) + Math.log(Math::PI)/2) / log10).ceil
  BigDecimal(1).sub(erf(x, high_prec), prec)
end

def 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::Internal::EXTRA_PREC + cnt
  x = x._decimal_shift(-cnt)
  # Decimal form of bit-burst algorithm
  # Calculate exp(x.xxxxxxxxxxxxxxxx) as
  # exp(x.xx) * exp(0.00xx) * exp(0.0000xxxx) * exp(0.00000000xxxxxxxx)
  x = x.mult(1, prec2)
  n = 2
  y = BigDecimal(1)
  BigDecimal.save_limit do
    BigDecimal.limit(0)
    while x != 0 do
      partial_x = x.truncate(n)
      x -= partial_x
      y = y.mult(_exp_binary_splitting(partial_x, prec2), prec2)
      n *= 2
    end
  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 expm1(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :expm1)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :expm1)
  return BigDecimal(-1) if x.infinite? == -1
  exp_prec = prec
  if x < -1
    # log10(exp(x)) = x * log10(e)
    lg_e = 0.4342944819032518
    exp_prec = prec + (lg_e * x).ceil + BigDecimal::Internal::EXTRA_PREC
  elsif x < 1
    exp_prec = prec - x.exponent + BigDecimal::Internal::EXTRA_PREC
  else
    exp_prec = prec
  end
  return BigDecimal(-1) if exp_prec <= 0
  exp(x, exp_prec).sub(1, prec)
end

def frexp(x)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, 0, :frexp)
  return [x, 0] unless x.finite?
  exponent = x.exponent
  [x._decimal_shift(-exponent), exponent]
end

def gamma(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma)
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  if x < 0.5
    raise Math::DomainError, 'Numerical argument is out of domain - gamma' if x.frac.zero?
    # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z)
    pi = PI(prec2)
    sin = _sinpix(x, pi, prec2)
    return pi.div(gamma(1 - x, prec2).mult(sin, prec2), prec)
  elsif x.frac.zero? && x < 1000 * prec
    return _gamma_positive_integer(x, prec2).mult(1, prec)
  end
  a, sum = _gamma_spouge_sum_part(x, prec2)
  (x + (a - 1)).power(x - 0.5, prec2).mult(BigMath.exp(1 - x, prec2), prec2).mult(sum, prec)
end

def hypot(x, y, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :hypot)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :hypot)
  y = BigDecimal::Internal.coerce_to_bigdecimal(y, prec, :hypot)
  return BigDecimal::Internal.nan_computation_result if x.nan? || y.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite? || y.infinite?
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  sqrt(x.mult(x, prec2) + y.mult(y, prec2), prec)
end

def ldexp(x, exponent)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, 0, :ldexp)
  x.finite? ? x._decimal_shift(exponent) : x
end

def lgamma(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma)
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  if x < 0.5
    return [BigDecimal::INFINITY, 1] if x.frac.zero?
    loop do
      # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z)
      pi = PI(prec2)
      sin = _sinpix(x, pi, prec2)
      log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec2).first + BigMath.log(sin.abs, prec2), prec)
      return [log_gamma, sin > 0 ? 1 : -1] if prec2 + log_gamma.exponent > prec + BigDecimal::Internal::EXTRA_PREC
      # Retry with higher precision if loss of significance is too large
      prec2 = prec2 * 3 / 2
    end
  elsif x.frac.zero? && x < 1000 * prec
    log_gamma = BigMath.log(_gamma_positive_integer(x, prec2), prec)
    [log_gamma, 1]
  else
    # if x is close to 1 or 2, increase precision to reduce loss of significance
    diff1_exponent = (x - 1).exponent
    diff2_exponent = (x - 2).exponent
    extremely_near_one = diff1_exponent < -prec2
    extremely_near_two = diff2_exponent < -prec2
    if extremely_near_one || extremely_near_two
      # If x is extreamely close to base = 1 or 2, linear interpolation is accurate enough.
      # Taylor expansion at x = base is: (x - base) * digamma(base) + (x - base) ** 2 * trigamma(base) / 2 + ...
      # And we can ignore (x - base) ** 2 and higher order terms.
      base = extremely_near_one ? 1 : 2
      d = BigDecimal(1)._decimal_shift(1 - prec2)
      log_gamma_d, sign = lgamma(base + d, prec2)
      return [log_gamma_d.mult(x - base, prec2).div(d, prec), sign]
    end
    prec2 += [-diff1_exponent, -diff2_exponent, 0].max
    a, sum = _gamma_spouge_sum_part(x, prec2)
    log_gamma = BigMath.log(sum, prec2).add((x - 0.5).mult(BigMath.log(x.add(a - 1, prec2), prec2), prec2) + 1 - x, prec)
    [log_gamma, 1]
  end
end

def 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::Internal::EXTRA_PREC
  # Reduce x to near 1
  if x > 1.01 || x < 0.99
    # log(x) = log(x/exp(logx_approx)) + logx_approx
    logx_approx = BigDecimal(BigDecimal::Internal.float_log(x), 0)
    x = x.div(exp(logx_approx, prec2), prec2)
  else
    logx_approx = BigDecimal(0)
  end
  # Solve exp(y) - x = 0 with Newton's method
  # Repeat: y -= (exp(y) - x) / exp(y)
  y = BigDecimal(BigDecimal::Internal.float_log(x), 0)
  exp_additional_prec = [-(x - 1).exponent, 0].max
  BigDecimal::Internal.newton_loop(prec2) do |p|
    expy = exp(y, p + exp_additional_prec)
    y = y.sub(expy.sub(x, p).div(expy, p), p)
  end
  y.add(logx_approx, prec)
end

def log10(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :log10)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log10)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2
  v = log(x, prec2).div(log(BigDecimal(10), prec2), prec2)
  # Perform half-up rounding to calculate log10(10**n)==n correctly in every rounding mode
  v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
  v.mult(1, prec)
end

def log1p(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :log1p)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log1p)
  raise Math::DomainError, 'Out of domain argument for log1p' if x < -1
  return log(x + 1, prec)
end

def log2(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :log2)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log2)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2
  v = log(x, prec2).div(log(BigDecimal(2), prec2), prec2)
  # Perform half-up rounding to calculate log2(2**n)==n correctly in every rounding mode
  v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP)
  v.mult(1, prec)
end

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::Internal::EXTRA_PREC
  sign, x = _sin_periodic_reduction(x, n)
  _sin_around_zero(x, n).mult(sign, prec)
end

def sinh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :sinh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sinh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite?
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC
  prec2 -= x.exponent if x.exponent < 0
  e = exp(x, prec2)
  (e - BigDecimal(1).div(e, prec2)).div(2, prec)
end

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

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

def tanh(x, prec)
  prec = BigDecimal::Internal.coerce_validate_prec(prec, :tanh)
  x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :tanh)
  return BigDecimal::Internal.nan_computation_result if x.nan?
  return BigDecimal(x.infinite?) if x.infinite?
  prec2 = prec + BigDecimal::Internal::EXTRA_PREC + [-x.exponent, 0].max
  e = exp(x, prec2)
  einv = BigDecimal(1).div(e, prec2)
  (e - einv).div(e + einv, prec)
end