Barton,
So is nfloat(I(42),[],20) = 8.7541921419809776312b-3 the right answer? It is certainly faster than setting fpprec:1000
and calling bfloat(I(42)). That is the brute force way which I really don't like.
Rich
--------------------------------------------------
From: "Barton Willis" <willisb at unk.edu>
Sent: Friday, March 26, 2010 3:51 PM
To: "Maxima List" <maxima at math.utexas.edu>
Subject: unstable recursions & nfloat
> It's a standard example, but I thought somebody might enjoy the Maxima
> way.
> We need nfloat, so load hypergeometric.
>
> (%i1) load(hypergeometric)$
>
> (%i2) i(n) := if n = 0 then 1-1/%e else n * i(n-1) - 1/%e$
>
> (%i17) [float(i(42)),
> float(expand(i(42))),
> block([fpprec : 52], bfloat(i(42))), block([fpprec : 56],
> bfloat(i(42))),
> nfloat(i(42),[],20)];
>
> (%o17) [-4.7467928326479132*10^34,
> 0.0,
> 2.732067241502428018740764014347935258990086011897504b-2,
> 8.7536981302551482132924825234864659298743235134456890016b-3,
> 8.7541921419809776312b-3]
>
> Solving the recursion for I(42) is identical to solving a linear
> system with coefficient matrix
>
> (%i19) m : genmatrix(lambda([i,j], if i-j = 0 then 1 else if i - j = 1
> then -i else 0), 42,42)$
>
> And the condition number of m is big--and I mean big
>
> (%i22) mat_norm(m,inf) * mat_norm(m^^-1, inf);
> (%o22) 103810448683368212798606246049077583803159875197431095
>
>
>
> Barton
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>