diff --git a/include/boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp
index 22491a70a..2064caa51 100644
--- a/include/boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp
+++ b/include/boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp
@@ -225,7 +225,7 @@ cardinal_cubic_b_spline_imp<Real>::cardinal_cubic_b_spline_imp(BidiIterator f, B
     // mapsto
     // 1 0 -1 | r0
     // 0 1 1/2| (r1 - r0)/4
-    super_diagonal[1] = 0.5;
+    super_diagonal[1] = static_cast<Real>(0.5);
     rhs[1] = (rhs[1] - rhs[0])/4;
 
     // Now do a tridiagonal row reduction the standard way, until just before the last row:
diff --git a/test/cardinal_cubic_b_spline_test.cpp b/test/cardinal_cubic_b_spline_test.cpp
index 9dcef343b..90452d6e5 100644
--- a/test/cardinal_cubic_b_spline_test.cpp
+++ b/test/cardinal_cubic_b_spline_test.cpp
@@ -14,8 +14,13 @@
 #include <boost/test/tools/floating_point_comparison.hpp>
 #include <boost/math/interpolators/cardinal_cubic_b_spline.hpp>
 #include <boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp>
+#include <boost/math/constants/constants.hpp>
 #include <boost/multiprecision/cpp_bin_float.hpp>
 
+#if __has_include(<stdfloat>)
+#  include <stdfloat>
+#endif
+
 using boost::multiprecision::cpp_bin_float_50;
 using boost::math::constants::third;
 using boost::math::constants::half;
@@ -26,12 +31,12 @@ void test_b3_spline()
     std::cout << "Testing evaluation of spline basis functions on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
     // Outside the support:
     Real eps = std::numeric_limits<Real>::epsilon();
-    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline<Real>(2.5), (Real) 0);
-    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline<Real>(-2.5), (Real) 0);
-    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_prime<Real>(2.5), (Real) 0);
-    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_prime<Real>(-2.5), (Real) 0);
-    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_double_prime<Real>(2.5), (Real) 0);
-    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_double_prime<Real>(-2.5), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline<Real>(Real(2.5)), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline<Real>(Real(-2.5)), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_prime<Real>(Real(2.5)), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_prime<Real>(Real(-2.5)), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_double_prime<Real>(Real(2.5)), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::interpolators::detail::b3_spline_double_prime<Real>(Real(-2.5)), (Real) 0);
 
 
     // On the boundary of support:
@@ -52,7 +57,7 @@ void test_b3_spline()
     // Properties: B3 is an even function, B3' is an odd function.
     for (size_t i = 1; i < 200; ++i)
     {
-        Real arg = i*0.01;
+        Real arg = i*Real(0.01);
         BOOST_CHECK_CLOSE(boost::math::interpolators::detail::b3_spline<Real>(arg), boost::math::interpolators::detail::b3_spline<Real>(arg), eps);
         BOOST_CHECK_CLOSE(boost::math::interpolators::detail::b3_spline_prime<Real>(-arg), -boost::math::interpolators::detail::b3_spline_prime<Real>(arg), eps);
         BOOST_CHECK_CLOSE(boost::math::interpolators::detail::b3_spline_double_prime<Real>(-arg), boost::math::interpolators::detail::b3_spline_double_prime<Real>(arg), eps);
@@ -76,8 +81,8 @@ void test_interpolation_condition()
         v[i] = dis(gen);
     }
 
-    Real step = 0.01;
-    Real a = 5;
+    Real step = static_cast<Real>(0.01);
+    Real a = Real(5);
     boost::math::interpolators::cardinal_cubic_b_spline<Real> spline(v.data(), v.size(), a, step);
 
     for (size_t i = 0; i < v.size(); ++i)
@@ -99,21 +104,21 @@ void test_constant_function()
     Real constant = 50.2;
     for (size_t i = 0; i < v.size(); ++i)
     {
-        v[i] = 50.2;
+        v[i] = Real(50.2);
     }
 
