1#!/usr/bin/env perl
2use strict;
3use warnings;
4
5use Test::More;
6use Math::Prime::Util qw/primes nth_prime nth_twin_prime
7                         nth_prime_lower nth_prime_upper
8                         nth_prime_approx nth_twin_prime_approx
9                         nth_semiprime is_semiprime
10                         inverse_li/;
11
12my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
13my $usexs = Math::Prime::Util::prime_get_config->{'xs'};
14my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
15my $broken64 = (18446744073709550592 == ~0);
16
17my $nsmallprimes = 1000;
18my $nth_small_prime = 7919;  # nth_prime(1000)
19
20my %pivals32 = (
21                  1 => 0,
22                 10 => 4,
23                100 => 25,
24               1000 => 168,
25              10000 => 1229,
26             100000 => 9592,
27            1000000 => 78498,
28);
29
30#  Powers of 10: http://oeis.org/A006988/b006988.txt
31#  Powers of  2: http://oeis.org/A033844/b033844.txt
32my %nthprimes32 = (
33                  1 => 2,
34                 10 => 29,
35                100 => 541,
36               1000 => 7919,
37              10000 => 104729,
38             100000 => 1299709,
39            1000000 => 15485863,
40           10000000 => 179424673,
41          100000000 => 2038074743,
42  # Some values that estimate right around the value
43            6305537 => 110040379,
44            6305538 => 110040383,
45            6305539 => 110040391,
46            6305540 => 110040407,
47            6305541 => 110040467,
48            6305542 => 110040499,
49            6305543 => 110040503,
50);
51my %nthprimes64 = (
52         1000000000 => 22801763489,
53        10000000000 => 252097800623,
54       100000000000 => 2760727302517,
55      1000000000000 => 29996224275833,
56     10000000000000 => 323780508946331,
57    100000000000000 => 3475385758524527,
58   1000000000000000 => 37124508045065437,
59  10000000000000000 => 394906913903735329,
60 100000000000000000 => 4185296581467695669,
61);
62my %nthprimes_small = map { $_ => $nthprimes32{$_} }
63                      grep { ($_ <= 10_000_000) || $extra }
64                      keys %nthprimes32;
65
66my @small_primes = (undef, @{primes($nth_small_prime)});
67
68my %ntpcs = (
69             5 =>                   29,
70            50 =>                 1487,
71           500 =>                32411,
72          5000 =>               557519,
73         50000 =>              8264957,
74        500000 =>            115438667,
75       5000000 =>           1523975909,
76      50000000 =>          19358093939,
77     500000000 =>         239211160649,
78);
79
80my %nthsemi = (
81          1234 =>          4497,
82         12345 =>         51019,
83        123456 =>        573355,
84       1234567 =>       6365389,
85);
86$nthsemi{12345678}    =     69914722  if $usexs || $extra;
87$nthsemi{123456789}   =    760797011  if $usexs && $extra;
88$nthsemi{1234567890}  =   8214915893  if $usexs && $extra && $use64;
89$nthsemi{8589934592}  =  60662588879  if $usexs && $extra && $use64;
90$nthsemi{17179869184} = 123806899739  if $usexs && $extra && $use64;
91
92plan tests => 0 + 2*scalar(keys %pivals32)
93                + 1
94                + 3*scalar(keys %nthprimes32)
95                + scalar(keys %nthprimes_small)
96                + $use64 * 3 * scalar(keys %nthprimes64)
97                + 3   # nth_prime_lower with max index
98                + 3   # nth_twin_prime
99                + 3   # inverse_li
100                + scalar(keys %ntpcs)   # nth_twin_prime_approx
101                + 2 + scalar(keys %nthsemi)   # nth_semiprime
102                + (($extra && $use64 && $usexs) ? 1 : 0);
103
104
105while (my($n, $pin) = each (%pivals32)) {
106  my $next = $pin+1;
107  cmp_ok( $pin ? nth_prime($pin) : 0, '<=', $n, "nth_prime($pin) <= $n");
108  cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
109}
110
111{
112  my @nth_primes = map { nth_prime($_) }  (0 .. $#small_primes);
113  is_deeply( \@nth_primes, \@small_primes, "nth_prime for primes 0 .. $#small_primes" );
114}
115
116while (my($n, $nth) = each (%nthprimes32)) {
117  cmp_ok( nth_prime_upper($n), '>=', $nth, "nth_prime($n) <= upper estimate" );
118  cmp_ok( nth_prime_lower($n), '<=', $nth, "nth_prime($n) >= lower estimate" );
119
120  my $approx = nth_prime_approx($n);
121  my $percent_limit = ($n >= 775) ? 1 : 2;
122  cmp_ok( abs($nth - $approx) / $nth, '<=', $percent_limit/100.0, "nth_prime_approx($n) = $approx within $percent_limit\% of $nth");
123}
124while (my($n, $nth) = each (%nthprimes_small)) {
125  is( nth_prime($n), $nth, "nth_prime($n) = $nth" );
126}
127
128if ($use64) {
129  while (my($n, $nth) = each (%nthprimes64)) {
130    cmp_ok( nth_prime_upper($n), '>=', $nth, "nth_prime($n) <= upper estimate" );
131    cmp_ok( nth_prime_lower($n), '<=', $nth, "nth_prime($n) >= lower estimate" );
132
133    my $approx = "" . nth_prime_approx($n);   # ensure not a bigint
134    my $percent_limit = 0.001;
135    cmp_ok( abs($nth - $approx) / $nth, '<=', $percent_limit/100.0, "nth_prime_approx($n) = $approx within $percent_limit\% of $nth");
136  }
137}
138
139my $maxindex   = $use64 ? '425656284035217743' : '203280221';
140my $maxindexp1 = $use64 ? '425656284035217744' : '203280222';
141my $maxprime   = $use64 ? '18446744073709551557' : '4294967291';
142cmp_ok( nth_prime_lower($maxindex), '<=', $maxprime, "nth_prime_lower(maxindex) <= maxprime");
143cmp_ok( nth_prime_upper($maxindex), '>=', $maxprime, "nth_prime_upper(maxindex) >= maxprime");
144cmp_ok( nth_prime_lower($maxindexp1), '>=', nth_prime_lower($maxindex), "nth_prime_lower(maxindex+1) >= nth_prime_lower(maxindex)");
145
146my $overindex = ($broken64) ? 425656284035217843 : $maxindexp1;
147
148if ($extra && $use64 && $usexs) {
149  # Test an nth prime value that uses the binary-search-on-R(n) algorithm
150  is( nth_prime(21234567890), 551990503367, "nth_prime(21234567890)" );
151}
152
153####################################3
154
155is( nth_twin_prime(0), undef, "nth_twin_prime(0) = undef" );
156is( nth_twin_prime(17), 239, "239 = 17th twin prime" );
157is( nth_twin_prime(1234), 101207, "101207 = 1234'th twin prime" );
158
159while (my($n, $nthtpc) = each (%ntpcs)) {
160  my $approx = nth_twin_prime_approx($n);
161  my $errorp = 100 * abs($nthtpc - $approx) / $nthtpc;
162  my $estr = sprintf "%8.6f%%", $errorp;
163  cmp_ok( $errorp, '<=', 2, "nth_twin_prime_approx($n) is $estr (got $approx, expected ~$nthtpc)");
164}
165
166####################################3
167
168is( nth_semiprime(0), undef, "nth_semiprime(0) = undef" );
169{
170  my $range = $extra ? 10000 : 500;
171  my @semiprimes = grep { is_semiprime($_) } 0 .. $range;
172  my $nsmall = scalar(@semiprimes);
173  my @nth_semis = map { nth_semiprime($_) } 1 .. $nsmall;
174  is_deeply(\@nth_semis, \@semiprimes, "nth_semiprime(1 .. $nsmall)");
175}
176while (my($n, $nthsemi) = each (%nthsemi)) {
177  is( nth_semiprime($n), $nthsemi, "nth_semiprime($n) = $nthsemi" );
178}
179
180####################################3
181
182is_deeply(
183  [ map { inverse_li($_) } 0 .. 50 ],
184  [qw/0 2 3 5 6 8 10 12 15 18 21 24 27 30 34 37 41 45 49 53 57 61 65 69 73 78 82 86 91 95 100 105 109 114 119 123 128 133 138 143 148 153 158 163 168 173 179 184 189 194 199/],
185  "inverse_li: Li^-1(0..50)"
186);
187# Allow +/- 2 for floating point differences in LogarithmicIntegral
188like(inverse_li(1000000000), qr/^2280162741[34567]$/, "inverse_li(1e9)");
189like(inverse_li(1100000000000), qr/^3310443690704[01234]$/, "inverse_li(11e11)");
190