Use portable expression for constant 2/sqrt(pi)
M_2_SQRTPI is not part of standard C++, replace it with an
equivalent expression. gcc immediately folds this and thus
generates the same code even at -O0.
Change-Id: I5d38e67a4508c1d43a1a8e5f1639f2243cdb143e
diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index da49f32..3264e2e 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -575,19 +575,28 @@
return y < x ? y : x;
}
-// erf is defined as an integral that cannot be expressed analyticaly
+// erf is defined as an integral that cannot be expressed analytically
// 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));
+ // We evaluate the constant as follows:
+ // 2 / sqrt(pi) = 1 / sqrt(atan(1.))
+ // On POSIX sytems it is defined as M_2_SQRTPI, but this is not
+ // portable and the type may not be T. The above expression
+ // evaluates to full precision with IEEE arithmetic and, since it's
+ // constant, the compiler can generate exactly the same code. gcc
+ // does so even at -O0.
+ return Jet<T, N>(erf(x.a), x.v * exp(-x.a * x.a) * (T(1) / sqrt(atan(T(1)))));
}
// 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));
+ // See in erf() above for the evaluation of the constant in the derivative.
+ return Jet<T, N>(erfc(x.a),
+ -x.v * exp(-x.a * x.a) * (T(1) / sqrt(atan(T(1)))));
}
// Bessel functions of the first kind with integer order equal to 0, 1, n.
