Imported Upstream version 1.51.0upstream/1.51.0

28 files changed, 2135 insertions, 87 deletions

diff --git a/boost/math/special_functions/detail/bessel_i0.hpp b/boost/math/special_functions/detail/bessel_i0.hpp
index 2c129facc7..7dc65d1a1b 100644
--- a/boost/math/special_functions/detail/bessel_i0.hpp
+++ b/boost/math/special_functions/detail/bessel_i0.hpp
@@ -21,8 +21,38 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T>
+T bessel_i0(T x);
+
+template <class T>
+struct bessel_i0_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         bessel_i0(T(1));
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T>
+const typename bessel_i0_initializer<T>::init bessel_i0_initializer<T>::initializer;
+
+template <typename T>
 T bessel_i0(T x)
 {
+    bessel_i0_initializer<T>::force_instantiate();
+
     static const T P1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375249e+15)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5050369673018427753e+14)),
diff --git a/boost/math/special_functions/detail/bessel_i1.hpp b/boost/math/special_functions/detail/bessel_i1.hpp
index aa4596cfd4..47f1b79883 100644
--- a/boost/math/special_functions/detail/bessel_i1.hpp
+++ b/boost/math/special_functions/detail/bessel_i1.hpp
@@ -21,8 +21,39 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T>
+T bessel_i1(T x);
+
+template <class T>
+struct bessel_i1_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         bessel_i1(T(1));
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T>
+const typename bessel_i1_initializer<T>::init bessel_i1_initializer<T>::initializer;
+
+template <typename T>
 T bessel_i1(T x)
 {
+
+    bessel_i1_initializer<T>::force_instantiate();
+
     static const T P1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4577180278143463643e+15)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7732037840791591320e+14)),
diff --git a/boost/math/special_functions/detail/bessel_j0.hpp b/boost/math/special_functions/detail/bessel_j0.hpp
index ee25d46f61..a07052d73e 100644
--- a/boost/math/special_functions/detail/bessel_j0.hpp
+++ b/boost/math/special_functions/detail/bessel_j0.hpp
@@ -22,8 +22,38 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T>
+T bessel_j0(T x);
+
+template <class T>
+struct bessel_j0_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         bessel_j0(T(1));
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T>
+const typename bessel_j0_initializer<T>::init bessel_j0_initializer<T>::initializer;
+
+template <typename T>
 T bessel_j0(T x)
 {
+    bessel_j0_initializer<T>::force_instantiate();
+    
     static const T P1[] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.1298668500990866786e+11)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7282507878605942706e+10)),
diff --git a/boost/math/special_functions/detail/bessel_j1.hpp b/boost/math/special_functions/detail/bessel_j1.hpp
index 3db2503ff6..09d862c240 100644
--- a/boost/math/special_functions/detail/bessel_j1.hpp
+++ b/boost/math/special_functions/detail/bessel_j1.hpp
@@ -22,8 +22,38 @@
 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/boost/math/special_functions/detail/bessel_jy.hpp b/boost/math/special_functions/detail/bessel_jy.hpp
index 19f951ab70..d60dda2d41 100644
--- a/boost/math/special_functions/detail/bessel_jy.hpp
+++ b/boost/math/special_functions/detail/bessel_jy.hpp
@@ -215,8 +215,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
    Cr = br + a / temp;
    Dr = br;
    Di = bi;
-   //std::cout << "C = " << Cr << " " << Ci << std::endl;
-   //std::cout << "D = " << Dr << " " << Di << std::endl;
    if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
    if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
    temp = Dr * Dr + Di * Di;
@@ -227,7 +225,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
    temp = fr;
    fr = temp * delta_r - fi * delta_i;
    fi = temp * delta_i + fi * delta_r;
-   //std::cout << fr << " " << fi << std::endl;
    for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
    {
       a = k - 0.5f;
@@ -239,8 +236,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
       Ci = bi - a * Ci / temp;
       Dr = br + a * Dr;
       Di = bi + a * Di;
-      //std::cout << "C = " << Cr << " " << Ci << std::endl;
-      //std::cout << "D = " << Dr << " " << Di << std::endl;
       if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
       if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
       temp = Dr * Dr + Di * Di;
@@ -253,7 +248,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
       fi = temp * delta_i + fi * delta_r;
       if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
          break;
-      //std::cout << fr << " " << fi << std::endl;
    }
    policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
    *p = fr;
@@ -491,6 +485,16 @@ int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
            CF2_jy(u, x, &p, &q, pol);                  // continued fraction CF2_jy
            T t = u / x - fu;                   // t = J'/J
            gamma = (p - t) / q;
+           //
+           // We can't allow gamma to cancel out to zero competely as it messes up
+           // the subsequent logic.  So pretend that one bit didn't cancel out
+           // and set to a suitably small value.  The only test case we've been able to
+           // find for this, is when v = 8.5 and x = 4*PI.
+           //
+           if(gamma == 0)
+           {
+              gamma = u * tools::epsilon<T>() / x;
+           }
            Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
 
            Jv = Ju * ratio;                    // normalization
diff --git a/boost/math/special_functions/detail/bessel_k0.hpp b/boost/math/special_functions/detail/bessel_k0.hpp
index 81407dab10..e209168e87 100644
--- a/boost/math/special_functions/detail/bessel_k0.hpp
+++ b/boost/math/special_functions/detail/bessel_k0.hpp
@@ -22,10 +22,40 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T, typename Policy>
+T bessel_k0(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_k0_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         bessel_k0(T(1), Policy());
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename bessel_k0_initializer<T, Policy>::init bessel_k0_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
 T bessel_k0(T x, const Policy& pol)
 {
     BOOST_MATH_INSTRUMENT_CODE(x);
 
+    bessel_k0_initializer<T, Policy>::force_instantiate();
+
     static const T P1[] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.4708152720399552679e+03)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9169059852270512312e+03)),
diff --git a/boost/math/special_functions/detail/bessel_k1.hpp b/boost/math/special_functions/detail/bessel_k1.hpp
index 225209f7ba..0d17cd3057 100644
--- a/boost/math/special_functions/detail/bessel_k1.hpp
+++ b/boost/math/special_functions/detail/bessel_k1.hpp
@@ -22,8 +22,38 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T, typename Policy>
+T bessel_k1(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_k1_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         bessel_k1(T(1), Policy());
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename bessel_k1_initializer<T, Policy>::init bessel_k1_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
 T bessel_k1(T x, const Policy& pol)
 {
+    bessel_k1_initializer<T, Policy>::force_instantiate();
+
     static const T P1[] = {
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2149374878243304548e+06)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.1938920065420586101e+05)),
diff --git a/boost/math/special_functions/detail/bessel_y0.hpp b/boost/math/special_functions/detail/bessel_y0.hpp
index e23f861bf0..289bda5f18 100644
--- a/boost/math/special_functions/detail/bessel_y0.hpp
+++ b/boost/math/special_functions/detail/bessel_y0.hpp
@@ -24,8 +24,38 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T, typename Policy>
+T bessel_y0(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_y0_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         bessel_y0(T(1), Policy());
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename bessel_y0_initializer<T, Policy>::init bessel_y0_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
 T bessel_y0(T x, const Policy& pol)
 {
+    bessel_y0_initializer<T, Policy>::force_instantiate();
+
     static const T P1[] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0723538782003176831e+11)),
         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.3716255451260504098e+09)),
diff --git a/boost/math/special_functions/detail/bessel_y1.hpp b/boost/math/special_functions/detail/bessel_y1.hpp
index b85e7011ea..caf09ffd26 100644
--- a/boost/math/special_functions/detail/bessel_y1.hpp
+++ b/boost/math/special_functions/detail/bessel_y1.hpp
@@ -24,8 +24,38 @@
 namespace boost { namespace math { namespace detail{
 
 template <typename T, typename Policy>
+T bessel_y1(T x, const Policy&);
+
+template <class T, class Policy>
+struct bessel_y1_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         bessel_y1(T(1), Policy());
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename bessel_y1_initializer<T, Policy>::init bessel_y1_initializer<T, Policy>::initializer;
+
+template <typename T, typename Policy>
 T bessel_y1(T x, const Policy& pol)
 {
+    bessel_y1_initializer<T, Policy>::force_instantiate();
+
     static const T P1[] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)),
diff --git a/boost/math/special_functions/detail/erf_inv.hpp b/boost/math/special_functions/detail/erf_inv.hpp
index f2f625f991..d51db9d52f 100644
--- a/boost/math/special_functions/detail/erf_inv.hpp
+++ b/boost/math/special_functions/detail/erf_inv.hpp
@@ -322,12 +322,47 @@ T erf_inv_imp(const T& p, const T& q, const Policy& pol, const boost::mpl::int_<
    return result;
 }
 
+template <class T, class Policy>
+struct erf_inv_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init();
+      }
+      static void do_init()
+      {
+         boost::math::erf_inv(static_cast<T>(0.25), Policy());
+         boost::math::erf_inv(static_cast<T>(0.55), Policy());
+         boost::math::erf_inv(static_cast<T>(0.95), Policy());
+         boost::math::erfc_inv(static_cast<T>(1e-15), Policy());
+         if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)) != 0)
+            boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)), Policy());
+         if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)) != 0)
+            boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)), Policy());
+         if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)) != 0)
+            boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)), Policy());
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename erf_inv_initializer<T, Policy>::init erf_inv_initializer<T, Policy>::initializer;
+
 } // namespace detail
 
 template <class T, class Policy>
 typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
 {
    typedef typename tools::promote_args<T>::type result_type;
+
    //
    // Begin by testing for domain errors, and other special cases:
    //
@@ -378,6 +413,8 @@ typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
+   detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
+
    //
    // And get the result, negating where required:
    //
@@ -389,6 +426,7 @@ template <class T, class Policy>
 typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
 {
    typedef typename tools::promote_args<T>::type result_type;
+
    //
    // Begin by testing for domain errors, and other special cases:
    //
@@ -445,6 +483,8 @@ typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
    // precision internally if it's appropriate:
    //
    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
+
+   detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
    //
    // And get the result, negating where required:
    //
diff --git a/boost/math/special_functions/detail/iconv.hpp b/boost/math/special_functions/detail/iconv.hpp
index 8916eaed1d..4256ffcc88 100644
--- a/boost/math/special_functions/detail/iconv.hpp
+++ b/boost/math/special_functions/detail/iconv.hpp
@@ -25,7 +25,7 @@ template <class T, class Policy>
 inline int iconv_imp(T v, Policy const& pol, mpl::false_ const&)
 {
    BOOST_MATH_STD_USING
-   return iround(v);
+   return iround(v, pol);
 }
 
 template <class T, class Policy>
diff --git a/boost/math/special_functions/detail/igamma_large.hpp b/boost/math/special_functions/detail/igamma_large.hpp
index f9a810c489..eb3d4ba93e 100644
--- a/boost/math/special_functions/detail/igamma_large.hpp
+++ b/boost/math/special_functions/detail/igamma_large.hpp
@@ -759,7 +759,6 @@ T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<113> const *)
    return result;
 }
 
