Skip to content

Implement BigMath.gamma and BigMath.lgamma #451

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Implement BigMath.gamma and BigMath.lgamma #451

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
110 changes: 110 additions & 0 deletions lib/bigdecimal/math.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@
# log10(x, prec)
# log1p(x, prec)
# expm1(x, prec)
# gamma(x, prec)
# lgamma(x, prec)
# PI (prec)
# E (prec) == exp(1.0,prec)
#
Expand Down Expand Up @@ -568,6 +570,114 @@ def expm1(x, prec)
exp_prec > 0 ? BigMath.exp(x, exp_prec).sub(1, prec) : BigDecimal(-1)
end

# call-seq:
# BigMath.gamma(decimal, numeric) -> BigDecimal
#
# Computes the gamma function of +decimal+ to the specified number of
# digits of precision, +numeric+.
#
# BigMath.gamma(BigDecimal('0.5'), 32).to_s
# #=> "0.17724538509055160272981674833411e1"
#
def gamma(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma)
prec2 = prec + BigDecimal.double_fig
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, prec).mult(sin, prec2), prec)
elsif x.frac.zero? && x < 10000
return BigDecimal((1..x - 1).inject(1, :*)).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

# call-seq:
# BigMath.lgamma(decimal, numeric) -> [BigDecimal, Integer]
#
# Computes the natural logarithm of the absolute value of the gamma function
# of +decimal+ to the specified number of digits of precision, +numeric+ and its sign.
#
# BigMath.lgamma(BigDecimal('0.5'), 32)
# #=> [0.57236494292470008707171367567653e0, 1]
#
def lgamma(x, prec)
prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma)
x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma)
prec2 = prec + BigDecimal.double_fig
if x < 0.5
return [BigDecimal::INFINITY, 1] 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)
log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec).first + BigMath.log(sin.abs, prec2), prec)
[log_gamma, sin > 0 ? 1 : -1]
elsif x.frac.zero? && x < 10000
log_gamma = BigMath.log(BigDecimal((1..x - 1).inject(1, :*)), prec)
[log_gamma, 1]
else
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

# Returns sum part: sqrt(2*pi) and c[k]/(x+k) terms of Spouge's approximation
private_class_method def _gamma_spouge_sum_part(x, prec) # :nodoc:
x -= 1
# Spouge's approximation
# x! = (x + a)**(x + 0.5) * exp(-x - a) * (sqrt(2 * pi) + (1..a - 1).sum{|k| c[k] / (x + k) } + epsilon)
# where c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(a - k) / (k - 1)!
# and epsilon is bounded by a**(-0.5) * (2 * pi) ** (-a - 0.5)

# Estimate required a for given precision
a = (prec / Math.log10(2 * Math::PI)).ceil

# Calculate exponent of c[k] in low precision to estimate required precision
low_prec = 16
log10f = Math.log(10)
x_low_prec = x.mult(1, low_prec)
loggamma_k = 0
ck_exponents = (1..a-1).map do |k|
loggamma_k += Math.log10(k - 1) if k > 1
-loggamma_k - k / log10f + (k - 0.5) * Math.log10(a - k) - BigMath.log10(x_low_prec.add(k, low_prec), low_prec)
end

# Estimate exponent of sum by Stirling's approximation
approx_sum_exponent = 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))

# Determine required precision of c[k]
prec2 = [ck_exponents.max.ceil - approx_sum_exponent.floor, 0].max + prec

einv = BigMath.exp(-1, prec2)
sum = (PI(prec) * 2).sqrt(prec).mult(BigMath.exp(-a, prec), prec)
y = BigDecimal(1)
(1..a - 1).each do |k|
# c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(-k) / (k-1)! / (x + k)
y = y.div(1 - k, prec2) if k > 1
y = y.mult(einv, prec2)
z = y.mult(BigDecimal((a - k) ** k), prec2).div(BigDecimal(a - k).sqrt(prec2).mult(x.add(k, prec2), prec2), prec2)
# sum += c[k] / (x + k)
sum = sum.add(z, prec2)
end
[a, sum]
end

