Closed GoogleCodeExporter closed 8 years ago
Christof,
You've stumbled across a corner case in the GMP library. The GMP library offers
minimal support for floating point. You might want to look at mpmath to get more
predictable results.
Below is a technical description of what is happening.
I'll assume a 32-bit system. The precision of a GMP mpf is just a "minimum"
precision
and GMP guarantees that "at least" that many bits are calculated. When gmpy
requests
a default precision of 53 bits, GMP rounds the precision up to a multiple of 32
(to
64), and then adds one more word at each end. The value of the mantissa is not
normalized so the most significant word (or limb, as GMP calls it) may contain
anywhere between 1 bit and 32 bits. And then depending on the underlying
algorithm,
the least significant limb may either be 0 or contain the next 32 bits of a
result.
Even though gmpy requested 53 bits of precision, GMP returns results with
anywhere
between 65 and 128 bits of precision. gmpy normally rounds the least
significant limb
away to reduce the result to have between 65 and 96 bits of precision.
The number of bits requested in fround are the number of bits that GMP will
guarantee
in a result. If we use the above paragraph with a requested precision of 0, we
would
end in a situation where the actual precision is somewhere between 1 and 32
bits (the
non-normalized most significant limb). Since this isn't useful, gmpy enforces a
minimum precision of 53 bits and that's the value you are getting.
gmpy may eventually get better floating point support but I would use mpmath
right now.
If this doesn't make sense, feel free to reply with more questions.
Thanks for using gmpy!
casevh
Original comment by casevh
on 6 Feb 2010 at 4:15
Hi Casevh,
Looks like I misunderstood what fround is supposed to do. What I wanted was
rounding
with a given digit precision of n, not with a bit precision.
Thanks for the interesting explanation anyway. I will have a look at mpmath
Christof
Original comment by christof...@googlemail.com
on 6 Feb 2010 at 8:51
Closing.
Original comment by casevh
on 7 Feb 2010 at 7:57
Original issue reported on code.google.com by
christof...@googlemail.com
on 5 Feb 2010 at 9:47