#!/usr/bin/env perl
use strict;
use warnings;

use Test::More;
use Math::Prime::Util qw/prime_count twin_prime_count
                         prime_count_lower prime_count_upper
                         prime_count_approx twin_prime_count_approx
                         ramanujan_prime_count/;

my $isxs  = Math::Prime::Util::prime_get_config->{'xs'};
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $extra = 0+(defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING});

#  Powers of 2:  http://oeis.org/A007053/b007053.txt
#  Powers of 10: http://oeis.org/A006880/b006880.txt
my %pivals32 = (
                   1 => 0,
                  10 => 4,
                 100 => 25,
                1000 => 168,
               10000 => 1229,
              100000 => 9592,
             1000000 => 78498,
            10000000 => 664579,
           100000000 => 5761455,
          1000000000 => 50847534,
               30239 => 3269,
               30249 => 3270,
               60067 => 6062,
               65535 => 6542,
            16777215 => 1077871,
          2147483647 => 105097565,
          4294967295 => 203280221,
);
my %pivals64 = (
         10000000000 => 455052511,
        100000000000 => 4118054813,
       1000000000000 => 37607912018,
      10000000000000 => 346065536839,
     100000000000000 => 3204941750802,
    1000000000000000 => 29844570422669,
   10000000000000000 => 279238341033925,
  100000000000000000 => 2623557157654233,
 1000000000000000000 => 24739954287740860,
10000000000000000000 => 234057667276344607,
         68719476735 => 2874398515,
       1099511627775 => 41203088796,
      17592186044415 => 597116381732,
     281474976710655 => 8731188863470,
    4503599627370495 => 128625503610475,
   72057594037927935 => 1906879381028850,
 1152921504606846975 => 28423094496953330,
18446744073709551615 => 425656284035217743,
);
my %pivals_small = map { $_ => $pivals32{$_} }
                   grep { ($_ <= 2000000) || $extra }
                   keys %pivals32;

#  ./primesieve 1e10 -o2**32 -c1
#  ./primesieve 24689 7973249 -c1
my %intervals = (
  "868396 to 9478505" => 563275,
  "1118105 to 9961674" => 575195,
  "24689 to 7973249" => 535368,
  "1e10 +2**16" => 2821,
  "17 to 13"    => 0,
  "0 to 1"      => 0,
  "0 to 2"      => 1,
  "1 to 3"      => 2,
  "3 to 17"     => 6,
  "4 to 17"     => 5,
  "4 to 16"     => 4,
  "191912783 +248" => 2,
  "191912784 +247" => 1,
  "191912783 +247" => 1,
  "191912784 +246" => 0,

  "1e14 +2**16" => 1973,
  "127976334671 +468" => 2,
  "127976334672 +467" => 1,
  "127976334671 +467" => 1,
  "127976334672 +466" => 0,
);
delete @intervals{ grep { (parse_range($_))[1] > ~0 } keys %intervals };

my %tpcs = (
              5000 =>           126,
            500000 =>          4565,
          50000000 =>        239101,
        5000000000 =>      14618166,
      500000000000 =>     986222314,
    50000000000000 =>   71018282471,
  5000000000000000 => 5357875276068,
);

my %rpcs = (
              5000 =>           302,
             50000 =>          2371,
            500000 =>         19492,
           5000000 =>        165440,
            135791 =>          5888,
             65536 =>          3030,
);

plan tests => 0 + 1
                + 3*scalar(keys %pivals32)
                + scalar(keys %pivals_small)
                + $use64 * 3 * scalar(keys %pivals64)
                + scalar(keys %intervals)
                + 1
                + 5 + 2*$extra # prime count specific methods
                + 3 + (($isxs && $use64) ? 1+2*scalar(keys %tpcs) : 0) # twin pc
                + 2 + (($isxs && $use64) ? 2+1*scalar(keys %rpcs) : 0) # ram pc
                + 0;

ok( eval { prime_count(13); 1; }, "prime_count in void context");

while (my($n, $pin) = each (%pivals32)) {
  cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
  cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
  my $approx_range = abs($pin - prime_count_approx($n));
  my $range_limit = ($n <= 100000000) ? 100 : 500;
  cmp_ok( $approx_range, '<=', $range_limit, "prime_count_approx($n) within $range_limit");
}
while (my($n, $pin) = each (%pivals_small)) {
  is( prime_count($n), $pin, "Pi($n) = $pin" );
}
if ($use64) {
  while (my($n, $pin) = each (%pivals64)) {
    cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
    cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
    my $approx = prime_count_approx($n);
    my $percent_limit = 0.0005;
    # This is the test we want:
    #cmp_ok( abs($pin - $approx) / $pin, '<=', $percent_limit/100.0, "prime_count_approx($n) within $percent_limit\% of Pi($n)");
    # Math rearranged so we don't lose all precision.
    cmp_ok( abs($pin - $approx) * (100.0 / $percent_limit), '<=', $pin, "prime_count_approx($n) within $percent_limit\% of Pi($n)");
  }
}