-
 }  // namespace detail
 }  // namespace math
 }  // namespace math
diff --git a/boost/math/special_functions/detail/lanczos_sse2.hpp b/boost/math/special_functions/detail/lanczos_sse2.hpp
index 6a3f3e5347..f8846bf376 100644
--- a/boost/math/special_functions/detail/lanczos_sse2.hpp
+++ b/boost/math/special_functions/detail/lanczos_sse2.hpp
@@ -13,7 +13,7 @@
 #include <emmintrin.h>
 
 #if defined(__GNUC__) || defined(__PGI)
-#define ALIGN16 __attribute__((aligned(16)))
+#define ALIGN16 __attribute__((__aligned__(16)))
 #else
 #define ALIGN16 __declspec(align(16))
 #endif
@@ -194,8 +194,11 @@ inline double lanczos13m53::lanczos_sum_expG_scaled<double>(const double& x)
 } // namespace math
 } // namespace boost
 
+#undef ALIGN16
+
 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_LANCZOS
 
 
 
 
+
diff --git a/boost/math/special_functions/detail/lgamma_small.hpp b/boost/math/special_functions/detail/lgamma_small.hpp
index 526a573583..ec28ed2adf 100644
--- a/boost/math/special_functions/detail/lgamma_small.hpp
+++ b/boost/math/special_functions/detail/lgamma_small.hpp
@@ -15,6 +15,14 @@
 namespace boost{ namespace math{ namespace detail{
 
 //
+// These need forward declaring to keep GCC happy:
+//
+template <class T, class Policy, class Lanczos>
+T gamma_imp(T z, const Policy& pol, const Lanczos& l);
+template <class T, class Policy>
+T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l);
+
+//
 // lgamma for small arguments:
 //
 template <class T, class Policy, class Lanczos>
diff --git a/boost/math/special_functions/digamma.hpp b/boost/math/special_functions/digamma.hpp
index 0dbddc77e2..1268b64dc9 100644
--- a/boost/math/special_functions/digamma.hpp
+++ b/boost/math/special_functions/digamma.hpp
@@ -407,6 +407,31 @@ T digamma_imp(T x, const Tag* t, const Policy& pol)
    return result;
 }
 
+//
+// Initializer: ensure all our constants are initialized prior to the first call of main:
+//
+template <class T, class Policy>
+struct digamma_initializer
+{
+   struct init
+   {
+      init()
+      {
+         boost::math::digamma(T(1.5), Policy());
+         boost::math::digamma(T(500), Policy());
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Policy>::initializer;
+
 } // namespace detail
 
 template <class T, class Policy>
@@ -433,6 +458,9 @@ inline typename tools::promote_args<T>::type
       >::type
    >::type tag_type;
 
+   // Force initialization of constants:
+   detail::digamma_initializer<result_type, Policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
       static_cast<value_type>(x),
       static_cast<const tag_type*>(0), pol), "boost::math::digamma<%1%>(%1%)");
diff --git a/boost/math/special_functions/erf.hpp b/boost/math/special_functions/erf.hpp
index 1abb59177f..e67332a61a 100644
--- a/boost/math/special_functions/erf.hpp
+++ b/boost/math/special_functions/erf.hpp
@@ -978,6 +978,59 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
    return result;
 } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<113>& t)
 
