Use C++17 Bessel functions
Move Bessel functions availability checks from configuration time to
inclusion time to be more robust and allow the use of ABI compatible
compilers (e.g., Ceres is compiled using Clang but is used in a project
compiled using GCC.)
Since libc++ does not yet implement special math functions, we fallback
to their POSIX implementation if available. However, then only the
deprecated BesselJ{0,1,n} are provided.
Fixes #814
Change-Id: Ic3e62452b36e90cb22644cc8e553e3dd1881193f
diff --git a/cmake/config.h.in b/cmake/config.h.in
index 0e7fd3a..1566795 100644
--- a/cmake/config.h.in
+++ b/cmake/config.h.in
@@ -76,10 +76,6 @@
// If defined Ceres was compiled without support for METIS via Eigen.
@CERES_NO_EIGEN_METIS@
-// If defined, Ceres was compiled with a version MSVC >= 2005 which
-// deprecated the standard POSIX names for bessel functions, replacing them
-// with underscore prefixed versions (e.g. j0() -> _j0()).
-@CERES_MSVC_USE_UNDERSCORE_PREFIXED_BESSEL_FUNCTIONS@
// CERES_NO_SPARSE should be automatically defined by config.h if Ceres was
// compiled without any sparse back-end. Verify that it has not subsequently
diff --git a/include/ceres/internal/port.h b/include/ceres/internal/port.h
index b8cb0ff..d78ed51 100644
--- a/include/ceres/internal/port.h
+++ b/include/ceres/internal/port.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2023 Google Inc. All rights reserved.
+// Copyright 2024 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -31,6 +31,8 @@
#ifndef CERES_PUBLIC_INTERNAL_PORT_H_
#define CERES_PUBLIC_INTERNAL_PORT_H_
+#include <cmath> // Necessary for __cpp_lib_math_special_functions feature test
+
// A macro to mark a function/variable/class as deprecated.
// We use compiler specific attributes rather than the c++
// attribute because they do not mix well with each other.
@@ -88,4 +90,22 @@
#define CERES_RESTORE_DEPRECATED_WARNING
#endif // defined(_MSC_VER)
+#if defined(__cpp_lib_math_special_functions) && \
+ ((__cpp_lib_math_special_functions >= 201603L) || \
+ defined(__STDCPP_MATH_SPEC_FUNCS__) && \
+ (__STDCPP_MATH_SPEC_FUNCS__ >= 201003L))
+// If defined, indicates whether C++17 Bessel functions (of the first kind) are
+// available. Some standard library implementations, such as libc++ (Android
+// NDK, Apple, Clang) do not yet provide these functions. Implementations that
+// do not support C++17, but support ISO 29124:2010, provide the functions if
+// __STDCPP_MATH_SPEC_FUNCS__ is defined by the implementation to a value at
+// least 201003L and if the user defines __STDCPP_WANT_MATH_SPEC_FUNCS__ before
+// including any standard library headers. Standard library Bessel functions are
+// preferred over any other implementation.
+#define CERES_HAS_CPP17_BESSEL_FUNCTIONS
+#elif defined(_SVID_SOURCE) || defined(_BSD_SOURCE) || defined(_XOPEN_SOURCE)
+// If defined, indicates that j0, j1, and jn from <math.h> are available.
+#define CERES_HAS_POSIX_BESSEL_FUNCTIONS
+#endif
+
#endif // CERES_PUBLIC_INTERNAL_PORT_H_
diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index f279ba3..3b7f23f 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2023 Google Inc. All rights reserved.
+// Copyright 2024 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -430,6 +430,9 @@
using std::copysign;
using std::cos;
using std::cosh;
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+using std::cyl_bessel_j;
+#endif // CERES_HAS_CPP17_BESSEL_FUNCTIONS
using std::erf;
using std::erfc;
using std::exp;
@@ -867,6 +870,9 @@
-x.v * exp(-x.a * x.a) * (T(1) / sqrt(atan(T(1)))));
}
+#if defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS) || \
+ defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS)
+
// 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
@@ -874,19 +880,33 @@
// function errors in client code (the specific warning is suppressed when
// Ceres itself is built).
inline double BesselJ0(double x) {
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+ return cyl_bessel_j(0, x);
+#else
CERES_DISABLE_DEPRECATED_WARNING
return j0(x);
CERES_RESTORE_DEPRECATED_WARNING
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
+
inline double BesselJ1(double x) {
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+ return cyl_bessel_j(1, x);
+#else
CERES_DISABLE_DEPRECATED_WARNING
return j1(x);
CERES_RESTORE_DEPRECATED_WARNING
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
+
inline double BesselJn(int n, double x) {
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+ return cyl_bessel_j(static_cast<double>(n), x);
+#else
CERES_DISABLE_DEPRECATED_WARNING
return jn(n, x);
CERES_RESTORE_DEPRECATED_WARNING
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
// For the formulae of the derivatives of the Bessel functions see the book:
@@ -899,26 +919,60 @@
// j0(a + h) ~= j0(a) - j1(a) h
template <typename T, int N>
inline Jet<T, N> BesselJ0(const Jet<T, N>& f) {
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+ return cyl_bessel_j(0, f);
+#else
return Jet<T, N>(BesselJ0(f.a), -BesselJ1(f.a) * f.v);
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
// See formula http://dlmf.nist.gov/10.6#E1
// j1(a + h) ~= j1(a) + 0.5 ( j0(a) - j2(a) ) h
template <typename T, int N>
inline Jet<T, N> BesselJ1(const Jet<T, N>& f) {
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+ return cyl_bessel_j(1, f);
+#else
return Jet<T, N>(BesselJ1(f.a),
T(0.5) * (BesselJ0(f.a) - BesselJn(2, f.a)) * f.v);
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
// See formula http://dlmf.nist.gov/10.6#E1
// j_n(a + h) ~= j_n(a) + 0.5 ( j_{n-1}(a) - j_{n+1}(a) ) h
template <typename T, int N>
inline Jet<T, N> BesselJn(int n, const Jet<T, N>& f) {
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+ return cyl_bessel_j(n, f);
+#else
return Jet<T, N>(
BesselJn(n, f.a),
T(0.5) * (BesselJn(n - 1, f.a) - BesselJn(n + 1, f.a)) * f.v);
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS) ||
+ // defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS)
+
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+
+// See formula http://dlmf.nist.gov/10.6#E1
+// j_n(a + h) ~= j_n(a) + 0.5 ( j_{n-1}(a) - j_{n+1}(a) ) h
+template <typename T, int N>
+inline Jet<T, N> cyl_bessel_j(double v, const Jet<T, N>& f) {
+ // See formula http://dlmf.nist.gov/10.6#E3
+ // j0(a + h) ~= j0(a) - j1(a) h
+ if (fpclassify(v) == FP_ZERO) {
+ return Jet<T, N>(cyl_bessel_j(0, f.a), -cyl_bessel_j(1, f.a) * f.v);
+ }
+
+ return Jet<T, N>(
+ cyl_bessel_j(v, f.a),
+ T(0.5) * (cyl_bessel_j(v - 1, f.a) - cyl_bessel_j(v + 1, f.a)) * f.v);
+}
+
+#endif // CERES_HAS_CPP17_BESSEL_FUNCTIONS
+
// Classification and comparison functionality referencing only the scalar part
// of a Jet. To classify the derivatives (e.g., for sanity checks), the dual
// part should be referenced explicitly. For instance, to check whether the
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index a006648..a920799 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2023 Google Inc. All rights reserved.
+// Copyright 2024 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
@@ -34,6 +34,7 @@
#include <algorithm>
#include <cfenv>
#include <cmath>
+#include <stdexcept>
#include "ceres/stringprintf.h"
#include "ceres/test_util.h"
@@ -75,7 +76,7 @@
return z;
}
-double const kTolerance = 1e-13;
+constexpr double kTolerance = 1e-13;
// Stores the floating-point environment containing active floating-point
// exceptions, rounding mode, etc., and restores it upon destruction.
@@ -115,16 +116,25 @@
return std::islessequal(relative_difference, max_abs_relative_difference);
}
-MATCHER_P(IsAlmostEqualTo, y, "") {
- const bool result = (AreAlmostEqual(arg.a, y.a, kTolerance) &&
- AreAlmostEqual(arg.v[0], y.v[0], kTolerance) &&
- AreAlmostEqual(arg.v[1], y.v[1], kTolerance));
+MATCHER_P2(IsAlmostEqualToWithTolerance,
+ y,
+ tolerance,
+ "is almost equal to " + testing::PrintToString(y) +
+ " with tolerance " + testing::PrintToString(tolerance)) {
+ const bool result = (AreAlmostEqual(arg.a, y.a, tolerance) &&
+ AreAlmostEqual(arg.v[0], y.v[0], tolerance) &&
+ AreAlmostEqual(arg.v[1], y.v[1], tolerance));
if (!result) {
*result_listener << "\nexpected - actual : " << y - arg;
}
return result;
}
+MATCHER_P(IsAlmostEqualTo, y, "") {
+ return ExplainMatchResult(
+ IsAlmostEqualToWithTolerance(y, kTolerance), arg, result_listener);
+}
+
const double kStep = 1e-8;
const double kNumericalTolerance = 1e-6; // Numeric derivation is quite inexact
@@ -269,23 +279,52 @@
}
}
+#if defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS) || \
+ defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
TEST(Jet, Bessel) {
J zero = J(0.0);
+ J z = MakeJet(0.1, -2.7, 1e-3);
+#ifdef CERES_HAS_POSIX_BESSEL_FUNCTIONS
EXPECT_THAT(BesselJ0(zero), IsAlmostEqualTo(J(1.0)));
EXPECT_THAT(BesselJ1(zero), IsAlmostEqualTo(zero));
EXPECT_THAT(BesselJn(2, zero), IsAlmostEqualTo(zero));
EXPECT_THAT(BesselJn(3, zero), IsAlmostEqualTo(zero));
- J z = MakeJet(0.1, -2.7, 1e-3);
-
EXPECT_THAT(BesselJ0(z), IsAlmostEqualTo(BesselJn(0, z)));
EXPECT_THAT(BesselJ1(z), IsAlmostEqualTo(BesselJn(1, z)));
- // See formula http://dlmf.nist.gov/10.6.E1
+ // See formula http://dlmf.nist.gov/10.6.E1
EXPECT_THAT(BesselJ0(z) + BesselJn(2, z),
IsAlmostEqualTo((2.0 / z) * BesselJ1(z)));
+#endif // CERES_HAS_POSIX_BESSEL_FUNCTIONS
+
+#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
+ EXPECT_THAT(cyl_bessel_j(0, zero), IsAlmostEqualTo(J(1.0)));
+ EXPECT_THAT(cyl_bessel_j(1, zero), IsAlmostEqualTo(zero));
+ EXPECT_THAT(cyl_bessel_j(2, zero), IsAlmostEqualTo(zero));
+ EXPECT_THAT(cyl_bessel_j(3, zero), IsAlmostEqualTo(zero));
+
+ EXPECT_THAT(cyl_bessel_j(0, z), IsAlmostEqualTo(BesselJn(0, z)));
+ EXPECT_THAT(cyl_bessel_j(1, z), IsAlmostEqualTo(BesselJn(1, z)));
+
+ // MSVC Bessel functions and their derivatives produce errors slightly above
+ // kTolerance. Provide an alternative variant with a relaxed threshold.
+ constexpr double kRelaxedTolerance = 10 * kTolerance;
+
+ // See formula http://dlmf.nist.gov/10.6.E1
+ EXPECT_THAT(cyl_bessel_j(0, z) + cyl_bessel_j(2, z),
+ IsAlmostEqualToWithTolerance((2.0 / z) * cyl_bessel_j(1, z),
+ kRelaxedTolerance));
+
+ // MSVC does not throw an exception on invalid first argument
+#ifndef _MSC_VER
+ EXPECT_THROW(cyl_bessel_j(-1, zero), std::domain_error);
+#endif // defined(_MSC_VER)
+#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
+#endif // defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS) ||
+ // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
TEST(Jet, Floor) {
{ // floor of a positive number works.
