support copysign jet
libc++ (default on macOS) uses the copysign function for std::complex
multiplication which causes compilation errors such as these if
std::complex is combined with jets:
/Applications/Xcode_12.4.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/../include/c++/v1/complex:607:19: error: no matching function for call to 'copysign'
__a = copysign(__libcpp_isinf_or_builtin(__a) ? _Tp(1) : _Tp(0), __a);
^~~~~~~~
foo.hpp:218:39: note: in instantiation of function template specialization 'std::__1::operator*<ceres::Jet<double, 3> >' requested here
std::complex<Scalar> result = lhs * rhs;
copysign is also useful in other situations where the sign is required
allowing it to compute without branching.
Change-Id: I84b2d23374f1bf3bae32833016bb686d97d74054
diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index be24617..c97e621 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -426,6 +426,42 @@
return Jet<T, N>(abs(f.a), copysign(T(1), f.a) * f.v);
}
+// copysign(a, b) composes a float with the magnitude of a and the sign of b.
+// Therefore, the function can be formally defined as
+//
+// copysign(a, b) = sgn(b)|a|
+//
+// where
+//
+// d/dx |x| = sgn(x)
+// d/dx sgn(x) = 2δ(x)
+//
+// sgn(x) being the signum function. Differentiating copysign(a, b) with respect
+// to a and b gives:
+//
+// d/da sgn(b)|a| = sgn(a) sgn(b)
+// d/db sgn(b)|a| = 2|a|δ(b)
+//
+// with the dual representation given by
+//
+// copysign(a + da, b + db) ~= sgn(b)|a| + (sgn(a)sgn(b) da + 2|a|δ(b) db)
+//
+// where δ(b) is the Dirac delta function.
+template <typename T, int N>
+inline Jet<T, N> copysign(const Jet<T, N>& f, const Jet<T, N> g) {
+ // The Dirac delta function δ(b) is undefined at b=0 (here it's
+ // infinite) and 0 everywhere else.
+ T d = g.a == T(0) ? std::numeric_limits<T>::infinity() : T(0);
+ T sa = copysign(T(1), f.a); // sgn(a)
+ T sb = copysign(T(1), g.a); // sgn(b)
+ // The second part of the infinitesimal is 2|a|δ(b) which is either infinity
+ // or 0 unless a or any of the values of the b infinitesimal are 0. In the
+ // latter case, the corresponding values become NaNs (multiplying 0 by
+ // infinity gives NaN). We drop the constant factor 2 since it does not change
+ // the result (its values will still be either 0, infinity or NaN).
+ return Jet<T, N>(copysign(f.a, g.a), sa * sb * f.v + abs(f.a) * d * g.v);
+}
+
// log(a + h) ~= log(a) + h / a
template <typename T, int N>
inline Jet<T, N> log(const Jet<T, N>& f) {
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index 4718a20..5631312 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -782,6 +782,103 @@
VL << "z = " << z;
ExpectJetsClose(J{y.a}, z);
}
+ { // copysign(x, +1)
+ J z = copysign(x, J{+1});
+ VL << "z = " << z;
+ ExpectJetsClose(x, z);
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(x, -1)
+ J z = copysign(x, J{-1});
+ VL << "z = " << z;
+ ExpectJetsClose(-x, z);
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(-x, +1)
+
+ J z = copysign(-x, J{+1});
+ VL << "z = " << z;
+ ExpectJetsClose(x, z);
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(-x, -1)
+ J z = copysign(-x, J{-1});
+ VL << "z = " << z;
+ ExpectJetsClose(-x, z);
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(-0, +1)
+ J z = copysign(MakeJet(-0, 1, 2), J{+1});
+ VL << "z = " << z;
+ ExpectJetsClose(MakeJet(+0, 1, 2), z);
+ EXPECT_FALSE(std::signbit(z.a));
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(-0, -1)
+ J z = copysign(MakeJet(-0, 1, 2), J{-1});
+ VL << "z = " << z;
+ ExpectJetsClose(MakeJet(-0, -1, -2), z);
+ EXPECT_TRUE(std::signbit(z.a));
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(+0, -1)
+ J z = copysign(MakeJet(+0, 1, 2), J{-1});
+ VL << "z = " << z;
+ ExpectJetsClose(MakeJet(-0, -1, -2), z);
+ EXPECT_TRUE(std::signbit(z.a));
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(+0, +1)
+ J z = copysign(MakeJet(+0, 1, 2), J{+1});
+ VL << "z = " << z;
+ ExpectJetsClose(MakeJet(+0, 1, 2), z);
+ EXPECT_FALSE(std::signbit(z.a));
+ EXPECT_TRUE(IsFinite(z.v[0]));
+ EXPECT_TRUE(IsFinite(z.v[1]));
+ }
+ { // copysign(+0, +0)
+ J z = copysign(MakeJet(+0, 1, 2), J{+0});
+ VL << "z = " << z;
+ EXPECT_FALSE(std::signbit(z.a));
+ EXPECT_TRUE(std::signbit(z.v[0]));
+ EXPECT_TRUE(std::signbit(z.v[1]));
+ EXPECT_TRUE(IsNaN(z.v[0]));
+ EXPECT_TRUE(IsNaN(z.v[1]));
+ }
+ { // copysign(+0, -0)
+ J z = copysign(MakeJet(+0, 1, 2), J{-0});
+ VL << "z = " << z;
+ EXPECT_FALSE(std::signbit(z.a));
+ EXPECT_TRUE(std::signbit(z.v[0]));
+ EXPECT_TRUE(std::signbit(z.v[1]));
+ EXPECT_TRUE(IsNaN(z.v[0]));
+ EXPECT_TRUE(IsNaN(z.v[1]));
+ }
+ { // copysign(-0, +0)
+ J z = copysign(MakeJet(-0, 1, 2), J{+0});
+ VL << "z = " << z;
+ EXPECT_FALSE(std::signbit(z.a));
+ EXPECT_TRUE(std::signbit(z.v[0]));
+ EXPECT_TRUE(std::signbit(z.v[1]));
+ EXPECT_TRUE(IsNaN(z.v[0]));
+ EXPECT_TRUE(IsNaN(z.v[1]));
+ }
+ { // copysign(-0, -0)
+ J z = copysign(MakeJet(-0, 1, 2), J{-0});
+ VL << "z = " << z;
+ EXPECT_FALSE(std::signbit(z.a));
+ EXPECT_TRUE(std::signbit(z.v[0]));
+ EXPECT_TRUE(std::signbit(z.v[1]));
+ EXPECT_TRUE(IsNaN(z.v[0]));
+ EXPECT_TRUE(IsNaN(z.v[1]));
+ }
}
TEST(Jet, JetsInEigenMatrices) {
