Re: Continued fraction problem
- To: mathgroup at smc.vnet.net
- Subject: [mg17120] Re: [mg17039] Continued fraction problem
- From: "Wolf, Hartmut" <hwolf at debis.com>
- Date: Sat, 17 Apr 1999 03:35:12 -0400
- Organization: debis Systemhaus
- References: <199904140611.CAA21070@smc.vnet.net.>
- Sender: owner-wri-mathgroup at wolfram.com
Alan W.Hopper schrieb:
>
> A quite artificial but nevertheless interesting construction
> in recreational maths is the Champernowne Constant, defined
> in Eric Weisstein's Concise Encyc. of Math.,
> as the decimal 0.1234567891011... , obtained by concatenating
> the positive integers.
>
> The continued fraction expansion of this number is stated to
> be ; (0,8,9,1,149083,1,1,1,4,1,1,1,3,4,1,1,1,15, a 166 digit #,...,
> (position 41)a 2504 digit #,...,(position 163)a 33102 digit #,...
>
> And using this code for the continued fraction period in list form ,
>
> In[1]:= cf[x_Real, n_]:=Module[{ip, fp=x, result = {}},
> Do[ip = Floor[fp];
> AppendTo[result, ip];
> fp = 1/(fp-ip), {n}];
> result];
>
> In[2]:= cf[0.12345678910111213141516171819202122,20]
>
> Out[3]=
> (Error messages detailing Indeterminate,
> Infinite expressions and ComplexInfinity.)
>
> {0,8,9,1,149083,1,1,1,4,1,1,1,3,4,1,1,1,15,1787142709274,
> ComplexInfinity}
>
> All correct up to where the 166 digit term should appear.
> (Inserting , N[....., 200] in the code does not help).
>
> The result above is better however than with ;
>
> <<NumberTheory`ContinuedFractions`
>
> ContinuedFraction[0.123456789101112131415161718192021, 20]
>
> which gives 150571 as the fifth term.
>
> I suppose that rounding errors and inadequate precision are coming
> into play with this situation as it is now.
>
> So I wonder if anyone knows a workable method of obtaining the correct
> continued fraction period of a difficult case like the Champernowne con=
st.
> with Mathematica, say up to that enormous 163 rd term, if it is at all
> feasible?
>
Dear Alan,
I think so. What we have to do is to carefully secure for sufficient prec=
ision for all calculation steps, and indeed there are some pitfalls.
My first presumption was: the problem is not the algorithm, but the preci=
sion of the calculation. This can best be seen if you don't use Real, but=
Rational arithmetic. Roughly speaking, each aN takes away its own length=
from the precision of the input. So to get your 2504 digit number at pos=
tion 41, you have to prepare your input at a precision of at least 3000 d=
ecimal places!
I took your (the straight forward) algorithm
In[20]:= continuedFraction[x_Rational,n_Integer]:=
Module[{aN,xN=x,r={}},
Do[aN=Floor[xN];
AppendTo[r,aN];
xN=(xN-aN)^-1,
{n}];
r]
To prepare the Champernowne Constant (i.e. an approximation to it) I firs=
t tried the algorithm:
f[x_,y_Integer]:=(ff*=10^-Ceiling[Log[10,y+1]];x+ ff*y)
ff=1;x=Fold[f,0,Range[13]]
12345678910111213/100000000000000000
This works fine say for y=400 (giving more than a thousand digits in 10=
seconds) but to tease your 33102 monster at position 163 I estimated a r=
untime of more than a year (on my weak machine AMD x486 at 133 MHz, so pe=
rhaps you can do it in a quarter).
So I came back to the most simple idea to prepare the Champernowne Consta=
nt through string concatenation. That works:
In[12]:= zStringList=ToString/@Range[30000];//Timing
Out[12]= {13.42 Second,Null}
In[13]:= zString=StringJoin@@zStringList;//Timing
Out[13]= {0.18 Second,Null}
In[14]:= StringLength[zString]
Out[14]= 138894
In[15]:= 9+90*2+900*3+9000*4+20001*5
Out[15]= 138894
In[16]:= xNenner=10^%;//Timing
Out[16]= {3.886 Second,Null}
In[17]:= xZ=Ehler=ToExpression[zString];//Timing
Out[17]= {785.329 Second,Null}
In[18]:= x=xZ=Ehler/xNenner;//Timing
Out[18]= {386.826 Second,Null}
Now we calculate the continued fraction
In[23]:= continuedFraction[x,500]//Timing
Out[23]= {28.371 Second,
{0,8,9,1,149083,1,1,1,4,1,1,<..many, many numbers..>,1,1,12,3,1,1,1,9,1,=
2,1}}
In[24]:= Ceiling[Log[10,#+0.5]]&/@ %[[2]]
Out[24]= {0,1,1,1,6,1,1,1,1,1,1,1,1,1,1,1,1,2,166,1,1,1,2,1,1,
1,1,1,1,1,1,1,1,2,1,3,1,2,1,2,2504,1,1,1,1,1,1,2,1,
1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,
1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
140,1,2,1,1,2,1,1,1,1,1,1,3,1,1,2,1,1,1,1,1,1,1,1,2,
1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,2,1,1,
1,1,1,1,1,1,1,1,1,1,33102,2,1,2,1,1,1,1,1,1,1,1,1,2,
1,1,2,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,2,1,1,1,1,2,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,2,1,1,1,1,1,1,
1,1,1,1,1,3,2,1,1,1,1,1,1,1,1,2,2,1,109,1,1,1,1,1,2,
3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,
2,2,1,2,1,2,1,1,1,1,1,1,2,1,1,1,1,1,2,1,1,1,2,1,1,2,
1,2,1,2,1,1,1,2,1,1,1,1,1,1,1,1,1,1,2,2,1,2,1,1,1,1,
1,1,1,1,1,1,1,2,1,2,2,1,1,1,1,1,2,2,1,1,1,1,1,1,1,1,
2468,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,2,1,1,1,2,
1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,2,1,1,1,3,1,1,1,1,1,1,1,2,1,1,1,2,1,1,
1,1,3,1,1,1,1,2,1,1,1,2,1,2,1,1,1,1,1,1,1,1,1,1,1,1,
136,1,1,1,2,1,1,1,1,2,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1}
In[25]:= Position[%,_?(#>3&)]
Out[25]= {{5},{19},{41},{102},{163},{247},{358},{460}}
In[26]:= Extract[%%,%]
Out[26]= {6,166,2504,140,33102,109,2468,136}
In[27]:= Plus@@%%%
Out[27]= 39198
In[28]:= Extract[%%%%%[[2]],Extract[%%%,{{1},{8}}]]
Out[28]= {149083,
1282409532859743054814016581451672169597214608073139006012029770452718\
223534299625694581836197983168908272712969697509897267073018852743}
Well, I'm not at all shure for that 136 digit number at position 460, but=
I indeed found a 33102 digit monster at position 163 as you promised! (i=
t's in my notebook at mailto:hwolf at debis.com if you want to see it). It t=
ook nearly 20 minutes to get Champerowne's Constant into my machine, but =
less than 30 seconds to spot that monster.
Well coming back, let's now try to do it with real arithmetic:
In[80]:= continuedFraction[x_Real,n_Integer]:=Module[{aN,xN=x,r={ }},
Do[aN=Floor[xN];
AppendTo[r,aN];
xN=(xN-aN)^-1,
{n}];
r]
In[87]:= f[x_,y_]:=(ff*=10^-Ceiling[Log[10,y+1]];x+ N[ff,30]*y)
In[90]:= ff=1; x=N[Fold[f,N[0.,30],N[Range[13],30]],30]
Out[90]= 0.1234567891011121
Not quite! The precision of this calculation is totally unsufficient! Why=
does it come only with
In[91]:= $MachinePrecision
Out[91]= 16
? I don't know (can someone tell me?) Looking at the formula: ff is a pur=
e rational, as is an Integer the Ceiling of the logarithm, which should b=
e precise enough. x and y (through Range[13]) are entered into the calcul=
ation x + ff * y with precision 30, so where does it get lost?
Anyways we can prepare the Champerowne's as a rational, and convert that =
to real
In[57]:= f[x_,y_]:=(ff*=10^-Ceiling[Log[10,y+1]];x+ ff*y)
In[72]:= ff=1; xRational=Fold[f,0,Range[200]];
In[70]:= xReal=N[xRational,500]
Out[70]=
0.12345678910111213141516171819202122232425262728293031323334353637383940=
41424\
3444546474849505152535455565758596061626364656667686970717273747576777879=
80818\
2838485868788899091929394959697989910010110210310410510610710810911011111=
21131\
1411511611711811912012112212312412512612712812913013113213313413513613713=
81391\
4014114214314414514614714814915015115215315415515615715815916016116216316=
41651\
6616716816917017117217317417517617717817918018118218318418518618718818919=
01911\
9219319419519619719819920000000000
In[74]:= continuedFraction[xReal,20]
Out[74]= {0,8,9,1,149083,1,1,1,4,1,1,1,3,4,1,1,1,15,
45754011139103107648364662824295611859960393971045755500066200439309026=
26592\
5631493795320774712865631386412093755035520946071830899845758014698631488=
33592\
141783010987,6}
In[75]:= Ceiling[Log[10,#+0.5]]&/@%
Out[75]= {0,1,1,1,6,1,1,1,1,1,1,1,1,1,1,1,1,2,166,1}
There were no error messages (but you will get some if you don't enter th=
e calculation with sufficient precision) and the value of that 166 digit =
number is the same as I get with pure rational arithmetic. The 5th number=
confirms to your testimony.
Kind regards
Hartmut Wolf
- References:
- Continued fraction problem
- From: "Alan W.Hopper" <awhopper@hermes.net.au>
- Continued fraction problem