The function nfloat uses a running error to approximately bound the
rounding error. When the
bound is too large, nfloat declares a do-over. Unlike interval methods
that properly set the
fpu rounding mode, it's possible for the running error bound to be too
small. See
http://www.amazon.com/Accuracy-Stability-Numerical-Algorithms-Nicholas/dp/0898715210
Example: (in %i28, the sum is quoted)
(%i9) load(hypergeometric)$
(%i28) nfloat('(sum(x^k/k!,k,0,999)),[x = -142.0], 42);
(%o28) 2.13886596489953912812261349689714081777545b-62
The sum is too ill-conditioned to sum it accurately using 42 decimal
digits:
(%i29) sum((-142.0b0)^k/k!,k,0,999), fpprec : 42;
(%o29) -3.08915603322448457558053810090231199675679b16
The exact value of the sum is (compare to %o28)
(%i30) bfloat(sum((-142)^k/k!,k,0,999)), fpprec : 42;
(%o30) 2.13886596489953912812261349689714081777545b-62
--Barton
maxima-bounces at math.utexas.edu wrote on 08/26/2010 10:16:04 AM:
> [image removed]
>
> Re: [Maxima] last bit function ?
>
> Richard Fateman
>
> to:
>
> Sheldon Newhouse
>
> 08/26/2010 10:16 AM
>
> Sent by:
>
> maxima-bounces at math.utexas.edu
>
> Cc:
>
> "maxima at math.utexas.edu"
>
> There are lots of interval arithmetic packages written in Lisp, at
> least some of which can be altered
> to be used from Maxima. (No bigfloats, but that could be added.) See,
> for example, the generic arithmetic package
> I wrote.
>
> There is a lot more activity in "Reliable Computation" that you can read
> about, including different approaches to
> improving accuracy, e.g. "taylor" methods, and methods that reduce
> imprecision by rearrangements or by maintaining
> dependencies.
>
> evaluating a polynomial of degree 30 at an interval point should
> definitely NOT be done by doing the "+" and "*" by
> interval operations. if x=[-1,1] then x*x is [-1,1]. but x^2 is
> [0,1].
>
> See http://www.cs.berkeley.edu/~fateman/papers/interval.pdf for some
> info and links into literature.
>
> There may even be interval packages for Maxima.
>
> None of this means you cannot write your own interval stuff for the fun
> of it.
>
>
> The short answer to your questions about last bits ... look at the lisp
> functions. For Maxima, try..
>
> ?integer\-decode\-float(1.2d5)
>
> For bfloats, try
> ?cadr (1.2b5)
>
> Note that if you set fpprec:16, then the number of BITS in the bfloat
> fraction (bfloat arithmetic is BINARY not decimal)
> can be viewed by looking at
> ?fpprec
>
> which is, in this case, 56.
>
> So if you care only about getting same accuracy, probably you will get
> about the same results if you abandon
> your "power of 10" calculations and do "power of 2" calculations with
> ?fpprec.
>
> You will still be rather inefficient, but using directed roundings in
> IEEE arithmetic, the "right" way to do some of this,
> is inhibited by some hardware which makes it inconvenient, and most
> software systems which compound the difficulty
> by not giving you access to the hardware.
>
>
>
>
> RJF
>
>
>
>
> On 8/26/2010 7:27 AM, Sheldon Newhouse wrote:
> > Hello,
> > Is there a fast implementation in maxima giving the last bit of a
> > float and bfloat? I am interested in using this to develop a fast and
> > accurate interval package. I have some naive ways to do this, but
> > they are very inefficient and give vast overestimation when compared
> > to interval packages which exist in C++ and fortran (at least in the
> > double precision versions).
> ....
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima