#!/usr/bin/env perl
use strict;
use warnings;
use Math::Prime::Util qw/prime_count prime_count_approx prime_count_lower prime_count_upper LogarithmicIntegral RiemannR/;
use Math::BigFloat;
$| = 1; # fast pipes
my %pivals = (
10 => 4,
100 => 25,
1000 => 168,
10000 => 1229,
100000 => 9592,
1000000 => 78498,
10000000 => 664579,
100000000 => 5761455,
1000000000 => 50847534,
10000000000 => 455052511,
100000000000 => 4118054813,
1000000000000 => 37607912018,
10000000000000 => 346065536839,
100000000000000 => 3204941750802,
'1000000000000000' => 29844570422669,
'10000000000000000' => 279238341033925,
'100000000000000000' => 2623557157654233,
'1000000000000000000' => 24739954287740860,
'10000000000000000000' => 234057667276344607,
);
printf(" N %12s %12s %12s %12s\n", "pc_approx", "Li", "LiCor", "R");
printf("----- %12s %12s %12s %12s\n", '-'x12,'-'x12,'-'x12,'-'x12);
foreach my $n (sort {$a<=>$b} keys %pivals) {
my $pin = $pivals{$n};
my $pca = prime_count_approx($n);
my $Lisub = sub { my $x = shift; return ($x < 2) ? 0 : (LogarithmicIntegral($x)-LogarithmicIntegral(2)+0.5); };
my $pcli = int($Lisub->($n));
my $pclicor = int( $Lisub->($n) - ($Lisub->(sqrt($n)) / 2) );
my $r = int(RiemannR($n)+0.5);
printf "10^%2d %12d %12d %12d %12d\n", length($n)-1,
abs($pca-$pin), abs($pcli-$pin), abs($pclicor-$pin), abs($r-$pin);
}
# Also see http://empslocal.ex.ac.uk/people/staff/mrwatkin/zeta/encoding1.htm
# for some ideas one how this could be made even more accurate.
print "\n";
print "Lower / Upper bounds. Percentages.\n";
print "\n";
printf(" N %12s %12s %12s %12s\n", "lower", "upper", "SchoenfeldL", "SchoenfeldU");
printf("----- %12s %12s %12s %12s\n", '-'x12,'-'x12,'-'x12,'-'x12);
foreach my $n (sort {$a<=>$b} keys %pivals) {
my ($pin, $pcl, $pcu, $scl, $scu) =
map { Math::BigFloat->new($_) }
($pivals{$n}, prime_count_lower($n), prime_count_upper($n), stoll($n));
#printf "10^%2d %12d %12d\n", length($n)-1, $pin-$pcl, $pcu-$pin;
printf "10^%2d %12.7f %12.7f %12.7f %12.7f\n",
length($n)-1, 100*($pin-$pcl)/$pin, 100*($pcu-$pin)/$pin,
100*($pin-$scl)/$pin, 100*($scu-$pin)/$pin;
}
sub schoenfeld {
my $x = shift;
my $lix = LogarithmicIntegral($x);
my $bound = (sqrt($x)*log($x)) / 8*3.1415926535;
($lix-$bound,$lix+$bound);
}
# http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html
sub stoll {
my $x = shift;
my $lix = LogarithmicIntegral($x);
my $bound = sqrt($x) * (log(log(log($x))) + exp(1) + 1) / (exp(1)*log($x));
($lix-$bound,$lix+$bound);
}