From 6b3e7122703e4c68a89a19a4e4f2df6b7b948def Mon Sep 17 00:00:00 2001 From: tompng Date: Mon, 2 Jun 2025 00:54:29 +0900 Subject: [PATCH] Implement BigMath.gamma and BigMath.lgamma Implement gamma and lgamma with Spouge's approximation --- lib/bigdecimal/math.rb | 110 ++++++++++++++++++++++++++++++++ test/bigdecimal/test_bigmath.rb | 43 +++++++++++++ 2 files changed, 153 insertions(+) diff --git a/lib/bigdecimal/math.rb b/lib/bigdecimal/math.rb index 16e095df..7f895eb4 100644 --- a/lib/bigdecimal/math.rb +++ b/lib/bigdecimal/math.rb @@ -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) # @@ -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 diff --git a/test/bigdecimal/test_bigmath.rb b/test/bigdecimal/test_bigmath.rb index 3fbf24c2..c1815a97 100644 --- a/test/bigdecimal/test_bigmath.rb +++ b/test/bigdecimal/test_bigmath.rb @@ -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) } @@ -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 @@ -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