+template <class T, class Policy, class tag>
+struct erf_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init(tag());
+      }
+      static void do_init(const mpl::int_<0>&){}
+      static void do_init(const mpl::int_<53>&)
+      {
+         boost::math::erf(static_cast<T>(1e-12), Policy());
+         boost::math::erf(static_cast<T>(0.25), Policy());
+         boost::math::erf(static_cast<T>(1.25), Policy());
+         boost::math::erf(static_cast<T>(2.25), Policy());
+         boost::math::erf(static_cast<T>(4.25), Policy());
+         boost::math::erf(static_cast<T>(5.25), Policy());
+      }
+      static void do_init(const mpl::int_<64>&)
+      {
+         boost::math::erf(static_cast<T>(1e-12), Policy());
+         boost::math::erf(static_cast<T>(0.25), Policy());
+         boost::math::erf(static_cast<T>(1.25), Policy());
+         boost::math::erf(static_cast<T>(2.25), Policy());
+         boost::math::erf(static_cast<T>(4.25), Policy());
+         boost::math::erf(static_cast<T>(5.25), Policy());
+      }
+      static void do_init(const mpl::int_<113>&)
+      {
+         boost::math::erf(static_cast<T>(1e-22), Policy());
+         boost::math::erf(static_cast<T>(0.25), Policy());
+         boost::math::erf(static_cast<T>(1.25), Policy());
+         boost::math::erf(static_cast<T>(2.125), Policy());
+         boost::math::erf(static_cast<T>(2.75), Policy());
+         boost::math::erf(static_cast<T>(3.25), Policy());
+         boost::math::erf(static_cast<T>(5.25), Policy());
+         boost::math::erf(static_cast<T>(7.25), Policy());
+         boost::math::erf(static_cast<T>(11.25), Policy());
+         boost::math::erf(static_cast<T>(12.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 erf_initializer<T, Policy, tag>::init erf_initializer<T, Policy, tag>::initializer;
+
 } // namespace detail
 
 template <class T, class Policy>
@@ -1017,6 +1070,8 @@ inline typename tools::promote_args<T>::type erf(T z, const Policy& /* pol */)
 
    BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());
 
+   detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp(
       static_cast<value_type>(z),
       false,
@@ -1061,6 +1116,8 @@ inline typename tools::promote_args<T>::type erfc(T z, const Policy& /* pol */)
 
    BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());
 
+   detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp(
       static_cast<value_type>(z),
       true,
diff --git a/boost/math/special_functions/expint.hpp b/boost/math/special_functions/expint.hpp
index 06bd78ff48..1c86d282fa 100644
--- a/boost/math/special_functions/expint.hpp
+++ b/boost/math/special_functions/expint.hpp
@@ -1473,6 +1473,94 @@ T expint_i_imp(T z, const Policy& pol, const mpl::int_<113>& tag)
    return result;
 }
 
+template <class T, class Policy, class tag>
+struct expint_i_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init(tag());
+      }
+      static void do_init(const mpl::int_<0>&){}
+      static void do_init(const mpl::int_<53>&)
+      {
+         boost::math::expint(T(5));
+         boost::math::expint(T(7));
+         boost::math::expint(T(18));
+         boost::math::expint(T(38));
+         boost::math::expint(T(45));
+      }
+      static void do_init(const mpl::int_<64>&)
+      {
+         boost::math::expint(T(5));
+         boost::math::expint(T(7));
+         boost::math::expint(T(18));
+         boost::math::expint(T(38));
+         boost::math::expint(T(45));
+      }
+      static void do_init(const mpl::int_<113>&)
+      {
+         boost::math::expint(T(5));
+         boost::math::expint(T(7));
+         boost::math::expint(T(17));
+         boost::math::expint(T(25));
+         boost::math::expint(T(40));
+         boost::math::expint(T(50));
+         boost::math::expint(T(80));
+         boost::math::expint(T(200));
+         boost::math::expint(T(220));
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy, class tag>
+const typename expint_i_initializer<T, Policy, tag>::init expint_i_initializer<T, Policy, tag>::initializer;
+
+template <class T, class Policy, class tag>
+struct expint_1_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init(tag());
+      }
+      static void do_init(const mpl::int_<0>&){}
+      static void do_init(const mpl::int_<53>&)
+      {
+         boost::math::expint(1, T(0.5));
+         boost::math::expint(1, T(2));
+      }
+      static void do_init(const mpl::int_<64>&)
+      {
+         boost::math::expint(1, T(0.5));
+         boost::math::expint(1, T(2));
+      }
+      static void do_init(const mpl::int_<113>&)
+      {
+         boost::math::expint(1, T(0.5));
+         boost::math::expint(1, T(2));
+         boost::math::expint(1, T(6));
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy, class tag>
+const typename expint_1_initializer<T, Policy, tag>::init expint_1_initializer<T, Policy, tag>::initializer;
+
 template <class T, class Policy>
 inline typename tools::promote_args<T>::type
    expint_forwarder(T z, const Policy& /*pol*/, mpl::true_ const&)
@@ -1504,6 +1592,8 @@ inline typename tools::promote_args<T>::type
       >::type
    >::type tag_type;
 
+   expint_i_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_i_imp(
       static_cast<value_type>(z),
       forwarding_policy(),
@@ -1550,6 +1640,8 @@ inline typename tools::promote_args<T>::type
       >::type
    >::type tag_type;
 
+   detail::expint_1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_imp(
       n,
       static_cast<value_type>(z),
diff --git a/boost/math/special_functions/expm1.hpp b/boost/math/special_functions/expm1.hpp
index 345220fcee..9ff2541fb1 100644
--- a/boost/math/special_functions/expm1.hpp
+++ b/boost/math/special_functions/expm1.hpp
@@ -65,34 +65,37 @@ namespace detail
      expm1_series& operator=(const expm1_series&);
   };
 
-template <class T, bool b = boost::is_pod<T>::value>
-struct expm1_init_on_startup
+template <class T, class Policy, class tag>
+struct expm1_initializer
 {
    struct init
    {
       init()
       {
-         boost::math::expm1(T(0.5f));
+         do_init(tag());
       }
-      void do_nothing()const{}
+      template <int N>
+      static void do_init(const mpl::int_<N>&){}
+      static void do_init(const mpl::int_<64>&)
+      {
+         expm1(T(0.5));
+      }
+      static void do_init(const mpl::int_<113>&)
+      {
+         expm1(T(0.5));
+      }
+      void force_instantiate()const{}
    };
-
-   static void do_nothing()
+   static const init initializer;
+   static void force_instantiate()
    {
-      initializer.do_nothing();
+      initializer.force_instantiate();
    }
-
-   static const init initializer;
 };
 
-template <class T, bool b>
-const typename expm1_init_on_startup<T, b>::init expm1_init_on_startup<T, b>::initializer;
+template <class T, class Policy, class tag>
+const typename expm1_initializer<T, Policy, tag>::init expm1_initializer<T, Policy, tag>::initializer;
 
-template <class T>
-struct expm1_init_on_startup<T, true>
-{
-   static void do_nothing(){}
-};
 //
 // Algorithm expm1 is part of C99, but is not yet provided by many compilers.
 //
@@ -133,8 +136,6 @@ T expm1_imp(T x, const mpl::int_<53>&, const P& pol)
 {
    BOOST_MATH_STD_USING
 
-   expm1_init_on_startup<T>::do_nothing();
-
    T a = fabs(x);
    if(a > T(0.5L))
    {
@@ -162,8 +163,6 @@ T expm1_imp(T x, const mpl::int_<64>&, const P& pol)
 {
    BOOST_MATH_STD_USING
 
-   expm1_init_on_startup<T>::do_nothing();
-
    T a = fabs(x);
    if(a > T(0.5L))
    {
@@ -207,8 +206,6 @@ T expm1_imp(T x, const mpl::int_<113>&, const P& pol)
 {
    BOOST_MATH_STD_USING
 
-   expm1_init_on_startup<T>::do_nothing();
-
    T a = fabs(x);
    if(a > T(0.5L))
    {
@@ -287,6 +284,8 @@ inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */)
       >::type
    >::type tag_type;
 
+   detail::expm1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
+   
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp(
       static_cast<value_type>(x),
       tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)");
diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp
index 1ae965f18c..86d15b7f2a 100644
--- a/boost/math/special_functions/gamma.hpp
+++ b/boost/math/special_functions/gamma.hpp
@@ -1258,6 +1258,101 @@ inline typename tools::promote_args<T>::type
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)");
 }
 
+template <class T, class Policy>
+struct igamma_initializer
+{
+   struct init
+   {
+      init()
+      {
+         typedef typename policies::precision<T, Policy>::type precision_type;
+
+         typedef typename mpl::if_<
+            mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
+            mpl::greater<precision_type, mpl::int_<113> > >,
+            mpl::int_<0>,
+            typename mpl::if_<
+               mpl::less_equal<precision_type, mpl::int_<53> >,
+               mpl::int_<53>,
+               typename mpl::if_<
+                  mpl::less_equal<precision_type, mpl::int_<64> >,
+                  mpl::int_<64>,
+                  mpl::int_<113>
+               >::type
+            >::type
+         >::type tag_type;
+
+         do_init(tag_type());
+      }
+      template <int N>
+      static void do_init(const mpl::int_<N>&)
+      {
+         boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
+      }
+      static void do_init(const mpl::int_<53>&){}
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
+
+template <class T, class Policy>
+struct lgamma_initializer
+{
+   struct init
+   {
+      init()
+      {
+         typedef typename policies::precision<T, Policy>::type precision_type;
+         typedef typename mpl::if_<
+            mpl::and_<
+               mpl::less_equal<precision_type, mpl::int_<64> >, 
+               mpl::greater<precision_type, mpl::int_<0> > 
+            >,
+            mpl::int_<64>,
+            typename mpl::if_<
+               mpl::and_<
+                  mpl::less_equal<precision_type, mpl::int_<113> >,
+                  mpl::greater<precision_type, mpl::int_<0> > 
+               >,
+               mpl::int_<113>, mpl::int_<0> >::type
+             >::type tag_type;
+         do_init(tag_type());
+      }
+      static void do_init(const mpl::int_<64>&)
+      {
+         boost::math::lgamma(static_cast<T>(2.5), Policy());
+         boost::math::lgamma(static_cast<T>(1.25), Policy());
+         boost::math::lgamma(static_cast<T>(1.75), Policy());
+      }
+      static void do_init(const mpl::int_<113>&)
+      {
+         boost::math::lgamma(static_cast<T>(2.5), Policy());
+         boost::math::lgamma(static_cast<T>(1.25), Policy());
+         boost::math::lgamma(static_cast<T>(1.5), Policy());
+         boost::math::lgamma(static_cast<T>(1.75), Policy());
+      }
+      static void do_init(const mpl::int_<0>&)
+      {
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
+
 template <class T1, class T2, class Policy>
 inline typename tools::promote_args<T1, T2>::type
    tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
@@ -1272,6 +1367,9 @@ inline typename tools::promote_args<T1, T2>::type
       policies::promote_double<false>, 
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
+
+   igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), false, true,
@@ -1285,6 +1383,7 @@ inline typename tools::promote_args<T1, T2>::type
    return tgamma(a, z, policies::policy<>(), tag);
 }
 
+
 } // namespace detail
 
 template <class T>
@@ -1308,6 +1407,9 @@ inline typename tools::promote_args<T>::type
       policies::promote_double<false>, 
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
+
+   detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)");
 }
 
@@ -1395,6 +1497,8 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
+   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), false, false,
@@ -1424,6 +1528,8 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
+   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), true, true,
@@ -1453,6 +1559,8 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
+   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), true, false,
diff --git a/boost/math/special_functions/hankel.hpp b/boost/math/special_functions/hankel.hpp
new file mode 100644
index 0000000000..bc3fc2d742
--- /dev/null
+++ b/boost/math/special_functions/hankel.hpp
@@ -0,0 +1,178 @@
+// Copyright John Maddock 2012.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_HANKEL_HPP
+#define BOOST_MATH_HANKEL_HPP
+
+#include <boost/math/special_functions/bessel.hpp>
+
+namespace boost{ namespace math{
+
+namespace detail{
+
+template <class T, class Policy>
+std::complex<T> hankel_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol, int sign)
+{
+   BOOST_MATH_STD_USING
+   static const char* function = "boost::math::cyl_hankel_1<%1%>(%1%,%1%)";
+
+   if(x < 0)
+   {
+      bool isint_v = floor(v) == v;
+      T j, y;
+      bessel_jy(v, -x, &j, &y, need_j | need_y, pol);
+      std::complex<T> cx(x), cv(v);
+      std::complex<T> j_result, y_result;
+      if(isint_v)
+      {
+         int s = (iround(j) & 1) ? -1 : 1;
+         j_result = j * s;
+         y_result = T(s) * (y - (2 / constants::pi<T>()) * (log(-x) - log(cx)) * j);
+      }
+      else
+      {
+         j_result = pow(cx, v) * pow(-cx, -v) * j;
+         T p1 = pow(-x, v);
+         std::complex<T> p2 = pow(cx, v);
+         y_result = p1 * y / p2
+            + (p2 / p1 - p1 / p2) * j / tan(constants::pi<T>() * v);
+      }
+      // multiply y_result by i:
+      y_result = std::complex<T>(-sign * y_result.imag(), sign * y_result.real());
+      return j_result + y_result;
+   }
+
+   if(x == 0)
+   {
+      if(v == 0)
+      {
+         // J is 1, Y is -INF
+         return std::complex<T>(1, sign * -policies::raise_overflow_error<T>(function, 0, pol));
+      }
+      else
+      {
+         // At least one of J and Y is complex infinity:
+         return std::complex<T>(policies::raise_overflow_error<T>(function, 0, pol), sign * policies::raise_overflow_error<T>(function, 0, pol));
+      }
+   }
+
+   T j, y;
+   bessel_jy(v, x, &j, &y, need_j | need_y, pol);
+   return std::complex<T>(j, sign * y);
+}
+
+template <class T, class Policy>
+std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign);
+
+template <class T, class Policy>
+inline std::complex<T> hankel_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol, int sign)
+{
+   BOOST_MATH_STD_USING  // ADL of std names.
+   int ival = detail::iconv(v, pol);
+   if(0 == v - ival)
+   {
+      return hankel_imp(ival, x, bessel_int_tag(), pol, sign);
+   }
+   return hankel_imp(v, x, bessel_no_int_tag(), pol, sign);
+}
+
+template <class T, class Policy>
+inline std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign)
+{
+   BOOST_MATH_STD_USING
+   if((std::abs(v < 200)) && (x > 0))
+      return std::complex<T>(bessel_jn(v, x, pol), sign * bessel_yn(v, x, pol));
+   return hankel_imp(static_cast<T>(v), x, bessel_no_int_tag(), pol, sign);
+}
+
+template <class T, class Policy>
+inline std::complex<T> sph_hankel_imp(T v, T x, const Policy& pol, int sign)
+{
+   BOOST_MATH_STD_USING
+   return constants::root_half_pi<T>() * hankel_imp(v + 0.5f, x, bessel_no_int_tag(), pol, sign) / sqrt(std::complex<T>(x));
+}
+
+} // namespace detail
+
+template <class T1, class T2, class Policy>
+inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
+   typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, 1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)");
+}
+
+template <class T1, class T2>
+inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x)
+{
+   return cyl_hankel_1(v, x, policies::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_2(T1 v, T2 x, const Policy& pol)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
+   typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, -1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)");
+}
+
+template <class T1, class T2>
+inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_2(T1 v, T2 x)
+{
+   return cyl_hankel_2(v, x, policies::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_1(T1 v, T2 x, const Policy& pol)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy, 
+      policies::promote_float<false>, 
+      policies::promote_double<false>, 
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), 1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)");
+}
+
+template <class T1, class T2>
+inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_1(T1 v, T2 x)
+{
+   return sph_hankel_1(v, x, policies::policy<>());
+}
+
+template <class T1, class T2, class Policy>
+inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_2(T1 v, T2 x, const Policy& pol)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy, 
+      policies::promote_float<false>, 
+      policies::promote_double<false>, 
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), -1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)");
+}
+
+template <class T1, class T2>
+inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_2(T1 v, T2 x)
+{
+   return sph_hankel_2(v, x, policies::policy<>());
+}
+
+}} // namespaces
+
+#endif // BOOST_MATH_HANKEL_HPP
\ No newline at end of file
diff --git a/boost/math/special_functions/lanczos.hpp b/boost/math/special_functions/lanczos.hpp
index 20ff969359..ed891549f1 100644
--- a/boost/math/special_functions/lanczos.hpp
+++ b/boost/math/special_functions/lanczos.hpp
@@ -34,6 +34,35 @@ namespace boost{ namespace math{ namespace lanczos{
 // http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at
 // http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision.
 //
+// Begin with a small helper to force initialization of constants prior
+// to main.  This makes the constant initialization thread safe, even
+// when called with a user-defined number type.
+//
+template <class Lanczos, class T>
+struct lanczos_initializer
+{
+   struct init
+   {
+      init()
+      {
+         T t(1);
+         Lanczos::lanczos_sum(t);
+         Lanczos::lanczos_sum_expG_scaled(t);
+         Lanczos::lanczos_sum_near_1(t);
+         Lanczos::lanczos_sum_near_2(t);
+         Lanczos::g();
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+template <class Lanczos, class T>
+typename lanczos_initializer<Lanczos, T>::init const lanczos_initializer<Lanczos, T>::initializer;
+//
 // Lanczos Coefficients for N=6 G=5.581
 // Max experimental error (with arbitary precision arithmetic) 9.516e-12
 // Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
@@ -47,6 +76,7 @@ struct lanczos6 : public mpl::int_<35>
    template <class T>
    static T lanczos_sum(const T& z)
    {
+      lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[6] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 8706.349592549009182288174442774377925882)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 8523.650341121874633477483696775067709735)),
@@ -69,6 +99,7 @@ struct lanczos6 : public mpl::int_<35>
    template <class T>
    static T lanczos_sum_expG_scaled(const T& z)
    {
+      lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[6] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 32.81244541029783471623665933780748627823)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 32.12388941444332003446077108933558534361)),
@@ -92,6 +123,7 @@ struct lanczos6 : public mpl::int_<35>
    template<class T>
    static T lanczos_sum_near_1(const T& dz)
    {
+      lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[5] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 2.044879010930422922760429926121241330235)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, -2.751366405578505366591317846728753993668)),
@@ -110,6 +142,7 @@ struct lanczos6 : public mpl::int_<35>
    template<class T>
    static T lanczos_sum_near_2(const T& dz)
    {
+      lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[5] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 5.748142489536043490764289256167080091892)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, -7.734074268282457156081021756682138251825)),
@@ -143,6 +176,7 @@ struct lanczos11 : public mpl::int_<60>
    template <class T>
    static T lanczos_sum(const T& z)
    {
+      lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[11] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 38474670393.31776828316099004518914832218)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 36857665043.51950660081971227404959150474)),
@@ -175,6 +209,7 @@ struct lanczos11 : public mpl::int_<60>
    template <class T>
    static T lanczos_sum_expG_scaled(const T& z)
    {
+      lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[11] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 709811.662581657956893540610814842699825)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 679979.847415722640161734319823103390728)),
@@ -208,6 +243,7 @@ struct lanczos11 : public mpl::int_<60>
    template<class T>
    static T lanczos_sum_near_1(const T& dz)
    {
+      lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[10] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 4.005853070677940377969080796551266387954)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, -13.17044315127646469834125159673527183164)),
@@ -231,6 +267,7 @@ struct lanczos11 : public mpl::int_<60>
    template<class T>
    static T lanczos_sum_near_2(const T& dz)
    {
+      lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[10] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 19.05889633808148715159575716844556056056)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, -62.66183664701721716960978577959655644762)),
@@ -269,6 +306,7 @@ struct lanczos13 : public mpl::int_<72>
    template <class T>
    static T lanczos_sum(const T& z)
    {
+      lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[13] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 44012138428004.60895436261759919070125699)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 41590453358593.20051581730723108131357995)),
@@ -305,6 +343,7 @@ struct lanczos13 : public mpl::int_<72>
    template <class T>
    static T lanczos_sum_expG_scaled(const T& z)
    {
+      lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[13] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 86091529.53418537217994842267760536134841)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 81354505.17858011242874285785316135398567)),
@@ -342,6 +381,7 @@ struct lanczos13 : public mpl::int_<72>
    template<class T>
    static T lanczos_sum_near_1(const T& dz)
    {
+      lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[12] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 4.832115561461656947793029596285626840312)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, -19.86441536140337740383120735104359034688)),
@@ -367,6 +407,7 @@ struct lanczos13 : public mpl::int_<72>
    template<class T>
    static T lanczos_sum_near_2(const T& dz)
    {
+      lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[12] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 26.96979819614830698367887026728396466395)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, -110.8705424709385114023884328797900204863)),
@@ -382,7 +423,7 @@ struct lanczos13 : public mpl::int_<72>
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, -0.9685385411006641478305219367315965391289e-9)),
       };
       T result = 0;
