Implement BigMath.gamma and BigMath.lgamma

Implement gamma and lgamma with Spouge's approximation

Author: tompng

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,123 @@ def expm1(x, prec)
     exp_prec > 0 ? 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 < 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
+
+  # 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 < 1000 * prec
+      log_gamma = BigMath.log(_gamma_positive_integer(x, prec2), 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
+
+  private_class_method def _gamma_positive_integer(x, prec) # :nodoc:
+    return x if x == 1
+    numbers = (1..x - 1).map {|i| BigDecimal(i) }
+    while numbers.size > 1
+      numbers = numbers.each_slice(2).map {|a, b| b ? a.mult(b, prec) : a }
+    end
+    numbers.first
+  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

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
@@ -469,4 +473,47 @@ def test_expm1
     assert_in_exact_precision(exp(BigDecimal("1.23e-10"), 120) - 1, expm1(BigDecimal("1.23e-10"), 100), 100)
     assert_in_exact_precision(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_equal(
+      BigDecimal('0.28242294079603478742934215780245355184774949260912e456569'),
+      BigMath.gamma(100000, 50)
+    )
+    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