Re^4: porting C code to Perl

[syphilis]

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