python - Stably computing large quantites through recursion -

i have 2 quantities & b defined recursion , through reference list of values x = [ x_1, x_2, ... x_n ], input program. program iterate on values in x , update & b according to:

for n in range(1,n)     a[n] = a[n-1] * exp(+x[n]) + b[n-1] * exp(-x[n])       b[n] = b[n-1] * exp(+x[n]) + a[n-1] * exp(-x[n])   

and starting values

a[0] = exp(+x[0]) b[0] = exp(-x[0]) 

the values in x not big numbers (always <10) n in hundreds, , because of exponentials final values of & b large. i'm worried because of form of recursion multiplying exponentially large numbers exponentially small ones , adding them scheme become quite numerically unstable.

ideally calculate log(a) , log(b) instead stop values becoming large. because of recursion scheme that's not possible, unless compute messier like

log_a[n] = x[n] + log_a[n-1] + log( 1 + exp(-2*x[n] + log_b[n-1]-log_a[n-1] ) ) 

is numerical stability right concerned here? , if log based scheme stabilise it?

we can rewrite first as:

for n in range(1,n)     a[n] = exp(log(a[n-1]) + x[n]) + exp(log(b[n-1]) - x[n])     b[n] = exp(log(b[n-1]) + x[n]) + exp(log(a[n-1]) - x[n])) 

then change our iteration variables:

for n in range(1,n)     log_a[n] = log(exp(log_a[n-1] + x[n]) + exp(log_b[n-1] - x[n]))     log_b[n] = log(exp(log_b[n-1] + x[n]) + exp(log_a[n-1] - x[n])) 

which can computed more stably using np.logaddexp:

for n in range(1,n)     log_a[n] = np.logaddexp(log_a[n-1] + x[n], log_b[n-1] - x[n])     log_b[n] = np.logaddexp(log_b[n-1] + x[n], log_a[n-1] - x[n]) 

the implementation of logaddexp can seen here