libs/math/example/continued_fractions.cpp

   1 //  (C) Copyright John Maddock 2018.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include <boost/math/tools/fraction.hpp>
   7 #include <iostream>
   8 #include <complex>
   9 #include <boost/multiprecision/cpp_complex.hpp>
  10
  11 //[golden_ratio_1
  12 template <class T>
  13 struct golden_ratio_fraction
  14 {
  15    typedef T result_type;
  16
  17    result_type operator()()
  18    {
  19       return 1;
  20    }
  21 };
  22 //]
  23
  24 //[cf_tan_fraction
  25 template <class T>
  26 struct tan_fraction
  27 {
  28 private:
  29    T a, b;
  30 public:
  31    tan_fraction(T v)
  32       : a(-v * v), b(-1)
  33    {}
  34
  35    typedef std::pair<T, T> result_type;
  36
  37    std::pair<T, T> operator()()
  38    {
  39       b += 2;
  40       return std::make_pair(a, b);
  41    }
  42 };
  43 //]
  44 //[cf_tan
  45 template <class T>
  46 T tan(T a)
  47 {
  48    tan_fraction<T> fract(a);
  49    return a / continued_fraction_b(fract, std::numeric_limits<T>::epsilon());
  50 }
  51 //]
  52 //[cf_expint_fraction
  53 template <class T>
  54 struct expint_fraction
  55 {
  56    typedef std::pair<T, T> result_type;
  57    expint_fraction(unsigned n_, T z_) : b(z_ + T(n_)), i(-1), n(n_) {}
  58    std::pair<T, T> operator()()
  59    {
  60       std::pair<T, T> result = std::make_pair(-static_cast<T>((i + 1) * (n + i)), b);
  61       b += 2;
  62       ++i;
  63       return result;
  64    }
  65 private:
  66    T b;
  67    int i;
  68    unsigned n;
  69 };
  70 //]
  71 //[cf_expint
  72 template <class T>
  73 inline std::complex<T> expint_as_fraction(unsigned n, std::complex<T> const& z)
  74 {
  75    boost::uintmax_t max_iter = 1000;
  76    expint_fraction<std::complex<T> > f(n, z);
  77    std::complex<T> result = boost::math::tools::continued_fraction_b(
  78       f,
  79       std::complex<T>(std::numeric_limits<T>::epsilon()),
  80       max_iter);
  81    result = exp(-z) / result;
  82    return result;
  83 }
  84 //]
  85 //[cf_upper_gamma_fraction
  86 template <class T>
  87 struct upper_incomplete_gamma_fract
  88 {
  89 private:
  90    typedef typename T::value_type scalar_type;
  91    T z, a;
  92    int k;
  93 public:
  94    typedef std::pair<T, T> result_type;
  95
  96    upper_incomplete_gamma_fract(T a1, T z1)
  97       : z(z1 - a1 + scalar_type(1)), a(a1), k(0)
  98    {
  99    }
 100
 101    result_type operator()()
 102    {
 103       ++k;
 104       z += scalar_type(2);
 105       return result_type(scalar_type(k) * (a - scalar_type(k)), z);
 106    }
 107 };
 108 //]
 109 //[cf_gamma_Q
 110 template <class T>
 111 inline std::complex<T> gamma_Q_as_fraction(const std::complex<T>& a, const std::complex<T>& z)
 112 {
 113    upper_incomplete_gamma_fract<std::complex<T> > f(a, z);
 114    std::complex<T> eps(std::numeric_limits<T>::epsilon());
 115    return pow(z, a) / (exp(z) *(z - a + T(1) + boost::math::tools::continued_fraction_a(f, eps)));
 116 }
 117 //]
 118 inline boost::multiprecision::cpp_complex_50 gamma_Q_as_fraction(const boost::multiprecision::cpp_complex_50& a, const boost::multiprecision::cpp_complex_50& z)
 119 {
 120    upper_incomplete_gamma_fract<boost::multiprecision::cpp_complex_50> f(a, z);
 121    boost::multiprecision::cpp_complex_50 eps(std::numeric_limits<boost::multiprecision::cpp_complex_50::value_type>::epsilon());
 122    return pow(z, a) / (exp(z) * (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps)));
 123 }
 124
 125
 126 int main()
 127 {
 128    using namespace boost::math::tools;
 129
 130    //[cf_gr
 131    golden_ratio_fraction<double> func;
 132    double gr = continued_fraction_a(
 133       func,
 134       std::numeric_limits<double>::epsilon());
 135    std::cout << "The golden ratio is: " << gr << std::endl;
 136    //]
 137
 138    std::cout << tan(0.5) << std::endl;
 139
 140    std::complex<double> arg(3, 2);
 141    std::cout << expint_as_fraction(5, arg) << std::endl;
 142
 143    std::complex<double> a(3, 3), z(3, 2);
 144    std::cout << gamma_Q_as_fraction(a, z) << std::endl;
 145
 146    boost::multiprecision::cpp_complex_50 am(3, 3), zm(3, 2);
 147    std::cout << gamma_Q_as_fraction(am, zm) << std::endl;
 148
 149    return 0;
 150 }