Beefy Boxes and Bandwidth Generously Provided by pair Networks
We don't bite newbies here... much

Re^3: porting C code to Perl

by danaj (Friar)
on Oct 26, 2017 at 02:50 UTC ( [id://1202042] : note . print w/replies, xml ) Need Help??

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

GMP is installed on many systems now, partly because GCC started requiring it. MPFR is much less common. GMP is a big help for lots of things, but MPFR much less so, and I used it via Math::MPFR so you'd need that installed also. It looks like it is used in prime_count_lower and prime_count_upper for large enough values, and the floating point routines harmreal, ExponentialIntegral, LogarithmicIntegral, RiemannZeta, RiemannR, and Pi. For Pi it's a very small savings vs.the standard GMP back end. But the floating point things like RiemannZeta are just hard to do in Perl or GMP's rudimentary mpf layer (for me at least). Using Math::MPFR (thanks syphilis!) is thousands of times faster than Math::BigFloat even using the fast GMP backend. The MPFR implementors are very sharp and also better practicing numerical analysts than me.

For AWS EC2 machines I usually just install the packages with apt-get or yum. I also do git and clone/build/install the latest versions of my modules but that's because I want the absolute latest. cpan is fine and as long as you have GMP you should just be able to cpan Math::Prime::Util Math::BigInt::GMP. But make sure it did find GMP (try cpan Math::Prime::Util::GMP to check).

14 seconds for Pi to 50k digits... Ouch. Moar fastar! At least it's faster than the default Math::BigFloat (which doesn't have the luxury of dropping into custom C code).

As Discipulus says I should try use integer as well to test the Perl spigot. I usually don't like that pragma because it's signed which is never what I want, but that's not relevent here. Now I'm pondering a digression into how to implement use uinteger. At first glance it doesn't look terrible hard to do in the core, but using a whole new hint bit is unlikely to fly. Maybe a test just to see if it would work...

Replies are listed 'Best First'.
Re^4: porting C code to Perl
by syphilis (Archbishop) on Nov 07, 2017 at 12:34 UTC
    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).

Re^4: porting C code to Perl
by marioroy (Prior) on Oct 26, 2017 at 06:48 UTC

    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.