...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

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

namespace boost{ namespace math{ namespace tools{ template <class Recurrence, class T> T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter); template <class Recurrence, class T> T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter); template <class NextCoefs, class T> T apply_recurrence_relation_forward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0); template <class T, class NextCoefs> T apply_recurrence_relation_backward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0); template <class Recurrence> struct forward_recurrence_iterator; template <class Recurrence> struct backward_recurrence_iterator; }}} // namespaces

All of the tools in this header require a description of the recurrence relation: this takes the form of a functor that returns a tuple containing the 3 coefficients, specifically, given a recurrence relation:

And a functor `F`

then the
expression:

F(n);

Returns a tuple containing { a_{n}, b_{n}, c_{n} }.

For example, the recurrence relation for the Bessel J and Y functions when written in this form is:

Therefore, given local variables *x* and *v*
of type `double`

the recurrence
relation for Bessel J and Y can be encoded in a lambda expression like this:

auto recurrence_functor_jy = [&](int n) { return std::make_tuple(1.0, -2 * (v + n) / x, 1.0); };

Similarly, the Bessel I and K recurrence relation differs just by the sign of the final term:

And this could be encoded as:

auto recurrence_functor_ik = [&](int n) { return std::make_tuple(1.0, -2 * (v + n) / x, -1.0); };

The tools are then as follows:

template <class Recurrence, class T> T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter);

Given a functor `r`

which encodes
the recurrence relation for function `F`

at some location *n*, then returns the ratio:

This calculation is stable only if recurrence is stable in the backwards direction. Further the ratio calculated is for the dominant solution (in the backwards direction) of the recurrence relation, if there are multiple solutions, then there is no guarantee that this will find the one you want or expect.

Argument *factor* is the tolerance required for convergence
of the continued fraction associated with the recurrence relation, and should
be no smaller than machine epsilon. Argument *max_iter*
sets the maximum number of permitted iterations in the associated continued
fraction.

template <class Recurrence, class T> T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter);

Given a functor `r`

which encodes
the recurrence relation for function F at some location *n*,
then returns the ratio:

This calculation is stable only if recurrence is stable in the forwards direction. Further the ratio calculated is for the dominant solution (in the forwards direction) of the recurrence relation, if there are multiple solutions, then there is no guarantee that this will find the one you want or expect.

Argument *factor* is the tolerance required for convergence
of the continued fraction associated with the recurrence relation, and should
be no smaller than machine epsilon. Argument *max_iter*
sets the maximum number of permitted iterations in the associated continued
fraction.

template <class NextCoefs, class T> T apply_recurrence_relation_forward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0);

Applies a recurrence relation in a stable forward direction, starting with
the values F_{n-1} and F_{n}.

- get_coefs
Functor that returns the coefficients of the recurrence relation. The coefficients should be centered on position

*second*.- number_of_steps
The number of steps to apply the recurrence relation onwards from

*second*.- first
The value of F

_{n-1}- second
The value of F

_{n}- log_scaling
When provided, the recurrence relations may be rescaled internally to avoid over/underflow issues. The result should be multiplied by

`exp(*log_scaling)`

to get the true value of the result.- previous
When provided, is set to the value of F

_{n + number_of_steps - 1}

Returns F_{n + number_of_steps}.

template <class NextCoefs, class T> T apply_recurrence_relation_backward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0);

Applies a recurrence relation in a stable backward direction, starting with
the values F_{n+1} and F_{n}.

- get_coefs
Functor that returns the coefficients of the recurrence relation. The coefficients should be centered on position

*second*.- number_of_steps
The number of steps to apply the recurrence relation backwards from

*second*.- first
The value of F

_{n+1}- second
The value of F

_{n}- log_scaling
When provided, the recurrence relations may be rescaled internally to avoid over/underflow issues. The result should be multiplied by

`exp(*log_scaling)`

to get the true value of the result.- previous
When provided, is set to the value of F

_{n - number_of_steps + 1}

Returns F_{n - number_of_steps}.

template <class Recurrence> struct forward_recurrence_iterator { typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type; forward_recurrence_iterator(const Recurrence& r, value_type f_n_minus_1, value_type f_n); forward_recurrence_iterator(const Recurrence& r, value_type f_n); /* Operators omitted for clarity */ };

Type `forward_recurrence_iterator`

defines a forward-iterator for a recurrence relation stable in the forward
direction. The constructors take the recurrence relation, plus either one
or two values: if only one value is provided, then the second is computed
by using the recurrence relation to calculate the function ratio.

template <class Recurrence> struct backward_recurrence_iterator { typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type; backward_recurrence_iterator(const Recurrence& r, value_type f_n_plus_1, value_type f_n); backward_recurrence_iterator(const Recurrence& r, value_type f_n); /* Operators omitted for clarity */ };

Type `backward_recurrence_iterator`

defines a forward-iterator for a recurrence relation stable in the backward
direction. The constructors take the recurrence relation, plus either one
or two values: if only one value is provided, then the second is computed
by using the recurrence relation to calculate the function ratio.

Note that *incrementing* this iterator moves the value
returned successively to F_{n-1}, F_{n-2} etc.