Pull in John D. Cook's Special Functions C# module and use it to impl…

…ement Real versions of gamma and expm1.
sorear · Oct 11, 2012 · bbf7bbd · bbf7bbd
1 parent 7df6137
commit bbf7bbd
Show file tree

Hide file tree

Showing 3 changed files with 910 additions and 0 deletions.
diff --git a/lib/Builtins.cs b/lib/Builtins.cs
@@ -8,6 +8,7 @@
 using System.Threading;
 using System.Text;
 using System.Net.Sockets;
+using SpecialMathFunctions;
 
 namespace Niecza {
     // IForeignInterpreter is used for runtime loading of the Perl 5
@@ -1011,6 +1012,44 @@ public partial class Builtins {
         }
     }
 
+    static readonly Func<Constants,Variable,Variable> gamma_d = gamma;
+    [ImplicitConsts] public static Variable gamma(Constants c, Variable a1) {
+        P6any o1 = a1.Fetch();
+        int r1;
+        if (!o1.mo.is_any)
+            return HandleSpecial1(c, a1,o1, gamma_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.Gamma(v1));
+        }
+    }
+
+    static readonly Func<Constants,Variable,Variable> expm1_d = expm1;
+    [ImplicitConsts] public static Variable expm1(Constants c, Variable a1) {
+        P6any o1 = a1.Fetch();
+        int r1;
+        if (!o1.mo.is_any)
+            return HandleSpecial1(c, a1,o1, expm1_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.ExpMinusOne(v1));
+        }
+    }
+
     static readonly Func<Constants,Variable,Variable,Variable> atan2_d = atan2;
     [ImplicitConsts] public static Variable atan2(Constants c, Variable a1, Variable a2) {
         P6any o1 = a1.Fetch();

diff --git a/lib/CORE.setting b/lib/CORE.setting
@@ -491,6 +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 split($matcher, $limit?, :$all?) {
         my $matchrx = (($matcher ~~ Regex) ?? $matcher !! /$matcher/);
         my $str = ~self;
@@ -3853,12 +3856,16 @@ sub acosech($x) is pure { $x.acosech }
 sub cotanh($x) is pure { $x.cotanh }
 sub acotanh($x) is pure { $x.acotanh }
 sub atan2($y, $x = 1) is pure { $y.atan2($x) }
+
 sub hypot($a, $b) is pure {
     my $max = $a.abs max $b.abs;
     my $min = $a.abs min $b.abs;
     my $r = $min / $max;
     $max * sqrt(1 + $r * $r);
 }
+sub gamma($x) is pure { $x.gamma }
+sub expm1($x) is pure { $x.expm1 }
+
 sub expmod($exp, $power, $mod) is pure { Q:CgOp { (expmod {$exp.Int} {$power.Int} {$mod.Int}) } }
 sub is-prime($candidate, $tries = 100) {
     # Miller-Rabin via TimToady