-      T z = z = 2;
+      T z = dz + 2;
       for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
       {
          result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
@@ -407,6 +448,7 @@ struct lanczos22 : public mpl::int_<120>
    template <class T>
    static T lanczos_sum(const T& z)
    {
+      lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[22] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 46198410803245094237463011094.12173081986)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 43735859291852324413622037436.321513777)),
@@ -461,6 +503,7 @@ struct lanczos22 : public mpl::int_<120>
    template <class T>
    static T lanczos_sum_expG_scaled(const T& z)
    {
+      lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[22] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 6939996264376682180.277485395074954356211)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 6570067992110214451.87201438870245659384)),
@@ -516,6 +559,7 @@ struct lanczos22 : public mpl::int_<120>
    template<class T>
    static T lanczos_sum_near_1(const T& dz)
    {
+      lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[21] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 8.318998691953337183034781139546384476554)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, -63.15415991415959158214140353299240638675)),
@@ -550,6 +594,7 @@ struct lanczos22 : public mpl::int_<120>
    template<class T>
    static T lanczos_sum_near_2(const T& dz)
    {
+      lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[21] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 75.39272007105208086018421070699575462226)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, -572.3481967049935412452681346759966390319)),
@@ -830,6 +875,7 @@ struct lanczos17m64 : public mpl::int_<64>
    template <class T>
    static T lanczos_sum(const T& z)
    {
+      lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[17] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 553681095419291969.2230556393350368550504)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 731918863887667017.2511276782146694632234)),
@@ -874,6 +920,7 @@ struct lanczos17m64 : public mpl::int_<64>
    template <class T>
    static T lanczos_sum_expG_scaled(const T& z)
    {
+      lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[17] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2715894658327.717377557655133124376674911)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3590179526097.912105038525528721129550434)),
@@ -919,6 +966,7 @@ struct lanczos17m64 : public mpl::int_<64>
    template<class T>
    static T lanczos_sum_near_1(const T& dz)
    {
+      lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[16] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.493645054286536365763334986866616581265)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -16.95716370392468543800733966378143997694)),
@@ -948,6 +996,7 @@ struct lanczos17m64 : public mpl::int_<64>
    template<class T>
    static T lanczos_sum_near_2(const T& dz)
    {
+      lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[16] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 23.56409085052261327114594781581930373708)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -88.92116338946308797946237246006238652361)),
@@ -991,6 +1040,7 @@ struct lanczos24m113 : public mpl::int_<113>
    template <class T>
    static T lanczos_sum(const T& z)
    {
+      lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[24] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 2029889364934367661624137213253.22102954656825019111612712252027267955023987678816620961507)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 2338599599286656537526273232565.2727349714338768161421882478417543004440597874814359063158)),
@@ -1049,6 +1099,7 @@ struct lanczos24m113 : public mpl::int_<113>
    template <class T>
    static T lanczos_sum_expG_scaled(const T& z)
    {
+      lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T num[24] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 3035162425359883494754.02878223286972654682199012688209026810841953293372712802258398358538)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 3496756894406430103600.16057175075063458536101374170860226963245118484234495645518505519827)),
@@ -1108,6 +1159,7 @@ struct lanczos24m113 : public mpl::int_<113>
    template<class T>
    static T lanczos_sum_near_1(const T& dz)
    {
+      lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[23] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 7.4734083002469026177867421609938203388868806387315406134072298925733950040583068760685908)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, -50.4225805042247530267317342133388132970816607563062253708655085754357843064134941138154171)),
@@ -1144,6 +1196,7 @@ struct lanczos24m113 : public mpl::int_<113>
    template<class T>
    static T lanczos_sum_near_2(const T& dz)
    {
+      lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main()
       static const T d[23] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 61.4165001061101455341808888883960361969557848005400286332291451422461117307237198559485365)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, -414.372973678657049667308134761613915623353625332248315105320470271523320700386200587519147)),
diff --git a/boost/math/special_functions/log1p.hpp b/boost/math/special_functions/log1p.hpp
index 9bae7165e4..989bdc21b6 100644
--- a/boost/math/special_functions/log1p.hpp
+++ b/boost/math/special_functions/log1p.hpp
@@ -258,6 +258,34 @@ T log1p_imp(T const& x, const Policy& pol, const mpl::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 mpl::int_<N>&){}
+      static void do_init(const mpl::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>
@@ -286,6 +314,9 @@ inline typename tools::promote_args<T>::type log1p(T x, const Policy&)
          >::type
       >::type
    >::type 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/boost/math/special_functions/math_fwd.hpp b/boost/math/special_functions/math_fwd.hpp
index 14364a3d5c..982cdf7ca3 100644
--- a/boost/math/special_functions/math_fwd.hpp
+++ b/boost/math/special_functions/math_fwd.hpp
@@ -28,6 +28,7 @@
 #include <boost/math/policies/policy.hpp>
 #include <boost/mpl/comparison.hpp>
 #include <boost/config/no_tr1/complex.hpp>
+#include <complex>
 
 #define BOOST_NO_MACRO_EXPAND /**/
 
@@ -614,6 +615,30 @@ namespace boost
    template <class T>
    typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann(unsigned v, T x);
 
+   template <class T1, class T2, class Policy>
+   std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol);
+
+   template <class T1, class T2>
+   std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x);
+
+   template <class T1, class T2, class Policy>
+   std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_2(T1 v, T2 x, const Policy& pol);
+
+   template <class T1, class T2>
+   std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_2(T1 v, T2 x);
+
+   template <class T1, class T2, class Policy>
+   std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_1(T1 v, T2 x, const Policy& pol);
+
+   template <class T1, class T2>
+   std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_1(T1 v, T2 x);
+
+   template <class T1, class T2, class Policy>
+   std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_2(T1 v, T2 x, const Policy& pol);
+
+   template <class T1, class T2>
+   std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_2(T1 v, T2 x);
+
    template <class T, class Policy>
    typename tools::promote_args<T>::type sin_pi(T x, const Policy&);
 
@@ -681,6 +706,13 @@ namespace boost
    template <class T, class Policy>
    typename tools::promote_args<T>::type zeta(T s, const Policy&);
 
+   // Owen's T function:
+   template <class T1, class T2, class Policy>
+   typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol);
+
+   template <class T1, class T2>
+   typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a);
+
    template <class T>
    typename tools::promote_args<T>::type zeta(T s);
 
@@ -1063,6 +1095,26 @@ namespace boost
    template <class T> T float_next(const T& a){ return boost::math::float_next(a, Policy()); }\
    template <class T> T float_prior(const T& a){ return boost::math::float_prior(a, Policy()); }\
    template <class T> T float_distance(const T& a, const T& b){ return boost::math::float_distance(a, b, Policy()); }\
