Add erf and erfc to jet.h, including tests in jet_test.cc
erf is necessary for evaluating Gaussian functions.
erfc was added because it is so similar to erf.
Change-Id: I5e470dbe013cc938fabb87cde3b0ebf26a90fff4
diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index 7aafaa0..c551d51 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -388,6 +388,8 @@
using std::ceil;
using std::cos;
using std::cosh;
+using std::erf;
+using std::erfc;
using std::exp;
using std::exp2;
using std::floor;
@@ -573,6 +575,21 @@
return y < x ? y : x;
}
+// erf is defined as an integral that cannot be expressed analyticaly
+// however, the derivative is trivial to compute
+// erf(x + h) = erf(x) + h * 2*exp(-x^2)/sqrt(pi)
+template <typename T, int N>
+inline Jet<T, N> erf(const Jet<T, N>& x) {
+ return Jet<T, N>(erf(x.a), x.v * M_2_SQRTPI * exp(-x.a*x.a));
+}
+
+// erfc(x) = 1-erf(x)
+// erfc(x + h) = erfc(x) + h * (-2*exp(-x^2)/sqrt(pi))
+template <typename T, int N>
+inline Jet<T, N> erfc(const Jet<T, N>& x) {
+ return Jet<T, N>(erfc(x.a), -x.v * M_2_SQRTPI * exp(-x.a*x.a));
+}
+
// Bessel functions of the first kind with integer order equal to 0, 1, n.
//
// Microsoft has deprecated the j[0,1,n]() POSIX Bessel functions in favour of
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index a6cad62..f1e9586 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -571,6 +571,28 @@
EXPECT_EQ(expected, b);
}
+ { // Check that erf works.
+ J a = MakeJet(10.123, -2.7, 1e-3);
+ J b = erf(a);
+ J expected = MakeJet(erf(a.a), 0.0, 0.0);
+ EXPECT_EQ(expected, b);
+ }
+ NumericalTest("erf", erf<double, 2>, -1.0);
+ NumericalTest("erf", erf<double, 2>, 1e-5);
+ NumericalTest("erf", erf<double, 2>, 0.5);
+ NumericalTest("erf", erf<double, 2>, 100.0);
+
+ { // Check that erfc works.
+ J a = MakeJet(10.123, -2.7, 1e-3);
+ J b = erfc(a);
+ J expected = MakeJet(erfc(a.a), 0.0, 0.0);
+ EXPECT_EQ(expected, b);
+ }
+ NumericalTest("erfc", erfc<double, 2>, -1.0);
+ NumericalTest("erfc", erfc<double, 2>, 1e-5);
+ NumericalTest("erfc", erfc<double, 2>, 0.5);
+ NumericalTest("erfc", erfc<double, 2>, 100.0);
+
{ // Check that cbrt(x * x * x) == x.
J z = x * x * x;
J w = cbrt(z);
