Beefy Boxes and Bandwidth Generously Provided by pair Networks
Welcome to the Monastery

Re^4: porting C code to Perl

by syphilis (Chancellor)
on Nov 07, 2017 at 12:34 UTC ( #1202921=note: print w/replies, xml ) Need Help??

in reply to Re^3: porting C code to Perl
in thread porting C code to Perl

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 ( 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:

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.

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).


Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://1202921]
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others taking refuge in the Monastery: (5)
As of 2018-06-20 19:16 GMT
Find Nodes?
    Voting Booth?
    Should cpanminus be part of the standard Perl release?

    Results (117 votes). Check out past polls.