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.