""" Infinite continued fractions via Lentz's method.
TODO More documentation.
"""
from __future__ import division, print_function, absolute_import
import numpy as np
# TODO some documentation here, better documentation throughout
# Note that if we jitted `lentz`, it would only be able to run on jitted functions.
[docs]def lentz(a, b, tol=1.e-10, N_min=0, N_max=np.Inf, tiny=1.e-30, args=None):
"""Compute a continued fraction via modified Lentz's method.
This implementation is by the book [1]_. The value to compute is:
b_0 + a_1/( b_1 + a_2/( b_2 + a_3/( b_3 + ...)))
where a_n = a(n, *args) and b_n = b(n, *args).
Parameters
----------
a: callable returning numeric.
b: callable returning numeric.
tol: float [default: 1.e-10]
Tolerance for termination of evaluation.
N_min: int [default: 0]
Minimum number of iterations to evaluate.
N_max: int or comparable [default: np.Inf]
Maximum number of iterations to evaluate.
tiny: float [default: 1.e-30]
Very small number to control convergence of Lentz's method when
there is cancellation in a denominator.
args: tuple [default: None]
Additional arguments to pass to the user-defined functions a, b.
If given, the additional arguments are passed to all
user-defined functions. So if, for example, `a` has the
signature `a(n, x, y)`, then `b` must have the same signature,
and args must be a tuple of length 2, `args=(x,y)`.
Returns
-------
(float, float, int)
The first element of the tuple is the value of the continued
fraction. The second element is the estimated error. The third
element is the number of iterations.
References
----------
.. [1] WH Press, SA Teukolsky, WT Vetterling, BP Flannery,
"Numerical Recipes," 3rd Ed., Cambridge University Press 2007,
ISBN 0521880688, 9780521880688 .
Examples
--------
Compute the square root of two using continued fractions:
>>> from qnm.contfrac import lentz
>>> def rt2b(n):
... if (n==0):
... return 1
... return 2
...
>>> def rt2a(n): return 1
>>> lentz(rt2a, rt2b)
(1.4142135623638004, 4.488287519421874e-11, 14)
Compute phi:
>>> phia = rt2a
>>> phib = rt2a
>>> lentz(phia, phib)
(1.6180339887802424, 6.785971784495359e-11, 25)
Compute pi:
>>> def pia(n):
... if (n==1):
... return 4.
... return (n-1.)*(n-1.)
...
>>> def pib(n):
... if (n==0):
... return 0.
... return 2*n-1.
...
>>> lentz(pia, pib, tol=1.e-15)
(3.1415926535897922, 8.881784197001252e-16, 21)
Compute e:
>>> def e_a(n):
... if (n==1):
... return 1.
... return (n-1.)
...
>>> def e_b(n):
... if (n==0):
... return 2.
... return n
...
>>> lentz(e_a, e_b, tol=1.e-15)
(2.7182818284590464, 3.3306690738754696e-16, 16)
Compute cotan(1):
>>> def cot1_a(n):
... return -1.
...
>>> def cot1_b(n):
... return 2.*n+1.
...
>>> lentz(cot1_a, cot1_b, tol=1.e-15)
(0.6420926159343306, 1.1102230246251565e-16, 9)
Compute tan(x):
>>> def tanx_a(n, x):
... if (n==1):
... return x
... else:
... return -x*x
...
>>> def tanx_b(n, x):
... if (n==0):
... return 0.
... else:
... return 2*n-1
...
>>> lentz(tanx_a, tanx_b, args=(1.,), tol=1.e-15)
(1.5574077246549018, 2.220446049250313e-16, 10)
"""
if args is None:
args = ()
if type(args) is not tuple:
raise ValueError("args={} must be of type tuple".format(args))
f_old = b(0, *args)
if (f_old == 0):
f_old = tiny
C_old = f_old
D_old = 0.
conv = False
j = 1
while ((not conv) and (j < N_max)):
aj, bj = a(j, *args), b(j, *args)
D_new = bj + aj * D_old
if (D_new == 0):
D_new = tiny
C_new = bj + aj / C_old
if (C_new == 0):
C_new = tiny
D_new = 1./D_new
Delta = C_new * D_new
f_new = f_old * Delta
if ((j > N_min) and (np.abs(Delta - 1.) < tol)): # converged
conv = True
# Set up for next iter
j = j + 1
D_old = D_new
C_old = C_new
f_old = f_new
# Success or failure can be assessed by the user
return f_new, np.abs(Delta - 1.), j-1
[docs]def lentz_gen(a, b, tol=1.e-10, N_min=0, N_max=np.Inf, tiny=1.e-30):
""" Compute a continued fraction via modified Lentz's method,
using generators rather than functions.
This implementation is by the book [1]_.
Parameters
----------
a: generator yielding numeric.
b: generator yielding numeric.
tol: float [default: 1.e-10]
Tolerance for termination of evaluation.
N_min: int [default: 0]
Minimum number of iterations to evaluate.
N_max: int or comparable [default: np.Inf]
Maximum number of iterations to evaluate.
tiny: float [default: 1.e-30]
Very small number to control convergence of Lentz's method when
there is cancellation in a denominator.
Returns
-------
(float, float, int)
The first element of the tuple is the value of the continued
fraction. The second element is the estimated error. The third
element is the number of iterations.
References
----------
.. [1] WH Press, SA Teukolsky, WT Vetterling, BP Flannery,
"Numerical Recipes," 3rd Ed., Cambridge University Press 2007,
ISBN 0521880688, 9780521880688 .
Examples
--------
Use generators to compute the square root of 2:
>>> from qnm.contfrac import lentz_gen
>>> import itertools
>>> def rt2b_g():
... yield 1
... for x in itertools.repeat(2):
... yield x
...
>>> def rt2a_g():
... for x in itertools.repeat(1):
... yield x
...
>>> lentz_gen(rt2a_g(), rt2b_g())
(1.4142135623638004, 4.488287519421874e-11, 14)
See the documentation for :meth:`lentz` for more examples.
"""
f_old = next(b) # 0
if (f_old == 0):
f_old = tiny
C_old = f_old
D_old = 0.
conv = False
j = 1
while ((not conv) and (j < N_max)):
aj, bj = next(a), next(b) # j and j. 'a' started with 1.
D_new = bj + aj * D_old
if (D_new == 0):
D_new = tiny
C_new = bj + aj / C_old
if (C_new == 0):
C_new = tiny
D_new = 1./D_new
Delta = C_new * D_new
f_new = f_old * Delta
if ((j > N_min) and (np.abs(Delta - 1.) < tol)): # converged
conv = True
# Set up for next iter
j = j + 1
D_old = D_new
C_old = C_new
f_old = f_new
# Success or failure can be assessed by the user
return f_new, np.abs(Delta - 1.), j-1