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...
| [reply] [d/l] [select] |
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] |
$ 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] |
Update: Added before/after results for Pi(1e9) to emphasize the performance increase.
I did an experiment on an i7-Broadwell machine running CentOS 7.3. The system gmp-devel libraries were built to support the least common denominator. Thus, leaving out optimizations for newer CPU instructions.
# default package build
MPN_PATH=" x86_64/k8 x86_64 generic"
# manually on the broadwell machine, includes extra optimizations
MPN_PATH=" x86_64/coreibwl x86_64/coreihwl x86_64/coreisbr x86_64/core
+inhm x86_64/core2 x86_64 generic"
I have everything point to the /opt/perl-5.26.1 folder to not impact the system in any way. Then, re-installed GMP modules including Math::Prime::Util.
cd gmp-6.1.2
./configure --prefix=/opt/perl-5.26.1
make -j 4
make check
make install
cd Math-GMP-2.15
/opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include"
+ LIBS="-L/opt/perl-5.26.1/lib -lgmp"
make
make test
make install
cd Math-BigInt-GMP-1.6004
/opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include"
+ LIBS="-L/opt/perl-5.26.1/lib -lgmp"
make
make test
make install
cd Math-Prime-Util-GMP-0.46
/opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include"
+ LIBS="-L/opt/perl-5.26.1/lib -lgmp"
make
make test
make install
cd Math-Prime-Util-master
/opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include"
+ LIBS="-L/opt/perl-5.26.1/lib -lgmp"
make
make test
make install
Beyond 20% gain is possible by simply building GMP manually. That will set MPN_PATH to utilize the full capabilities of the processor.
time /opt/perl-5.26.1/bin/perl -Mntheory=:all -E 'say Pi(1e6)' | wc -c
1000002
0m 1.524s using system gmp libraries
0m 1.168s using gmp libraries from /opt/perl-5.26.1/lib
time /opt/perl-5.26.1/bin/perl -Mntheory=:all -E 'say Pi(1e9)' | wc -c
1000000002
91m 7.581s using system gmp libraries
67m 40.951s using gmp libraries from /opt/perl-5.26.1/lib
Regards, Mario | [reply] [d/l] [select] |