diff --git a/include/boost/math/special_functions/detail/bessel_ik.hpp b/include/boost/math/special_functions/detail/bessel_ik.hpp
index 0c653b475..7a9903708 100644
--- a/include/boost/math/special_functions/detail/bessel_ik.hpp
+++ b/include/boost/math/special_functions/detail/bessel_ik.hpp
@@ -326,6 +326,9 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
     T scale = 1;
     T scale_sign = 1;
 
+    n = iround(v, pol);
+    u = v - n;                              // -1/2 <= u < 1/2
+
     if (((kind & need_i) == 0) && (fabs(4 * v * v - 25) / (8 * x) < tools::forth_root_epsilon<T>()))
     {
        // A&S 9.7.2
@@ -337,8 +340,6 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
     }
     else
     {
-       n = iround(v, pol);
-       u = v - n;                              // -1/2 <= u < 1/2
        BOOST_MATH_INSTRUMENT_VARIABLE(n);
        BOOST_MATH_INSTRUMENT_VARIABLE(u);
 
@@ -412,7 +413,7 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
        else
           Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do
     }
-    if (reflect)
+    if (reflect && (kind & need_i))
     {
         T z = (u + n % 2);
         T fact = (2 / pi<T>()) * (boost::math::sin_pi(z, pol) * Kv);
diff --git a/include/boost/math/special_functions/detail/bessel_j1.hpp b/include/boost/math/special_functions/detail/bessel_j1.hpp
index 6d354dcce..90f26bcc8 100644
--- a/include/boost/math/special_functions/detail/bessel_j1.hpp
+++ b/include/boost/math/special_functions/detail/bessel_j1.hpp
@@ -34,36 +34,9 @@ namespace boost { namespace math{  namespace detail{
 template <typename T>
 T bessel_j1(T x);
 
-template <class T>
-struct bessel_j1_initializer
-{
-   struct init
-   {
-      init()
-      {
-         do_init();
-      }
-      static void do_init()
-      {
-         bessel_j1(T(1));
-      }
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T>
-const typename bessel_j1_initializer<T>::init bessel_j1_initializer<T>::initializer;
-
 template <typename T>
 T bessel_j1(T x)
 {
-    bessel_j1_initializer<T>::force_instantiate();
-
     static const T P1[] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)),
diff --git a/include/boost/math/special_functions/detail/bessel_k0.hpp b/include/boost/math/special_functions/detail/bessel_k0.hpp
index f29ffa75c..69a48b779 100644
--- a/include/boost/math/special_functions/detail/bessel_k0.hpp
+++ b/include/boost/math/special_functions/detail/bessel_k0.hpp
@@ -46,40 +46,6 @@ namespace boost { namespace math { namespace detail{
 template <typename T>
 T bessel_k0(const T& x);
 
-template <class T, class tag>
-struct bessel_k0_initializer
-{
-   struct init
-   {
-      init()
-      {
-         do_init(tag());
-      }
-      static void do_init(const std::integral_constant<int, 113>&)
-      {
-         bessel_k0(T(0.5));
-         bessel_k0(T(1.5));
-      }
-      static void do_init(const std::integral_constant<int, 64>&)
-      {
-         bessel_k0(T(0.5));
-         bessel_k0(T(1.5));
-      }
-      template <class U>
-      static void do_init(const U&){}
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T, class tag>
-const typename bessel_k0_initializer<T, tag>::init bessel_k0_initializer<T, tag>::initializer;
-
-
 template <typename T, int N>
 T bessel_k0_imp(const T&, const std::integral_constant<int, N>&)
 {
@@ -505,7 +471,6 @@ inline T bessel_k0(const T& x)
       113 : -1
    > tag_type;
 
-   bessel_k0_initializer<T, tag_type>::force_instantiate();
    return bessel_k0_imp(x, tag_type());
 }
 
diff --git a/include/boost/math/special_functions/detail/bessel_k1.hpp b/include/boost/math/special_functions/detail/bessel_k1.hpp
index bd37f9021..b7259e467 100644
--- a/include/boost/math/special_functions/detail/bessel_k1.hpp
+++ b/include/boost/math/special_functions/detail/bessel_k1.hpp
@@ -46,41 +46,6 @@ namespace boost { namespace math { namespace detail{
    template <typename T>
    T bessel_k1(const T&);
 
-   template <class T, class tag>
-   struct bessel_k1_initializer
-   {
-      struct init
-      {
-         init()
-         {
-            do_init(tag());
-         }
-         static void do_init(const std::integral_constant<int, 113>&)
-         {
-            bessel_k1(T(0.5));
-            bessel_k1(T(2));
-            bessel_k1(T(6));
-         }
-         static void do_init(const std::integral_constant<int, 64>&)
-         {
-            bessel_k1(T(0.5));
-            bessel_k1(T(6));
-         }
-         template <class U>
-         static void do_init(const U&) {}
-         void force_instantiate()const {}
-      };
-      static const init initializer;
-      static void force_instantiate()
-      {
-         initializer.force_instantiate();
-      }
-   };
-
-   template <class T, class tag>
-   const typename bessel_k1_initializer<T, tag>::init bessel_k1_initializer<T, tag>::initializer;
-
-
    template <typename T, int N>
    inline T bessel_k1_imp(const T&, const std::integral_constant<int, N>&)
    {
@@ -547,7 +512,6 @@ namespace boost { namespace math { namespace detail{
          113 : -1
       > tag_type;
 
-      bessel_k1_initializer<T, tag_type>::force_instantiate();
       return bessel_k1_imp(x, tag_type());
    }
 
diff --git a/include/boost/math/special_functions/detail/polygamma.hpp b/include/boost/math/special_functions/detail/polygamma.hpp
index 8ca366314..a608f22c7 100644
--- a/include/boost/math/special_functions/detail/polygamma.hpp
+++ b/include/boost/math/special_functions/detail/polygamma.hpp
@@ -470,34 +470,12 @@ namespace boost { namespace math { namespace detail{
      return exp(power_terms) * ((s < 0) && ((n + 1) & 1) ? -1 : 1) * boost::math::sign(sum);
   }
 
-  template <class T, class Policy>
-  struct polygamma_initializer
-  {
-     struct init
-     {
-        init()
-        {
-           // Forces initialization of our table of coefficients and mutex:
-           boost::math::polygamma(30, T(-2.5f), Policy());
-        }
-        void force_instantiate()const{}
-     };
-     static const init initializer;
-     static void force_instantiate()
-     {
-        initializer.force_instantiate();
-     }
-  };
-
-  template <class T, class Policy>
-  const typename polygamma_initializer<T, Policy>::init polygamma_initializer<T, Policy>::initializer;
-
   template<class T, class Policy>
   inline T polygamma_imp(const int n, T x, const Policy &pol)
   {
     BOOST_MATH_STD_USING
     static const char* function = "boost::math::polygamma<%1%>(int, %1%)";
-    polygamma_initializer<T, Policy>::initializer.force_instantiate();
+
     if(n < 0)
        return policies::raise_domain_error<T>(function, "Order must be >= 0, but got %1%", static_cast<T>(n), pol);
     if(x < 0)
diff --git a/include/boost/math/special_functions/detail/unchecked_factorial.hpp b/include/boost/math/special_functions/detail/unchecked_factorial.hpp
index b528a24fe..8428f7400 100644
--- a/include/boost/math/special_functions/detail/unchecked_factorial.hpp
+++ b/include/boost/math/special_functions/detail/unchecked_factorial.hpp
@@ -1271,28 +1271,6 @@ struct max_factorial<BOOST_MATH_FLOAT128_TYPE>
 
 #endif
 
-template <class T>
-struct unchecked_factorial_initializer
-{
-   struct init
-   {
-      init()
-      {
-         boost::math::unchecked_factorial<T>(3);
-      }
-      void force_instantiate()const {}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T>
-const typename unchecked_factorial_initializer<T>::init unchecked_factorial_initializer<T>::initializer;
-
-
 template <class T, int N>
 inline T unchecked_factorial_imp(unsigned i, const std::integral_constant<int, N>&)
 {
diff --git a/include/boost/math/special_functions/log1p.hpp b/include/boost/math/special_functions/log1p.hpp
index 9b8a8e0eb..1fdd6dc52 100644
--- a/include/boost/math/special_functions/log1p.hpp
+++ b/include/boost/math/special_functions/log1p.hpp
@@ -263,34 +263,6 @@ T log1p_imp(T const& x, const Policy& pol, const std::integral_constant<int, 24>
    return result;
 }
 
-template <class T, class Policy, class tag>
-struct log1p_initializer
-{
-   struct init
-   {
-      init()
-      {
-         do_init(tag());
-      }
-      template <int N>
-      static void do_init(const std::integral_constant<int, N>&){}
-      static void do_init(const std::integral_constant<int, 64>&)
-      {
-         boost::math::log1p(static_cast<T>(0.25), Policy());
-      }
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T, class Policy, class tag>
-const typename log1p_initializer<T, Policy, tag>::init log1p_initializer<T, Policy, tag>::initializer;
-
-
 } // namespace detail
 
 template <class T, class Policy>
@@ -312,8 +284,6 @@ inline typename tools::promote_args<T>::type log1p(T x, const Policy&)
       precision_type::value <= 64 ? 64 : 0
    > tag_type;
 
-   detail::log1p_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
-
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::log1p_imp(static_cast<value_type>(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)");
 }
diff --git a/include/boost/math/special_functions/next.hpp b/include/boost/math/special_functions/next.hpp
index 02a208e4e..ee22ce78c 100644
--- a/include/boost/math/special_functions/next.hpp
+++ b/include/boost/math/special_functions/next.hpp
@@ -121,31 +121,6 @@ inline bool has_denorm_now() {
 template <class T>
 T get_min_shift_value();
 
-template <class T>
-struct min_shift_initializer
-{
-   struct init
-   {
-      init()
-      {
-         do_init();
-      }
-      static void do_init()
-      {
-         get_min_shift_value<T>();
-      }
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T>
-const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
-
 template <class T>
 inline T calc_min_shifted(const std::true_type&)
 {
@@ -166,8 +141,6 @@ template <class T>
 inline T get_min_shift_value()
 {
    static const T val = calc_min_shifted<T>(std::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
-   min_shift_initializer<T>::force_instantiate();
-
    return val;
 }
 
diff --git a/include/boost/math/special_functions/owens_t.hpp b/include/boost/math/special_functions/owens_t.hpp
index c794f0cac..e09608b39 100644
--- a/include/boost/math/special_functions/owens_t.hpp
+++ b/include/boost/math/special_functions/owens_t.hpp
@@ -1026,34 +1026,6 @@ namespace boost
             return val;
          } // RealType owens_t(RealType h, RealType a)
 
-         template <class T, class Policy, class tag>
-         struct owens_t_initializer
-         {
-            struct init
-            {
-               init()
-               {
-                  do_init(tag());
-               }
-               template <int N>
-               static void do_init(const std::integral_constant<int, N>&){}
-               static void do_init(const std::integral_constant<int, 64>&)
-               {
-                  boost::math::owens_t(static_cast<T>(7), static_cast<T>(0.96875), Policy());
-                  boost::math::owens_t(static_cast<T>(2), static_cast<T>(0.5), Policy());
-               }
-               void force_instantiate()const{}
-            };
-            static const init initializer;
-            static void force_instantiate()
-            {
-               initializer.force_instantiate();
-            }
-         };
-
-         template <class T, class Policy, class tag>
-         const typename owens_t_initializer<T, Policy, tag>::init owens_t_initializer<T, Policy, tag>::initializer;
-
       } // namespace detail
 
       template <class T1, class T2, class Policy>
@@ -1067,8 +1039,6 @@ namespace boost
             precision_type::value <= 64 ? 64 : 65
          > tag_type;
 
-         detail::owens_t_initializer<result_type, Policy, tag_type>::force_instantiate();
-            
          return policies::checked_narrowing_cast<result_type, Policy>(detail::owens_t(static_cast<value_type>(h), static_cast<value_type>(a), pol), "boost::math::owens_t<%1%>(%1%,%1%)");
       }
 
diff --git a/include/boost/math/special_functions/trigamma.hpp b/include/boost/math/special_functions/trigamma.hpp
index f74b43db1..d40be39cb 100644
--- a/include/boost/math/special_functions/trigamma.hpp
+++ b/include/boost/math/special_functions/trigamma.hpp
@@ -394,35 +394,6 @@ T trigamma_imp(T x, const std::integral_constant<int, 0>*, const Policy& pol)
 {
    return polygamma_imp(1, x, pol);
 }
-//
-// Initializer: ensure all our constants are initialized prior to the first call of main:
-//
-template <class T, class Policy>
-struct trigamma_initializer
-{
-   struct init
-   {
-      init()
-      {
-         typedef typename policies::precision<T, Policy>::type precision_type;
-         do_init(std::integral_constant<bool, precision_type::value && (precision_type::value <= 113)>());
-      }
-      void do_init(const std::true_type&)
-      {
-         boost::math::trigamma(T(2.5), Policy());
-      }
-      void do_init(const std::false_type&){}
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T, class Policy>
-const typename trigamma_initializer<T, Policy>::init trigamma_initializer<T, Policy>::initializer;
 
 } // namespace detail
 
@@ -446,9 +417,6 @@ inline typename tools::promote_args<T>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   // Force initialization of constants:
-   detail::trigamma_initializer<value_type, forwarding_policy>::force_instantiate();
-
    return policies::checked_narrowing_cast<result_type, Policy>(detail::trigamma_imp(
       static_cast<value_type>(x),
       static_cast<const tag_type*>(nullptr), forwarding_policy()), "boost::math::trigamma<%1%>(%1%)");
diff --git a/include/boost/math/special_functions/zeta.hpp b/include/boost/math/special_functions/zeta.hpp
index 1a3bba748..de74d3d05 100644
--- a/include/boost/math/special_functions/zeta.hpp
+++ b/include/boost/math/special_functions/zeta.hpp
@@ -1010,53 +1010,6 @@ T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
    return result;
 }
 
-template <class T, class Policy, class tag>
-struct zeta_initializer
-{
-   struct init
-   {
-      init()
-      {
-         do_init(tag());
-      }
-      static void do_init(const std::integral_constant<int, 0>&){ boost::math::zeta(static_cast<T>(5), Policy()); }
-      static void do_init(const std::integral_constant<int, 53>&){ boost::math::zeta(static_cast<T>(5), Policy()); }
-      static void do_init(const std::integral_constant<int, 64>&)
-      {
-         boost::math::zeta(static_cast<T>(0.5), Policy());
-         boost::math::zeta(static_cast<T>(1.5), Policy());
-         boost::math::zeta(static_cast<T>(3.5), Policy());
-         boost::math::zeta(static_cast<T>(6.5), Policy());
-         boost::math::zeta(static_cast<T>(14.5), Policy());
-         boost::math::zeta(static_cast<T>(40.5), Policy());
-
-         boost::math::zeta(static_cast<T>(5), Policy());
-      }
-      static void do_init(const std::integral_constant<int, 113>&)
-      {
-         boost::math::zeta(static_cast<T>(0.5), Policy());
-         boost::math::zeta(static_cast<T>(1.5), Policy());
-         boost::math::zeta(static_cast<T>(3.5), Policy());
-         boost::math::zeta(static_cast<T>(5.5), Policy());
-         boost::math::zeta(static_cast<T>(9.5), Policy());
-         boost::math::zeta(static_cast<T>(16.5), Policy());
-         boost::math::zeta(static_cast<T>(25.5), Policy());
-         boost::math::zeta(static_cast<T>(70.5), Policy());
-
-         boost::math::zeta(static_cast<T>(5), Policy());
-      }
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T, class Policy, class tag>
-const typename zeta_initializer<T, Policy, tag>::init zeta_initializer<T, Policy, tag>::initializer;
-
 } // detail
 
 template <class T, class Policy>
@@ -1078,8 +1031,6 @@ inline typename tools::promote_args<T>::type zeta(T s, const Policy&)
       precision_type::value <= 113 ? 113 : 0
    > tag_type;
 
-   detail::zeta_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
-
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp(
       static_cast<value_type>(s),
       static_cast<value_type>(1 - static_cast<value_type>(s)),