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 <d...@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 _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits