Permalink
Browse files

Import Cook's implementation of LogGamma.

  • Loading branch information...
1 parent 7ae412a commit 664c87c05cc1de808b07dae419353145a923465c @colomon colomon committed Oct 13, 2012
Showing with 62 additions and 2 deletions.
  1. +19 −0 lib/Builtins.cs
  2. +4 −2 lib/CORE.setting
  3. +39 −0 lib/SpecialMathFunctions.cs
View
@@ -1031,6 +1031,25 @@ public partial class Builtins {
}
}
+ static readonly Func<Constants,Variable,Variable> lgamma_d = lgamma;
+ [ImplicitConsts] public static Variable lgamma(Constants c, Variable a1) {
+ P6any o1 = a1.Fetch();
+ int r1;
+ if (!o1.mo.is_any)
+ return HandleSpecial1(c, a1,o1, lgamma_d);
+ P6any n1 = GetNumber(a1, o1, out r1);
+
+// if (r1 == NR_COMPLEX) {
+// Complex v1 = PromoteToComplex(r1, n1);
+// return MakeComplex(Math.Exp(v1.re) * Math.Cos(v1.im),
+// Math.Exp(v1.re) * Math.Sin(v1.im));
+// }
+ {
+ double v1 = PromoteToFloat(r1, n1);
+ return MakeFloat(SpecialMathFunctions.SpecialFunctions.LogGamma(v1));
+ }
+ }
+
static readonly Func<Constants,Variable,Variable> expm1_d = expm1;
[ImplicitConsts] public static Variable expm1(Constants c, Variable a1) {
P6any o1 = a1.Fetch();
View
@@ -491,8 +491,9 @@ my class Cool {
method rand() { self * rand; }
method roots($n) { roots(self, $n); }
- method gamma() { Q:CgOp { (gamma {self}) } }
- method expm1() { Q:CgOp { (expm1 {self}) } }
+ method gamma() { Q:CgOp { (gamma {self}) } }
+ method lgamma() { Q:CgOp { (lgamma {self}) } }
+ method expm1() { Q:CgOp { (expm1 {self}) } }
method split($matcher, $limit?, :$all?) {
my $matchrx = (($matcher ~~ Regex) ?? $matcher !! /$matcher/);
@@ -3864,6 +3865,7 @@ sub hypot($a, $b) is pure {
$max * sqrt(1 + $r * $r);
}
sub gamma($x) is pure { $x.gamma }
+sub lgamma($x) is pure { $x.lgamma }
sub expm1($x) is pure { $x.expm1 }
sub expmod($exp, $power, $mod) is pure { Q:CgOp { (expmod {$exp.Int} {$power.Int} {$mod.Int}) } }
@@ -360,6 +360,45 @@ public static double Gamma(double x)
return result;
}
+ public static double LogGamma(double x)
+ {
+ Validate.IsNonNegative(x);
+
+ if (x < 12.0)
+ {
+ return Math.Log(Math.Abs(Gamma(x)));
+ }
+
+ // Abramowitz and Stegun 6.1.41
+ // Asymptotic series should be good to at least 11 or 12 figures
+ // For error analysis, see Whittiker and Watson
+ // A Course in Modern Analysis (1927), page 252
+
+ double[] c =
+ {
+ 1.0/12.0,
+ -1.0/360.0,
+ 1.0/1260.0,
+ -1.0/1680.0,
+ 1.0/1188.0,
+ -691.0/360360.0,
+ 1.0/156.0,
+ -3617.0/122400.0
+ };
+ double z = 1.0/(x*x);
+ double sum = c[7];
+ for (int i=6; i >= 0; i--)
+ {
+ sum *= z;
+ sum += c[i];
+ }
+ double series = sum/x;
+
+ double halfLogTwoPi = 0.91893853320467274178032973640562;
+ double logGamma = (x - 0.5)*Math.Log(x) - x + halfLogTwoPi + series;
+ return logGamma;
+ }
+
/// <summary>
/// Factorial. Is exact for return values that fit into an int. Uses table look-up.
/// </summary>

0 comments on commit 664c87c

Please sign in to comment.