# Returns sin(pi * x), for gamma reflection formula calculation
private_class_method def _sinpix(x, pi, prec) # :nodoc:
x = x % 2
sign = x > 1 ? -1 : 1
x %= 1
x = 1 - x if x > 0.5 # to avoid sin(pi*x) loss of precision for x close to 1
sign * sin(x.mult(pi, prec), prec)
end

# call-seq:
# PI(numeric) -> BigDecimal
Expand Down
43 changes: 43 additions & 0 deletions test/bigdecimal/test_bigmath.rb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,8 @@ def test_consistent_precision_acceptance
assert_consistent_precision_acceptance {|prec| BigMath.log10(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.log1p(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.expm1(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.gamma(x, prec) }
assert_consistent_precision_acceptance {|prec| BigMath.lgamma(x, prec) }

assert_consistent_precision_acceptance {|prec| BigMath.E(prec) }
assert_consistent_precision_acceptance {|prec| BigMath.PI(prec) }
Expand Down Expand Up @@ -112,6 +114,8 @@ def test_coerce_argument
assert_equal(log10(bd, N), log10(f, N))
assert_equal(log1p(bd, N), log1p(f, N))
assert_equal(expm1(bd, N), expm1(f, N))
assert_equal(gamma(bd, N), gamma(f, N))
assert_equal(lgamma(bd, N), lgamma(f, N))
end

def test_sqrt
Expand Down Expand Up @@ -471,4 +475,43 @@ def test_expm1
assert_in_exact_precision(BigMath.exp(BigDecimal("1.23e-10"), 120) - 1, expm1(BigDecimal("1.23e-10"), 100), 100)
assert_in_exact_precision(BigMath.exp(123, 120) - 1, expm1(BigDecimal("123"), 100), 100)
end

def test_gamma
[-1.8, -0.7, 0.6, 1.5, 2.4].each do |x|
assert_in_epsilon(Math.gamma(x), gamma(BigDecimal(x.to_s), N))
end
[1, 2, 3, 10, 16].each do |x|
assert_equal(Math.gamma(x).round, gamma(BigDecimal(x), N))
end
sqrt_pi = PI(120).sqrt(120)
assert_equal(sqrt_pi.mult(1, 100), gamma(BigDecimal("0.5"), 100))
assert_equal((sqrt_pi * 4).div(3, 100), gamma(BigDecimal("-1.5"), 100))
assert_converge_in_precision {|n| gamma(BigDecimal("0.3"), n) }
assert_converge_in_precision {|n| gamma(BigDecimal("-1.9" + "9" * 30), n) }
assert_converge_in_precision {|n| gamma(BigDecimal("1234.56789"), n) }
assert_converge_in_precision {|n| gamma(BigDecimal("-987.654321"), n) }
end

def test_lgamma
[-2, -1, 0].each do |x|
l, sign = lgamma(BigDecimal(x), N)
assert(l.infinite?)
assert_equal(1, sign)
end
[-1.8, -0.7, 0.6, 1, 1.5, 2, 2.4, 3, 1e+300].each do |x|
l, sign = Math.lgamma(x)
bigl, bigsign = lgamma(BigDecimal(x.to_s), N)
assert_in_epsilon(l, bigl)
assert_equal(sign, bigsign)
end
assert_equal([BigMath.log(PI(120).sqrt(120), 100), 1], lgamma(BigDecimal("0.5"), 100))
assert_converge_in_precision {|n| lgamma(BigDecimal("-1." + "9" * 30), n).first }
assert_converge_in_precision {|n| lgamma(BigDecimal("-3." + "0" * 30 + "1"), n).first }
assert_converge_in_precision {|n| lgamma(BigDecimal("10"), n).first }
assert_converge_in_precision {|n| lgamma(BigDecimal("0.3"), n).first }
assert_converge_in_precision {|n| lgamma(BigDecimal("-1.9" + "9" * 30), n).first }
assert_converge_in_precision {|n| lgamma(BigDecimal("987.65421"), n).first }
assert_converge_in_precision {|n| lgamma(BigDecimal("-1234.56789"), n).first }
assert_converge_in_precision {|n| lgamma(BigDecimal("1e+400"), n).first }
end
end