1#!/usr/bin/env perl 2use strict; 3use warnings; 4 5# This is a subset of our tests. You really should run the whole test suite 6# on the PP code. What this will do is basic regression testing. 7my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING}; 8my $use64 = ~0 > 4294967295 && ~0 != 18446744073709550592; 9 10use Test::More; 11my @small_primes = qw/ 122 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 1373 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 14179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 15283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 16419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 17547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 18661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 19811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 20947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 1051 211061 1063 1069 22/; # next prime is 1087 23 24my @primes = qw/ 251129 1327 9551 15683 19609 31397 155921 265 11 29 97 127 541 907 1151 1361 9587 15727 19661 31469 156007 360749 27370373 492227 1349651 1357333 2010881 4652507 17051887 20831533 47326913 28122164969 189695893 191913031 29/; 30 31my @composites = qw/ 320 4 6 8 9 10 12 14 15 16 18 20 21 22 339 2047 1373653 25326001 3215031751 34561 1105 1729 2465 2821 6601 8911 10585 15841 29341 41041 46657 52633 3562745 63973 75361 101101 340561 488881 852841 1857241 6733693 369439201 17236801 23382529 34657141 56052361 146843929 37341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 384681 5461 6601 7957 8321 52633 88357 3966066 173645446 7500135 115501463 40/; 41 42# pseudoprimes to various small prime bases 43my %pseudoprimes = ( 44 2 => [ qw/2047 42799 4335241 1078467589/ ], 45 3 => [ qw/121 44287 4252381 1075490821/ ], 46 5 => [ qw/781 38081 4265257 1075156291/ ], 47 7 => [ qw/25 325 29857 4411681 1074439981/ ], 48 11 => [ qw/133 43213 4224533 1076929261/ ], 49 13 => [ qw/85 35371 4336879 1079159203/ ], 50 17 => [ qw/9 91 71071 4224533 1076237119/ ], 51 19 => [ qw/9 49 49771 4384693 1074718783/ ], 52 23 => [ qw/169 25201 4219129 1079063371/ ], 53 29 => [ qw/15 91 48133 4219129 1075151447/ ], 54 31 => [ qw/15 49 29341 4270657 1073833843/ ], 55 37 => [ qw/9 451 59563 4287817 1075430539/ ], 56 61 => [ qw/217 79381 4219129 1079326249/ ], 57 73 => [ qw/205 34219 4321153 1074220489/ ], 58 psp2 => [ qw/341 561 29341 4259905 1073823745/ ], 59 psp3 => [ qw/91 121 44287 4252381 1073827147/ ], 60 lucas => [ qw/9179 10877 44099 4259789 1074039119/ ], 61 slucas => [ qw/5459 5777 75077 4309631 1080085439/ ], 62 eslucas => [ qw/989 3239 5777 72389 4226777 1076503199/ ], 63 aeslucas1 => [ qw/989 10469 39059 4269341 1076503199/ ], 64 aeslucas2 => [ qw/4531 12209 62479 4403027 1074695441/ ], 65); 66# Test a pseudoprime larger than 2^32. 67push @{$pseudoprimes{2}}, 75792980677 if $use64; 68my $num_pseudoprimes = 0; 69foreach my $ppref (values %pseudoprimes) { 70 push @composites, @$ppref; 71 $num_pseudoprimes += scalar @$ppref; 72} 73{ 74 my %uniq; 75 $uniq{$_}++ for (@composites); 76 @composites = sort {$a<=>$b} keys %uniq; 77} 78 79my %small_single = ( 80 0 => [], 81 1 => [], 82 2 => [2], 83 3 => [2, 3], 84 4 => [2, 3], 85 5 => [2, 3, 5], 86 6 => [2, 3, 5], 87 7 => [2, 3, 5, 7], 88 11 => [2, 3, 5, 7, 11], 89 18 => [2, 3, 5, 7, 11, 13, 17], 90 19 => [2, 3, 5, 7, 11, 13, 17, 19], 91 20 => [2, 3, 5, 7, 11, 13, 17, 19], 92); 93 94my %small_range = ( 95 "3 to 9" => [3,5,7], 96 "2 to 20" => [2,3,5,7,11,13,17,19], 97 "30 to 70" => [31,37,41,43,47,53,59,61,67], 98 "70 to 30" => [], 99 "20 to 2" => [], 100 "1 to 1" => [], 101 "2 to 2" => [2], 102 "3 to 3" => [3], 103 "2 to 3" => [2,3], 104 "2 to 5" => [2,3,5], 105 "3 to 6" => [3,5], 106 "3 to 7" => [3,5,7], 107 "4 to 8" => [5,7], 108 "2010733 to 2010881" => [2010733,2010881], 109 "2010734 to 2010880" => [], 110 "3088 to 3164" => [3089,3109,3119,3121,3137,3163], 111 "3089 to 3163" => [3089,3109,3119,3121,3137,3163], 112 "3090 to 3162" => [3109,3119,3121,3137], 113 "3842610773 to 3842611109" => [3842610773,3842611109], 114 "3842610774 to 3842611108" => [], 115); 116 117my %primegaps = ( 118 19609 => 52, 119 360653 => 96, 120 2010733 => 148, 121); 122 123my %pivals32 = ( 124 1 => 0, 125 10 => 4, 126 100 => 25, 127 1000 => 168, 128 10000 => 1229, 129 100000 => 9592, 130 1000000 => 78498, 131 10000000 => 664579, 132 100000000 => 5761455, 133 1000000000 => 50847534, 134 60067 => 6062, 135 65535 => 6542, 136 16777215 => 1077871, 137 2147483647 => 105097565, 138 4294967295 => 203280221, 139); 140my %pivals_small = map { $_ => $pivals32{$_} } 141 grep {$_ <= 80000} 142 keys %pivals32; 143 144my %pi_intervals = ( 145 "1e9 +2**14" => 785, 146 "17 to 13" => 0, 147 "3 to 17" => 6, 148 "4 to 17" => 5, 149 "4 to 16" => 4, 150 "191912783 +248" => 2, 151 "191912784 +247" => 1, 152 "191912783 +247" => 1, 153 "191912784 +246" => 0, 154); 155my %extra_pi_intervals = ( 156 "868396 to 9478505" => 563275, 157 "1118105 to 9961674" => 575195, 158 "24689 to 7973249" => 535368, 159); 160# Add extra intervals to pi_intervals if we're doing release testing 161@pi_intervals{keys %extra_pi_intervals} = values %extra_pi_intervals if $extra; 162 163# Remove any entries where the high value is too large for us 164# ikegami++ for the delete from a hash slice idea 165delete @pi_intervals{ grep { (parse_range($_))[1] > ~0 } keys %pi_intervals }; 166 167my %nthprimes32 = ( 168 1 => 2, 169 10 => 29, 170 100 => 541, 171 1000 => 7919, 172 10000 => 104729, 173 100000 => 1299709, 174 1000000 => 15485863, 175 10000000 => 179424673, 176 100000000 => 2038074743, 177); 178my %nthprimes_small = map { $_ => $nthprimes32{$_} } 179 grep { $extra ? ($_ <= 2_000_000) : ($_ <= 5_000) } 180 keys %nthprimes32; 181 182my %eivals = ( 183 -10 => -0.00000415696892968532438, 184 -0.5 => -0.55977359477616, 185 -0.1 => -1.8229239584193906660809, 186 -0.001 => -6.33153936413615, 187 -0.00001 => -10.9357198000436956, 188 -0.00000001 => -17.843465089050832587, 189 0.693147180559945 => 1.0451637801174927848446, # log2 190 1 => 1.8951178163559367554665, 191 1.5 => 3.3012854491297978379574, 192 2 => 4.9542343560018901633795, 193 5 => 40.185275355803177455091, 194 10 => 2492.2289762418777591384, 195 12 => 14959.532666397528852292, 196 20 => 25615652.664056588820481, 197 40 => 6039718263611241.5783592, 198 41 => 16006649143245041.110700, 199); 200my %livals = ( 201 0 => 0, 202 1.01 => -4.0229586739299358695031, 203 2 => 1.0451637801174927848446, 204 10 => 6.1655995047872979375230, 205 24 => 11.200315795232698830550, 206 1000 => 177.60965799015222668764, 207 100000 => 9629.8090010507982050343, 208 100000000 => 5762209.3754480314675691, 209 4294967295 => 203284081.95454158906409, 210 10000000000 => 455055614.58662307560953, 211 100000000000 => 4118066400.6216115150394, 212); 213my %rvals = ( 214 1.01 => 1.0060697180622924796117, 215 2 => 1.5410090161871318832885, 216 10 => 4.5645831410050902398658, 217 1000 => 168.35944628116734806491, 218 1000000 => 78527.399429127704858870, 219 10000000 => 664667.44756474776798535, 220 4294967295 => 203280697.51326064541983, 221 10000000000 => 455050683.30684692446315, 22218446744073709551615 => 4.25656284014012122706963685602e17, 223); 224my %rzvals = ( 225 2 => 0.6449340668482264364724151666, 226 2.5 => 0.3414872572509171797567696934, 227 4.5 => 0.0547075107614542640229672890, 228 7 => 0.0083492773819228268397975498, 229 8.5 => 0.0028592508824156277133439825, 230 20.6 => 0.0000006293391573578212882457, 231 80 => 8.27180612553034e-25, 232 180 => 6.52530446799852e-55, 233); 234my %ipp = ( 235 5 => 2, 236 10 => 0, 237 49 => 0, 238 347 => 2, 239 697 => 0, 240 7080233 => 2, 241 7080249 => 0, 242 17471059 => 2, 243 17471061 => 0, 244 36010357 => 2, 245 36010359 => 0, 246); 247 248plan tests => 2 + 249 3 + 250 3 + scalar(keys %small_single) + scalar(keys %small_range) + 251 2*scalar(keys %primegaps) + 8 + 1 + 1 + 1 + 252 scalar(keys %pivals_small) + scalar(keys %pi_intervals) + 253 6 + # PC, pc approx 254 2*scalar(keys %pivals_small) + scalar(keys %nthprimes_small) + 255 4 + scalar(keys %pseudoprimes) + 256 scalar(keys %eivals) + scalar(keys %livals) + scalar(keys %rvals) + scalar(keys %rzvals) + 257 ($extra ? 4 : 0) + # Bigfloat RiemannZeta 258 1 + 1 + # factor 259 10 + 7*3 + # factoring subs 260 1 + # HOLF 261 ($extra ? 3 : 0) + # HOLF extra 262 ($extra ? 3 : 0) + # factor stage 2 263 10 + # AKS 264 ($use64 ? 3 : 2) + # Lucas and BLS75 primality proofs 265 6 + # M-R and Lucas on bigint 266 2 + # PC and NP approx 267 65 + # Misc util.pm functions 268 ($extra ? 1 : 0) + # twin prime count approx 269 scalar(keys %ipp) + # is_prob_prime 270 1; 271 272use Math::Prime::Util qw/primes 273 prime_count_approx nth_prime_approx 274 prime_get_config prime_set_config 275 consecutive_integer_lcm 276 primorial pn_primorial partitions miller_rabin_random 277 is_prob_prime 278 /; 279use Math::BigInt; 280use Math::BigFloat; 281require_ok 'Math::Prime::Util::PP'; 282require_ok 'Math::Prime::Util::PrimalityProving'; 283 284 # This function skips some setup 285 undef *primes; 286 *primes = \&Math::Prime::Util::PP::primes; 287 288 *prime_count = \&Math::Prime::Util::PP::prime_count; 289 *prime_count_lower = \&Math::Prime::Util::PP::prime_count_lower; 290 *prime_count_upper = \&Math::Prime::Util::PP::prime_count_upper; 291 *nth_prime = \&Math::Prime::Util::PP::nth_prime; 292 undef *prime_count_approx; 293 undef *nth_prime_approx; 294 *prime_count_approx = \&Math::Prime::Util::PP::prime_count_approx; 295 *nth_prime_approx = \&Math::Prime::Util::PP::nth_prime_approx; 296 297 *twin_prime_count = \&Math::Prime::Util::PP::twin_prime_count; 298 *nth_twin_prime = \&Math::Prime::Util::PP::nth_twin_prime; 299 *twin_prime_count_approx = \&Math::Prime::Util::PP::twin_prime_count_approx; 300 301 *is_prime = \&Math::Prime::Util::PP::is_prime; 302 *next_prime = \&Math::Prime::Util::PP::next_prime; 303 *prev_prime = \&Math::Prime::Util::PP::prev_prime; 304 305 *is_pseudoprime = \&Math::Prime::Util::PP::is_pseudoprime; 306 *is_strong_pseudoprime = \&Math::Prime::Util::PP::is_strong_pseudoprime; 307 *is_lucas_pseudoprime = \&Math::Prime::Util::PP::is_lucas_pseudoprime; 308 *is_strong_lucas_pseudoprime = \&Math::Prime::Util::PP::is_strong_lucas_pseudoprime; 309 *is_extra_strong_lucas_pseudoprime = \&Math::Prime::Util::PP::is_extra_strong_lucas_pseudoprime; 310 *is_almost_extra_strong_lucas_pseudoprime = \&Math::Prime::Util::PP::is_almost_extra_strong_lucas_pseudoprime; 311 *is_frobenius_underwood_pseudoprime = \&Math::Prime::Util::PP::is_frobenius_underwood_pseudoprime; 312 *is_perrin_pseudoprime = \&Math::Prime::Util::PP::is_perrin_pseudoprime; 313 *is_frobenius_pseudoprime = \&Math::Prime::Util::PP::is_frobenius_pseudoprime; 314 *is_aks_prime = \&Math::Prime::Util::PP::is_aks_prime; 315 316 *factor = \&Math::Prime::Util::PP::factor; 317 318 *gcd = \&Math::Prime::Util::PP::gcd; 319 *lcm = \&Math::Prime::Util::PP::lcm; 320 321 *moebius = \&Math::Prime::Util::PP::moebius; 322 *euler_phi = \&Math::Prime::Util::PP::euler_phi; 323 *jordan_totient = \&Math::Prime::Util::PP::jordan_totient; 324 *mertens = \&Math::Prime::Util::PP::mertens; 325 *exp_mangoldt = \&Math::Prime::Util::PP::exp_mangoldt; 326 *chebyshev_theta= \&Math::Prime::Util::PP::chebyshev_theta; 327 *chebyshev_psi = \&Math::Prime::Util::PP::chebyshev_psi; 328 329 *znprimroot = \&Math::Prime::Util::PP::znprimroot; 330 *znorder = \&Math::Prime::Util::PP::znorder; 331 *znlog = \&Math::Prime::Util::PP::znlog; 332 *binomial = \&Math::Prime::Util::PP::binomial; 333 *stirling = \&Math::Prime::Util::PP::stirling; 334 *bernfrac = \&Math::Prime::Util::PP::bernfrac; 335 *valuation = \&Math::Prime::Util::PP::valuation; 336 *gcdext = \&Math::Prime::Util::PP::gcdext; 337 *invmod = \&Math::Prime::Util::PP::invmod; 338 *vecmin = \&Math::Prime::Util::PP::vecmin; 339 *vecmax = \&Math::Prime::Util::PP::vecmax; 340 *vecsum = \&Math::Prime::Util::PP::vecsum; 341 *vecprod = \&Math::Prime::Util::PP::vecprod; 342 *liouville = \&Math::Prime::Util::PP::liouville; 343 *carmichael_lambda = \&Math::Prime::Util::PP::carmichael_lambda; 344 *forperm = \&Math::Prime::Util::PP::forperm; 345 *forcomb = \&Math::Prime::Util::PP::forcomb; 346 *forpart = \&Math::Prime::Util::PP::forpart; 347 *Pi = \&Math::Prime::Util::PP::Pi; 348 349 *RiemannR = \&Math::Prime::Util::PP::RiemannR; 350 *RiemannZeta = \&Math::Prime::Util::PP::RiemannZeta; 351 *LogarithmicIntegral = \&Math::Prime::Util::PP::LogarithmicIntegral; 352 *ExponentialIntegral = \&Math::Prime::Util::PP::ExponentialIntegral; 353 *LambertW = \&Math::Prime::Util::PP::LambertW; 354 355# Turn off use of BRS - ECM tries to use this. 356# prime_set_config( irand => sub { int(rand(4294967296.0)) } ); 357 358############################################################################### 359 360$_ = 'this should not change'; 361 362{ 363 my %small_primes = map { $_ => 1 } @small_primes; 364 my @isprime = map { is_prime($_) } (0 .. 1086); 365 my @exprime = map { $small_primes{$_} ? 2 : 0 } (0 .. 1086); 366 is_deeply( \@isprime, \@exprime, "is_prime 0 .. 1086" ); 367} 368{ 369 my @isprime = map { is_prime($_) ? "$_ is prime" : "$_ is composite" } 370 @primes, @composites; 371 my @exprime = map { "$_ is prime" } @primes; 372 push @exprime, map { "$_ is composite" } @composites; 373 is_deeply( \@isprime, \@exprime, "is_prime for selected numbers" ); 374} 375 376is_deeply( Math::Prime::Util::PP::trial_primes(80), 377 [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79], 378 "Trial primes 2-80" ); 379 380############################################################################### 381 382is_deeply( primes(1069), \@small_primes, "Primes between 0 and 1069" ); 383is_deeply( primes(1070), \@small_primes, "Primes between 0 and 1070" ); 384is_deeply( primes(1086), \@small_primes, "Primes between 0 and 1086" ); 385 386while (my($high, $expect) = each (%small_single)) { 387 is_deeply( primes($high), $expect, "primes($high) should return [@{$expect}]") 388; 389} 390 391while (my($range, $expect) = each (%small_range)) { 392 my($low,$high) = $range =~ /(\d+) to (\d+)/; 393 is_deeply( primes($low, $high), $expect, "primes($low,$high) should return [@{$expect}]"); 394} 395 396############################################################################### 397 398while (my($base, $range) = each (%primegaps)) { 399 is( next_prime($base), $base+$range, "next prime of $base is $base+$range" ); 400 is( prev_prime($base+$range), $base, "prev prime of $base+$range is $base" ); 401} 402 403is( next_prime(19608), 19609, "next prime of 19608 is 19609" ); 404is( next_prime(19610), 19661, "next prime of 19610 is 19661" ); 405is( next_prime(19660), 19661, "next prime of 19660 is 19661" ); 406is( prev_prime(19662), 19661, "prev prime of 19662 is 19661" ); 407is( prev_prime(19660), 19609, "prev prime of 19660 is 19609" ); 408is( prev_prime(19610), 19609, "prev prime of 19610 is 19609" ); 409 410is( prev_prime(2), undef, "Previous prime of 2 returns undef" ); 411if ($use64) { 412 is( next_prime(18446744073709551611), "18446744073709551629", "Next prime of ~0-4 returns bigint next prime" ); 413} else { 414 is( next_prime(4294967291), "4294967311", "Next prime of ~0-4 returns bigint next prime" ); 415} 416 417{ 418 my @exprime = map { "next_prime($_) == 2010881" } (2010733..2010880); 419 my @isprime = map { "next_prime($_) == ".next_prime($_) } (2010733..2010880); 420 is_deeply(\@isprime, \@exprime, "next_prime for 148 primes before primegap end 2010881"); 421} 422{ 423 my @exprime = map { "prev_prime($_) == 2010733" } (2010734..2010881); 424 my @isprime = map { "prev_prime($_) == ".prev_prime($_) } (2010734..2010881); 425 is_deeply(\@isprime, \@exprime, "prev_prime for 148 primes before primegap start 2010733"); 426} 427# Similar test case to 2010870, where m=0 and next_prime is at m=1 428is(next_prime(1234567890), 1234567891, "next_prime(1234567890) == 1234567891)"); 429 430############################################################################### 431 432while (my($n, $pin) = each (%pivals_small)) { 433 is( prime_count($n), $pin, "Pi($n) = $pin" ); 434} 435 436while (my($range, $expect) = each (%pi_intervals)) { 437 my($low,$high) = parse_range($range); 438 is( prime_count($low,$high), $expect, "prime_count($range) = $expect"); 439} 440 441# These are small enough they should be exact. 442is( prime_count_lower(450), 87, "prime_count_lower(450)" ); 443is( prime_count_upper(450), 87, "prime_count_upper(450)" ); 444# Make sure these are about right 445cmp_closeto( prime_count_lower(1234567), 95360, 60, "prime_count_lower(1234567) in range" ); 446cmp_closeto( prime_count_upper(1234567), 95360, 60, "prime_count_upper(1234567) in range" ); 447cmp_closeto( prime_count_lower(412345678), 21958997, 1500, "prime_count_lower(412345678) in range" ); 448cmp_closeto( prime_count_upper(412345678), 21958997, 1500, "prime_count_upper(412345678) in range" ); 449 450############################################################################### 451 452while (my($n, $pin) = each (%pivals_small)) { 453 my $next = $pin+1; 454 cmp_ok( $pin ? nth_prime($pin) : 0, '<=', $n, "nth_prime($pin) <= $n"); 455 cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n"); 456} 457 458while (my($n, $nth) = each (%nthprimes_small)) { 459 is( nth_prime($n), $nth, "nth_prime($n) = $nth" ); 460} 461 462############################################################################### 463 464is( is_strong_pseudoprime(0, 2), 0, "MR with 0 shortcut composite"); 465is( is_strong_pseudoprime(1, 2), 0, "MR with 0 shortcut composite"); 466is( is_strong_pseudoprime(2, 2), 1, "MR with 2 shortcut prime"); 467is( is_strong_pseudoprime(3, 2), 1, "MR with 3 shortcut prime"); 468 469while (my($base, $ppref) = each (%pseudoprimes)) { 470 my $npseudos = scalar @$ppref; 471 my @expmr = map { 1 } @$ppref; 472 my @gotmr; 473 if ($base =~ /^psp(\d+)/) { 474 my $pbase = $1; 475 @gotmr = map { is_pseudoprime($_, $pbase) } @$ppref; 476 } elsif ($base =~ /^aeslucas(\d+)/) { 477 my $inc = $1; 478 @gotmr = map { is_almost_extra_strong_lucas_pseudoprime($_, $inc) } @$ppref; 479 } elsif ($base eq 'eslucas') { 480 @gotmr = map { is_extra_strong_lucas_pseudoprime($_) } @$ppref; 481 } elsif ($base eq 'slucas') { 482 @gotmr = map { is_strong_lucas_pseudoprime($_) } @$ppref; 483 } elsif ($base eq 'lucas') { 484 @gotmr = map { is_lucas_pseudoprime($_) } @$ppref; 485 } else { 486 @gotmr = map { is_strong_pseudoprime($_, $base) } @$ppref; 487 } 488 is_deeply(\@gotmr, \@expmr, "$npseudos pseudoprimes (base $base)"); 489} 490 491############################################################################### 492 493while (my($n, $ein) = each (%eivals)) { 494 cmp_closeto( ExponentialIntegral($n), $ein, 0.00000001 * abs($ein), "Ei($n) ~= $ein"); 495} 496while (my($n, $lin) = each (%livals)) { 497 cmp_closeto( LogarithmicIntegral($n), $lin, 0.00000001 * abs($lin), "li($n) ~= $lin"); 498} 499while (my($n, $rin) = each (%rvals)) { 500 cmp_closeto( RiemannR($n), $rin, 0.00000001 * abs($rin), "R($n) ~= $rin"); 501} 502while (my($n, $zin) = each (%rzvals)) { 503 cmp_closeto( RiemannZeta($n), $zin, 0.00000001 * abs($zin), "Zeta($n) ~= $zin"); 504} 505cmp_closeto( LambertW(6588), 6.86636957140619, 0.000000001, "LambertW(6588)"); 506if ($extra) { 507 my ($n, $zin); 508 ($n, $zin) = (4.5, $rzvals{4.5}); 509 cmp_closeto( RiemannZeta(Math::BigFloat->new($n)), $zin, 0.00000001 * abs($zin), "Zeta($n) ~= $zin"); 510 ($n, $zin) = (20.6, $rzvals{20.6}); 511 cmp_closeto( RiemannZeta(Math::BigFloat->new($n)), $zin, 0.00000001 * abs($zin), "Zeta($n) ~= $zin"); 512 ($n, $zin) = (80, $rzvals{80}); 513 cmp_closeto( RiemannZeta(Math::BigFloat->new($n)), $zin, 0.00000001 * abs($zin), "Zeta($n) ~= $zin"); 514 ($n, $zin) = (180, $rzvals{180}); 515 cmp_closeto( RiemannZeta(Math::BigFloat->new($n)), $zin, 0.00000001 * abs($zin), "Zeta($n) ~= $zin"); 516} 517 518############################################################################### 519 520#foreach my $n (@primes) { 521# my @f = factor($n); 522# is_deeply( \@f, [$n], "factor prime $n yields $n" ); 523#} 524{ 525 my $ntests = scalar @primes; 526 my @expfactor = map { "$_" } @primes; 527 my @gotfactor = map { join(' * ', factor($_)) } @primes; 528 is_deeply( \@gotfactor, \@expfactor, "test factoring for $ntests primes"); 529} 530{ 531 my $ntests = scalar @composites; 532 my @expfactor = map { "$_ factored correctly" } @composites; 533 my @gotfactor; 534 535 foreach my $n (@composites) { 536 my @f = factor($n); 537 my $facstring = join(' * ', @f); 538 539 if ($n < 2) { 540 push @gotfactor, (@f == 1 && $f[0] == $n) 541 ? "$n factored correctly" 542 : "$n not correct: $facstring"; 543 next; 544 } 545 my $product = 1; $product = int($product * $_) for @f; 546 my $allprime = 1; $allprime *= is_prime($_) for @f; 547 if (@f >= 2 && $product == $n && $allprime) { 548 push @gotfactor, "$n factored correctly"; 549 } else { 550 push @gotfactor, "$n not correct: $facstring"; 551 } 552 } 553 is_deeply( \@gotfactor, \@expfactor, "test factoring for $ntests composites"); 554} 555 556# The PP factor code does small trials, then loops doing 64k rounds of HOLF 557# if the composite is less than a half word, followed by 64k rounds each of 558# prho with a = {3,5,7,11,13}. Most numbers are handled by these. The ones 559# that aren't end up being too slow for us to put in a test. So we'll try 560# running the various factoring methods manually. 561{ 562 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::holf_factor(403) ], 563 [ 13, 31 ], 564 "holf(403)" ); 565 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::fermat_factor(403) ], 566 [ 13, 31 ], 567 "fermat(403)" ); 568 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::prho_factor(403) ], 569 [ 13, 31 ], 570 "prho(403)" ); 571 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor(403) ], 572 [ 13, 31 ], 573 "pbrent(403)" ); 574 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pminus1_factor(403) ], 575 [ 13, 31 ], 576 "pminus1(403)" ); 577 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::prho_factor(851981) ], 578 [ 13, 65537 ], 579 "prho(851981)" ); 580 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor(851981) ], 581 [ 13, 65537 ], 582 "pbrent(851981)" ); 583 #is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::ecm_factor(851981) ], 584 # [ 13, 65537 ], 585 # "ecm(851981)" ); 586 # Try to force using stage 2. 587 SKIP: { 588 skip "Skipping ecm stage 2 tests", 1 if defined $Math::Prime::Util::GMP::VERSION && $Math::Prime::Util::GMP::VERSION < 0.20; 589 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::ecm_factor(101303039, 5, 100000,100) ], 590 [ 1013, 100003 ], 591 "ecm(101303039)" ); 592 } 593 my $n64 = $use64 ? 55834573561 : Math::BigInt->new("55834573561"); 594 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::prho_factor($n64) ], 595 [ 13, 4294967197 ], 596 "prho(55834573561)" ); 597 is_deeply( [ sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor($n64) ], 598 [ 13, 4294967197 ], 599 "pbrent(55834573561)" ); 600 # 1013 4294967197 4294967291 601 my $nbig = Math::BigInt->new("18686551294184381720251"); 602 my @nfac; 603 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::prho_factor($nbig); 604 is(scalar @nfac, 2, "prho finds a factor of 18686551294184381720251"); 605 is($nfac[0] * $nfac[1], $nbig, "prho found a correct factor"); 606 ok($nfac[0] != 1 && $nfac[1] != 1, "prho didn't return a degenerate factor"); 607 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::pbrent_factor($nbig); 608 is(scalar @nfac, 2, "pbrent finds a factor of 18686551294184381720251"); 609 is($nfac[0] * $nfac[1], $nbig, "pbrent found a correct factor"); 610 ok($nfac[0] != 1 && $nfac[1] != 1, "pbrent didn't return a degenerate factor"); 611 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::pminus1_factor($nbig); 612 is(scalar @nfac, 2, "pminus1 finds a factor of 18686551294184381720251"); 613 is($nfac[0] * $nfac[1], $nbig, "pminus1 found a correct factor"); 614 ok($nfac[0] != 1 && $nfac[1] != 1, "pminus1 didn't return a degenerate factor"); 615 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::ecm_factor($nbig); 616 SKIP: { 617 skip "Skipping ecm test", 3 if defined $Math::Prime::Util::GMP::VERSION && $Math::Prime::Util::GMP::VERSION < 0.20; 618 is(scalar @nfac, 2, "ecm finds a factor of 18686551294184381720251"); 619 is($nfac[0] * $nfac[1], $nbig, "ecm found a correct factor"); 620 ok($nfac[0] != 1 && $nfac[1] != 1, "ecm didn't return a degenerate factor"); 621 } 622 623 $nbig = Math::BigInt->new("73786976930493367637"); 624 # Check stage 2 p-1. Fast with Math::BigInt::GMP, slow without. 625 SKIP: { 626 skip "Skipping p-1 stage 2 tests", 3 unless $extra; 627 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::pminus1_factor($nbig, 27000, 35000); 628 is(scalar @nfac, 2, "pminus1 finds a factor of 73786976930493367637"); 629 is($nfac[0] * $nfac[1], $nbig, "pminus1 found a correct factor"); 630 ok($nfac[0] != 1 && $nfac[1] != 1, "pminus1 didn't return a degenerate factor"); 631 } 632 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::fermat_factor($nbig); 633 is(scalar @nfac, 2, "fermat finds a factor of 73786976930493367637"); 634 is($nfac[0] * $nfac[1], $nbig, "fermat found a correct factor"); 635 ok($nfac[0] != 1 && $nfac[1] != 1, "fermat didn't return a degenerate factor"); 636 if ($extra) { 637 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::holf_factor($nbig); 638 is(scalar @nfac, 2, "holf finds a factor of 18686551294184381720251"); 639 is($nfac[0] * $nfac[1], $nbig, "holf found a correct factor"); 640 ok($nfac[0] != 1 && $nfac[1] != 1, "holf didn't return a degenerate factor"); 641 } 642 { 643 $nbig = Math::BigInt->new("99999999999979999998975857"); 644 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::holf_factor($nbig); 645 is_deeply(\@nfac, [9999999998987,10000000001011], "holf correctly factors 99999999999979999998975857"); 646 } 647 SKIP: { 648 # Unfortunately we can't guarantee this isn't found in stage 1. 649 skip "ecm stage 2", 3 unless $extra; 650 $nbig = Math::BigInt->new("14270401808568703916861"); 651 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::ecm_factor($nbig, 5, 2000, 40); 652 is(scalar @nfac, 2, "ecm(5,2000) finds a factor of 14270401808568703916861"); 653 is($nfac[0] * $nfac[1], $nbig, "ecm(5,2000) found a correct factor"); 654 ok($nfac[0] != 1 && $nfac[1] != 1, "ecm(5,2000) didn't return a degenerate factor"); 655 } 656} 657 658##### Some numbers that go to stage 2 of tests 659if ($extra) { 660 my $nbig = Math::BigInt->new("9087500560545072247139"); 661 my @nfac; 662 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::pminus1_factor($nbig,1000,10000); 663 is_deeply( [@nfac], ["24133","376559091722747783"], "p-1 stage 2 finds factors of $nbig" ); 664 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::trial_factor($nbig, 50000); 665 is_deeply( [@nfac], ["24133","376559091722747783"], "trial factor finds factors of $nbig" ); 666 @nfac = sort {$a<=>$b} Math::Prime::Util::PP::ecm_factor($nbig, 10,1000,100); 667 is_deeply( [@nfac], ["24133","376559091722747783"], "ecm factor finds factors of $nbig" ); 668} 669 670##### AKS primality test. Be very careful with performance. 671is( is_aks_prime(1), 0, "AKS: 1 is composite (less than 2)" ); 672is( is_aks_prime(2), 1, "AKS: 2 is prime" ); 673is( is_aks_prime(3), 1, "AKS: 3 is prime" ); 674is( is_aks_prime(4), 0, "AKS: 4 is composite" ); 675is( is_aks_prime(64), 0, "AKS: 64 is composite (perfect power)" ); 676is( is_aks_prime(65), 0, "AKS: 65 is composite (caught in trial)" ); 677is( is_aks_prime(23), 1, "AKS: 23 is prime (r >= n)" ); 678is( is_aks_prime(70747), 0, "AKS: 70747 is composite (n mod r)" ); 679SKIP: { 680 skip "Skipping PP AKS test without EXTENDED_TESTING", 2 unless $extra; 681 diag "32-bit Perl will be very slow for AKS" unless $use64; 682 is( is_aks_prime(1009), 1, "AKS: 1009 is prime (passed anr test)" ); 683 is( is_aks_prime(74513), 0, "AKS: 74513 is composite (failed anr test)" ); 684} 685 686is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_lucas(100003)], 687 [2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN 100003\n\nType Lucas\nN 100003\nQ[1] 2\nQ[2] 3\nQ[3] 7\nQ[4] 2381\nA 2\n"], 688 "primality_proof_lucas(100003)" ); 689# Had to reduce these to make borked up Perl 5.6.2 work. 690#is_deeply( [Math::Prime::Util::PP::primality_proof_bls75("210596120454733723")], 691# [2, ["210596120454733723", "n-1", [2, 3, 82651, "47185492693"], [2, 2, 2, 2]]], 692# "primality_proof_bls75(210596120454733723)" ); 693is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_bls75(1490266103)], 694 [2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN 1490266103\n\nType BLS5\nN 1490266103\nQ[1] 13\nQ[2] 19\nQ[3] 1597\nQ[4] 1889\nA[0] 5\n----\n"], 695 "primality_proof_bls75(1490266103)" ); 696if ($use64) { 697 is_deeply( [Math::Prime::Util::PrimalityProving::primality_proof_bls75(27141057803)], 698 [2, "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN 27141057803\n\nType BLS5\nN 27141057803\nQ[1] 47533\nQ[2] 285497\n----\n"], 699 "primality_proof_bls75(27141057803)" ); 700} 701 702{ 703 my $n = Math::BigInt->new("168790877523676911809192454171451"); 704 is( is_strong_pseudoprime( $n, 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47), 1, "168790877523676911809192454171451 looks prime with bases 2..52" ); 705 is( is_strong_pseudoprime( $n, 53), 0, "168790877523676911809192454171451 found composite with base 53" ); 706 is( is_strong_lucas_pseudoprime($n), 0, "168790877523676911809192454171451 is not a strong Lucas pseudoprime" ); 707 SKIP: { 708 skip "Old Perl+bigint segfaults in F-U code", 1 if $] < 5.008; 709 is( is_frobenius_underwood_pseudoprime($n), 0, "168790877523676911809192454171451 is not a Frobenius pseudoprime" ); 710 } 711 #is( is_perrin_pseudoprime($n), 0, "168790877523676911809192454171451 is not a Perrin pseudoprime" ); 712 is(is_perrin_pseudoprime(517697641), 1, "517697641 is a Perrin pseudoprime"); 713 is(is_frobenius_pseudoprime(517697641), 0, "517697641 is not a Frobenius pseudoprime"); 714} 715 716{ 717 my $ntha = nth_prime_approx(1287248); 718 ok( $ntha >= 20274907 && $ntha <= 20284058, "nth_prime_approx(1287248) in range" ); 719 my $pca = prime_count_approx(128722248); 720 ok( $pca >= 7309252 && $pca <= 7310044, "prime_count_approx(128722248) in range" ); 721} 722 723{ 724 # Test some functions usually not tested in Util.pm 725 my $xs = prime_get_config->{'xs'}; 726 my $gmp = prime_get_config->{'gmp'}; 727 my $verbose = prime_get_config->{'verbose'}; 728 prime_set_config(xs=>0, gmp=>0); 729 730 is( consecutive_integer_lcm(13), 360360, "consecutive_integer_lcm(13)" ); 731 is( consecutive_integer_lcm(52), Math::BigInt->new("3099044504245996706400"), "consecutive_integer_lcm(52)" ); 732 733 is_deeply( [moebius(513,537)], 734 [qw/0 1 1 0 1 -1 1 0 -1 0 -1 0 0 1 1 0 0 -1 0 0 1 -1 1 0 1/], 735 "moebius(513,537)" ); 736 is( moebius(42199), 1, "moebius(42199)" ); 737 is( liouville(444456), 1, "liouville(444456)" ); 738 is( liouville(562894), -1, "liouville(562894)" ); 739 740 is( mertens(4219), -13, "mertens(4219)" ); 741 742 is_deeply( [euler_phi(1513,1537)], 743 [qw/1408 756 800 756 1440 440 1260 576 936 760 1522 504 1200 648 1016 760 1380 384 1530 764 864 696 1224 512 1456/], 744 "euler_phi(1513,1537)" ); 745 is( euler_phi(324234), 108072, "euler_phi(324234)" ); 746 is( jordan_totient(4, 899), "653187225600", "jordan_totient(4, 899)" ); 747 is( carmichael_lambda(324234), 18012, "carmichael_lambda(324234)" ); 748 749 is( exp_mangoldt(16), 2, "exp_mangoldt of power of 2 = 2" ); 750 is( exp_mangoldt(14), 1, "exp_mangoldt of even = 1" ); 751 is( exp_mangoldt(21), 1, "exp_mangoldt of 21 = 1" ); 752 is( exp_mangoldt(23), 23, "exp_mangoldt of 23 = 23" ); 753 is( exp_mangoldt(27), 3, "exp_mangoldt of 27 (3^3) = 3" ); 754 755 is_deeply( [map { scalar znprimroot($_) } (-11, 0, 8, 3, 1729, 10, 5109721)], 756 [2, undef, undef, 2, undef, 3, 94], 757 "znprimroot" ); 758 759 is(znorder(2,35), 12, "znorder(2,35) = 12"); 760 is(znorder(7,35), undef, "znorder(7,35) = undef"); 761 is(znorder(67,999999749), 30612237, "znorder(67,999999749) = 30612237"); 762 763 is(znlog(5678, 5, 10007), 8620, "znlog(5678, 5, 10007)"); 764 765 is(binomial(35,16), 4059928950, "binomial(35,16)"); 766 is(binomial(228,12), "30689926618143230620", "binomial(228,12)"); 767 is(binomial(-23,-26), -2300, "binomial(-23,-26) should be -2300"); 768 769 is(stirling(12,4,2), '611501', "S(12,4)" ); 770 is(stirling(12,4,1), '105258076', "s(12,4)" ); 771 772 is_deeply( [bernfrac(0)], [1,1], "bernfrac(0)" ); 773 is_deeply( [bernfrac(1)], [1,2], "bernfrac(1)" ); 774 is_deeply( [bernfrac(2)], [1,6], "bernfrac(2)" ); 775 is_deeply( [bernfrac(3)], [0,1], "bernfrac(3)" ); 776 is_deeply( [bernfrac(12)], [-691,2730], "bernfrac(12)" ); 777 is_deeply( [bernfrac(13)], [0,1], "bernfrac(12)" ); 778 779 is_deeply( [gcdext(23948236,3498248)], [2263, -15492, 52], "gcdext(23948236,3498248)" ); 780 781 is( valuation(1879048192,2), 28, "valuation(1879048192,2)"); 782 is( valuation(96552,6), 3, "valuation(96552,6)"); 783 784 is(invmod(45,59), 21, "invmod(45,59)"); 785 is(invmod(14,28474), undef, "invmod(14,28474)"); 786 is(invmod(42,-2017), 1969, "invmod(42,-2017)"); 787 788 is(vecsum(15, 30, 45), 90, "vecsum(15,30,45)"); 789 is(vecsum(4294966296,4294965296,4294964296), "12884895888", "vecsum(2^32-1000,2^32-2000,2^32-3000)"); 790 is(vecprod(15, 30, 45), 20250, "vecprod(15,30,45)"); 791 is(vecprod(4294966296,4294965296,4294964296), "79228051833847139970490254336", "vecprod(2^32-1000,2^32-2000,2^32-3000)"); 792 is(vecmin(4294966296,4294965296,4294964296), 4294964296, "vecmin(2^32-1000,2^32-2000,2^32-3000)"); 793 is(vecmax(4294966296,4294965296,4294964296), 4294966296, "vecmax(2^32-1000,2^32-2000,2^32-3000)"); 794 795 cmp_closeto( chebyshev_theta(7001), 6929.27483821865062, 0.006929, "chebyshev_theta(7001) =~ 6929.2748"); 796 cmp_closeto( chebyshev_psi(6588), 6597.07452996633704, 0.006597, "chebyshev_psi(6588) =~ 6597.07453"); 797 798 while (my($n, $isp) = each (%ipp)) { 799 is( is_prob_prime($n), $isp, "is_prob_prime($n) should be $isp" ); 800 } 801 802 is( primorial(24), 223092870, "primorial(24)" ); 803 is( primorial(118), "31610054640417607788145206291543662493274686990", "primorial(118)" ); 804 is( pn_primorial(7), 510510, "pn_primorial(7)" ); 805 is( partitions(74), 7089500, "partitions(74)" ); 806 is( miller_rabin_random(4294967281, 20), "0", "Miller-Rabin random 40 on composite" ); 807 808 { my @t; 809 Math::Prime::Util::_generic_forprimes(sub {push @t,$_}, 2387234,2387303); 810 is_deeply( [@t], [2387237,2387243,2387249,2387269,2387291,2387299,2387303], 811 "generic forprimes 2387234,2387303" ); 812 } 813 { my @t; 814 Math::Prime::Util::_generic_forcomposites(sub {push @t,$_}, 15202630,15202641); 815 is_deeply( [@t], [15202630,15202632,15202634,15202635,15202636,15202638,15202640,15202641], "generic forcomposites 15202630,15202641" ); 816 } 817 { my @t; 818 Math::Prime::Util::_generic_foroddcomposites(sub {push @t,$_}, 15202630,15202641); 819 is_deeply( [@t], [15202635,15202641], "generic foroddcomposites 15202630,15202641" ); 820 } 821 { my $k = 0; 822 Math::Prime::Util::_generic_fordivisors(sub {$k += $_+int(sqrt($_))},92834); 823 is( $k, 168921, "generic fordivisors: d|92834: k+=d+int(sqrt(d))" ); 824 } 825 826 { my @p; 827 forcomb(sub { push @p, [@_] }, 3, 2); 828 is_deeply( \@p, [ [0,1], [0,2], [1,2] ], "forcomb(3,2)" ); 829 } 830 { my @p; 831 forperm(sub { push @p, [@_] }, 3); 832 is_deeply( \@p, [ [0,1,2], [0,2,1], [1,0,2], [1,2,0], [2,0,1], [2,1,0] ], "forperm(3)" ); 833 } 834 { my @p; 835 forpart(sub { push @p, [@_] }, 4); 836 is_deeply( \@p, [ [1,1,1,1],[1,1,2],[1,3],[2,2],[4] ], "forpart(4)" ); 837 } 838 839 is( Pi(82), "3.141592653589793238462643383279502884197169399375105820974944592307816406286208999", "Pi(82)" ); 840 841 is( gcd(-30,-90,90), 30, "gcd(-30,-90,90) = 30" ); 842 is( lcm(11926,78001,2211), 2790719778, "lcm(11926,78001,2211) = 2790719778" ); 843 844 is( twin_prime_count(4321), 114, "twin_prime_count(4321)" ); 845 cmp_closeto( twin_prime_count_approx(Math::BigInt->new("4123456784123")), "6950213327", 14937 * 2, "twin_prime_count_approx(4123456784123)" ); 846 if ($extra) { 847 cmp_closeto( twin_prime_count_approx(Math::BigInt->new("412345678412345678412345678")), "149939117920176008847283", 1e11, "twin_prime_count_approx(412345678412345678412345678)" ); 848 } 849 is( nth_twin_prime(249), 13217, "nth_twin_prime(249)" ); 850 851 prime_set_config(xs=>$xs, gmp=>$gmp, verbose=>$verbose); 852} 853 854is( $_, 'this should not change', "Nobody clobbered \$_" ); 855 856############################################################################### 857 858sub parse_range { 859 my($range) = @_; 860 my($low,$high); 861 my $fixnum = sub { 862 my $nstr = shift; 863 $nstr =~ s/^(\d+)e(\d+)$/$1*(10**$2)/e; 864 $nstr =~ s/^(\d+)\*\*(\d+)$/$1**$2/e; 865 die "Unknown string in test" unless $nstr =~ /^\d+$/; 866 $nstr; 867 }; 868 if ($range =~ /(\S+)\s+to\s+(\S+)/) { 869 $low = $fixnum->($1); 870 $high = $fixnum->($2); 871 } elsif ($range =~ /(\S+)\s*\+\s*(\S+)/) { 872 $low = $fixnum->($1); 873 $high = $low + $fixnum->($2); 874 } else { 875 die "Can't parse test data"; 876 } 877 ($low,$high); 878} 879 880sub cmp_closeto { 881 my $got = shift; 882 my $expect = shift; 883 my $tolerance = shift; 884 my $message = shift; 885 cmp_ok( abs($got - $expect), '<=', $tolerance, $message ); 886} 887