Re: [Gcl-devel] Apology (was: Bug in mpz_get_d)

[Enrique Perez-Terron]

On Mon, 2009-12-14 at 15:11 +0100, Torbjorn Granlund wrote:

> Please explain what you get and what you expect.  *One* example of
> incorrect+correct data is what we need.

I aplolgise that I did not check the list of recipients more carefully.
I meant to follow up to the Gcl-devel list, and was not aware that
the parent post already had added gmp-bugs to the recipients.  Reporting
a bug to Gmp was mentioned as an option later in the thread, and I am
confused to find that the original poster had already done so.

> Enrique Perez-Terron <address@hidden> writes:
> 
>   If stability is desired, mpfr is the one to use. GMP does not aim at
>   stability, as far as I understand.
> 
> You misunderstand.  GMP very much aims at stability.

You are right.  Mpfr adds the option to specify the rounding policy.
That does not mean that Gmp's rounding is unpredictable.  My thinking
and my choice of words were inexact.

Gmp seems to round toward zero, rather than to nearest.  I appreciate
that users of a package like Gmp normally should set the required
precision to suits their needs with a reasonable margin, and not get
hurt if the rounding is another than the canonical one.

Different application areas require different rounding policies, and
rounding to zero is not unreasonable, especially if it is faster.

Now lisp is an application area independent tool, and it has this
(float) command that converts to a hardware-supported floating point
type.  I suppose rounding to nearest is the least surprising to users of
such general software. I suppose there is some cost to this rounding,
but suspect this natural trade-off tips in favor of agreeing with the
canonical rounding.

For the benefit of the CC-ed list readers, I add some explanatory
remarks.

In the test data used by the OP, the value of "a" happens to be

   a == 2^1000  - 1.

In a binary representation, this is a 1000-bit field with all ones.
Adding one to this number produces a 1001-bit field with a leading one
and a thousand zeros. One can expect rounding policy transitions to
happen here.

The "standard" floating-point format splits a 64-bit field in three
areas (s,e,m) of 1, 11, and 52 bits respectively.  If these fields are
thought as containing unsigned integers s,e,m, the value of the number
represented is normally

  (-1)^s  *  2^(e - 1023) * ( 1 + m/2^52 ).

This gives a choice between representing "a" as

  2^1000     (with an "m" part that is all zeros)

or

  2^1000 - 2^947    (with an "m" part that is all ones)

It appears that Gmp truncates the '1's in the binary representation of
"a" that do not fit in the "m" field. This is rounding toward zero, and
chooses the second alternative above.

The leading 20 digits of the number 2^1000 - 2^947 in a decimal
representation are

  1071 5086 0718 6267 2019

The leading 20 digits of the number 2^1000 in a decimal representation
are

  1071 5086 0718 6267 3209

These are the digits that show in the printout from the OP's post and
from my little test program. This explains the differences in the
outputs as effects of rounding.

The Mpfr library has an additional parameter to its functions, that
specifies the desired rounding policy. In the test program I used
rounding to nearest and got for 2^1000 - 1, 2^1000 rather than 2^1000 -
2^947.

-Enrique