diff --git a/include/boost/math/differentiation/finite_difference.hpp b/include/boost/math/differentiation/finite_difference.hpp
index 5e4433e00..1375bac7a 100644
--- a/include/boost/math/differentiation/finite_difference.hpp
+++ b/include/boost/math/differentiation/finite_difference.hpp
@@ -153,7 +153,7 @@ namespace detail {
 
       const Real eps = (numeric_limits<Real>::epsilon)();
       // Error bound ~eps^4/5
-      Real h = pow(11.25*eps, static_cast<Real>(1) / static_cast<Real>(5));
+      Real h = pow(Real(11.25)*eps, static_cast<Real>(1) / static_cast<Real>(5));
       h = detail::make_xph_representable(x, h);
       Real ymth = f(x - 2 * h);
       Real yth = f(x + 2 * h);
@@ -222,7 +222,7 @@ namespace detail {
       // Mathematica code to get the error:
       // Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}]
       // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h.
-      Real h = pow(551.25*eps, static_cast<Real>(1) / static_cast<Real>(9));
+      Real h = pow(Real(551.25)*eps, static_cast<Real>(1) / static_cast<Real>(9));
       h = detail::make_xph_representable(x, h);
 
       Real yh = f(x + h);
diff --git a/test/test_numerical_differentiation.cpp b/test/test_numerical_differentiation.cpp
index 1b2f53962..0f2489a61 100644
--- a/test/test_numerical_differentiation.cpp
+++ b/test/test_numerical_differentiation.cpp
@@ -16,6 +16,10 @@
 #include <boost/math/special_functions/next.hpp>
 #include <boost/math/differentiation/finite_difference.hpp>
 
+#if __has_include(<stdfloat>)
+#  include <stdfloat>
+#endif
+
 using std::abs;
 using std::pow;
 using boost::math::differentiation::finite_difference_derivative;
@@ -31,7 +35,7 @@ void test_order(size_t points_to_test)
     std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
     //std::cout << std::fixed << std::scientific;
     auto f = [](Real t) { return boost::math::cyl_bessel_j<Real>(1, t); };
-    Real min = -100000.0;
+    Real min = Real(-100000.0);
     Real max = -min;
     Real x = min;
     Real max_error = 0;
@@ -95,7 +99,7 @@ void test_bessel()
 
     Real computed = finite_difference_derivative<decltype(f), Real, 1>(f, x);
     Real expected = cyl_bessel_j_prime(12, x);
-    Real error_estimate = 4*abs(f(x))*sqrt(eps);
+    Real error_estimate = Real(4*abs(f(x))*sqrt(eps));
     //std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
     //std::cout << "cyl_bessel_j_prime: " << expected << std::endl;
     //std::cout << "First order fd    : " << computed << std::endl;
@@ -218,28 +222,41 @@ void test_complex_step()
 
 BOOST_AUTO_TEST_CASE(numerical_differentiation_test)
 {
+    constexpr size_t points_to_test = 1000;
+    
+    #ifdef __STDCPP_FLOAT32_T__
+    test_complex_step<std::float32_t>();
+    test_bessel<std::float32_t>();
+    test_order<std::float32_t, 1>(points_to_test);
+    test_order<std::float32_t, 2>(points_to_test);
+    test_order<std::float32_t, 4>(points_to_test);
+    test_order<std::float32_t, 6>(points_to_test);
+    test_order<std::float32_t, 8>(points_to_test);
+    #else
     test_complex_step<float>();
-    test_complex_step<double>();
-
     test_bessel<float>();
-    test_bessel<double>();
-
-
-    size_t points_to_test = 1000;
     test_order<float, 1>(points_to_test);
-    test_order<double, 1>(points_to_test);
-
-
     test_order<float, 2>(points_to_test);
-    test_order<double, 2>(points_to_test);
-
     test_order<float, 4>(points_to_test);
-    test_order<double, 4>(points_to_test);
-
     test_order<float, 6>(points_to_test);
-    test_order<double, 6>(points_to_test);
-
     test_order<float, 8>(points_to_test);
-    test_order<double, 8>(points_to_test);
+    #endif
 
+    #ifdef __STDCPP_FLOAT64_T__
+    test_complex_step<std::float64_t>();
+    test_bessel<std::float64_t>();
+    test_order<std::float64_t, 1>(points_to_test);
+    test_order<std::float64_t, 2>(points_to_test);
+    test_order<std::float64_t, 4>(points_to_test);
+    test_order<std::float64_t, 6>(points_to_test);
+    test_order<std::float64_t, 8>(points_to_test);
+    #else
+    test_complex_step<double>();
+    test_bessel<double>();
+    test_order<double, 1>(points_to_test);   
+    test_order<double, 2>(points_to_test);
+    test_order<double, 4>(points_to_test);
+    test_order<double, 6>(points_to_test);
+    test_order<double, 8>(points_to_test);
+    #endif
 }