The continued logarithm of π
The continued logarithm of π
Briefly, a continued logarithm is an arbitrarily long bit string approximating a real number arbitrarily well, and supports arbitrarily precise bit-at-a-time algorithms for rational functions of these numbers. (See http://www.tweedledum.com/rwg/cfup.htm , p 47+)
A six-bits-per-term π series:
In[]:=
Product[MatrixForm@{{512,0},{32(-37+42k),}},{k,∞}]MatrixForm[{{ooπ,0},{oo,"?"}}]
3
k
3
k
3
(-1+2k)
Out[]=
∞
∏
k
512 3 k | 0 |
32 3 k | 3 (2k-1) |
ooπ | 0 |
oo | ? |
Where ∏ means matrix product, not Mathematica product, and oo is some quantity that blows up with the number of product terms, the same way you compute continued fractions with 2×2 matrices. An incorrect form of this series is derived in https://dspace.mit.edu/handle/1721.1/6088 .
We initialize the work matrix m to the first term of the matrix product:
We initialize the work matrix m to the first term of the matrix product:
In[]:=
MatrixForm[m={{512,0},{32(-37+42k),}}/.k1]
3
k
3
k
3
(-1+2k)
Out[]//MatrixForm=
512 | 0 |
160 | 1 |
This represents the function
In[]:=
Divide@@(m.{t,1})
Out[]=
512t
1+160t
where t is the tail of the series, starting with k=2 rather than 1. It should be easy to show that, in general,
In[]:=
3/8/k<t<3/8/(k+1)
Out[]=
3
8k
3
8(1+k)
giving bounds
m/.{{t3/8/1.},{t3/8/2.}}
Out[]=
{3.14754,3.09677}
Since these both exceed 2, we commence output with
In[]:=
Style[cl@π={1},22]
Out[]=
{1}
and left multiply m by the divide-by-2 matrix:
In[]:=
MatrixForm[m={{1,0},{0,2}}.m]
Out[]//MatrixForm=
512 | 0 |
320 | 2 |
It costs almost nothing to remove the common power of 2:
In[]:=
MatrixForm[m=m/2]
Out[]//MatrixForm=
256 | 0 |
160 | 1 |
representing the function
In[]:=
Divide@@(m.{t,1})
Out[]=
256t
1+160t
We still are on input term k=1, and can still use
In[]:=
%/.{{t3/8/1.},{t3/8/2.}}
Out[]=
{1.57377,1.54839}
which is smack between 1 and 2, which we celebrate with
In[]:=
Style[AppendTo[cl@π,0],22]
Out[]=
{1,0}
and left multiply m by the subtract-1-and-reciprocate matrix:
In[]:=
(m={{0,1},{1,-1}}.m)//MatrixForm
Out[]//MatrixForm=
160 | 1 |
96 | -1 |
representing
In[]:=
Divide@@(m.{t,1})
Out[]=
1+160t
-1+96t
Again using
In[]:=
%/.{{t3/8/1.},{t3/8/2.}}
Out[]=
{1.74286,1.82353}
which dictates another
In[]:=
Style[AppendTo[cl@π,0],22]
Out[]=
{1,0,0}
and
In[]:=
(m={{0,1},{1,-1}}.m)//MatrixForm
Out[]//MatrixForm=
96 | -1 |
64 | 2 |
representing
In[]:=
Divide@@(m.{t,1})
Out[]=
-1+96t
2+64t
In[]:=
%/.{{t3/8/1.},{t3/8/2.}}
Out[]=
{1.34615,1.21429}
In[]:=
Style[AppendTo[cl@π,0],22]
Out[]=
{1,0,0,0}
In[]:=
(m={{0,1},{1,-1}}.m)//MatrixForm
Out[]//MatrixForm=
64 | 2 |
32 | -3 |
In[]:=
Divide@@(m.{t,1})
Out[]=
2+64t
-3+32t
In[]:=
%/.{{t3/8/1.},{t3/8/2.}}
Out[]=
{2.88889,4.66667}
This unambiguously exceeds 2, so
In[]:=
Style[AppendTo[cl@π,1],22]
Out[]=
{1,0,0,0,1}
In[]:=
MatrixForm[m={{1,0},{0,2}}.m]
Out[]//MatrixForm=
64 | 2 |
64 | -6 |
In[]:=
Divide@@(m.{t,1})
Out[]=
2+64t
-6+64t
In[]:=
%/.{{t3/8/1.},{t3/8/2.}}
Out[]=
{1.44444,2.33333}
Our interval of uncertainty contains 2! At last we gobble (right multiply) the k=2 term.
In[]:=
MatrixForm[m=m.{{512,0},{32(-37+42k),}}/.k2]
3
k
3
k
3
(-1+2k)
Out[]//MatrixForm=
286208 | 54 |
189952 | -162 |
In[]:=
Divide@@(m.{t,1})
Out[]=
54+286208t
-162+189952t
(Remembering to bump k !)
In[]:=
%/.{{t3/8/2.},{t3/8/3.}}
Out[]=
{1.51515,1.51938}
(Giving six more bits of precision!)
In[]:=
Style[AppendTo[cl@π,0],22]
Out[]=
{1,0,0,0,1,0}
In[]:=
(m={{0,1},{1,-1}}.m)//MatrixForm
Whoa, a burst of three 1s!
(Since we knew it was three ones and not four, we could have output 1,1,1,0 and left multiplied m by
This process may be continued indefinitely, or until the integers in m overflow.