erf.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 1,085 行 · 第 1/3 页

//  (C) Copyright John Maddock 2006.
//  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_SPECIAL_ERF_HPP
#define BOOST_MATH_SPECIAL_ERF_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/tools/config.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/math/policies/error_handling.hpp>

namespace boost{ namespace math{

namespace detail
{

//
// Asymptotic series for large z:
//
template <class T>
struct erf_asympt_series_t
{
   erf_asympt_series_t(T z) : xx(2 * -z * z), tk(1)
   {
      BOOST_MATH_STD_USING
      result = -exp(-z * z) / sqrt(boost::math::constants::pi<T>());
      result /= z;
   }

   typedef T result_type;

   T operator()()
   {
      BOOST_MATH_STD_USING
      T r = result;
      result *= tk / xx;
      tk += 2;
      if( fabs(r) < fabs(result))
         result = 0;
      return r;
   }
private:
   T result;
   T xx;
   int tk;
};
//
// How large z has to be in order to ensure that the series converges:
//
template <class T>
inline float erf_asymptotic_limit_N(const T&)
{
   return (std::numeric_limits<float>::max)();
}
inline float erf_asymptotic_limit_N(const mpl::int_<24>&)
{
   return 2.8F;
}
inline float erf_asymptotic_limit_N(const mpl::int_<53>&)
{
   return 4.3F;
}
inline float erf_asymptotic_limit_N(const mpl::int_<64>&)
{
   return 4.8F;
}
inline float erf_asymptotic_limit_N(const mpl::int_<106>&)
{
   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(-z, invert, pol, t);
      else
         return 1 + erf_imp(-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::digits<T, Policy>(), max_iter, 1);
      policies::check_series_iterations("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::digits<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(-z, invert, pol, t);
      else if(z < -0.5)
         return 2 - erf_imp(-z, invert, pol, t);
      else
         return 1 + erf_imp(-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)
   {
      //
      // We're going to calculate erf:
      //
      if(z == 0)
      {
         result = T(0);
      }
      else if(z < 1e-10)
      {
         result = static_cast<T>(z * 1.125f + z * 0.003379167095512573896158903121545171688L);
      }
      else
      {
         // Maximum Deviation Found:                     1.561e-17
         // Expected Error Term:                         1.561e-17
         // Maximum Relative Change in Control Points:   1.155e-04
         // Max Error found at double precision =        2.961182e-17

         static const T Y = 1.044948577880859375f;
         static const T P[] = {    
            0.0834305892146531832907L,
            -0.338165134459360935041L,
            -0.0509990735146777432841L,
            -0.00772758345802133288487L,
            -0.000322780120964605683831L,
         };
         static const T Q[] = {    
            1L,
            0.455004033050794024546L,
            0.0875222600142252549554L,
            0.00858571925074406212772L,
            0.000370900071787748000569L,
         };
         result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
      }
   }
   else if((z < 14) || ((z < 28) && invert))
   {
      //
      // We'll be calculating erfc:
      //
      invert = !invert;
      if(z < 1.5f)
      {
         // Maximum Deviation Found:                     3.702e-17
         // Expected Error Term:                         3.702e-17
         // Maximum Relative Change in Control Points:   2.845e-04
         // Max Error found at double precision =        4.841816e-17
         static const T Y = 0.405935764312744140625f;
         static const T P[] = {    
            -0.098090592216281240205L,
            0.178114665841120341155L,
            0.191003695796775433986L,
            0.0888900368967884466578L,
            0.0195049001251218801359L,
            0.00180424538297014223957L,
         };
         static const T Q[] = {    
            1L,
            1.84759070983002217845L,
            1.42628004845511324508L,
            0.578052804889902404909L,
            0.12385097467900864233L,
            0.0113385233577001411017L,
            0.337511472483094676155e-5L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 0.5) / tools::evaluate_polynomial(Q, z - 0.5);
         result *= exp(-z * z) / z;
      }
      else if(z < 2.5f)
      {
         // Max Error found at double precision =        6.599585e-18
         // Maximum Deviation Found:                     3.909e-18
         // Expected Error Term:                         3.909e-18
         // Maximum Relative Change in Control Points:   9.886e-05
         static const T Y = 0.50672817230224609375f;
         static const T P[] = {    
            -0.0243500476207698441272L,
            0.0386540375035707201728L,
            0.04394818964209516296L,
            0.0175679436311802092299L,
            0.00323962406290842133584L,
            0.000235839115596880717416L,
         };
         static const T Q[] = {    
            1L,
            1.53991494948552447182L,
            0.982403709157920235114L,
            0.325732924782444448493L,
            0.0563921837420478160373L,
            0.00410369723978904575884L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 1.5) / tools::evaluate_polynomial(Q, z - 1.5);
         result *= exp(-z * z) / z;
      }
      else if(z < 4.5f)
      {
         // Maximum Deviation Found:                     1.512e-17
         // Expected Error Term:                         1.512e-17
         // Maximum Relative Change in Control Points:   2.222e-04
         // Max Error found at double precision =        2.062515e-17
         static const T Y = 0.5405750274658203125f;
         static const T P[] = {    
            0.00295276716530971662634L,
            0.0137384425896355332126L,
            0.00840807615555585383007L,
            0.00212825620914618649141L,
            0.000250269961544794627958L,
            0.113212406648847561139e-4L,
         };
         static const T Q[] = {    
            1L,
            1.04217814166938418171L,
            0.442597659481563127003L,
            0.0958492726301061423444L,
            0.0105982906484876531489L,
            0.000479411269521714493907L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 3.5) / tools::evaluate_polynomial(Q, z - 3.5);
         result *= exp(-z * z) / z;
      }
      else
      {
         // Max Error found at double precision =        2.997958e-17
         // Maximum Deviation Found:                     2.860e-17
         // Expected Error Term:                         2.859e-17
         // Maximum Relative Change in Control Points:   1.357e-05
         static const T Y = 0.5579090118408203125f;
         static const T P[] = {    
            0.00628057170626964891937L,
            0.0175389834052493308818L,
            -0.212652252872804219852L,
            -0.687717681153649930619L,
            -2.5518551727311523996L,
            -3.22729451764143718517L,
            -2.8175401114513378771L,
         };
         static const T Q[] = {    
            1L,
            2.79257750980575282228L,
            11.0567237927800161565L,
            15.930646027911794143L,
            22.9367376522880577224L,
            13.5064170191802889145L,
            5.48409182238641741584L,
         };
         result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
         result *= exp(-z * z) / z;
      }
   }
   else
   {
      //
      // Any value of z larger than 28 will underflow to zero:
      //
      result = 0;
      invert = !invert;
   }

   if(invert)
   {
      result = 1 - result;
   }

   return result;
} // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<53>& t)