Henry,
Here is some code but the product continued fractions are not pretty.
Everyone else,
This is almost a small article.
To conform with Abramowitz and Stegun conventions, continued fractions
are always encoded
b0 , a1 , b1 , a2 , b2 , a3 , b3 , ...
which represents
b0 + a1
-------
b1 + a2
-------
b2 + a3
-------
b3 + ...
The part shown is the third convergent. The first convergent ends with
b1, the second with b2, and so on. Yes, the 0th convergent ends with b0.
cv =: (1 0 $~ #) # +`%/\ NB. Calculates convergents
pi4cf NB. Continued fraction for pi%4 (A and S 4.4.43 page 81)
0 1 1 1 3 4 5 9 7 16 9 25 11 36 13 49 15 64 17 81 19 100 21
9!:11 [ 4 NB. Print precision 4 figures
,. 7 {. cv pi4cf NB. Show first seven convergents
0
1
0.75
0.7917
0.7843
0.7856
0.7854
o. 1%4 NB. pi%4
0.7854
AB is my name for an array
(A0 , B0) ,. (A1 , B1) ,. (A2 , B2) ,. ...
which has numerators and denominators for the k^th convergent Ak % Bk .
The array is defined by
buildAB =: 3 : 0
NB. y is b0,a1,b1,a2,b2, ...
b0 =. {. y
ab =. _2 ]\ }. y
AB =. 1 0 ,. b0 , 1
while. # ab do.
AB =. AB ,. (_2 {."1 AB) +/ . * {. ab
ab =. }. ab
end.
}."1 AB
)
7 {."1 pi4AB =: buildAB pi4cf NB. Show columns through ,. A6,B6
0 1 3 19 160 1744 2.318e4
1 1 4 24 204 2220 2.952e4
,. 7 {. %/ pi4AB NB. Show first seven convergents from Ak % Bk
(compare above)
0
1
0.75
0.7917
0.7843
0.7856
0.7854
buildcf =: 3 : 0 NB. Builds continued fraction from table AB
NB. y is (A0,B0),.(A1,B1),. ...
b0 =. (<0 0) { y
AB =. 1 0 ,. y
ab =. ''
while. 2 < {: $ AB do.
ab =. ab , (2 {"1 AB) %. 2 {."1 AB
AB =. }."1 AB
end.
b0 , , ab
)
buildcf x: pi4AB NB. compare to pi4cf above
0 1 1 1 3 4 5 9 7 16 9 25 11 36 13 49 15 64 17 81 19 100 21
ln2cf NB. Continued fraction for ^. 2 (ln 2) (A and S 4.1.39 page
68)
0 1 1 1 2 1 3 4 4 4 5 9 6 9 7 16 8 16 9 25 10 25 11
ln2AB =: buildAB ln2cf
pi4ln2AB =: pi4AB * ln2AB NB. AB for product (o. 1%4)*(^. 2)
7 {."1 pi4ln2AB
0 1 6 133 5760 3.628e5 3.645e7
1 1 12 240 1.061e4 6.66e5 6.695e7
pi4ln2cf =: buildcf x: pi4ln2AB NB. PRODUCT CONTINUED FRACTION
13 {. pi4ln2cf NB. product not pretty
0 1 1 6 6 26 107r6 2372r13 456r13 248192r593 31615r593 660555r554 45171r554
,. 7 {. %/ pi4ln2AB NB. Show first seven convergents from Ak % Bk
0
1
0.5
0.5542
0.543
0.5447
0.5444
,. 7 {. x:^:_1 cv pi4ln2cf NB. Show first 7 convergents from product
0
1
0.5
0.5542
0.543
0.5447
0.5444
(o. 1%4)*(^.2)
0.5444
Kip Murray
Kip Murray wrote:
Henry,
I do not have code but can point you to formulas in the Handbook of
Mathematical Functions by Abramowitz and Stegun, Chapter 3 Elementary
Mathematical Methods, Section 3.10 Theorems on Continued Fractions
(page 19 in my edition). Be sure to read both columns of formulas.
There is no unique continued fraction product for two continued
fractions, but you could try the following to find _a_ product. First
given two infinite continued fractions f and g use the recursive
matrix multiplication given in formula (3) of Section 3.10 to find the
numerators and denominators of the n^th convergents
f_n = A_n/B_n
n = 1,2,3,...
g_n = C_n/D_n
(conventional notation, _n means subscript n). The convergents
converge to the value of the continued fraction.
Let
E_n = A_n * C_n
n = 1,2,3...
F_n = B_n * D_n
Then, to find a continued fraction h (this is your product) with
convergents
h_n = E_n/F_n n = 1,2,3...
solve the matrix recursion (3) for for the entries e_n and f_n in
h = f_0 + e_1
---------
f_1 + e_2
--------
f_2 + ...
(I hope that comes out right on your screen.)
The recursion (3) written for the continued fraction h is
(E_n , F_n) = ((E_(n-1) , F_(n-1)) ,. (E_(n-2) , F_(n-2))) o (f_n , e_n)
using a mix of standard notation and J, where o is +/ . *
Starting values are E_(-1) = 1, E_0 = f_0, F_(-1) = 0, F_0 = 1. You
can use f_0 = 0.
Please get the correct formulas from Abramowitz and Stegun -- in case
I mistyped something. And of course you have %. for solving for (f_n
, e_n).
Kip Murray
Henry Rich wrote:
Does anyone have J code for multiplying two numbers expressed as
simple continued fractions, producing a continued-fraction result?
Henry Rich
------------------------------------------------------------------------
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
------------------------------------------------------------------------
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