-    Real step = 0.02;
+    Real step = static_cast<Real>(0.02);
     Real a = 5;
     boost::math::interpolators::cardinal_cubic_b_spline<Real> spline(v.data(), v.size(), a, step);
 
     for (size_t i = 0; i < v.size(); ++i)
     {
         // Do not test at interpolation point; we already know it works there:
-        Real y = spline(i*step + a + 0.001);
+        Real y = spline(i*step + a + Real(0.001));
         BOOST_CHECK_CLOSE(y, constant, 10*std::numeric_limits<Real>::epsilon());
-        Real y_prime = spline.prime(i*step + a + 0.002);
+        Real y_prime = spline.prime(i*step + a + Real(0.002));
         BOOST_CHECK_SMALL(y_prime, 5000*std::numeric_limits<Real>::epsilon());
-        Real y_double_prime = spline.double_prime(i*step + a + 0.002);
+        Real y_double_prime = spline.double_prime(i*step + a + Real(0.002));
         BOOST_CHECK_SMALL(y_double_prime, 5000*std::numeric_limits<Real>::epsilon());
 
     }
@@ -123,9 +128,9 @@ void test_constant_function()
 
     for (size_t i = 0; i < v.size(); ++i)
     {
-        Real y = spline(i*step + a + 0.002);
+        Real y = spline(i*step + a + Real(0.002));
         BOOST_CHECK_CLOSE(y, constant, std::numeric_limits<Real>::epsilon());
-        Real y_prime = spline.prime(i*step + a + 0.002);
+        Real y_prime = spline.prime(i*step + a + Real(0.002));
         BOOST_CHECK_SMALL(y_prime, std::numeric_limits<Real>::epsilon());
     }
 
@@ -137,9 +142,9 @@ void test_constant_function()
     for (size_t i = 0; i < v.size(); ++i)
     {
        // Do not test at interpolation point; we already know it works there:
-       Real y = spline2(i*step + a + 0.001);
+       Real y = spline2(i*step + a + Real(0.001));
        BOOST_CHECK_CLOSE(y, constant, 10 * std::numeric_limits<Real>::epsilon());
-       Real y_prime = spline2.prime(i*step + a + 0.002);
+       Real y_prime = spline2.prime(i*step + a + Real(0.002));
        BOOST_CHECK_SMALL(y_prime, 5000 * std::numeric_limits<Real>::epsilon());
     }
 
@@ -148,9 +153,9 @@ void test_constant_function()
 
     for (size_t i = 0; i < v.size(); ++i)
     {
-       Real y = spline2(i*step + a + 0.002);
+       Real y = spline2(i*step + a + Real(0.002));
        BOOST_CHECK_CLOSE(y, constant, std::numeric_limits<Real>::epsilon());
-       Real y_prime = spline2.prime(i*step + a + 0.002);
+       Real y_prime = spline2.prime(i*step + a + Real(0.002));
        BOOST_CHECK_SMALL(y_prime, std::numeric_limits<Real>::epsilon());
     }
 }
@@ -164,7 +169,7 @@ void test_affine_function()
     std::vector<Real> v(500);
     Real a = 10;
     Real b = 8;
-    Real step = 0.005;
+    Real step = static_cast<Real>(0.005);
 
     auto f = [a, b](Real x) { return a*x + b; };
     for (size_t i = 0; i < v.size(); ++i)