+   \
+   template <class RT1, class RT2>\
+   inline typename boost::math::tools::promote_args<RT1, RT2>::type owens_t(RT1 a, RT2 z){ return boost::math::owens_t(a, z, Policy()); }\
+   \
+   template <class T1, class T2>\
+   inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> cyl_hankel_1(T1 v, T2 x)\
+   {  return boost::math::cyl_hankel_1(v, x, Policy()); }\
+   \
+   template <class T1, class T2>\
+   inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> cyl_hankel_2(T1 v, T2 x)\
+   { return boost::math::cyl_hankel_2(v, x, Policy()); }\
+   \
+   template <class T1, class T2>\
+   inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> sph_hankel_1(T1 v, T2 x)\
+   { return boost::math::sph_hankel_1(v, x, Policy()); }\
+   \
+   template <class T1, class T2>\
+   inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> sph_hankel_2(T1 v, T2 x)\
+   { return boost::math::sph_hankel_2(v, x, Policy()); }\
+
 
 
 #endif // BOOST_MATH_SPECIAL_MATH_FWD_HPP
diff --git a/boost/math/special_functions/nonfinite_num_facets.hpp b/boost/math/special_functions/nonfinite_num_facets.hpp
index 9fa61481b5..84d3f1070a 100644
--- a/boost/math/special_functions/nonfinite_num_facets.hpp
+++ b/boost/math/special_functions/nonfinite_num_facets.hpp
@@ -1,8 +1,9 @@
 #ifndef BOOST_MATH_NONFINITE_NUM_FACETS_HPP
 #define BOOST_MATH_NONFINITE_NUM_FACETS_HPP
 
-// Copyright (c) 2006 Johan Rade
-// Copyright 2011 Paul A. Bristow (comments)
+// Copyright 2006 Johan Rade
+// Copyright 2012 K R Walker
+// Copyright 2011, 2012 Paul A. Bristow 
 
 // Distributed under the Boost Software License, Version 1.0.
 // (See accompanying file LICENSE_1_0.txt
@@ -28,10 +29,9 @@
 #include <boost/math/special_functions/sign.hpp>
 
 #ifdef _MSC_VER
