#!/usr/bin/env perl
use strict;
use warnings;
use Math::PariInit qw( primes=10000000 stack=1e8 );
use Math::Pari qw/pari2iv/;
use Math::Prime::Util qw/:all/;
use Data::Dumper;
$|=1;

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 = 80_000;
print "Comparing for small inputs: 0 - $small\n";

foreach my $n (0 .. $small) {
  print '.' unless ($n+1) % int($small/80);
  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);

  for my $k (2,3,9,10) {
    die "valuation($n,$k)" unless Math::Pari::valuation($n,$k) == valuation($n,$k);
  }

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

  if ($n > 1) {
    for (1..10) {
      my $k;  do { $k = int(rand(50)) } while !($k % $n);
      die "binomial($n,$k)" unless Math::Pari::binomial($n,$k) == binomial($n,$k);
      my $negn = - ($n >> 1);
      die "binomial($negn,$k)" unless Math::Pari::binomial($negn,$k) == binomial($negn,$k);
    }
  }

  {
    my $d = $n+3;
    my @gmpu  = gcdext($n,$d);
    my $gpari = Math::Pari::bezout($n,$d);
    die "gcdext($n,$d)" unless $gmpu[0] == $gpari->[0] && $gmpu[1] == $gpari->[1] && $gmpu[2] == $gpari->[2];
  }

  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 $invmod = invmod($m, $n);
    if (defined $invmod) {
      die "invmod($m, $n)" unless Math::Pari::lift(PARI "Mod(1/$m,$n)") == $invmod;
    } else {
      eval { PARI "Mod(1/$m,$n)" };
      die "invmod($m, $n) defined in Pari" unless $@ =~ /impossible inverse/
        || ($m == 0 && $@ =~ /division by zero/);
    }
  }

  { 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) == $order
    } 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;

  {
  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 { [ 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);

  for my $k (2,3,9,10) {
    die "valuation($n,$k)" unless Math::Pari::valuation($n,$k) == valuation($n,$k);
  }

  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";
  # 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,$n);
  die "jordan_totient(3,$n)"
    unless Math::Pari::sumdiv($n,"d","d^3*moebius($n/d)")
        == jordan_totient(3,$n);

  if ($n > 2) {
    for (1..10) {
      my $k;  do { $k = int(rand(10)) } while !($k % $n);
      die "binomial($n,$k)" unless Math::Pari::binomial($n,$k) == binomial($n,$k);
      my $negn = - ($n >> 1);
      die "binomial($negn,$k)" unless Math::Pari::binomial($negn,$k) == binomial($negn,$k);
    }
  }

  # 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: primorial / pn_primorial

  # TODO: carmichael lambda?  Pari doesn't have it.

  { my $m = int(rand($n-1));
    my $invmod = invmod($m, $n);
    if (defined $invmod) {
      die "invmod($m, $n)" unless Math::Pari::lift(PARI "Mod(1/$m,$n)") == $invmod;
    } else {
      eval { PARI "Mod(1/$m,$n)" };
      die "invmod($m, $n) defined in Pari" unless $@ =~ /impossible inverse/
        || ($m == 0 && $@ =~ /division by zero/);
    }
  }

  { 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) == $order;
    } else {
      eval { Math::Pari::znorder($mn); };
      die "znorder($m, $n) defined in Pari" unless $@ =~ /not an element/;
    }
  }

  # TODO: znlog with reasonable values

  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 = $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, 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 = $small + int(rand(10**12));
  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)" );
}

# Pari's primepi in 2.1-2.5 is strangely lacking

{ my $a = (int(rand(2**32)) << 32) + int(rand(2**32));
  my $b = (int(rand(2**32)) << 32) + int(rand(2**32));
  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);
  $a >>= 1 if $a > 2**63;
  die "kronecker(-$a,$b)" unless Math::Pari::kronecker(-$a,$b) == kronecker(-$a,$b);
  $b >>= 1 if $b > 2**63;
  die "kronecker($a,-$b)" unless Math::Pari::kronecker($a,-$b) == kronecker($a,-$b);
  die "kronecker(-$a,-$b)" unless Math::Pari::kronecker(-$a,-$b) == kronecker(-$a,-$b);

  { my @gmpu  = gcdext($a,$b);   my $gpari = Math::Pari::bezout($a,$b);
    die "gcdext($a,$b)" unless $gmpu[0] == $gpari->[0] && $gmpu[1] == $gpari->[1] && $gmpu[2] == $gpari->[2]; }
}
{ my $a = int(rand(2**32));
  my $b = int(rand(2**32));
  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;
}




use Bytes::Random::Secure qw/random_string_from/;
sub ndigit_rand {
  my($digits, $howmany) = @_;
  die "digits must be > 0" if $digits < 1;
  $howmany = 1 unless defined $howmany;
  my @nums = map { random_string_from("123456789",1) . random_string_from("0123456789",$digits-1) } 1 .. $howmany;
  if (10**$digits > ~0) {  @nums = map { Math::BigInt->new($_) } @nums;  }
  else                  {  @nums = map { int($_) } @nums;                }
  return wantarray ? @nums : $nums[0];
}


sub array_compare {
  my($a1, $a2, $text) = @_;
  #eq_or_diff $a1, $a2, $text;
  die "$text wrong count ",scalar @$a1," ",scalar @$a2 unless @$a1 == @$a2;
  foreach my $i (0 .. $#$a1) {
    if (ref($a1->[$i])) {
      array_compare($a1->[$i],$a2->[$i], "> $text");
    } else {
    #print "a1: ", Dumper($a1), "\na2: ", Dumper($a2), "\n" unless $a1->[$i] == $a2->[$i];
      die "$text entry $i  $a1->[$i] != $a2->[$i]" unless $a1->[$i] == $a2->[$i];
    }
  }
}