@@ -176,7 +181,7 @@ void test_affine_function()
 
     for (size_t i = 0; i < v.size() - 1; ++i)
     {
-        Real arg = i*step + 0.0001;
+        Real arg = static_cast<Real>(i*step + Real(0.0001));
         Real y = spline(arg);
         BOOST_CHECK_CLOSE(y, f(arg), sqrt(std::numeric_limits<Real>::epsilon()));
         Real y_prime = spline.prime(arg);
@@ -188,7 +193,7 @@ void test_affine_function()
 
     for (size_t i = 0; i < v.size() - 1; ++i)
     {
-        Real arg = i*step + 0.0001;
+        Real arg = static_cast<Real>(i*step + Real(0.0001));
         Real y = spline(arg);
         BOOST_CHECK_CLOSE(y, f(arg), sqrt(std::numeric_limits<Real>::epsilon()));
         Real y_prime = spline.prime(arg);
@@ -203,10 +208,10 @@ void test_quadratic_function()
     using std::sqrt;
     std::cout << "Testing that quadratic functions are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
     std::vector<Real> v(500);
-    Real a = 1.2;
-    Real b = -3.4;
-    Real c = -8.6;
-    Real step = 0.01;
+    Real a = static_cast<Real>(1.2);
+    Real b = static_cast<Real>(-3.4);
+    Real c = static_cast<Real>(-8.6);
+    Real step = static_cast<Real>(0.01);
 
     auto f = [a, b, c](Real x) { return a*x*x + b*x + c; };
     for (size_t i = 0; i < v.size(); ++i)
@@ -218,11 +223,11 @@ void test_quadratic_function()
 
     for (size_t i = 0; i < v.size() -1; ++i)
     {
-        Real arg = i*step + 0.001;
+        Real arg = static_cast<Real>(i*step + Real(0.001));
         Real y = spline(arg);
-        BOOST_CHECK_CLOSE(y, f(arg), 0.1);
+        BOOST_CHECK_CLOSE(y, f(arg), Real(0.1));
         Real y_prime = spline.prime(arg);
-        BOOST_CHECK_CLOSE(y_prime, 2*a*arg + b, 2.0);
+        BOOST_CHECK_CLOSE(y_prime, 2*a*arg + b, Real(2.0));
     }
 }
 
@@ -233,9 +238,9 @@ void test_circ_conic_function()
     using std::sqrt;
     std::cout << "Testing that conic section of a circle is interpolated correctly by cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name() << '\n';
     std::vector<Real> v(500);
-    Real cv = 0.1;
-    Real w = 2.0;
-    Real step = 2 * w / (v.size() - 1);
+    Real cv = Real(0.1);
+    Real w = Real(2.0);
+    Real step = Real(2 * w / (v.size() - 1));
 
     auto f = [cv](Real x) { return cv * x * x / (1 + sqrt(1 - cv * cv * x * x)); };
     auto df = [cv](Real x) { return cv * x / sqrt(1 - cv * cv * x * x); };
@@ -256,11 +261,11 @@ void test_circ_conic_function()
 
     for (size_t i = 0; i < v.size() - 1; ++i)
     {
-        Real arg = -w + i * step + 0.001;
+        Real arg = -w + i * step + Real(0.001);
         Real y = spline(arg);
-        BOOST_CHECK_CLOSE(y, f(arg), 2.0);
+        BOOST_CHECK_CLOSE(y, f(arg), Real(2.0));
         Real y_prime = spline.prime(arg);
-        BOOST_CHECK_CLOSE(y_prime, df(arg), 1.0);
+        BOOST_CHECK_CLOSE(y_prime, df(arg), Real(1.0));
     }
 }
 
@@ -272,7 +277,7 @@ void test_trig_function()
     std::mt19937 gen;
     std::vector<Real> v(500);
     Real x0 = 1;
-    Real step = 0.125;
+    Real step = Real(0.125);
 
     for (size_t i = 0; i < v.size(); ++i)
     {
@@ -287,9 +292,9 @@ void test_trig_function()
     {
         Real x = abscissa(gen);
         Real y = spline(x);
-        BOOST_CHECK_CLOSE(y, sin(x), 1.0);
+        BOOST_CHECK_CLOSE(y, Real(sin(x)), Real(1.0));
         auto y_prime = spline.prime(x);
-        BOOST_CHECK_CLOSE(y_prime, cos(x), 2.0);
+        BOOST_CHECK_CLOSE(y_prime, Real(cos(x)), Real(2.0));
     }
 }
 
@@ -299,7 +304,8 @@ void test_copy_move()
     std::cout << "Testing that copy/move operation succeed on cubic b spline\n";
     std::vector<Real> v(500);
     Real x0 = 1;
-    Real step = 0.125;
+    Real step = static_cast<Real>(0.125);
+    constexpr Real tol = Real(0.01);
 
     for (size_t i = 0; i < v.size(); ++i)
     {
@@ -312,22 +318,22 @@ void test_copy_move()
     // Default constructor should compile so that splines can be member variables:
     boost::math::interpolators::cardinal_cubic_b_spline<Real> d;
     d = boost::math::interpolators::cardinal_cubic_b_spline<Real>(v.data(), v.size(), x0, step);
-    BOOST_CHECK_CLOSE(d(x0), sin(x0), 0.01);
+    BOOST_CHECK_CLOSE(d(x0), Real(sin(x0)), tol);
     // Passing to lambda should compile:
     auto f = [=](Real x) { return d(x); };
     // Make sure this variable is used.
-    BOOST_CHECK_CLOSE(f(x0), sin(x0), 0.01);
+    BOOST_CHECK_CLOSE(f(x0), Real(sin(x0)), tol);
 
     // Move operations should compile.
     auto s = std::move(spline);
 
     // Copy operations should compile:
     boost::math::interpolators::cardinal_cubic_b_spline<Real> c = d;
-    BOOST_CHECK_CLOSE(c(x0), sin(x0), 0.01);
+    BOOST_CHECK_CLOSE(c(x0), Real(sin(x0)), tol);
 
     // Test with std::bind:
-    auto h = std::bind(&boost::math::interpolators::cardinal_cubic_b_spline<double>::operator(), &s, std::placeholders::_1);
-    BOOST_CHECK_CLOSE(h(x0), sin(x0), 0.01);
+    auto h = std::bind(&boost::math::interpolators::cardinal_cubic_b_spline<Real>::operator(), &s, std::placeholders::_1);
+    BOOST_CHECK_CLOSE(h(x0), Real(sin(x0)), tol);
 }
 
 template<class Real>
