view release on metacpan or  search on metacpan

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

{
   return 6.5F;
}
inline float erf_asymptotic_limit_N(const mpl::int_<113>&)
{
   return 6.8F;
}

template <class T, class Policy>
inline T erf_asymptotic_limit()
{
   typedef typename policies::precision<T, Policy>::type precision_type;
   typedef typename mpl::if_<
      mpl::less_equal<precision_type, mpl::int_<24> >,
      typename mpl::if_<
         mpl::less_equal<precision_type, mpl::int_<0> >,
         mpl::int_<0>,
         mpl::int_<24>
      >::type,
      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>,
            typename mpl::if_<
               mpl::less_equal<precision_type, mpl::int_<106> >,
               mpl::int_<106>,
               typename mpl::if_<
                  mpl::less_equal<precision_type, mpl::int_<113> >,
                  mpl::int_<113>,
                  mpl::int_<0>
               >::type
            >::type
         >::type
      >::type
   >::type tag_type;
   return erf_asymptotic_limit_N(tag_type());
}

template <class T, class Policy, class Tag>
T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
{
   BOOST_MATH_STD_USING

   BOOST_MATH_INSTRUMENT_CODE("Generic erf_imp called");

   if(z < 0)
   {
      if(!invert)
         return -erf_imp(T(-z), invert, pol, t);
      else
         return 1 + erf_imp(T(-z), false, pol, t);
   }

   T result;

   if(!invert && (z > detail::erf_asymptotic_limit<T, Policy>()))
   {
      detail::erf_asympt_series_t<T> s(z);
      boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
      result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, 1);
      policies::check_series_iterations<T>("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol);
   }
   else
   {
      T x = z * z;
      if(x < 0.6)
      {
         // Compute P:
         result = z * exp(-x);
         result /= sqrt(boost::math::constants::pi<T>());
         if(result != 0)
            result *= 2 * detail::lower_gamma_series(T(0.5f), x, pol);
      }
      else if(x < 1.1f)
      {
         // Compute Q:
         invert = !invert;
         result = tgamma_small_upper_part(T(0.5f), x, pol);
         result /= sqrt(boost::math::constants::pi<T>());
      }
      else
      {
         // Compute Q:
         invert = !invert;
         result = z * exp(-x);
         result /= sqrt(boost::math::constants::pi<T>());
         result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon<T, Policy>());
      }
   }
   if(invert)
      result = 1 - result;
   return result;
}

template <class T, class Policy>
T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
{
   BOOST_MATH_STD_USING

   BOOST_MATH_INSTRUMENT_CODE("53-bit precision erf_imp called");

   if(z < 0)
   {
      if(!invert)
         return -erf_imp(T(-z), invert, pol, t);
      else if(z < -0.5)
         return 2 - erf_imp(T(-z), invert, pol, t);
      else
         return 1 + erf_imp(T(-z), false, pol, t);
   }

   T result;

   //
   // Big bunch of selection statements now to pick
   // which implementation to use,
   // try to put most likely options first:
   //
   if(z < 0.5)
   {
      //