while (my($range, $expect) = each (%intervals)) {
  my($low,$high) = parse_range($range);
  is( prime_count($low,$high), $expect, "prime_count($range) = $expect");
}

# Defect found in prime binary search
is( prime_count(130066574), 7381740, "prime_count(130066574) = 7381740");

sub parse_range {
  my($range) = @_;
  my($low,$high);
  my $fixnum = sub {
    my $nstr = shift;
    $nstr =~ s/^(\d+)e(\d+)$/$1*(10**$2)/e;
    $nstr =~ s/^(\d+)\*\*(\d+)$/$1**$2/e;
    die "Unknown string in test" unless $nstr =~ /^\d+$/;
    $nstr;
  };
  if ($range =~ /(\S+)\s+to\s+(\S+)/) {
    $low = $fixnum->($1);
    $high = $fixnum->($2);
  } elsif ($range =~ /(\S+)\s*\+\s*(\S+)/) {
    $low = $fixnum->($1);
    $high = $low + $fixnum->($2);
  } else {
    die "Can't parse test data";
  }
  ($low,$high);
}

# TODO: intervals.  From primesieve:
#    155428406, // prime count 2^32 interval starting at 10^12
#    143482916, // prime count 2^32 interval starting at 10^13
#    133235063, // prime count 2^32 interval starting at 10^14
#    124350420, // prime count 2^32 interval starting at 10^15
#    116578809, // prime count 2^32 interval starting at 10^16
#    109726486, // prime count 2^32 interval starting at 10^17
#    103626726, // prime count 2^32 interval starting at 10^18
#    98169972}; // prime count 2^32 interval starting at 10^19

# Make sure each specific algorithm isn't broken.
SKIP: {
  skip "Not XS -- skipping direct primecount tests", 2 unless $isxs;
  # This has to be above SIEVE_LIMIT in lehmer.c and lmo.c or nothing happens.
  #is(Math::Prime::Util::_XS_lehmer_pi  (66123456),3903023,"XS Lehmer count");
  #is(Math::Prime::Util::_XS_meissel_pi (66123456),3903023,"XS Meissel count");
  #is(Math::Prime::Util::_XS_legendre_pi(66123456),3903023,"XS Legendre count");
  #is(Math::Prime::Util::_XS_LMOS_pi    (66123456),3903023,"XS LMOS count");
  is(Math::Prime::Util::_XS_LMO_pi     (66123456), 3903023,"XS LMO count");
  is(Math::Prime::Util::_XS_segment_pi (66123456), 3903023,"XS segment count");
}

require_ok 'Math::Prime::Util::PP';
is(Math::Prime::Util::PP::_lehmer_pi   (1456789), 111119, "PP Lehmer count");
is(Math::Prime::Util::PP::_sieve_prime_count(145678), 13478, "PP sieve count");
if ($extra) {
  is(Math::Prime::Util::PP::_lehmer_pi   (3456789), 247352, "PP Lehmer count");
  is(Math::Prime::Util::PP::_sieve_prime_count(3456789), 247352, "PP sieve count");
}

####### Twin prime counts
is(twin_prime_count(13,31), 2, "twin prime count 13 to 31");
is(twin_prime_count(10**8,10**8+34587), 137, "twin prime count 10^8 to +34587");
is(twin_prime_count(654321), 5744, "twin prime count 654321");
if ($isxs && $use64) {
  is(twin_prime_count(1000000000123456), 1177209242446, "twin prime count 1000000000123456");
  while (my($n, $tpc) = each (%tpcs)) {
    is(twin_prime_count($n), $tpc, "twin prime count $n");
    my $errorp = 100 * abs($tpc - twin_prime_count_approx($n)) / $tpc;
    my $estr = sprintf "%8.6f%%", $errorp;
    cmp_ok( $errorp, '<=', 2, "twin_prime_count_approx($n) is $estr");
  }
}

####### Ramanujan prime counts
is(ramanujan_prime_count(13,31), 2, "Ramanujan prime count 13 to 31");
is(ramanujan_prime_count(1357), 94, "Ramanujan prime count 1357");
if ($isxs && $use64) {
  is(ramanujan_prime_count(10**8,10**8+34587), 927, "Ramanujan prime count 10^8 to +34587");
  is(ramanujan_prime_count(654321), 24973, "Ramanujan prime count 654321");
  while (my($n, $rpc) = each (%rpcs)) {
    is(ramanujan_prime_count($n), $rpc, "Ramanujan prime count $n");
    #my $errorp = 100 * abs($tpc - ramanujan_prime_count_approx($n)) / $tpc;
    #my $estr = sprintf "%8.6f%%", $errorp;
    #cmp_ok( $errorp, '<=', 2, "ramanujan_prime_count_approx($n) is $estr");
  }
}