GMP's rudimentary mpf layer
I think "rudimentary" is a fair description, and the GMP's mpf layer documentation even states that "new projects should consider using the GMP extension library MPFR (http://mpfr.org/) instead".
I personally think the same advice should be given wrt "old projects", too.
One limitation of the mpf layer is that the allocated precision is always a multiple of limb size - where limb size will commonly be either 32 or 64 bits.
Another limitation is the absence of both Inf and NaN.
Also, the rounding rules are rather odd.
The mpf_get_d function which converts the mpf_t to a double (which it then returns) rounds towards zero (ie truncates) - and that seems to be the rounding rule generally used throughout.
But if one assigns a higher precision mpf_t to a lower precision mpf_t one doesn't always get what one expects.
The following demo has two 128-bit mpf_t variables ($large_1 and $large_2) that have been assigned slightly different values.
However, using any sane rounding rules, both of those values should round to the same 64-bit value - as both values have the first 65 bits set && there are other set bits further along.
To test this, I assign $large_1 to the 64-bit $little_1, and I assign $large_2 to the 64-bit $little_2, expecting that $little_1 == $little_2. But Rmpf_cmp() reports that $little_1 != $little_2.
And Rmpf_get_str() also shows that $little_1 != $little_2.
Rmpf_out_str() outputs identical base 2 representations for $little_1 and $little_2 but displays 65 set bits ... which seems rather strange for a 64-bit precision mpf_t:
use strict;
use warnings;
use Math::GMPf qw(:mpf);
my $large_1 = Rmpf_init2(128);
my $large_2 = Rmpf_init2(128);
my $little_1 = Rmpf_init2(64);
my $little_2 = Rmpf_init2(64);
Rmpf_set_str($large_1, ('1' x 65) . ('0' x 2) . ('1' x 61), 2);
Rmpf_set_str($large_2, ('1' x 65) . ('0' x 42) . ('1' x 21), 2);
Rmpf_set($little_1, $large_1);
Rmpf_set($little_2, $large_2);
if(Rmpf_cmp($little_1, $little_2)) {print "WTF\n"}
print $little_1, "\n", $little_2, "\n";
Rmpf_out_str($little_1, 2, 0);
print "\n";
Rmpf_out_str($little_2, 2, 0);
print "\n";
print Rmpf_get_prec($little_1), "\n";
print Rmpf_get_prec($little_2), "\n";
__END__
OUTPUTS:
WTF
0.340282366920938463456e39
0.340282366920938463454e39
0.11111111111111111111111111111111111111111111111111111111111111111e12
+8
0.11111111111111111111111111111111111111111111111111111111111111111e12
+8
64
64
These are not Math::GMPf bugs, because the following C script outputs the same results (though the format is not identical).
#include <stdio.h>
#include <gmp.h>
int main(void) {
mp_exp_t exponent;
char *out_1, *out_2;
mpf_t large_1, large_2, little_1, little_2;
char *str_1 =
"1111111111111111111111111111111111111111111111111111111111111111100
+1111111111111111111111111111111111111111111111111111111111111";
char *str_2 =
"1111111111111111111111111111111111111111111111111111111111111111100
+0000000000000000000000000000000000000000111111111111111111111";
mpf_init2(large_1, 128);
mpf_init2(large_2, 128);
mpf_init2(little_1, 64);
mpf_init2(little_2, 64);
mpf_set_str(large_1, str_1, 2);
mpf_set_str(large_2, str_2, 2);
mpf_set(little_1, large_1);
mpf_set(little_2, large_2);
if(mpf_cmp(little_1, little_2)) printf("WTF\n");
out_1 = mpf_get_str(NULL, &exponent, 10, 0, little_1);
printf("%s %d\n", out_1, exponent);
out_2 = mpf_get_str(NULL, &exponent, 10, 0, little_2);
printf("%s %d\n", out_2, exponent);
mpf_out_str(stdout, 2, 0, little_1);
printf("\n");
mpf_out_str(stdout, 2, 0, little_2);
printf("\n");
printf("%d\n%d\n", mpf_get_prec(little_1),
mpf_get_prec(little_2));
return 0;
}
/*
OUTPUTS:
WTF
340282366920938463456 39
340282366920938463454 39
0.11111111111111111111111111111111111111111111111111111111111111111e12
+8
0.11111111111111111111111111111111111111111111111111111111111111111e12
+8
64
64
*/
The mpf documentation does state:
<quote>
Internally, GMP sometimes calculates with higher precision than that of the destination variable in order to limit errors. Final results are always truncated to the destination variable’s precision.
</quote>
I suspect that this might account for the above anomalies - though I don't see any obvious connection.
One avoids these types of anomalies with the MPFR library.
Doing the same exercise with MPFR, one finds that $little_1 == $little_2. And that value is 3.40282366920938463463e38 (rounding to nearest), and 3.40282366920938463445e38 (rounding down).
Cheers,
Rob | [reply]
[d/l]
[select] |
Out of curiosity, I looked at perllocal.pod. It looks like Math::MPFR installed automatically during the Math::Prime::Util cpan installation.
$ vi /opt/perl-5.26.1/lib/5.26.1/darwin-thread-multi-2level/perllocal.
+pod
Math::Prime::Util::GMP
Math::Prime::Util
-- auto-installed Math::MPFR
-- auto-installed Math::BigInt::GMP
Regarding use integer (thanks haukex), the spigot Perl code ran 1.36 times faster on my Mac. Math::Prime::Util is already incredibly fast. Now wonder in amazement on what you will do next.
| [reply]
[d/l]
[select] |