-#   pragma warning(push)
-#   pragma warning(disable : 4127) // conditional expression is constant.
-#   pragma warning(disable : 4706) // assignment within conditional expression.
-#   pragma warning(disable : 4224) // formal parameter 'version' was previously defined as a type.
+#  pragma warning(push)
+#  pragma warning(disable : 4127) // conditional expression is constant.
+#  pragma warning(disable : 4706) // assignment within conditional expression.
 #endif
 
 namespace boost {
@@ -65,16 +65,14 @@ namespace boost {
 
     protected:
       virtual OutputIterator do_put(
-        OutputIterator it, std::ios_base& iosb,
-        CharType fill, double val) const
+        OutputIterator it, std::ios_base& iosb, CharType fill, double val) const
       {
         put_and_reset_width(it, iosb, fill, val);
         return it;
       }
 
       virtual OutputIterator do_put(
-        OutputIterator it, std::ios_base& iosb,
-        CharType fill, long double val) const
+        OutputIterator it, std::ios_base& iosb,  CharType fill, long double val) const
       {
         put_and_reset_width(it, iosb, fill, val);
         return it;
@@ -93,86 +91,136 @@ namespace boost {
         OutputIterator& it, std::ios_base& iosb,
         CharType fill, ValType val) const
       {
-        switch((boost::math::fpclassify)(val)) {
+        static const CharType prefix_plus[2] = { '+', '\0' };
+        static const CharType prefix_minus[2] = { '-', '\0' };
+        static const CharType body_inf[4] = { 'i', 'n', 'f', '\0' };
+        static const CharType body_nan[4] = { 'n', 'a', 'n', '\0' };
+        static const CharType* null_string = 0;
+
+        switch((boost::math::fpclassify)(val))
+        {
 
         case FP_INFINITE:
           if(flags_ & trap_infinity)
+          {
             throw std::ios_base::failure("Infinity");
+          }
           else if((boost::math::signbit)(val))
-            put_num_and_fill(it, iosb, "-", "inf", fill);
+          { // negative infinity.
+            put_num_and_fill(it, iosb, prefix_minus, body_inf, fill, val);
+          }
           else if(iosb.flags() & std::ios_base::showpos)
-            put_num_and_fill(it, iosb, "+", "inf", fill);
+          { // Explicit "+inf" wanted.
+            put_num_and_fill(it, iosb, prefix_plus, body_inf, fill, val);
+          }
           else
-            put_num_and_fill(it, iosb, "", "inf", fill);
+          { // just "inf" wanted.
+            put_num_and_fill(it, iosb, null_string, body_inf, fill, val);
+          }
           break;
 
         case FP_NAN:
           if(flags_ & trap_nan)
+          {
             throw std::ios_base::failure("NaN");
+          }
           else if((boost::math::signbit)(val))
-            put_num_and_fill(it, iosb, "-", "nan", fill);
+          { // negative so "-nan".
+            put_num_and_fill(it, iosb, prefix_minus, body_nan, fill, val);
+          }
           else if(iosb.flags() & std::ios_base::showpos)
-            put_num_and_fill(it, iosb, "+", "nan", fill);
+          { // explicit "+nan" wanted.
+            put_num_and_fill(it, iosb, prefix_plus, body_nan, fill, val);
+          }
           else
-            put_num_and_fill(it, iosb, "", "nan", fill);
+          { // Just "nan".
+            put_num_and_fill(it, iosb, null_string, body_nan, fill, val);
+          }
           break;
 
         case FP_ZERO:
-          if(flags_ & signed_zero) {
-            if((boost::math::signbit)(val))
-              put_num_and_fill(it, iosb, "-", "0", fill);
-            else if(iosb.flags() & std::ios_base::showpos)
-              put_num_and_fill(it, iosb, "+", "0", fill);
-            else
-              put_num_and_fill(it, iosb, "", "0", fill);
+          if((flags_ & signed_zero) && ((boost::math::signbit)(val)))
+          { // Flag set to distinguish between positive and negative zero.
+            // But string "0" should have stuff after decimal point if setprecision and/or exp format. 
+
+            std::basic_ostringstream<CharType> zeros; // Needs to be CharType version.
+
+            // Copy flags, fill, width and precision.
+            zeros.flags(iosb.flags());
+            zeros.unsetf(std::ios::showpos); // Ignore showpos because must be negative.
+            zeros.precision(iosb.precision());
+            //zeros.width is set by put_num_and_fill
+            zeros.fill(static_cast<char>(fill));
+            zeros << ValType(0);
+            put_num_and_fill(it, iosb, prefix_minus, zeros.str().c_str(), fill, val);
           }
           else
-            put_num_and_fill(it, iosb, "", "0", fill);
+          { // Output the platform default for positive and negative zero.
+            put_num_and_fill(it, iosb, null_string, null_string, fill, val);
+          }
           break;
 
-        default:
-          it = std::num_put<CharType, OutputIterator>::do_put(
-            it, iosb, fill, val);
+        default:  // Normal non-zero finite value.
+          it = std::num_put<CharType, OutputIterator>::do_put(it, iosb, fill, val);
           break;
         }
       }
 
+      template<class ValType>
       void put_num_and_fill(
-        OutputIterator& it, std::ios_base& iosb, const char* prefix,
-        const char* body, CharType fill) const
+        OutputIterator& it, std::ios_base& iosb, const CharType* prefix,
+          const CharType* body, CharType fill, ValType val) const
       {
-        int width = (int)std::strlen(prefix) + (int)std::strlen(body);
-        std::ios_base::fmtflags adjust
-          = iosb.flags() & std::ios_base::adjustfield;
+        int prefix_length = prefix ? (int)std::char_traits<CharType>::length(prefix) : 0;
+        int body_length = body ? (int)std::char_traits<CharType>::length(body) : 0;
+        int width = prefix_length + body_length;
+        std::ios_base::fmtflags adjust = iosb.flags() & std::ios_base::adjustfield;
         const std::ctype<CharType>& ct
           = std::use_facet<std::ctype<CharType> >(iosb.getloc());
 
-        if(adjust != std::ios_base::internal && adjust != std::ios_base::left)
-          put_fill(it, iosb, fill, width);
+        if(body || prefix)
+        { // adjust == std::ios_base::right, so leading fill needed.
+          if(adjust != std::ios_base::internal && adjust != std::ios_base::left)
+            put_fill(it, iosb, fill, width);
+        }
 
-        while(*prefix)
-          *it = ct.widen(*(prefix++));
+        if(prefix)
+        { // Adjust width for prefix.
+          while(*prefix)
+            *it = *(prefix++);
+          iosb.width( iosb.width() - prefix_length );
+          width -= prefix_length;
+        }
 
-        if(adjust == std::ios_base::internal)
-          put_fill(it, iosb, fill, width);
+        if(body)
+        { // 
+          if(adjust == std::ios_base::internal)
+          { // Put fill between sign and digits.
+            put_fill(it, iosb, fill, width);
+          }
+          if(iosb.flags() & std::ios_base::uppercase)
+          {
+              while(*body)
+                *it = ct.toupper(*(body++));
+          }
+          else
+          {
+            while(*body)
+              *it = *(body++);
+          }
 
-        if(iosb.flags() & std::ios_base::uppercase) {
-          while(*body)
-            *it = ct.toupper(ct.widen(*(body++)));
+          if(adjust == std::ios_base::left)
+            put_fill(it, iosb, fill, width);
         }
-        else {
-          while(*body)
-            *it = ct.widen(*(body++));
+        else
+        {
+          it = std::num_put<CharType, OutputIterator>::do_put(it, iosb, fill, val);
         }
-
-        if(adjust == std::ios_base::left)
-          put_fill(it, iosb, fill, width);
       }
 
       void put_fill(
-        OutputIterator& it, std::ios_base& iosb,
-        CharType fill, int width) const
-      {
+        OutputIterator& it, std::ios_base& iosb, CharType fill, int width) const
+      { // Insert fill chars.
         for(std::streamsize i = iosb.width() - static_cast<std::streamsize>(width); i > 0; --i)
           *it = fill;
       }
@@ -540,4 +588,5 @@ namespace boost {
 #   pragma warning(pop)
 #endif
 
-#endif
+#endif // BOOST_MATH_NONFINITE_NUM_FACETS_HPP
+
diff --git a/boost/math/special_functions/owens_t.hpp b/boost/math/special_functions/owens_t.hpp
new file mode 100644
index 0000000000..98d6380c39
--- /dev/null
+++ b/boost/math/special_functions/owens_t.hpp
@@ -0,0 +1,1061 @@
+// (C) Benjamin Sobotta 2012
+
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_OWENS_T_HPP
+#define BOOST_OWENS_T_HPP
+
+// Reference:
+// Mike Patefield, David Tandy
+// FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION
+// Journal of Statistical Software, 5 (5), 1-25
+
+#ifdef _MSC_VER
+#  pragma once
+#endif
+
+#include <boost/config/no_tr1/cmath.hpp>
+#include <boost/math/special_functions/erf.hpp>
+#include <boost/math/special_functions/expm1.hpp>
+#include <boost/throw_exception.hpp>
+#include <boost/assert.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/tools/big_constant.hpp>
+
+#include <stdexcept>
+
+namespace boost
+{
+   namespace math
+   {
+      namespace detail
+      {
+         // owens_t_znorm1(x) = P(-oo<Z<=x)-0.5 with Z being normally distributed.
+         template<typename RealType>
+         inline RealType owens_t_znorm1(const RealType x)
+         {
+            using namespace boost::math::constants;
+            return erf(x*one_div_root_two<RealType>())*half<RealType>();
+         } // RealType owens_t_znorm1(const RealType x)
+
+         // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed.
+         template<typename RealType>
+         inline RealType owens_t_znorm2(const RealType x)
+         {
+            using namespace boost::math::constants;
+            return erfc(x*one_div_root_two<RealType>())*half<RealType>();
+         } // RealType owens_t_znorm2(const RealType x)
+
+         // Auxiliary function, it computes an array key that is used to determine
+         // the specific computation method for Owen's T and the order thereof
+         // used in owens_t_dispatch.
+         template<typename RealType>
+         inline unsigned short owens_t_compute_code(const RealType h, const RealType a)
+         {
+            static const RealType hrange[] =
+            {0.02, 0.06, 0.09, 0.125, 0.26, 0.4,  0.6,  1.6,  1.7,  2.33,  2.4,  3.36, 3.4,  4.8};
+
+            static const RealType arange[] = {0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999};
+            /*
+            original select array from paper:
+            1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9
+            1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9
+            2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10
+            2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10
+            2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11
+            2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12
+            2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12
+            2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12
+            */                  
+            // subtract one because the array is written in FORTRAN in mind - in C arrays start @ zero
+            static const unsigned short select[] =
+            {
+               0,    0 ,   1  , 12   ,12 ,  12  , 12  , 12 ,  12  , 12  , 12  , 15  , 15 ,  15  ,  8,
+               0  ,  1  ,  1   , 2 ,   2   , 4  ,  4  , 13 ,  13  , 14  , 14 ,  15  , 15  , 15  ,  8,
+               1  ,  1   , 2 ,   2  ,  2  ,  4   , 4  , 14  , 14 ,  14  , 14 ,  15  , 15 ,  15  ,  9,
+               1  ,  1   , 2 ,   4  ,  4  ,  4   , 4  ,  6  ,  6 ,  15  , 15 ,  15 ,  15 ,  15  ,  9,
+               1  ,  2   , 2  ,  4  ,  4  ,  5   , 5  ,  7  ,  7  , 16   ,16 ,  16 ,  11 ,  11 ,  10,
+               1  ,  2   , 4  ,  4   , 4  ,  5   , 5  ,  7  ,  7  , 16  , 16 ,  16 ,  11  , 11 ,  11,
+               1  ,  2   , 3  ,  3  ,  5  ,  5   , 7  ,  7  , 16 ,  16  , 16 ,  16 ,  16  , 11 ,  11,
+               1  ,  2   , 3   , 3   , 5  ,  5 ,  17  , 17  , 17 ,  17  , 16 ,  16 ,  16 ,  11 ,  11
+            };
+
+            unsigned short ihint = 14, iaint = 7;
+            for(unsigned short i = 0; i != 14; i++)
+            {
+               if( h <= hrange[i] )
+               {
+                  ihint = i;
+                  break;
+               }
+            } // for(unsigned short i = 0; i != 14; i++)
+
+            for(unsigned short i = 0; i != 7; i++)
+            {
+               if( a <= arange[i] )
+               {
+                  iaint = i;
+                  break;
+               }
+            } // for(unsigned short i = 0; i != 7; i++)
+
+            // interprete select array as 8x15 matrix
+            return select[iaint*15 + ihint];
+
+         } // unsigned short owens_t_compute_code(const RealType h, const RealType a)
+
+         template<typename RealType>
+         inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<53>&)
+         {
+            static const unsigned short ord[] = {2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 0, 4, 7, 8, 20, 0, 0}; // 18 entries
+
+            BOOST_ASSERT(icode<18);
+
+            return ord[icode];
+         } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<53> const&)
+
+         template<typename RealType>
+         inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<64>&)
+        {
+           // method ================>>>       {1, 1, 1, 1, 1,  1,  1,  1,  2,  2,  2,  3, 4,  4,  4,  4,  5, 6}
+           static const unsigned short ord[] = {3, 4, 5, 6, 8, 11, 13, 19, 10, 20, 30,  0, 7, 10, 11, 23,  0, 0}; // 18 entries
+
+          BOOST_ASSERT(icode<18);
+
+          return ord[icode];
+        } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<64> const&)
+
+         template<typename RealType, typename Policy>
+         inline unsigned short owens_t_get_order(const unsigned short icode, RealType r, const Policy&)
+         {
+            typedef typename policies::precision<RealType, Policy>::type precision_type;
+            typedef typename mpl::if_<
+               mpl::or_<
+                  mpl::less_equal<precision_type, mpl::int_<0> >,
+                  mpl::greater<precision_type, mpl::int_<53> >
+               >,
+               mpl::int_<64>,
+               mpl::int_<53>
+            >::type tag_type;
+
+            return owens_t_get_order_imp(icode, r, tag_type());
+         }
+
+         // compute the value of Owen's T function with method T1 from the reference paper
+         template<typename RealType>
+         inline RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)
+         {
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+            const RealType hs = -h*h*half<RealType>();
+            const RealType dhs = exp( hs );
+            const RealType as = a*a;
+
+            unsigned short j=1;
+            RealType jj = 1;
+            RealType aj = a * one_div_two_pi<RealType>();
+            RealType dj = expm1( hs );
+            RealType gj = hs*dhs;
+
+            RealType val = atan( a ) * one_div_two_pi<RealType>();
+
+            while( true )
+            {
+               val += dj*aj/jj;
+
+               if( m <= j )
+                  break;
+
+               j++;
+               jj += static_cast<RealType>(2);
+               aj *= as;
+               dj = gj - dj;
+               gj *= hs / static_cast<RealType>(j);
+            } // while( true )
+
+            return val;
+         } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)
+
+         // compute the value of Owen's T function with method T2 from the reference paper
+         template<typename RealType, class Policy>
+         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::false_&)
+         {
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+            const unsigned short maxii = m+m+1;
+            const RealType hs = h*h;
+            const RealType as = -a*a;
+            const RealType y = static_cast<RealType>(1) / hs;
+
+            unsigned short ii = 1;
+            RealType val = 0;
+            RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
+            RealType z = owens_t_znorm1(ah)/h;
+
+            while( true )
+            {
+               val += z;
+               if( maxii <= ii )
+               {
+                  val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
+                  break;
+               } // if( maxii <= ii )
+               z = y * ( vi - static_cast<RealType>(ii) * z );
+               vi *= as;
+               ii += 2;
+            } // while( true )
+
+            return val;
+         } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
+
+         // compute the value of Owen's T function with method T3 from the reference paper
+         template<typename RealType>
+         inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<53>&)
+         {
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+	    const unsigned short m = 20;
+
+            static const RealType c2[] =
+            {
+               0.99999999999999987510,
+               -0.99999999999988796462,      0.99999999998290743652,
+               -0.99999999896282500134,      0.99999996660459362918,
+               -0.99999933986272476760,      0.99999125611136965852,
+               -0.99991777624463387686,      0.99942835555870132569,
+               -0.99697311720723000295,      0.98751448037275303682,
+               -0.95915857980572882813,      0.89246305511006708555,
+               -0.76893425990463999675,      0.58893528468484693250,
+               -0.38380345160440256652,      0.20317601701045299653,
+               -0.82813631607004984866E-01,  0.24167984735759576523E-01,
+               -0.44676566663971825242E-02,  0.39141169402373836468E-03
+            };
+
+            const RealType as = a*a;
+            const RealType hs = h*h;
+            const RealType y = static_cast<RealType>(1)/hs;
+
+            RealType ii = 1;
+            unsigned short i = 0;
+            RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
+            RealType zi = owens_t_znorm1(ah)/h;
+            RealType val = 0;
+
+            while( true )
+            {
+               BOOST_ASSERT(i < 21);
+               val += zi*c2[i];
+               if( m <= i ) // if( m < i+1 )
+               {
+                  val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
+                  break;
+               } // if( m < i )
+               zi = y * (ii*zi - vi);
+               vi *= as;
+               ii += 2;
+               i++;
+            } // while( true )
+
+            return val;
+         } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
+
+        // compute the value of Owen's T function with method T3 from the reference paper
+        template<class RealType>
+        inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<64>&)
+        {
+          BOOST_MATH_STD_USING
+          using namespace boost::math::constants;
+          
+          const unsigned short m = 30;
+
+          static const RealType c2[] =
+          {
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142),
+            BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4),
+              BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6)
+          };
+
+          const RealType as = a*a;
+          const RealType hs = h*h;
+          const RealType y = 1 / hs;
+
+          RealType ii = 1;
+          unsigned short i = 0;
+          RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
+          RealType zi = owens_t_znorm1(ah)/h;
+          RealType val = 0;
+
+          while( true )
+          {
+              BOOST_ASSERT(i < 31);
+              val += zi*c2[i];
+              if( m <= i ) // if( m < i+1 )
+              {
+                val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
+                break;
+              } // if( m < i )
+              zi = y * (ii*zi - vi);
+              vi *= as;
+              ii += 2;
+              i++;
+          } // while( true )
+
+          return val;
+        } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
+
+        template<class RealType, class Policy>
+        inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy&)
+        {
+            typedef typename policies::precision<RealType, Policy>::type precision_type;
+            typedef typename mpl::if_<
+               mpl::or_<
+                  mpl::less_equal<precision_type, mpl::int_<0> >,
+                  mpl::greater<precision_type, mpl::int_<53> >
+               >,
+               mpl::int_<64>,
+               mpl::int_<53>
+            >::type tag_type;
+
+            return owens_t_T3_imp(h, a, ah, tag_type());
+        }
+
+         // compute the value of Owen's T function with method T4 from the reference paper
+         template<typename RealType>
+         inline RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
+         {
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+            const unsigned short maxii = m+m+1;
+            const RealType hs = h*h;
+            const RealType as = -a*a;
+
+            unsigned short ii = 1;
+            RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>();
+            RealType yi = 1;
+            RealType val = 0;
+
+            while( true )
+            {
+               val += ai*yi;
+               if( maxii <= ii )
+                  break;
+               ii += 2;
+               yi = (static_cast<RealType>(1)-hs*yi) / static_cast<RealType>(ii);
+               ai *= as;
+            } // while( true )
+
+            return val;
+         } // RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
+
+         // compute the value of Owen's T function with method T5 from the reference paper
+         template<typename RealType>
+         inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<53>&)
+         {
+            BOOST_MATH_STD_USING
+            /*
+               NOTICE:
+               - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
+                 polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
+                 quadrature, because T5(h,a,m) contains only x^2 terms.
+               - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
+                 of 1/(2*pi) according to T5(h,a,m).
+             */
+
+            const unsigned short m = 13;
+            static const RealType pts[] = {0.35082039676451715489E-02,
+               0.31279042338030753740E-01,  0.85266826283219451090E-01,
+               0.16245071730812277011,      0.25851196049125434828,
+               0.36807553840697533536,      0.48501092905604697475,
+               0.60277514152618576821,      0.71477884217753226516,
+               0.81475510988760098605,      0.89711029755948965867,
+               0.95723808085944261843,      0.99178832974629703586};
+            static const RealType wts[] = { 0.18831438115323502887E-01,
+               0.18567086243977649478E-01,  0.18042093461223385584E-01,
+               0.17263829606398753364E-01,  0.16243219975989856730E-01,
+               0.14994592034116704829E-01,  0.13535474469662088392E-01,
+               0.11886351605820165233E-01,  0.10070377242777431897E-01,
+               0.81130545742299586629E-02,  0.60419009528470238773E-02,
+               0.38862217010742057883E-02,  0.16793031084546090448E-02};
+
+            const RealType as = a*a;
+            const RealType hs = -h*h*boost::math::constants::half<RealType>();
+
+            RealType val = 0;
+            for(unsigned short i = 0; i < m; ++i)
+            {
+               BOOST_ASSERT(i < 13);
+               const RealType r = static_cast<RealType>(1) + as*pts[i];
+               val += wts[i] * exp( hs*r ) / r;
+            } // for(unsigned short i = 0; i < m; ++i)
+
+            return val*a;
+         } // RealType owens_t_T5(const RealType h, const RealType a)
+
+        // compute the value of Owen's T function with method T5 from the reference paper
+        template<typename RealType>
+        inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<64>&)
+        {
+          BOOST_MATH_STD_USING
+            /*
+              NOTICE:
+              - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
+              polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
+              quadrature, because T5(h,a,m) contains only x^2 terms.
+              - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
+              of 1/(2*pi) according to T5(h,a,m).
+            */
+
+          const unsigned short m = 19;
+          static const RealType pts[] = {
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321)
+          };
+          static const RealType wts[] = {
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947),
+               BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578)
+          };
+
+          const RealType as = a*a;
+          const RealType hs = -h*h*boost::math::constants::half<RealType>();
+
+          RealType val = 0;
+          for(unsigned short i = 0; i < m; ++i)
+            {
+              BOOST_ASSERT(i < 19);
+              const RealType r = 1 + as*pts[i];
+              val += wts[i] * exp( hs*r ) / r;
+            } // for(unsigned short i = 0; i < m; ++i)
+
+          return val*a;
+        } // RealType owens_t_T5(const RealType h, const RealType a)
+
+        template<class RealType, class Policy>
+        inline RealType owens_t_T5(const RealType h, const RealType a, const Policy&)
+        {
+            typedef typename policies::precision<RealType, Policy>::type precision_type;
+            typedef typename mpl::if_<
+               mpl::or_<
+                  mpl::less_equal<precision_type, mpl::int_<0> >,
+                  mpl::greater<precision_type, mpl::int_<53> >
+               >,
+               mpl::int_<64>,
+               mpl::int_<53>
+            >::type tag_type;
+
+            return owens_t_T5_imp(h, a, tag_type());
+        }
+
+
+         // compute the value of Owen's T function with method T6 from the reference paper
+         template<typename RealType>
+         inline RealType owens_t_T6(const RealType h, const RealType a)
+         {
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+            const RealType normh = owens_t_znorm2( h );
+            const RealType y = static_cast<RealType>(1) - a;
+            const RealType r = atan2(y, static_cast<RealType>(1 + a) );
+
+            RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>();
+
+            if( r != 0 )
+               val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>();
+
+            return val;
+         } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)
+
+         template <class T, class Policy>
+         std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)
+         {
+            //
+            // This is the same series as T1, but:
+            // * The Taylor series for atan has been combined with that for T1, 
+            //   reducing but not eliminating cancellation error.
+            // * The resulting alternating series is then accelerated using method 1
+            //   from H. Cohen, F. Rodriguez Villegas, D. Zagier, 
+            //   "Convergence acceleration of alternating series", Bonn, (1991).
+            //
+            BOOST_MATH_STD_USING
+            static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)";
+            T half_h_h = h * h / 2;
+            T a_pow = a;
+            T aa = a * a;
+            T exp_term = exp(-h * h / 2);
+            T one_minus_dj_sum = exp_term; 
+            T sum = a_pow * exp_term;
+            T dj_pow = exp_term;
+            T term = sum;
+            T abs_err;
+            int j = 1;
+
+            //
+            // Normally with this form of series acceleration we can calculate
+            // up front how many terms will be required - based on the assumption
+            // that each term decreases in size by a factor of 3.  However,
+            // that assumption does not apply here, as the underlying T1 series can 
+            // go quite strongly divergent in the early terms, before strongly
+            // converging later.  Various "guestimates" have been tried to take account
+            // of this, but they don't always work.... so instead set "n" to the 
+            // largest value that won't cause overflow later, and abort iteration
+            // when the last accelerated term was small enough...
+            //
+            int n;
+            try
+            {
+               n = itrunc(T(tools::log_max_value<T>() / 6));
+            }
+            catch(...)
+            {
+               n = (std::numeric_limits<int>::max)();
+            }
+            n = (std::min)(n, 1500);
+            T d = pow(3 + sqrt(T(8)), n);
+            d = (d + 1 / d) / 2;
+            T b = -1;
+            T c = -d;
+            c = b - c;
+            sum *= c;
+            b = -n * n * b * 2;
+            abs_err = ldexp(fabs(sum), -tools::digits<T>());
+
+            while(j < n)
+            {
+               a_pow *= aa;
+               dj_pow *= half_h_h / j;
+               one_minus_dj_sum += dj_pow;
+               term = one_minus_dj_sum * a_pow / (2 * j + 1);
+               c = b - c;
+               sum += c * term;
+               abs_err += ldexp(std::max(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>());
+               b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));
+               ++j;
+               //
+               // Include an escape route to prevent calculating too many terms:
+               //
+               if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term)))
+                  break;
+            }
+            abs_err += fabs(c * term);
+            if(sum < 0)  // sum must always be positive, if it's negative something really bad has happend:
+               policies::raise_evaluation_error(function, 0, T(0), pol);
+            return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum);
+         }
+
+         template<typename RealType, class Policy>
+         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::true_&)
+         {
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+            const unsigned short maxii = m+m+1;
+            const RealType hs = h*h;
+            const RealType as = -a*a;
+            const RealType y = static_cast<RealType>(1) / hs;
+
+            unsigned short ii = 1;
+            RealType val = 0;
+            RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
+            RealType z = owens_t_znorm1(ah)/h;
+            RealType last_z = fabs(z);
+            RealType lim = policies::get_epsilon<RealType, Policy>();
+
+            while( true )
+            {
+               val += z;
+               //
+               // This series stops converging after a while, so put a limit
+               // on how far we go before returning our best guess:
+               //
+               if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))
+               {
+                  val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
+                  break;
+               } // if( maxii <= ii )
+               last_z = fabs(z);
+               z = y * ( vi - static_cast<RealType>(ii) * z );
+               vi *= as;
+               ii += 2;
+            } // while( true )
+
+            return val;
+         } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
+
+         template<typename RealType, class Policy>
+         inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
+         {
+            //
+            // This is the same series as T2, but with acceleration applied.
+            // Note that we have to be *very* careful to check that nothing bad
+            // has happened during evaluation - this series will go divergent
+            // and/or fail to alternate at a drop of a hat! :-(
+            //
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+            const RealType hs = h*h;
+            const RealType as = -a*a;
+            const RealType y = static_cast<RealType>(1) / hs;
+
+            unsigned short ii = 1;
+            RealType val = 0;
+            RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
+            RealType z = boost::math::detail::owens_t_znorm1(ah)/h;
+            RealType last_z = fabs(z);
+
+            //
+            // Normally with this form of series acceleration we can calculate
+            // up front how many terms will be required - based on the assumption
+            // that each term decreases in size by a factor of 3.  However,
+            // that assumption does not apply here, as the underlying T1 series can 
+            // go quite strongly divergent in the early terms, before strongly
+            // converging later.  Various "guestimates" have been tried to take account
+            // of this, but they don't always work.... so instead set "n" to the 
+            // largest value that won't cause overflow later, and abort iteration
+            // when the last accelerated term was small enough...
+            //
+            int n;
+            try
+            {
+               n = itrunc(RealType(tools::log_max_value<RealType>() / 6));
+            }
+            catch(...)
+            {
+               n = (std::numeric_limits<int>::max)();
+            }
+            n = (std::min)(n, 1500);
+            RealType d = pow(3 + sqrt(RealType(8)), n);
+            d = (d + 1 / d) / 2;
+            RealType b = -1;
+            RealType c = -d;
+            int s = 1;
+
+            for(int k = 0; k < n; ++k)
+            {
+               //
+               // Check for both convergence and whether the series has gone bad:
+               //
+               if(
+                  (fabs(z) > last_z)     // Series has gone divergent, abort
+                  || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z))  // Convergence!
+                  || (z * s < 0)         // Series has stopped alternating - all bets are off - abort.
+                  )
+               {
+                  break;
+               }
+               c = b - c;
+               val += c * s * z;
+               b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1));
+               last_z = fabs(z);
+               s = -s;
+               z = y * ( vi - static_cast<RealType>(ii) * z );
+               vi *= as;
+               ii += 2;
+            } // while( true )
+            RealType err = fabs(c * z) / val;
+            return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err);
+         } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
+
+         template<typename RealType, typename Policy>
+         inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol)
+         {
+            BOOST_MATH_STD_USING
+            
+            const RealType hs = h*h;
+            const RealType as = -a*a;
+
+            unsigned short ii = 1;
+            RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) );
+            RealType yi = 1.0;
+            RealType val = 0.0;
+
+            RealType lim = boost::math::policies::get_epsilon<RealType, Policy>();
+
+            while( true )
+            {
+               RealType term = ai*yi;
+               val += term;
+               if((yi != 0) && (fabs(val * lim) > fabs(term)))
+                  break;
+               ii += 2;
+               yi = (1.0-hs*yi) / static_cast<RealType>(ii);
+               ai *= as;
+               if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))
+                  policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);
+            } // while( true )
+
+            return val;
+         } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)
+
+
+         // This routine dispatches the call to one of six subroutines, depending on the values
+         // of h and a.
+         // preconditions: h >= 0, 0<=a<=1, ah=a*h
+         //
+         // Note there are different versions for different precisions....
+         template<typename RealType, typename Policy>
+         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, mpl::int_<64> const&)
+         {
+            // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper:
+            BOOST_MATH_STD_USING
+            //
+            // Handle some special cases first, these are from
+            // page 1077 of Owen's original paper:
+            //
+            if(h == 0)
+            {
+               return atan(a) * constants::one_div_two_pi<RealType>();
+            }
+            if(a == 0)
+            {
+               return 0;
+            }
+            if(a == 1)
+            {
+               return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
+            }
+            if(a >= tools::max_value<RealType>())
+            {
+               return owens_t_znorm2(RealType(fabs(h)));
+            }
+            RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case
+            const unsigned short icode = owens_t_compute_code(h, a);
+            const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);
+            static const unsigned short meth[] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6}; // 18 entries
+
+            // determine the appropriate method, T1 ... T6
+            switch( meth[icode] )
+            {
+            case 1: // T1
+               val = owens_t_T1(h,a,m);
+               break;
+            case 2: // T2
+               typedef typename policies::precision<RealType, Policy>::type precision_type;
+               typedef mpl::bool_<(precision_type::value == 0) || (precision_type::value > 64)> tag_type;
+               val = owens_t_T2(h, a, m, ah, pol, tag_type());
+               break;
+            case 3: // T3
+               val = owens_t_T3(h,a,ah, pol);
+               break;
+            case 4: // T4
+               val = owens_t_T4(h,a,m);
+               break;
+            case 5: // T5
+               val = owens_t_T5(h,a, pol);
+               break;
+            case 6: // T6
+               val = owens_t_T6(h,a);
+               break;
+            default:
+               BOOST_THROW_EXCEPTION(std::logic_error("selection routine in Owen's T function failed"));
+            }
+            return val;
+         }
+
+         template<typename RealType, typename Policy>
+         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<65>&)
+         {
+            // Arbitrary precision version:
+            BOOST_MATH_STD_USING
+            //
+            // Handle some special cases first, these are from
+            // page 1077 of Owen's original paper:
+            //
+            if(h == 0)
+            {
+               return atan(a) * constants::one_div_two_pi<RealType>();
+            }
+            if(a == 0)
+            {
+               return 0;
+            }
+            if(a == 1)
+            {
+               return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
+            }
+            if(a >= tools::max_value<RealType>())
+            {
+               return owens_t_znorm2(RealType(fabs(h)));
+            }
+            // Attempt arbitrary precision code, this will throw if it goes wrong:
+            typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;
+            std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>());
+            RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000;
+            bool have_t1(false), have_t2(false);
+            if(ah < 3)
+            {
+               try
+               {
+                  have_t1 = true;
+                  p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
+                  if(p1.second < target_precision)
+                     return p1.first;
+               }
+               catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
+            }
+            if(ah > 1)
+            {
+               try
+               {
+                  have_t2 = true;
+                  p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
+                  if(p2.second < target_precision)
+                     return p2.first;
+               }
+               catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
+            }
+            //
+            // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
+            // is fairly low compared to T4.
+            //
+            if(!have_t1)
+            {
+               try
+               {
+                  have_t1 = true;
+                  p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
+                  if(p1.second < target_precision)
+                     return p1.first;
+               }
+               catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
+            }
+            //
+            // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
+            // is fairly low compared to T4.
+            //
+            if(!have_t2)
+            {
+               try
+               {
+                  have_t2 = true;
+                  p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
+                  if(p2.second < target_precision)
+                     return p2.first;
+               }
+               catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
+            }
+            //
+            // OK, nothing left to do but try the most expensive option which is T4,
+            // this is often slow to converge, but when it does converge it tends to
+            // be accurate:
+            try
+            {
+               return T4_mp(h, a, pol);
+            }
+            catch(const boost::math::evaluation_error&){}  // T4 may fail and throw, that's OK
+            //
+            // Now look back at the results from T1 and T2 and see if either gave better
+            // results than we could get from the 64-bit precision versions.
+            //
+            if((std::min)(p1.second, p2.second) < 1e-20)
+            {
+               return p1.second < p2.second ? p1.first : p2.first;
+            }
+            //
+            // We give up - no arbitrary precision versions succeeded!
+            //
+            return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
+         } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah)
+         template<typename RealType, typename Policy>
+         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<0>&)
+         {
+            // We don't know what the precision is until runtime:
+            if(tools::digits<RealType>() <= 64)
+               return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
+            return owens_t_dispatch(h, a, ah, pol, mpl::int_<65>());
+         }
+         template<typename RealType, typename Policy>
+         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol)
+         {
+            // Figure out the precision and forward to the correct version:
+            typedef typename policies::precision<RealType, Policy>::type precision_type;
+            typedef typename mpl::if_c<
+               precision_type::value == 0,
+               mpl::int_<0>,
+               typename mpl::if_c<
+                  precision_type::value <= 64,
+                  mpl::int_<64>,
+                  mpl::int_<65>
+               >::type
+            >::type tag_type;
+            return owens_t_dispatch(h, a, ah, pol, tag_type());
+         }
+         // compute Owen's T function, T(h,a), for arbitrary values of h and a
+         template<typename RealType, class Policy>
+         inline RealType owens_t(RealType h, RealType a, const Policy& pol)
+         {
+            BOOST_MATH_STD_USING
+            // exploit that T(-h,a) == T(h,a)
+            h = fabs(h);
+
+            // Use equation (2) in the paper to remap the arguments
+            // such that h>=0 and 0<=a<=1 for the call of the actual
+            // computation routine.
+
+            const RealType fabs_a = fabs(a);
+            const RealType fabs_ah = fabs_a*h;
+
+            RealType val = 0.0; // avoid compiler warnings, 0.0 will be overwritten in any case
+
+            if(fabs_a <= 1)
+            {
+               val = owens_t_dispatch(h, fabs_a, fabs_ah, pol);
+            } // if(fabs_a <= 1.0)
+            else 
+            {
+               if( h <= 0.67 )
+               {
+                  const RealType normh = owens_t_znorm1(h);
+                  const RealType normah = owens_t_znorm1(fabs_ah);
+                  val = static_cast<RealType>(1)/static_cast<RealType>(4) - normh*normah -
+                     owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
+               } // if( h <= 0.67 )
+               else
+               {
+                  const RealType normh = detail::owens_t_znorm2(h);
+                  const RealType normah = detail::owens_t_znorm2(fabs_ah);
+                  val = constants::half<RealType>()*(normh+normah) - normh*normah -
+                     owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
+               } // else [if( h <= 0.67 )]
+            } // else [if(fabs_a <= 1)]
+
+            // exploit that T(h,-a) == -T(h,a)
+            if(a < 0)
+            {
+               return -val;
+            } // if(a < 0)
+
+            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 mpl::int_<N>&){}
+               static void do_init(const mpl::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>
+      inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol)
+      {
+         typedef typename tools::promote_args<T1, T2>::type result_type;
+         typedef typename policies::evaluation<result_type, Policy>::type value_type;
+         typedef typename policies::precision<value_type, Policy>::type precision_type;
+         typedef typename mpl::if_c<
+               precision_type::value == 0,
+               mpl::int_<0>,
+               typename mpl::if_c<
+                  precision_type::value <= 64,
+                  mpl::int_<64>,
+                  mpl::int_<65>
+               >::type
+            >::type 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%)");
+      }
+
+      template <class T1, class T2>
+      inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a)
+      {
+         return owens_t(h, a, policies::policy<>());
+      }
+
+
+   } // namespace math
+} // namespace boost
+
+#endif
+// EOF
diff --git a/boost/math/special_functions/round.hpp b/boost/math/special_functions/round.hpp
index f11c6aeb1d..2b4497e198 100644
--- a/boost/math/special_functions/round.hpp
+++ b/boost/math/special_functions/round.hpp
@@ -76,7 +76,7 @@ inline boost::long_long_type llround(const T& v, const Policy& pol)
    BOOST_MATH_STD_USING
    T r = boost::math::round(v, pol);
    if((r > (std::numeric_limits<boost::long_long_type>::max)()) || (r < (std::numeric_limits<boost::long_long_type>::min)()))
