boost/math/special_functions/detail/hypergeometric_series.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2014 Anton Bikineev
   3 //  Copyright 2014 Christopher Kormanyos
   4 //  Copyright 2014 John Maddock
   5 //  Copyright 2014 Paul Bristow
   6 //  Distributed under the Boost
   7 //  Software License, Version 1.0. (See accompanying file
   8 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   9 //
  10 #ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP
  11 #define BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP
  12
  13 #include <boost/math/tools/series.hpp>
  14 #include <boost/math/special_functions/trunc.hpp>
  15 #include <boost/math/policies/error_handling.hpp>
  16
  17   namespace boost { namespace math { namespace detail {
  18
  19   // primary template for term of Taylor series
  20   template <class T, unsigned p, unsigned q>
  21   struct hypergeometric_pFq_generic_series_term;
  22
  23   // partial specialization for 0F1
  24   template <class T>
  25   struct hypergeometric_pFq_generic_series_term<T, 0u, 1u>
  26   {
  27     typedef T result_type;
  28
  29     hypergeometric_pFq_generic_series_term(const T& b, const T& z)
  30        : n(0), term(1), b(b), z(z)
  31     {
  32     }
  33
  34     T operator()()
  35     {
  36       BOOST_MATH_STD_USING
  37       const T r = term;
  38       term *= ((1 / ((b + n) * (n + 1))) * z);
  39       ++n;
  40       return r;
  41     }
  42
  43   private:
  44     unsigned n;
  45     T term;
  46     const T b, z;
  47   };
  48
  49   // partial specialization for 1F0
  50   template <class T>
  51   struct hypergeometric_pFq_generic_series_term<T, 1u, 0u>
  52   {
  53     typedef T result_type;
  54
  55     hypergeometric_pFq_generic_series_term(const T& a, const T& z)
  56        : n(0), term(1), a(a), z(z)
  57     {
  58     }
  59
  60     T operator()()
  61     {
  62       BOOST_MATH_STD_USING
  63       const T r = term;
  64       term *= (((a + n) / (n + 1)) * z);
  65       ++n;
  66       return r;
  67     }
  68
  69   private:
  70     unsigned n;
  71     T term;
  72     const T a, z;
  73   };
  74
  75   // partial specialization for 1F1
  76   template <class T>
  77   struct hypergeometric_pFq_generic_series_term<T, 1u, 1u>
  78   {
  79     typedef T result_type;
  80
  81     hypergeometric_pFq_generic_series_term(const T& a, const T& b, const T& z)
  82        : n(0), term(1), a(a), b(b), z(z)
  83     {
  84     }
  85
  86     T operator()()
  87     {
  88       BOOST_MATH_STD_USING
  89       const T r = term;
  90       term *= (((a + n) / ((b + n) * (n + 1))) * z);
  91       ++n;
  92       return r;
  93     }
  94
  95   private:
  96     unsigned n;
  97     T term;
  98     const T a, b, z;
  99   };
 100
 101   // partial specialization for 1F2
 102   template <class T>
 103   struct hypergeometric_pFq_generic_series_term<T, 1u, 2u>
 104   {
 105     typedef T result_type;
 106
 107     hypergeometric_pFq_generic_series_term(const T& a, const T& b1, const T& b2, const T& z)
 108        : n(0), term(1), a(a), b1(b1), b2(b2), z(z)
 109     {
 110     }
 111
 112     T operator()()
 113     {
 114       BOOST_MATH_STD_USING
 115       const T r = term;
 116       term *= (((a + n) / ((b1 + n) * (b2 + n) * (n + 1))) * z);
 117       ++n;
 118       return r;
 119     }
 120
 121   private:
 122     unsigned n;
 123     T term;
 124     const T a, b1, b2, z;
 125   };
 126
 127   // partial specialization for 2F0
 128   template <class T>
 129   struct hypergeometric_pFq_generic_series_term<T, 2u, 0u>
 130   {
 131     typedef T result_type;
 132
 133     hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& z)
 134        : n(0), term(1), a1(a1), a2(a2), z(z)
 135     {
 136     }
 137
 138     T operator()()
 139     {
 140       BOOST_MATH_STD_USING
 141       const T r = term;
 142       term *= (((a1 + n) * (a2 + n) / (n + 1)) * z);
 143       ++n;
 144       return r;
 145     }
 146
 147   private:
 148     unsigned n;
 149     T term;
 150     const T a1, a2, z;
 151   };
 152
 153   // partial specialization for 2F1
 154   template <class T>
 155   struct hypergeometric_pFq_generic_series_term<T, 2u, 1u>
 156   {
 157     typedef T result_type;
 158
 159     hypergeometric_pFq_generic_series_term(const T& a1, const T& a2, const T& b, const T& z)
 160        : n(0), term(1), a1(a1), a2(a2), b(b), z(z)
 161     {
 162     }
 163
 164     T operator()()
 165     {
 166       BOOST_MATH_STD_USING
 167       const T r = term;
 168       term *= (((a1 + n) * (a2 + n) / ((b + n) * (n + 1))) * z);
 169       ++n;
 170       return r;
 171     }
 172
 173   private:
 174     unsigned n;
 175     T term;
 176     const T a1, a2, b, z;
 177   };
 178
 179   // we don't need to define extra check and make a polinom from
 180   // series, when p(i) and q(i) are negative integers and p(i) >= q(i)
 181   // as described in functions.wolfram.alpha, because we always
 182   // stop summation when result (in this case numerator) is zero.
 183   template <class T, unsigned p, unsigned q, class Policy>
 184   inline T sum_pFq_series(detail::hypergeometric_pFq_generic_series_term<T, p, q>& term, const Policy& pol)
 185   {
 186     BOOST_MATH_STD_USING
 187     boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 188 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 189     const T zero = 0;
 190     const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
 191 #else
 192     const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 193 #endif
 194     policies::check_series_iterations<T>("boost::math::hypergeometric_pFq_generic_series<%1%>(%1%,%1%,%1%)", max_iter, pol);
 195     return result;
 196   }
 197
 198   template <class T, class Policy>
 199   inline T hypergeometric_0F1_generic_series(const T& b, const T& z, const Policy& pol)
 200   {
 201     detail::hypergeometric_pFq_generic_series_term<T, 0u, 1u> s(b, z);
 202     return detail::sum_pFq_series(s, pol);
 203   }
 204
 205   template <class T, class Policy>
 206   inline T hypergeometric_1F0_generic_series(const T& a, const T& z, const Policy& pol)
 207   {
 208     detail::hypergeometric_pFq_generic_series_term<T, 1u, 0u> s(a, z);
 209     return detail::sum_pFq_series(s, pol);
 210   }
 211
 212   template <class T, class Policy>
 213   inline T log_pochhammer(T z, unsigned n, const Policy pol, int* s = 0)
 214   {
 215      BOOST_MATH_STD_USING
 216 #if 0
 217      if (z < 0)
 218      {
 219         if (n < -z)
 220         {
 221            if(s)
 222             *s = (n & 1 ? -1 : 1);
 223            return log_pochhammer(T(-z + (1 - (int)n)), n, pol);
 224         }
 225         else
 226         {
 227            int cross = itrunc(ceil(-z));
 228            return log_pochhammer(T(-z + (1 - cross)), cross, pol, s) + log_pochhammer(T(cross + z), n - cross, pol);
 229         }
 230      }
 231      else
 232 #endif
 233      {
 234         if (z + n < 0)
 235         {
 236            T r = log_pochhammer(T(-z - n + 1), n, pol, s);
 237            if (s)
 238               *s *= (n & 1 ? -1 : 1);
 239            return r;
 240         }
 241         int s1, s2;
 242         T r = boost::math::lgamma(T(z + n), &s1, pol) - boost::math::lgamma(z, &s2, pol);
 243         if(s)
 244            *s = s1 * s2;
 245         return r;
 246      }
 247   }
 248
 249   template <class T, class Policy>
 250   inline T hypergeometric_1F1_generic_series(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling, const char* function)
 251   {
 252      BOOST_MATH_STD_USING
 253      T sum(0), term(1), upper_limit(sqrt(boost::math::tools::max_value<T>())), diff;
 254      T lower_limit(1 / upper_limit);
 255      unsigned n = 0;
 256      int log_scaling_factor = itrunc(boost::math::tools::log_max_value<T>()) - 2;
 257      T scaling_factor = exp(T(log_scaling_factor));
 258      T term_m1 = 0;
 259      int local_scaling = 0;
 260      //
 261      // When a is very small, then (a+n)/n => 1 faster than
 262      // z / (b+n) => 1, as a result the series starts off
 263      // converging, then at some unspecified time very gradually
 264      // starts to diverge, potentially resulting in some very large
 265      // values being missed.  As a result we need a check for small
 266      // a in the convergence critera.  Note that this issue occurs
 267      // even when all the terms are positive.
 268      //
 269      bool small_a = fabs(a) < 0.25;
 270
 271      unsigned summit_location = 0;
 272      bool have_minima = false;
 273      T sq = 4 * a * z + b * b - 2 * b * z + z * z;
 274      if (sq >= 0)
 275      {
 276         T t = (-sqrt(sq) - b + z) / 2;
 277         if (t > 1)  // Don't worry about a minima between 0 and 1.
 278            have_minima = true;
 279         t = (sqrt(sq) - b + z) / 2;
 280         if (t > 0)
 281            summit_location = itrunc(t);
 282      }
 283
 284      if (summit_location > boost::math::policies::get_max_series_iterations<Policy>() / 4)
 285      {
 286         //
 287         // Skip forward to the location of the largest term in the series and
 288         // evaluate outwards from there:
 289         //
 290         int s1, s2;
 291         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);
 292         //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;
 293         local_scaling = itrunc(term);
 294         log_scaling += local_scaling;
 295         term = s1 * s2 * exp(term - local_scaling);
 296         //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;
 297         n = summit_location;
 298      }
 299      else
 300         summit_location = 0;
 301
 302      T saved_term = term;
 303      int saved_scale = local_scaling;
 304
 305      do
 306      {
 307         sum += term;
 308         //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;
 309         if (fabs(sum) >= upper_limit)
 310         {
 311            sum /= scaling_factor;
 312            term /= scaling_factor;
 313            log_scaling += log_scaling_factor;
 314            local_scaling += log_scaling_factor;
 315         }
 316         if (fabs(sum) < lower_limit)
 317         {
 318            sum *= scaling_factor;
 319            term *= scaling_factor;
 320            log_scaling -= log_scaling_factor;
 321            local_scaling -= log_scaling_factor;
 322         }
 323         term_m1 = term;
 324         term *= (((a + n) / ((b + n) * (n + 1))) * z);
 325         if (n - summit_location > boost::math::policies::get_max_series_iterations<Policy>())
 326            return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
 327         ++n;
 328         diff = fabs(term / sum);
 329      } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)) || (small_a && n < 10));
 330
 331      //
 332      // See if we need to go backwards as well:
 333      //
 334      if (summit_location)
 335      {
 336         //
 337         // Backup state:
 338         //
 339         term = saved_term * exp(T(local_scaling - saved_scale));
 340         n = summit_location;
 341         term *= (b + (n - 1)) * n / ((a + (n - 1)) * z);
 342         --n;
 343
 344         do
 345         {
 346            sum += term;
 347            //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;
 348            if (n == 0)
 349               break;
 350            if (fabs(sum) >= upper_limit)
 351            {
 352               sum /= scaling_factor;
 353               term /= scaling_factor;
 354               log_scaling += log_scaling_factor;
 355               local_scaling += log_scaling_factor;
 356            }
 357            if (fabs(sum) < lower_limit)
 358            {
 359               sum *= scaling_factor;
 360               term *= scaling_factor;
 361               log_scaling -= log_scaling_factor;
 362               local_scaling -= log_scaling_factor;
 363            }
 364            term_m1 = term;
 365            term *= (b + (n - 1)) * n / ((a + (n - 1)) * z);
 366            if (summit_location - n > boost::math::policies::get_max_series_iterations<Policy>())
 367               return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
 368            --n;
 369            diff = fabs(term / sum);
 370         } while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)));
 371      }
 372
 373      if (have_minima && n && summit_location)
 374      {
 375         //
 376         // There are a few terms starting at n == 0 which
 377         // haven't been accounted for yet...
 378         //
 379         unsigned backstop = n;
 380         n = 0;
 381         term = exp(T(-local_scaling));
 382         do
 383         {
 384            sum += term;
 385            //std::cout << n << " " << term << " " << sum << std::endl;
 386            if (fabs(sum) >= upper_limit)
 387            {
 388               sum /= scaling_factor;
 389               term /= scaling_factor;
 390               log_scaling += log_scaling_factor;
 391            }
 392            if (fabs(sum) < lower_limit)
 393            {
 394               sum *= scaling_factor;
 395               term *= scaling_factor;
 396               log_scaling -= log_scaling_factor;
 397            }
 398            //term_m1 = term;
 399            term *= (((a + n) / ((b + n) * (n + 1))) * z);
 400            if (n > boost::math::policies::get_max_series_iterations<Policy>())
 401               return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
 402            if (++n == backstop)
 403               break; // we've caught up with ourselves.
 404            diff = fabs(term / sum);
 405         } while ((diff > boost::math::policies::get_epsilon<T, Policy>())/* || (fabs(term_m1) < fabs(term))*/);
 406      }
 407      //std::cout << sum << std::endl;
 408      return sum;
 409   }
 410
 411   template <class T, class Policy>
 412   inline T hypergeometric_1F2_generic_series(const T& a, const T& b1, const T& b2, const T& z, const Policy& pol)
 413   {
 414     detail::hypergeometric_pFq_generic_series_term<T, 1u, 2u> s(a, b1, b2, z);
 415     return detail::sum_pFq_series(s, pol);
 416   }
 417
 418   template <class T, class Policy>
 419   inline T hypergeometric_2F0_generic_series(const T& a1, const T& a2, const T& z, const Policy& pol)
 420   {
 421     detail::hypergeometric_pFq_generic_series_term<T, 2u, 0u> s(a1, a2, z);
 422     return detail::sum_pFq_series(s, pol);
 423   }
 424
 425   template <class T, class Policy>
 426   inline T hypergeometric_2F1_generic_series(const T& a1, const T& a2, const T& b, const T& z, const Policy& pol)
 427   {
 428     detail::hypergeometric_pFq_generic_series_term<T, 2u, 1u> s(a1, a2, b, z);
 429     return detail::sum_pFq_series(s, pol);
 430   }
 431
 432   } } } // namespaces
 433
 434 #endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_SERIES_HPP