From bb4dd8289b351fae6b55e303f189127a394a1edd Mon Sep 17 00:00:00 2001
From: Anas Nashif <anas.nashif@intel.com>
Date: Mon, 26 Aug 2013 08:15:55 -0400
Subject: Imported Upstream version 1.51.0

---
 boost/math/special_functions/detail/bessel_i0.hpp  | 30 ++++++++++++++++
 boost/math/special_functions/detail/bessel_i1.hpp  | 31 +++++++++++++++++
 boost/math/special_functions/detail/bessel_j0.hpp  | 30 ++++++++++++++++
 boost/math/special_functions/detail/bessel_j1.hpp  | 30 ++++++++++++++++
 boost/math/special_functions/detail/bessel_jy.hpp  | 16 +++++----
 boost/math/special_functions/detail/bessel_k0.hpp  | 30 ++++++++++++++++
 boost/math/special_functions/detail/bessel_k1.hpp  | 30 ++++++++++++++++
 boost/math/special_functions/detail/bessel_y0.hpp  | 30 ++++++++++++++++
 boost/math/special_functions/detail/bessel_y1.hpp  | 30 ++++++++++++++++
 boost/math/special_functions/detail/erf_inv.hpp    | 40 ++++++++++++++++++++++
 boost/math/special_functions/detail/iconv.hpp      |  2 +-
 .../math/special_functions/detail/igamma_large.hpp |  1 -
 .../math/special_functions/detail/lanczos_sse2.hpp |  5 ++-
 .../math/special_functions/detail/lgamma_small.hpp |  8 +++++
 14 files changed, 304 insertions(+), 9 deletions(-)

(limited to 'boost/math/special_functions/detail')

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
@@ -20,9 +20,39 @@
 
 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
@@ -20,9 +20,40 @@
 
 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
@@ -21,9 +21,39 @@
 
 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
@@ -21,9 +21,39 @@
 
 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
@@ -21,11 +21,41 @@
 
 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
@@ -21,9 +21,39 @@
 
 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
@@ -23,9 +23,39 @@
 
 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
@@ -23,9 +23,39 @@
 
 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
@@ -14,6 +14,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:
 //
-- 
cgit v1.2.3