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

[Kevin Ryde]

Frederik Schaffalitzky <address@hidden> writes:
> 
> holds exactly (as deemed by mpf_cmp()), for many numbers `a' there is no
> power of 1/2 which has the property that
>   a + epsilon == a

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.

I'm not exactly sure why the extras are held, I think it ends up being
more efficient to hold them than to strip them when it turns out
they're not needed for the target precision.

> I am not sure if this is strictly speaking a bug as GMP makes few
> guarantees about rounding behaviour, but on the other hand it is hard to
> imagine an addition algorithm which produces such results and which isn't
> rather careless about trailing digits. 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.

> [Whatever the merit of this practice, I have a numerical algorithm which
> assumes this sort of behaviour for its termination criterion.]

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.