view release on metacpan or  search on metacpan

src/boost/math/special_functions/detail/bessel_jy_asym.hpp  view on Meta::CPAN

         mpl::equal_to<precision_type, mpl::int_<0> >,
         mpl::greater<precision_type, mpl::int_<113> > >,
      mpl::int_<0>,
      typename mpl::if_<
         mpl::greater<precision_type, mpl::int_<64> >,
         mpl::int_<113>,
         typename mpl::if_<
            mpl::greater<precision_type, mpl::int_<53> >,
            mpl::int_<64>,
            mpl::int_<53>
         >::type
      >::type
   >::type type;
};

template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<0>&)
{
   // default case:
   BOOST_MATH_STD_USING
   T v2 = (std::max)(T(3), T(v * v));
   return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f));
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<53>&)
{
   // double case:
   T v2 = (std::max)(T(3), T(v * v));
   return v2 * 33 /*73*/;
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<64>&)
{
   // 80-bit extended-double case:
   T v2 = (std::max)(T(3), T(v * v));
   return v2 * 121 /*266*/;
}
template <class T>
inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<113>&)
{
   // 128-bit long double case:
   T v2 = (std::max)(T(3), T(v * v));
   return v2 * 39154 /*85700*/;
}

template <class T, class Policy>
void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
{
   T c = 1;
   T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol);
   T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
   T f = (p - q) / v;
   T g_prefix = boost::math::sin_pi(v / 2, pol);
   g_prefix *= g_prefix * 2 / v;
   T g = f + g_prefix * q;
   T h = p;
   T c_mult = -x * x / 4;

   T y(c * g), y1(c * h);

   for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
   {
      f = (k * f + p + q) / (k*k - v*v);
      p /= k - v;
      q /= k + v;
      c *= c_mult / k;
      T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol);
      g = f + g_prefix * q;
      h = -k * g + p;
      y += c * g;
      y1 += c * h;
      if(c * g / tools::epsilon<T>() < y)
         break;
   }

   *Y = -y;
   *Y1 = (-2 / x) * y1;
}

template <class T, class Policy>
T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
{
   BOOST_MATH_STD_USING  // ADL of std names
   T s = 1;
   T mu = 4 * v * v;
   T ex = 8 * x;
   T num = mu - 1;
   T denom = ex;

   s -= num / denom;

   num *= mu - 9;
   denom *= ex * 2;
   s += num / denom;

   num *= mu - 25;
   denom *= ex * 3;
   s -= num / denom;

   // Try and avoid overflow to the last minute:
   T e = exp(x/2);

   s = e * (e * s / sqrt(2 * x * constants::pi<T>()));

   return (boost::math::isfinite)(s) ? 
      s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol);
}

}}} // namespaces

#endif