Synopsis

#include <boost/math/tools/fraction.hpp>

namespace boost{ namespace math{ namespace tools{

template <class Gen, class U>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_b(Gen& g, const U& tolerance, boost::uintmax_t& max_terms)

template <class Gen, class U>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_b(Gen& g, const U& tolerance)

template <class Gen, class U>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_a(Gen& g, const U& tolerance, boost::uintmax_t& max_terms)

template <class Gen, class U>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_a(Gen& g, const U& tolerance)

//
// These interfaces are present for legacy reasons, and are now deprecated:
//
template <class Gen>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_b(Gen& g, int bits);

template <class Gen>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms);

template <class Gen>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_a(Gen& g, int bits);

template <class Gen>
typename detail::fraction_traits<Gen>::result_type
   continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms);

}}} // namespaces

Description

Continued fractions are a common method of approximation. These functions all evaluate the continued fraction described by the generator type argument. The functions with an "_a" suffix evaluate the fraction:

and those with a "_b" suffix evaluate the fraction:

This latter form is somewhat more natural in that it corresponds with the usual definition of a continued fraction, but note that the first a value returned by the generator is discarded. Further, often the first a and b values in a continued fraction have different defining equations to the remaining terms, which may make the "_a" suffixed form more appropriate.

The generator type should be a function object which supports the following operations:

In all the continued fraction evaluation functions the tolerance parameter is the precision desired in the result, evaluation of the fraction will continue until the last term evaluated leaves the relative error in the result less than tolerance. The deprecated interfaces take a number of digits precision here, internally they just convert this to a tolerance and forward call.

If the optional max_terms parameter is specified then no more than max_terms calls to the generator will be made, and on output, max_terms will be set to actual number of calls made. This facility is particularly useful when profiling a continued fraction for convergence.

Implementation

Internally these algorithms all use the modified Lentz algorithm: refer to Numeric Recipes in C++, W. H. Press et all, chapter 5, (especially 5.2 Evaluation of continued fractions, p 175 - 179) for more information, also Lentz, W.J. 1976, Applied Optics, vol. 15, pp. 668-671.

Examples

All of these examples are in continued_fractions.cpp.

The golden ratio phi = 1.618033989... can be computed from the simplest continued fraction of all:

We begin by defining a generator function:

template <class T>
struct golden_ratio_fraction
{
   typedef T result_type;

   result_type operator()()
   {
      return 1;
   }
};

The golden ratio can then be computed to double precision using:

golden_ratio_fraction<double> func;
double gr = continued_fraction_a(
   func,
   std::numeric_limits<double>::epsilon());
std::cout << "The golden ratio is: " << gr << std::endl;

It's more usual though to have to define both the a's and the b's when evaluating special functions by continued fractions, for example the tan function is defined by:

So its generator object would look like:

template <class T>
struct tan_fraction
{
private:
   T a, b;
public:
   tan_fraction(T v)
      : a(-v * v), b(-1)
   {}

   typedef std::pair<T, T> result_type;

   std::pair<T, T> operator()()
   {
      b += 2;
      return std::make_pair(a, b);
   }
};

Notice that if the continuant is subtracted from the b terms, as is the case here, then all the a terms returned by the generator will be negative. The tangent function can now be evaluated using:

template <class T>
T tan(T a)
{
   tan_fraction<T> fract(a);
   return a / continued_fraction_b(fract, std::numeric_limits<T>::epsilon());
}

Notice that this time we're using the "_b" suffixed version to evaluate the fraction: we're removing the leading a term during fraction evaluation as it's different from all the others.

Now we'll look at a couple of complex number examples, starting with the exponential integral which can be calculated via:

So our functor looks like this:

template <class T>
struct expint_fraction
{
   typedef std::pair<T, T> result_type;
   expint_fraction(unsigned n_, T z_) : b(z_ + T(n_)), i(-1), n(n_) {}
   std::pair<T, T> operator()()
   {
      std::pair<T, T> result = std::make_pair(-static_cast<T>((i + 1) * (n + i)), b);
      b += 2;
      ++i;
      return result;
   }
private:
   T b;
   int i;
   unsigned n;
};

We can finish the example by wrapping everything up in a function:

template <class T>
inline std::complex<T> expint_as_fraction(unsigned n, std::complex<T> const& z)
{
   boost::uintmax_t max_iter = 1000;
   expint_fraction<std::complex<T> > f(n, z);
   std::complex<T> result = boost::math::tools::continued_fraction_b(
      f,
      std::complex<T>(std::numeric_limits<T>::epsilon()),
      max_iter);
   result = exp(-z) / result;
   return result;
}

Notice how the termination condition is still expressed as a complex number, albeit one with zero imaginary part.

Our final example will use continued_fraction_a, in fact there is only one special function in our code which uses that variant, and it's the upper incomplete gamma function (Q), which can be calculated via:

In this case the first couple of terms are different from the rest, so our fraction will start with the first "regular" a term:

template <class T>
struct upper_incomplete_gamma_fract
{
private:
   typedef typename T::value_type scalar_type;
   T z, a;
   int k;
public:
   typedef std::pair<T, T> result_type;

   upper_incomplete_gamma_fract(T a1, T z1)
      : z(z1 - a1 + scalar_type(1)), a(a1), k(0)
   {
   }

   result_type operator()()
   {
      ++k;
      z += scalar_type(2);
      return result_type(scalar_type(k) * (a - scalar_type(k)), z);
   }
};

So now we can implement Q, this time using continued_fraction_a:

template <class T>
inline std::complex<T> gamma_Q_as_fraction(const std::complex<T>& a, const std::complex<T>& z)
{
   upper_incomplete_gamma_fract<std::complex<T> > f(a, z);
   std::complex<T> eps(std::numeric_limits<T>::epsilon());
   return pow(z, a) / (exp(z) *(z - a + T(1) + boost::math::tools::continued_fraction_a(f, eps)));
}