@@ -336,7 +342,7 @@ void test_outside_interval()
     std::cout << "Testing that the spline can be evaluated outside the interpolation interval\n";
     std::vector<Real> v(400);
     Real x0 = 1;
-    Real step = 0.125;
+    Real step = Real(0.125);
 
     for (size_t i = 0; i < v.size(); ++i)
     {
@@ -354,54 +360,70 @@ void test_outside_interval()
 
 BOOST_AUTO_TEST_CASE(test_cubic_b_spline)
 {
+    #ifdef __STDCPP_FLOAT32_T__
+
+    test_b3_spline<std::float32_t>();
+    test_interpolation_condition<std::float32_t>();
+    test_constant_function<std::float32_t>();
+    test_affine_function<std::float32_t>();
+    test_quadratic_function<std::float32_t>();
+    test_circ_conic_function<std::float32_t>();
+    test_trig_function<std::float32_t>();
+
+    #else
+    
     test_b3_spline<float>();
-    test_b3_spline<double>();
-#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-    test_b3_spline<long double>();
-#endif
-    test_b3_spline<cpp_bin_float_50>();
-
     test_interpolation_condition<float>();
-    test_interpolation_condition<double>();
-#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-    test_interpolation_condition<long double>();
-#endif
-    test_interpolation_condition<cpp_bin_float_50>();
-
     test_constant_function<float>();
-    test_constant_function<double>();
-#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-    test_constant_function<long double>();
-#endif
-    test_constant_function<cpp_bin_float_50>();
-
     test_affine_function<float>();
-    test_affine_function<double>();
-#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-    test_affine_function<long double>();
-#endif
-    test_affine_function<cpp_bin_float_50>();
-
     test_quadratic_function<float>();
-    test_quadratic_function<double>();
-#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-    test_quadratic_function<long double>();
-#endif
-    test_quadratic_function<cpp_bin_float_50>();
-
     test_circ_conic_function<float>();
-    test_circ_conic_function<double>();
-#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-    test_circ_conic_function<long double>();
-#endif
-
     test_trig_function<float>();
-    test_trig_function<double>();
-#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-    test_trig_function<long double>();
-#endif
-    test_trig_function<cpp_bin_float_50>();
 
+    #endif
+
+    #ifdef __STDCPP_FLOAT64_T__
+
+    test_b3_spline<std::float64_t>();
+    test_interpolation_condition<std::float64_t>();
+    test_constant_function<std::float64_t>();
+    test_affine_function<std::float64_t>();
+    test_quadratic_function<std::float64_t>();
+    test_circ_conic_function<std::float64_t>();
+    test_trig_function<std::float64_t>();
+    test_copy_move<std::float64_t>();
+    test_outside_interval<std::float64_t>();
+
+    #else
+
+    test_b3_spline<double>();
+    test_interpolation_condition<double>();
+    test_constant_function<double>();
+    test_affine_function<double>();
+    test_quadratic_function<double>();
+    test_circ_conic_function<double>();
+    test_trig_function<double>();
     test_copy_move<double>();
     test_outside_interval<double>();
+
+    #endif
+
+    #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+
+    test_b3_spline<long double>();
+    test_interpolation_condition<long double>();
+    test_constant_function<long double>();
+    test_affine_function<long double>();
+    test_quadratic_function<long double>();
+    test_circ_conic_function<long double>();
+    test_trig_function<long double>();
+
+    #endif
+
+    test_b3_spline<cpp_bin_float_50>();
+    test_interpolation_condition<cpp_bin_float_50>();
+    test_constant_function<cpp_bin_float_50>();
+    test_affine_function<cpp_bin_float_50>();
+    test_quadratic_function<cpp_bin_float_50>();
+    test_trig_function<cpp_bin_float_50>();
 }