1[section:recurrence Tools For 3-Term Recurrence Relations] 2 3[h4 Synopsis] 4 5`` 6#include <boost/math/tools/recurrence.hpp> 7`` 8 9 namespace boost{ namespace math{ namespace tools{ 10 11 template <class Recurrence, class T> 12 T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter); 13 14 template <class Recurrence, class T> 15 T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter); 16 17 template <class NextCoefs, class T> 18 T apply_recurrence_relation_forward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0); 19 20 template <class T, class NextCoefs> 21 T apply_recurrence_relation_backward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0); 22 23 template <class Recurrence> 24 struct forward_recurrence_iterator; 25 26 template <class Recurrence> 27 struct backward_recurrence_iterator; 28 29 }}} // namespaces 30 31[h4 Description] 32 33All of the tools in this header require a description of the recurrence relation: this takes the form of 34a functor that returns a tuple containing the 3 coefficients, specifically, given a recurrence relation: 35 36[/\Large $$ a_nF_{n-1} + b_nF_n + c_nF_{n+1} = 0 $$] 37[equation three_term_recurrence] 38 39And a functor `F` then the expression: 40 41[expression F(n);] 42 43Returns a tuple containing [role serif_italic { a[sub n], b[sub n], c[sub n] }]. 44 45For example, the recurrence relation for the Bessel J and Y functions when written in this form is: 46 47[/\Large $$ J_{v-1}(x) - \frac{2v}{x}J_v(x) + J_{v+1}(x)= 0 $$] 48[$../equations/three_term_recurrence_bessel_jy.svg] 49 50Therefore, given local variables /x/ and /v/ of type `double` the recurrence relation for Bessel J and Y can be encoded 51in a lambda expression like this: 52 53 auto recurrence_functor_jy = [&](int n) { return std::make_tuple(1.0, -2 * (v + n) / x, 1.0); }; 54 55Similarly, the Bessel I and K recurrence relation differs just by the sign of the final term: 56 57[/\Large $$ I_{v-1}(x) - \frac{2v}{x}I_v(x) - I_{v+1}(x)= 0 $$] 58[$../equations/three_term_recurrence_bessel_ik.svg] 59 60And this could be encoded as: 61 62 auto recurrence_functor_ik = [&](int n) { return std::make_tuple(1.0, -2 * (v + n) / x, -1.0); }; 63 64The tools are then as follows: 65 66 template <class Recurrence, class T> 67 T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter); 68 69Given a functor `r` which encodes the recurrence relation for function `F` at some location /n/, then returns the ratio: 70 71[/\Large $$ F_n / F_{n-1} $$] 72[$../equations/three_term_recurrence_backwards_ratio.svg] 73 74This calculation is stable only if recurrence is stable in the backwards direction. Further the ratio calculated 75is for the dominant solution (in the backwards direction) of the recurrence relation, if there are multiple solutions, 76then there is no guarantee that this will find the one you want or expect. 77 78Argument /factor/ is the tolerance required for convergence of the continued fraction associated with 79the recurrence relation, and should be no smaller than machine epsilon. Argument /max_iter/ sets 80the maximum number of permitted iterations in the associated continued fraction. 81 82 template <class Recurrence, class T> 83 T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter); 84 85Given a functor `r` which encodes the recurrence relation for function F at some location /n/, then returns the ratio: 86 87[/\Large $$ F_n / F_{n+1} $$] 88[$../equations/three_term_recurrence_forwards_ratio.svg] 89 90This calculation is stable only if recurrence is stable in the forwards direction. Further the ratio calculated 91is for the dominant solution (in the forwards direction) of the recurrence relation, if there are multiple solutions, 92then there is no guarantee that this will find the one you want or expect. 93 94Argument /factor/ is the tolerance required for convergence of the continued fraction associated with 95the recurrence relation, and should be no smaller than machine epsilon. Argument /max_iter/ sets 96the maximum number of permitted iterations in the associated continued fraction. 97 98 template <class NextCoefs, class T> 99 T apply_recurrence_relation_forward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0); 100 101Applies a recurrence relation in a stable forward direction, starting with the values F[sub n-1] and F[sub n]. 102 103[variablelist 104 [[get_coefs] [Functor that returns the coefficients of the recurrence relation. The coefficients should be centered on position /second/.]] 105 [[number_of_steps][The number of steps to apply the recurrence relation onwards from /second/.]] 106 [[first] [The value of F[sub n-1]]] 107 [[second] [The value of F[sub n]]] 108 [[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.]] 109 [[previous][When provided, is set to the value of F[sub n + number_of_steps - 1]] 110 ] 111] 112 113Returns F[sub n + number_of_steps]. 114 115 template <class NextCoefs, class T> 116 T apply_recurrence_relation_backward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0); 117 118Applies a recurrence relation in a stable backward direction, starting with the values F[sub n+1] and F[sub n]. 119 120[variablelist 121 [[get_coefs] [Functor that returns the coefficients of the recurrence relation. The coefficients should be centered on position /second/.]] 122 [[number_of_steps][The number of steps to apply the recurrence relation backwards from /second/.]] 123 [[first] [The value of F[sub n+1]]] 124 [[second] [The value of F[sub n]]] 125 [[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.]] 126 [[previous][When provided, is set to the value of F[sub n - number_of_steps + 1]] 127 ] 128] 129 130Returns F[sub n - number_of_steps]. 131 132 template <class Recurrence> 133 struct forward_recurrence_iterator 134 { 135 typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type; 136 137 forward_recurrence_iterator(const Recurrence& r, value_type f_n_minus_1, value_type f_n); 138 forward_recurrence_iterator(const Recurrence& r, value_type f_n); 139 /* Operators omitted for clarity */ 140 }; 141 142Type `forward_recurrence_iterator` defines a forward-iterator for a recurrence relation stable in the 143forward direction. The constructors take the recurrence relation, plus either one or two values: if 144only one value is provided, then the second is computed by using the recurrence relation to calculate the function ratio. 145 146 template <class Recurrence> 147 struct backward_recurrence_iterator 148 { 149 typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type; 150 151 backward_recurrence_iterator(const Recurrence& r, value_type f_n_plus_1, value_type f_n); 152 backward_recurrence_iterator(const Recurrence& r, value_type f_n); 153 /* Operators omitted for clarity */ 154 }; 155 156Type `backward_recurrence_iterator` defines a forward-iterator for a recurrence relation stable in the 157backward direction. The constructors take the recurrence relation, plus either one or two values: if 158only one value is provided, then the second is computed by using the recurrence relation to calculate the function ratio. 159 160Note that /incrementing/ this iterator moves the value returned successively to F[sub n-1], F[sub n-2] etc. 161 162[endsect] [/section:recurrence Tools For 3-Term Recurrence Relations] 163 164[/ 165 Copyright 2019 John Maddock. 166 Distributed under the Boost Software License, Version 1.0. 167 (See accompanying file LICENSE_1_0.txt or copy at 168 http://www.boost.org/LICENSE_1_0.txt). 169] 170 171