-      return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", 0, v, 0LL, pol));
+      return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", 0, v, static_cast<boost::long_long_type>(0), pol));
    return static_cast<boost::long_long_type>(r);
 }
 template <class T>
diff --git a/boost/math/special_functions/trunc.hpp b/boost/math/special_functions/trunc.hpp
index 520ae89f5d..7346afe6d1 100644
--- a/boost/math/special_functions/trunc.hpp
+++ b/boost/math/special_functions/trunc.hpp
@@ -76,7 +76,7 @@ inline boost::long_long_type lltrunc(const T& v, const Policy& pol)
    BOOST_MATH_STD_USING
    T r = boost::math::trunc(v, pol);
    if((r > (std::numeric_limits<boost::long_long_type>::max)()) || (r < (std::numeric_limits<boost::long_long_type>::min)()))
-      return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::lltrunc<%1%>(%1%)", 0, v, 0LL, pol));
+      return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::lltrunc<%1%>(%1%)", 0, v, static_cast<boost::long_long_type>(0), pol));
    return static_cast<boost::long_long_type>(r);
 }
 template <class T>
diff --git a/boost/math/special_functions/zeta.hpp b/boost/math/special_functions/zeta.hpp
index 4ed5f6a705..011182718e 100644
--- a/boost/math/special_functions/zeta.hpp
+++ b/boost/math/special_functions/zeta.hpp
@@ -542,7 +542,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<64>&)
 }
 
 template <class T, class Policy>
-T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<113>&)
+T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<113>&)
 {
    BOOST_MATH_STD_USING
    T result;
@@ -896,6 +896,49 @@ 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 mpl::int_<0>&){}
+      static void do_init(const mpl::int_<53>&){}
+      static void do_init(const mpl::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());
+      }
+      static void do_init(const mpl::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), Policy());
+         boost::math::zeta(static_cast<T>(70), 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>
@@ -929,6 +972,8 @@ inline typename tools::promote_args<T>::type zeta(T s, const Policy&)
    >::type tag_type;
    //typedef mpl::int_<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)),