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

include/ceres/jet.h

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

internal/ceres/jet_test.cc

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);