Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:40:06 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2014 Anton Bikineev
0003 //  Copyright 2014 Christopher Kormanyos
0004 //  Copyright 2014 John Maddock
0005 //  Copyright 2014 Paul Bristow
0006 //  Distributed under the Boost
0007 //  Software License, Version 1.0. (See accompanying file
0008 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0009 //
0010 #ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP
0011 #define BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP
0012 
0013 #include <cmath>
0014 #include <cstdint>
0015 #include <boost/math/tools/series.hpp>
0016 #include <boost/math/special_functions/gamma.hpp>
0017 #include <boost/math/special_functions/trunc.hpp>
0018 #include <boost/math/policies/error_handling.hpp>
0019 
0020   namespace boost { namespace math { namespace detail {
0021 
0022   // primary template for term of Taylor series
0023   template <class T, unsigned p, unsigned q>
0024   struct hypergeometric_pFq_generic_series_term;
0025 
0026   // partial specialization for 0F1
0027   template <class T>
0028   struct hypergeometric_pFq_generic_series_term<T, 0u, 1u>
0029   {
0030     typedef T result_type;
0031 
0032     hypergeometric_pFq_generic_series_term(const T& b, const T& z)
0033        : n(0), term(1), b(b), z(z)
0034     {
0035     }
0036 
0037     T operator()()
0038     {
0039       BOOST_MATH_STD_USING
0040       const T r = term;
0041       term *= ((1 / ((b + n) * (n + 1))) * z);
0042       ++n;
0043       return r;
0044     }
0045 
0046   private:
0047     unsigned n;
0048     T term;
0049     const T b, z;
0050   };
0051 
0052   // partial specialization for 1F0
0053   template <class T>
0054   struct hypergeometric_pFq_generic_series_term<T, 1u, 0u>
0055   {
0056     typedef T result_type;
0057 
0058     hypergeometric_pFq_generic_series_term(const T& a, const T& z)
0059        : n(0), term(1), a(a), z(z)
0060     {
0061     }
0062 
0063     T operator()()
0064     {
0065       BOOST_MATH_STD_USING
0066       const T r = term;
0067       term *= (((a + n) / (n + 1)) * z);
0068       ++n;
0069       return r;
0070     }
0071 
0072   private:
0073     unsigned n;
0074     T term;
0075     const T a, z;
0076   };
0077 
0078   // partial specialization for 1F1
0079   template <class T>
0080   struct hypergeometric_pFq_generic_series_term<T, 1u, 1u>
0081   {
0082     typedef T result_type;
0083 
0084     hypergeometric_pFq_generic_series_term(const T& a, const T& b, const T& z)
0085        : n(0), term(1), a(a), b(b), z(z)
0086     {
0087     }
0088 
0089     T operator()()
0090     {
0091       BOOST_MATH_STD_USING
0092       const T r = term;
0093       term *= (((a + n) / ((b + n) * (n + 1))) * z);
0094       ++n;
0095       return r;
0096     }
0097 
0098   private:
0099     unsigned n;
0100     T term;
0101     const T a, b, z;
0102   };
0103 
0104   // partial specialization for 1F2
0105   template <class T>
0106   struct hypergeometric_pFq_generic_series_term<T, 1u, 2u>
0107   {
0108     typedef T result_type;
0109 
0110     hypergeometric_pFq_generic_series_term(const T& a, const T& b1, const T& b2, const T& z)
0111        : n(0), term(1), a(a), b1(b1), b2(b2), z(z)
0112     {
0113     }
0114 
0115     T operator()()
0116     {
0117       BOOST_MATH_STD_USING
0118       const T r = term;
0119       term *= (((a + n) / ((b1 + n) * (b2 + n) * (n + 1))) * z);
0120       ++n;
0121       return r;
0122     }
0123 
0124   private:
0125     unsigned n;
0126     T term;
0127     const T a, b1, b2, z;
0128   };
0129 
0130   // partial specialization for 2F0
0131   template <class T>
0132   struct hypergeometric_pFq_generic_series_term<T, 2u, 0u>
0133   {
0134     typedef T result_type;
0135 
0136     hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& z)
0137        : n(0), term(1), a1(a1), a2(a2), z(z)
0138     {
0139     }
0140 
0141     T operator()()
0142     {
0143       BOOST_MATH_STD_USING
0144       const T r = term;
0145       term *= (((a1 + n) * (a2 + n) / (n + 1)) * z);
0146       ++n;
0147       return r;
0148     }
0149 
0150   private:
0151     unsigned n;
0152     T term;
0153     const T a1, a2, z;
0154   };
0155 
0156   // partial specialization for 2F1
0157   template <class T>
0158   struct hypergeometric_pFq_generic_series_term<T, 2u, 1u>
0159   {
0160     typedef T result_type;
0161 
0162     hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& b, const T& z)
0163        : n(0), term(1), a1(a1), a2(a2), b(b), z(z)
0164     {
0165     }
0166 
0167     T operator()()
0168     {
0169       BOOST_MATH_STD_USING
0170       const T r = term;
0171       term *= (((a1 + n) * (a2 + n) / ((b + n) * (n + 1))) * z);
0172       ++n;
0173       return r;
0174     }
0175 
0176   private:
0177     unsigned n;
0178     T term;
0179     const T a1, a2, b, z;
0180   };
0181 
0182   // we don't need to define extra check and make a polinom from
0183   // series, when p(i) and q(i) are negative integers and p(i) >= q(i)
0184   // as described in functions.wolfram.alpha, because we always
0185   // stop summation when result (in this case numerator) is zero.
0186   template <class T, unsigned p, unsigned q, class Policy>
0187   inline T sum_pFq_series(detail::hypergeometric_pFq_generic_series_term<T, p, q>& term, const Policy& pol)
0188   {
0189     BOOST_MATH_STD_USING
0190     std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0191 
0192     const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0193 
0194     policies::check_series_iterations<T>("boost::math::hypergeometric_pFq_generic_series<%1%>(%1%,%1%,%1%)", max_iter, pol);
0195     return result;
0196   }
0197 
0198   template <class T, class Policy>
0199   inline T hypergeometric_0F1_generic_series(const T& b, const T& z, const Policy& pol)
0200   {
0201     detail::hypergeometric_pFq_generic_series_term<T, 0u, 1u> s(b, z);
0202     return detail::sum_pFq_series(s, pol);
0203   }
0204 
0205   template <class T, class Policy>
0206   inline T hypergeometric_1F0_generic_series(const T& a, const T& z, const Policy& pol)
0207   {
0208     detail::hypergeometric_pFq_generic_series_term<T, 1u, 0u> s(a, z);
0209     return detail::sum_pFq_series(s, pol);
0210   }
0211 
0212   template <class T, class Policy>
0213   inline T log_pochhammer(T z, unsigned n, const Policy pol, int* s = nullptr)
0214   {
0215      BOOST_MATH_STD_USING
0216 #if 0
0217      if (z < 0)
0218      {
0219         if (n < -z)
0220         {
0221            if(s)
0222             *s = (n & 1 ? -1 : 1);
0223            return log_pochhammer(T(-z + (1 - (int)n)), n, pol);
0224         }
0225         else
0226         {
0227            int cross = itrunc(ceil(-z));
0228            return log_pochhammer(T(-z + (1 - cross)), cross, pol, s) + log_pochhammer(T(cross + z), n - cross, pol);
0229         }
0230      }
0231      else
0232 #endif
0233      {
0234         if (z + n < 0)
0235         {
0236            T r = log_pochhammer(T(-z - n + 1), n, pol, s);
0237            if (s)
0238               *s *= (n & 1 ? -1 : 1);
0239            return r;
0240         }
0241         int s1, s2;
0242         auto r = static_cast<T>(boost::math::lgamma(T(z + n), &s1, pol) - boost::math::lgamma(z, &s2, pol));
0243         if(s)
0244            *s = s1 * s2;
0245         return r;
0246      }
0247   }
0248 
0249   template <class T, class Policy>
0250   inline T hypergeometric_1F1_generic_series(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling, const char* function)
0251   {
0252      BOOST_MATH_STD_USING
0253      T sum(0), term(1), upper_limit(sqrt(boost::math::tools::max_value<T>())), diff;
0254      T lower_limit(1 / upper_limit);
0255      unsigned n = 0;
0256      long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<T>()) - 2;
0257      T scaling_factor = exp(T(log_scaling_factor));
0258      T term_m1 = 0;
0259      long long local_scaling = 0;
0260      //
0261      // When a is very small, then (a+n)/n => 1 faster than
0262      // z / (b+n) => 1, as a result the series starts off
0263      // converging, then at some unspecified time very gradually
0264      // starts to diverge, potentially resulting in some very large
0265      // values being missed.  As a result we need a check for small
0266      // a in the convergence criteria.  Note that this issue occurs
0267      // even when all the terms are positive.
0268      //
0269      bool small_a = fabs(a) < 0.25;
0270 
0271      unsigned summit_location = 0;
0272      bool have_minima = false;
0273      T sq = 4 * a * z + b * b - 2 * b * z + z * z;
0274      if (sq >= 0)
0275      {
0276         T t = (-sqrt(sq) - b + z) / 2;
0277         if (t > 1)  // Don't worry about a minima between 0 and 1.
0278            have_minima = true;
0279         t = (sqrt(sq) - b + z) / 2;
0280         if (t > 0)
0281            summit_location = itrunc(t);
0282      }
0283 
0284      if (summit_location > boost::math::policies::get_max_series_iterations<Policy>() / 4)
0285      {
0286         //
0287         // Skip forward to the location of the largest term in the series and
0288         // evaluate outwards from there:
0289         //
0290         int s1, s2;
0291         term = log_pochhammer(a, summit_location, pol, &s1) + summit_location * log(z) - log_pochhammer(b, summit_location, pol, &s2) - lgamma(T(summit_location + 1), pol);
0292         //std::cout << term << " " << log_pochhammer(boost::multiprecision::mpfr_float(a), summit_location, pol, &s1) + summit_location * log(boost::multiprecision::mpfr_float(z)) - log_pochhammer(boost::multiprecision::mpfr_float(b), summit_location, pol, &s2) - lgamma(boost::multiprecision::mpfr_float(summit_location + 1), pol) << std::endl;
0293         local_scaling = lltrunc(term);
0294         log_scaling += local_scaling;
0295         term = s1 * s2 * exp(term - local_scaling);
0296         //std::cout << term << " " << exp(log_pochhammer(boost::multiprecision::mpfr_float(a), summit_location, pol, &s1) + summit_location * log(boost::multiprecision::mpfr_float(z)) - log_pochhammer(boost::multiprecision::mpfr_float(b), summit_location, pol, &s2) - lgamma(boost::multiprecision::mpfr_float(summit_location + 1), pol) - local_scaling) << std::endl;
0297         n = summit_location;
0298      }
0299      else
0300         summit_location = 0;
0301 
0302      T saved_term = term;
0303      long long saved_scale = local_scaling;
0304 
0305      do
0306      {
0307         sum += term;
0308         //std::cout << n << " " << term * exp(boost::multiprecision::mpfr_float(local_scaling)) << " " << rising_factorial(boost::multiprecision::mpfr_float(a), n) * pow(boost::multiprecision::mpfr_float(z), n) / (rising_factorial(boost::multiprecision::mpfr_float(b), n) * factorial<boost::multiprecision::mpfr_float>(n)) << std::endl;
0309         if (fabs(sum) >= upper_limit)
0310         {
0311            sum /= scaling_factor;
0312            term /= scaling_factor;
0313            log_scaling += log_scaling_factor;
0314            local_scaling += log_scaling_factor;
0315         }
0316         if (fabs(sum) < lower_limit)
0317         {
0318            sum *= scaling_factor;
0319            term *= scaling_factor;
0320            log_scaling -= log_scaling_factor;
0321            local_scaling -= log_scaling_factor;
0322         }
0323         term_m1 = term;
0324         term *= (((a + n) / ((b + n) * (n + 1))) * z);
0325         if (n - summit_location > boost::math::policies::get_max_series_iterations<Policy>())
0326            return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
0327         ++n;
0328         diff = fabs(term / sum);
0329      } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)) || (small_a && n < 10));
0330 
0331      //
0332      // See if we need to go backwards as well:
0333      //
0334      if (summit_location)
0335      {
0336         //
0337         // Backup state:
0338         //
0339         term = saved_term * exp(T(local_scaling - saved_scale));
0340         n = summit_location;
0341         term *= (b + (n - 1)) * n / ((a + (n - 1)) * z);
0342         --n;
0343 
0344         do
0345         {
0346            sum += term;
0347            //std::cout << n << " " << term * exp(boost::multiprecision::mpfr_float(local_scaling)) << " " << rising_factorial(boost::multiprecision::mpfr_float(a), n) * pow(boost::multiprecision::mpfr_float(z), n) / (rising_factorial(boost::multiprecision::mpfr_float(b), n) * factorial<boost::multiprecision::mpfr_float>(n)) << std::endl;
0348            if (n == 0)
0349               break;
0350            if (fabs(sum) >= upper_limit)
0351            {
0352               sum /= scaling_factor;
0353               term /= scaling_factor;
0354               log_scaling += log_scaling_factor;
0355               local_scaling += log_scaling_factor;
0356            }
0357            if (fabs(sum) < lower_limit)
0358            {
0359               sum *= scaling_factor;
0360               term *= scaling_factor;
0361               log_scaling -= log_scaling_factor;
0362               local_scaling -= log_scaling_factor;
0363            }
0364            term_m1 = term;
0365            term *= (b + (n - 1)) * n / ((a + (n - 1)) * z);
0366            if (summit_location - n > boost::math::policies::get_max_series_iterations<Policy>())
0367               return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
0368            --n;
0369            diff = fabs(term / sum);
0370         } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)));
0371      }
0372 
0373      if (have_minima && n && summit_location)
0374      {
0375         //
0376         // There are a few terms starting at n == 0 which
0377         // haven't been accounted for yet...
0378         //
0379         unsigned backstop = n;
0380         n = 0;
0381         term = exp(T(-local_scaling));
0382         do
0383         {
0384            sum += term;
0385            //std::cout << n << " " << term << " " << sum << std::endl;
0386            if (fabs(sum) >= upper_limit)
0387            {
0388               sum /= scaling_factor;
0389               term /= scaling_factor;
0390               log_scaling += log_scaling_factor;
0391            }
0392            if (fabs(sum) < lower_limit)
0393            {
0394               sum *= scaling_factor;
0395               term *= scaling_factor;
0396               log_scaling -= log_scaling_factor;
0397            }
0398            //term_m1 = term;
0399            term *= (((a + n) / ((b + n) * (n + 1))) * z);
0400            if (n > boost::math::policies::get_max_series_iterations<Policy>())
0401               return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
0402            if (++n == backstop)
0403               break; // we've caught up with ourselves.
0404            diff = fabs(term / sum);
0405         } while ((diff > boost::math::policies::get_epsilon<T, Policy>())/* || (fabs(term_m1) < fabs(term))*/);
0406      }
0407      //std::cout << sum << std::endl;
0408      return sum;
0409   }
0410 
0411   template <class T, class Policy>
0412   inline T hypergeometric_1F2_generic_series(const T& a, const T& b1, const T& b2, const T& z, const Policy& pol)
0413   {
0414     detail::hypergeometric_pFq_generic_series_term<T, 1u, 2u> s(a, b1, b2, z);
0415     return detail::sum_pFq_series(s, pol);
0416   }
0417 
0418   template <class T, class Policy>
0419   inline T hypergeometric_2F0_generic_series(const T& a1, const T& a2, const T& z, const Policy& pol)
0420   {
0421     detail::hypergeometric_pFq_generic_series_term<T, 2u, 0u> s(a1, a2, z);
0422     return detail::sum_pFq_series(s, pol);
0423   }
0424 
0425   template <class T, class Policy>
0426   inline T hypergeometric_2F1_generic_series(const T& a1, const T& a2, const T& b, const T& z, const Policy& pol)
0427   {
0428     detail::hypergeometric_pFq_generic_series_term<T, 2u, 1u> s(a1, a2, b, z);
0429     return detail::sum_pFq_series(s, pol);
0430   }
0431 
0432   } } } // namespaces
0433 
0434 #endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!