# Worksheet for computing long reals> restart; Digits := 60:# First we define a procedure for doing the grunt work> longreal := proc(e)> local c,b, first, base, mant, i;> c := e;> b := simplify(ceil(ln(c)/ln(2)));> first := false;> base := 0;> mant := {};> for i from 0 to 52 while c>0 do>    if c - 1/(2^(-b+i)) >= 0>       then if not first then>          first := true;>          base := b-i>       end if;>       mant := mant union {i-b+base};>       c := c - 1/(2^(-b+i));>       if c = 0 then break end if>    end if> end do;> return([base,mant]);> end:# We feed in the real number we want to represent> a:= 6.023*10^23;> lr := longreal(a);# Now let's compute (to sufficient accuracy) what this representation's# exact real value happens to be.> ma := 0: for j in lr[2] do ma := evalf(ma + 1/(2^j)) end do: ma;> ac := evalf(ma*2^lr[1]);# With absolute error> abs(a - ac);# and relative error> abs((a-ac)/a);# Try it out with various values of a. Note that the relative error# always has about the same order of magnitude (unless the# representation is exact)> # Next, let's wrap the whole approximation into a proc so we can examine# computation problems# > approxlr := proc(e)> local mb, j, lrl, cb;> lrl := longreal(e):> mb := 0: for j in lrl[2] do >    mb := evalf(mb + 1/(2^j)) > end do:> cb := mb*2^lrl[1];> return(cb);> end:# Now use this to find the decimal version of the long real# approximations> f1 := 1/1000000;> f2 := 1/1000001;> r1 := approxlr(f1);> r2 := approxlr(f2);# Compute the relative errors in adding and subtracting these number:> appsum := approxlr(r1+r2);> exactsum := f1 + f2;> relerrsum := abs((exactsum - appsum)/exactsum);> appdiff := approxlr(r1-r2);> exactdiff := f1 - f2;> relerrdiff := abs((exactdiff - appdiff)/exactdiff);# Evaluate a function which has a subtraction problem> f:= x -> sqrt(x^2+1) - 1;> f(1.01);> ans := f(0.01);> > g := x -> x^2/(sqrt(x^2+1)+1);> g(0.01);> Digits := 10;> apprf := f(0.01);> apprg := g(0.01);> Digits := 60;> abs((ans - apprf)/ans);> abs((ans - apprg)/ans);>