Re: unexpected rounding behaviour of mpf_add(). bug?

[Frederik Schaffalitzky]

On 30 Apr 2001, Kevin Ryde wrote:

> Frederik Schaffalitzky <address@hidden> writes:

> Some routines can generate extra bits on the end of "a", above the
> requested precision.  Other routines, like mpf_add, end up truncating
> those, giving the results you find.

Yes, I found out by looking through the sources, although to me it looks
as if what is really happening is not excess precision on the end of "a"
but truncation during the addition. All my operands have precision
(_mp_prec) 3 but the square root has size (_mp_size) 4, i.e. it uses the
full mantissa available for that precision. In the addition routine, the
least significant limb is discarded before it is detected that the smaller
operand was completely cancelled when lining up the exponents so the
result ends up having only three limbs of precision (which is less than
the 64 bits that I asked for). It puzzles me that although the number of
limbs allocated is always _mp_prec + 1, often this full mantissa is not
used (being one limb less).


> > ............................. Is there some way I can make mpf_t
> > behave in the way that I desire?
> 
> Not easily, I think.  Maybe it suffices to use an assumed epsilon
> which is a/2^256 or some such.

Unfortunately my algorithm is (C++) templated over scalar type and I would
prefer not to put in per-type epsilons etc. However, I did find one way to
solve my problem, which was to zero out the least significant limb of any
mpf_t for which _mp_size == _mp_prec + 1. Since the addition routine does
that anyway, not much is lost... :)


> mpfr might give you more joy (see www.mpfr.org for the latest), it has
> guaranteed ieee style rounding.  Not sure there's an actual function
> to give an epsilon, but there's a vaguely internal-use
> mpfr_add_one_ulp for adding one.

That looks interesting (and it appears I already have the sources), thank
you. I don't really want to know what epsilon is, but just want to be able
to assume that the arithmetic behaves as if there is one.

Thanks for your help,

Frederik