Ref: Kahan (Jan 2007): How Futile are Mindless Assessments of Roundoff in Floating-Point Computation ?

I wanted to know, how Maple behaves on the example given in \2475 of that working paper (where other

systems are considered, mostly Matlab) and thought it might be of interest here in the forum.

> restart; Digits:=16; interface(version);

Consider the following recursion (Muller 1980, Kahan 2007)

> x[0] = 4, x[1] = 4.25, x[n+1] = (108-(815-1500/x[n-1])/x[n]); x:='x':

which can be easily coded in Maple as

> K0:=proc(n::nonnegint)
option remember; # to improve speed
if n = 0 then
return 4;
elif n = 1 then
return 4.25;
else
108-(815-1500/K0(n-2))/K0(n-1);
end if;
end proc:

That numerical routine gives

> 'K0(80)': '%'=%;

while coding with rationals gives a different result:

> K:=proc(n::nonnegint) option remember;
if n = 0 then return 4 end if;
if n = 1 then return 17/4 end if;
108-(815-1500/K(n-2))/K(n-1);
end proc:

'K(80)': '%'=evalf(%), '%'=%;

Since working with rationals is 'exact', the correct answer should be 5.

Indeed: Kahan shows the recursion is solved by

> ( a*3^(n+1) + b*5^(n+1)+ c*100^(n+1) ) / ( a*3^(n) + b*5^(n)+ c*100^(n) ):
eval(%,[a=1,b=1,c=0]):
X:=unapply(%,n);
'X(n)': '%'=%;

so 5 is the limit for large n.

Obviously there are rounding errors, Kahan talks about that.

I wanted to know, how Maple behaves. First convert to rationals in a lazy way:

> K1:=proc(n::nonnegint) option remember;
if n = 0 then
return 4;
elif n = 1 then
return 4.25;
else
108-(815-1500/K1(n-2))/K1(n-1);
end if;
return convert(%,rational);
end proc:
'K1(80)': '%'=evalf(%); #, '%'=%;

Still the same problem ...

One has to inrease precision considerably (I could not verify, that 0.7*N Digits are

enough for a good result, as stated in the paper):

> N:=80;
oldDigits:=Digits:
forget(K1): gc(); # clear tables for various N
Digits:=max(round(3/2*N),Digits); # increase precision
'K1(N)': evalf(%): '%%'= evalf[oldDigits](%);
Digits:=oldDigits:

The correct work with rationals is in *each* recursion step (which is K from above),

in K1 the 2nd step was not converted to rationals

> K2:=proc(n::nonnegint) option remember;
local result;
if n = 0 then
result:= 4;
elif n = 1 then
result:= 17/4.0;
else
result:=108-(815-1500/K2(n-2))/K2(n-1);
end if;
return convert(result,rational);
end proc:

> N:=80; forget(K2): gc();
'K2(N)': '%'= evalf(%);

>

This post was generated using the MaplePrimes File Manager

View 1_mindless.mw on MapleNet or Download 1_mindless.mw
View file details