Permalink
Browse files

Add Cook's Erf implementation.

  • Loading branch information...
1 parent 15c5b26 commit 960f3280663f19f3dc31b890472e1c8742d3e1e1 @colomon colomon committed Oct 13, 2012
Showing with 44 additions and 0 deletions.
  1. +19 −0 lib/Builtins.cs
  2. +2 −0 lib/CORE.setting
  3. +23 −0 lib/SpecialMathFunctions.cs
View
19 lib/Builtins.cs
@@ -1088,6 +1088,25 @@ public partial class Builtins {
}
}
+ static readonly Func<Constants,Variable,Variable> erf_d = erf;
+ [ImplicitConsts] public static Variable erf(Constants c, Variable a1) {
+ P6any o1 = a1.Fetch();
+ int r1;
+ if (!o1.mo.is_any)
+ return HandleSpecial1(c, a1,o1, erf_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.Erf(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();
View
2 lib/CORE.setting
@@ -495,6 +495,7 @@ my class Cool {
method lgamma() { Q:CgOp { (lgamma {self}) } }
method expm1() { Q:CgOp { (expm1 {self}) } }
method log1p() { Q:CgOp { (log1p {self}) } }
+ method erf() { Q:CgOp { (erf {self}) } }
method split($matcher, $limit?, :$all?) {
my $matchrx = (($matcher ~~ Regex) ?? $matcher !! /$matcher/);
@@ -3869,6 +3870,7 @@ sub gamma($x) is pure { $x.gamma }
sub lgamma($x) is pure { $x.lgamma }
sub expm1($x) is pure { $x.expm1 }
sub log1p($x) is pure { $x.log1p }
+sub erf($x) is pure { $x.erf }
sub expmod($exp, $power, $mod) is pure { Q:CgOp { (expmod {$exp.Int} {$power.Int} {$mod.Int}) } }
sub is-prime($candidate, $tries = 100) {
View
23 lib/SpecialMathFunctions.cs
@@ -873,6 +873,29 @@ public static double Square(double x)
return x*x;
}
+ public static double Erf(double x)
+ {
+ // constants
+ double a1 = 0.254829592;
+ double a2 = -0.284496736;
+ double a3 = 1.421413741;
+ double a4 = -1.453152027;
+ double a5 = 1.061405429;
+ double p = 0.3275911;
+
+ // Save the sign of x
+ int sign = 1;
+ if (x < 0)
+ sign = -1;
+ x = Math.Abs(x);
+
+ // A&S formula 7.1.26
+ double t = 1.0 / (1.0 + p*x);
+ double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*Math.Exp(-x*x);
+
+ return sign*y;
+ }
+
///////////////////////////////////////////////////////////////////////////////
////// Utility functions //////////////////////////////////////////////////////

0 comments on commit 960f328

Please sign in to comment.