[libmath-prime-util-perl] 176/181: Updated Pari compare test
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:19 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit 3e7296520ced51e703064c327bcb5831d736fc54
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Jan 13 18:21:38 2014 -0800
Updated Pari compare test
---
xt/pari-compare.pl | 205 ++++++++++++++++++++++++++++++++++++++++++++++-------
1 file changed, 179 insertions(+), 26 deletions(-)
diff --git a/xt/pari-compare.pl b/xt/pari-compare.pl
index a7d83af..666cf99 100755
--- a/xt/pari-compare.pl
+++ b/xt/pari-compare.pl
@@ -1,37 +1,174 @@
#!/usr/bin/env perl
use strict;
use warnings;
-use Math::Prime::Util qw/:all/;
use Math::PariInit qw( primes=10000000 stack=1e8 );
-use Math::Pari;
+use Math::Pari qw/pari2iv/;
+use Math::Prime::Util qw/:all/;
use Data::Dumper;
-use Test::Differences;
$|=1;
-# For these tests we're looking at small inputs.
+BEGIN {
+ use Config;
+ die "Tests have 64-bit assumptions" if $Config{uvsize} < 8;
+ die "Tests need double floats" if $Config{nvsize} < 8;
+ no Config;
+}
+
+my $small = 100000;
+print "Comparing for small inputs: 0 - $small\n";
+
+foreach my $n (0 .. $small) {
+ die "isprime($n)" unless Math::Pari::isprime($n) == !!is_prime($n);
+ die "is_prob_prime($n)" unless Math::Pari::isprime($n) == !!is_prob_prime($n);
+ die "next_prime($n)" unless Math::Pari::nextprime($n+1) == next_prime($n);
+ die "prev_prime($n)" unless Math::Pari::precprime($n-1) == prev_prime($n);
+ next if $n == 0;
+
+ my($pn,$pc) = @{Math::Pari::factorint($n)};
+ my @f1 = map { [ pari2iv($pn->[$_]), pari2iv($pc->[$_])] } 0 .. $#$pn;
+ array_compare( \@f1, [factor_exp($n)], "factor_exp($n)" );
+ @f1 = map { ($_->[0]) x $_->[1] } @f1;
+ array_compare( \@f1, [factor($n)], "factor($n)" );
+
+ array_compare( [map { pari2iv($_) } @{Math::Pari::divisors($n)}], [divisors($n)], "divisors($n)" );
+
+ die "omega($n)" unless Math::Pari::omega($n) == factor_exp($n);
+ die "bigomega($n)" unless Math::Pari::bigomega($n) == factor($n);
+ die "numdiv($n)" unless Math::Pari::numdiv($n) == divisors($n);
+
+ foreach my $k (0..4) {
+ die "sigma($n,$k)" unless Math::Pari::sigma($n,$k) == divisor_sum($n,$k);
+ }
+
+ die "moebius($n)" unless Math::Pari::moebius($n) == moebius($n);
+ die "euler_phi($n)" unless Math::Pari::eulerphi($n) == euler_phi($n);
+
+ my $d = PARI "d";
+ die "jordan_totient(2,$n)"
+ unless Math::Pari::sumdiv($n,"d","d^2*moebius($n/d)")
+ == jordan_totient(2,$n);
+ die "jordan_totient(3,$n)"
+ unless Math::Pari::sumdiv($n,"d","d^3*moebius($n/d)")
+ == jordan_totient(3,$n);
+
+ die "nth_prime($n)" unless Math::Pari::prime($n) == nth_prime($n);
+
+ # All the pari2iv calls are very time-consuming
+ if ($n < 1000) {
+ array_compare( [map { pari2iv($_) } @{Math::Pari::primes($n)}], primes(nth_prime($n)), "primes($n)" );
+ }
+
+ # Math Pari's forprime is super slow for some reason. Pari/gp isn't this slow.
+ if ($n < 1000) {
+ my $m = $n+int(rand(10**4));
+ PARI "s1=0";
+ PARI "forprime(X=$n,$m,s1=s1+X)";
+ my $s1 = PARI('s1');
+
+ my $s2 = 0; forprimes { $s2 += $_ } $n,$m;
+ die "forprimes($n,$m) $s1 != $s2" unless $s1 == $s2;
+ }
+ {
+ my $d = PARI "d";
+ my @a1; Math::Pari::fordiv($n, $d, sub { push @a1, pari2iv($d)});
+ my @a2; fordivisors { push @a2, $_ } $n;
+ array_compare( \@a1, \@a2, "fordivisors($n)" );
+ }
+
+ { my $m = int(rand($n-1));
+ my $mn = PARI "Mod($m,$n)";
+ my $order = znorder($m, $n);
+ if (defined $order) {
+ die "znorder($m, $n)" unless Math::Pari::znorder($mn) == znorder($m,$n);
+ } else {
+ eval { Math::Pari::znorder($mn); };
+ die "znorder($m, $n) defined in Pari" unless $@ =~ /not an element/;
+ }
+ }
+
+ # Pari's znprimroot is iffy for non-primes
+ if (is_prime($n)) {
+ my $g = znprimroot($n);
+ die "znprimroot($n)" unless Math::Pari::znprimroot($n) == $g;
+ my $a = 1 + int(rand($n-2));
+ my $gn = PARI "Mod($g,$n)";
+ my $log = znlog($a, $g, $n);
+ die "znlog($a, $g, $n) should be defined" unless defined $log;
+ die "znlog($a, $g, $n)" unless Math::Pari::znlog($a,$gn) == $log;
+ }
+
+ if ($n < 100) {
+ foreach my $d (0 .. 9) {
+ my $arg = $n + $d/10;
+ next if $arg < 0.1;
+ my $e1 = -Math::Pari::eint1(-$arg);
+ my $e2 = ExponentialIntegral($arg);
+ die "ExponentialIntegral($arg) $e1 != $e2" if abs($e1 - $e2) > $e1*1e-14;
+ }
+ }
+ if ($n > 1) {
+ my $arg = $n;
+ my $e1 = -Math::Pari::eint1(-log($arg));
+ my $e2 = LogarithmicIntegral($arg);
+ die "LogarithmicIntegral($arg) $e1 != $e2" if abs($e1 - $e2) > $e1*1e-14;
+ }
+
+ {
+ my $s = 50.0/$small;
+ if ($s != 1.0) {
+ my $zeta1 = Math::Pari::zeta($s) - 1;
+ my $zeta2 = RiemannZeta($s);
+ die "zeta($s) $zeta1 != $zeta2" if abs($zeta1 - $zeta2) > abs($zeta1) * 1e-14;
+ }
+ }
+
+ print "." unless $n % 1250;
+}
+
+print "\nkronecker, gcd, and lcm for small values\n";
+foreach my $a (-400 .. 400) {
+ foreach my $b (-400 .. 400) {
+ # Pari 2.1's gcd doesn't work right for 0,-x and -x,0. Pari 2.2.3 fixed.
+ if ($a != 0 && $b != 0) {
+ die "gcd($a,$b)" unless Math::Pari::gcd($a,$b) == gcd($a,$b);
+ }
+ die "kronecker($a,$b)" unless Math::Pari::kronecker($a,$b) == kronecker($a,$b);
+ die "lcm($a,$b)" unless Math::Pari::lcm($a,$b) == lcm($a,$b);
+ }
+ print "." unless (400+$a) % 20;
+}
+
+print "\nloop forever with random values\n";
+
+# forcomposites in Pari 2.6, not Math::Pari's 2.1
my $loops = 0;
while (1) {
+ my $n;
-{ my $n = (int(rand(2**32)) << 32) + int(rand(2**32));
+ {
+ do { $n = (int(rand(2**32)) << 32) + int(rand(2**32)) } while $n < $small;
die "isprime($n)" unless Math::Pari::isprime($n) == !!is_prime($n);
die "is_prob_prime($n)" unless Math::Pari::isprime($n) == !!is_prob_prime($n);
die "next_prime($n)" unless Math::Pari::nextprime($n+1) == next_prime($n);
die "prev_prime($n)" unless Math::Pari::precprime($n-1) == prev_prime($n);
my($pn,$pc) = @{Math::Pari::factorint($n)};
- my @f1 = map { [ int("$pn->[$_]"), int("$pc->[$_]")] } 0 .. $#$pn;
+ my @f1 = map { [ pari2iv($pn->[$_]), pari2iv($pc->[$_])] } 0 .. $#$pn;
array_compare( \@f1, [factor_exp($n)], "factor_exp($n)" );
@f1 = map { ($_->[0]) x $_->[1] } @f1;
array_compare( \@f1, [factor($n)], "factor($n)" );
- my @d1 = map { int("$_") } @{Math::Pari::divisors($n)};
- array_compare( \@d1, [divisors($n)], "divisors($n)" );
+ array_compare( [map { pari2iv($_) } @{Math::Pari::divisors($n)}], [divisors($n)], "divisors($n)" );
die "omega($n)" unless Math::Pari::omega($n) == factor_exp($n);
die "bigomega($n)" unless Math::Pari::bigomega($n) == factor($n);
die "numdiv($n)" unless Math::Pari::numdiv($n) == divisors($n);
+ foreach my $k (0..4) {
+ die "sigma($n,$k)" unless Math::Pari::sigma($n,$k) == divisor_sum($n,$k);
+ }
+
die "moebius($n)" unless Math::Pari::moebius($n) == moebius($n);
die "euler_phi($n)" unless Math::Pari::eulerphi($n) == euler_phi($n);
@@ -39,51 +176,63 @@ while (1) {
# TODO: our jordan_totient should auto-bigint
die "jordan_totient(2,$n)"
unless Math::Pari::sumdiv($n,"d","d^2*moebius($n/d)")
- == jordan_totient(2,Math::BigInt->new("$n"));
+ == jordan_totient(2,$n);
die "jordan_totient(3,$n)"
unless Math::Pari::sumdiv($n,"d","d^3*moebius($n/d)")
- == jordan_totient(3,Math::BigInt->new("$n"));
+ == jordan_totient(3,$n);
# TODO: exp_mangoldt:
# Lambda(n)={
# v=factor(n);
# if(matsize(v)[1]!=1,return(0),return(log(v[1,1])));
# };
+ # TODO: chebyshev_theta, chebyshev_psi
# Chebyshev Psi(x)=sum(n=2,floor(x),Lambda(n));
# TODO: partitions. new Pari has this as numbpart.
# See OEIS A000041 for some alternate Pari functions
-# TODO:
-# chebyshev_theta chebyshev_psi
-# divisor_sum
-# carmichael_lambda znorder znprimroot znlog legendre_phi
-# ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
+ # TODO: primorial / pn_primorial
+ # TODO: carmichael lambda? Pari doesn't have it.
-}
+ { my $m = int(rand($n-1));
+ my $mn = PARI "Mod($m,$n)";
+ my $order = znorder($m, $n);
+ if (defined $order) {
+ die "znorder($m, $n)" unless Math::Pari::znorder($mn) == znorder($m,$n);
+ } else {
+ eval { Math::Pari::znorder($mn); };
+ die "znorder($m, $n) defined in Pari" unless $@ =~ /not an element/;
+ }
+ }
-{ my $n = 1+int(rand(10**5)); # Pari doesn't like 0
- die "nth_prime($n)" unless Math::Pari::prime($n) == nth_prime($n);
-}
+ # TODO: znlog with reasonable values
-{ my $n = int(rand(10**4));
- array_compare( [map { int("$_") } @{Math::Pari::primes($n)}], primes(nth_prime($n)), "primes($n)" );
-}
+ if ($n > 1) {
+ my $arg = $n;
+ my $e1 = -Math::Pari::eint1(-log($arg));
+ my $e2 = LogarithmicIntegral($arg);
+ die "LogarithmicIntegral($arg) $e1 != $e2" if abs($e1 - $e2) > $e1*1e-12;
+ }
+ # TODO: RiemannZeta
+ }
-{ my $a = int(rand(10**6));
+
+
+{ my $a = $small + int(rand(10**6));
my $b = $a+int(rand(10**4));
my $x = PARI "x";
- my @a1; Math::Pari::forprime($x,$a,$b,sub { push @a1, int("$x") });
+ my @a1; Math::Pari::forprime($x,$a,$b,sub { push @a1, pari2iv($x) });
my @a2; forprimes { push @a2, $_ } $a,$b;
array_compare( \@a1, \@a2, "forprimes($a,$b)" );
}
# forcomposites in Pari 2.6, not Math::Pari's 2.1
-{ my $n = int(rand(10**12));
+{ my $n = $small + int(rand(10**12));
my $d = PARI "d";
- my @a1; Math::Pari::fordiv($n, $d, sub { push @a1, int("$d") });
+ my @a1; Math::Pari::fordiv($n, $d, sub { push @a1, pari2iv($d) });
my @a2; fordivisors { push @a2, $_ } $n;
array_compare( \@a1, \@a2, "fordivisors($n)" );
}
@@ -106,6 +255,10 @@ while (1) {
die "lcm($a,$b)" unless Math::Pari::lcm($a,$b) == lcm($a,$b);
}
+{ my $n = random_prime(10000,~0);
+ die "znprimroot($n)" unless Math::Pari::znprimroot($n) == znprimroot($n);
+}
+
$loops++;
print "." unless $loops % 100;
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list