BrianGladman / mpir

Multiple Precision Integers and Rationals
GNU General Public License v3.0
75 stars 36 forks source link

Fix for mpf_get_2exp_d on Win64. MPIR.Net updates for that API and for mpf_cmp_z. #10

Closed adyache closed 5 years ago

adyache commented 5 years ago

The point of the mpf_get_2exp_d update is to fix returned values on Win64, where the mpf exponent is an mp_bit_cnt, which is a 32-bit int, and the exponent value (stored as # of limbs) would exceed 32 bits when recalculated as # of bits. See included test.

SanderBouwhuis commented 5 years ago

@adyache Does your fix also influence C++? I see you talking about .Net only.

adyache commented 5 years ago

Yes, mpf_get_2exp_d is a C++ API. When I created a .Net wrapper for it and started testing, I realized that it was returning incorrect values for large exponents, so this includes my proposed fix.

BrianGladman commented 5 years ago

I have done the merge now

BrianGladman commented 5 years ago

The original GMP functions: double mpz_get_d_2exp(signed long exp2, mpz_srcptr src) double mpf_get_d_2exp(signed long exp2, mpf_srcptr src)

both return a signed long value through a pointer, which can lead to problems because GMP is designed to use 64 bit integers whereas Windows defines longs to be 32 bits.

This can mean that using a long pointer will result in a 64-bit value being written to a pointer address for a 32-bit value and this will result in corruption of the memory just above this value. .

The revised 'Windows' versions reverse the define so that it is the double returned via a pointer which avoids the potentially disastrous results of memory corruption since doubles are the same length on Linux and Windows (the function return value might still be truncated and cause errors but at least it won't now corrupt memory).

SanderBouwhuis commented 5 years ago

So, only if the exponent goes above 2^32 - 1 we would run into trouble?

BrianGladman commented 5 years ago

So, only if the exponent goes above 2^32 - 1 we would run into trouble?

Assuming that you are asking about the new Windows versions of these functions:

mpir_si mpz_get_2exp_d (double r, mpz_srcptr src) mpir_si mpf_get_2exp_d (double r, mpz_srcptr src)

their return values are valid in the range -2^63 .. 2^63 - 1.

If you place the return values in a Windows signed long value the result will be valid in the range -2^31 .. 2^31 - 1.

If you place the result in an unsigned long the result will be valid in the range 0 .. 2^31 - 1.

SanderBouwhuis commented 5 years ago

So you are not changing the original functions, but adding new ones? Wouldn't int64* have been better than returning a double? It's not at all intuitive that the returned double has the range of a 64-bit signed int.

BrianGladman commented 5 years ago

The old functions are still there on Windows x64 since they are required for backwards compatibility.

You have misunderstood how the new functions work. The old and the new functions are not API compatible with the old ones and require a change in any user code that calls them.

Both functions return two values - a signed 64-bit integer and a floating point double. Both the old and the new functions return two identical values, the difference being that the new functions return the double through a pointer and the signed 64-bit integer through the function's return value. The old functions do the reverse of this, returning the 64-bit integer through a pointer and the double via the function return value.

SanderBouwhuis commented 5 years ago

Aaaaaaaaaaaaah... now I get it! Ok, that makes a lot more sense. Thank you for your continued effort in supporting this math lib. Just this week I permanently replaced my other math lib with your MPIR / MPFR library.

adyache commented 5 years ago

Just my 2c as a clarification: the exponent of an mpf number is stored in an mp_bit_cnt, which on Win64 is a signed int32, a.k.a. long a.k.a int, not a signed int64 a.k.a. long long as on linux. It is however expressed in limbs, which means it's capable of expressing binary exponents in the range -2^37..2^37-1, the binary exponent being the exponent in limbs multiplied by 64, plus change. The fix in this PR corrected this calculation of exp in bits from exp in limbs, so it no longer overflows and allows the entire range above to be returnable. The overflow in question affected both new and old APIs, because one of them calls the other.