Codebase list libmath-prime-util-perl / c5a026d
Extend prime count bounds based on new paper Dana Jacobsen 8 years ago
6 changed file(s) with 79 addition(s) and 60 deletion(s). Raw diff Collapse all Expand all
+2
-0
Changes less more
1010 - Faster, nonrecursive divisors_from_factors routine.
1111
1212 - gcdext(0,0) returns (0,0,0) to match GMP and Pari/GP.
13
14 - Use better prime count lower/upper bounds from Büthe 2015.
1315
1416
1517 0.55 2015-10-19
+50
-28
lib/Math/Prime/Util/PP.pm less more
14941494 # Dusart 2010 x/logx*(1+1/logx+2.0/logxlogx) x >= 88783
14951495 # Axler 2014 (1.2) ""+... x >= 1332450001
14961496 # Axler 2014 (1.2) x/(logx-1-1/logx-...) x >= 1332479531
1497 # Büthe 2015 (1.9) li(x)-(sqrtx/logx)*(...) x <= 10^19
14971498 # Büthe 2014 Th 2 li(x)-logx*sqrtx/8Pi x > 2657, x <= 1.4*10^25
14981499
1499 if ($x >= 5589603006000 && # Schoenfeld / Büthe 2014 Th 7.4
1500 ($x < 1.4e25 || Math::Prime::Util::prime_get_config()->{'assume_rh'}) ) {
1500 if ($x < 599) { # Decent for small numbers
1501 $result = $x / ($fl1 - 0.7);
1502 } elsif ($x < 52600000) { # Dusart 2010 tweaked
1503 if ($x < 2700) { $a = 0.30; }
1504 elsif ($x < 5500) { $a = 0.90; }
1505 elsif ($x < 19400) { $a = 1.30; }
1506 elsif ($x < 32299) { $a = 1.60; }
1507 elsif ($x < 88783) { $a = 1.83; }
1508 elsif ($x < 176000) { $a = 1.99; }
1509 elsif ($x < 315000) { $a = 2.11; }
1510 elsif ($x < 1100000) { $a = 2.19; }
1511 elsif ($x < 4500000) { $a = 2.31; }
1512 else { $a = 2.35; }
1513 $result = ($x/$fl1) * ($one + $one/$fl1 + $a/$fl2);
1514 } elsif ($x < 1.4e25 || Math::Prime::Util::prime_get_config()->{'assume_rh'}){
1515 # Büthe 2014/2015
15011516 if (_MPFR_available()) {
15021517 my $wantbf = (defined $bignum::VERSION || ref($x) =~ /^Math::Big/);
15031518 my $xdigits = length($x);
15141529 Math::MPFR::Rmpfr_set_prec($sqx, $bit_precision);
15151530 Math::MPFR::Rmpfr_sqrt($sqx, $rx, $bit_precision);
15161531 Math::MPFR::Rmpfr_log($rx, $rx, $rnd);
1517 Math::MPFR::Rmpfr_mul($rx, $rx, $sqx, $rnd);
1518 Math::MPFR::Rmpfr_const_pi($sqx, $rnd);
1519 Math::MPFR::Rmpfr_mul_ui($sqx, $sqx, 8, $rnd);
1520 Math::MPFR::Rmpfr_div($rx, $rx, $sqx, $rnd);
1532 # rx = log(x) lix = li(x) sqx = sqrt(x)
1533 if ($x < 1e19) { # Büthe 2015 1.9
1534 #Math::MPFR::Rmpfr_div($sqx, $sqx, $rx, $rnd);
1535 #$rx = 1.94 + 3.88/$rx + 27.57/($rx*$rx);
1536 my $tmp = Math::MPFR->new();
1537 Math::MPFR::Rmpfr_set_prec($tmp, $bit_precision);
1538 my $acc = Math::MPFR->new();
1539 Math::MPFR::Rmpfr_set_prec($acc, $bit_precision);
1540 Math::MPFR::Rmpfr_set_d($acc, 1.94, $rnd);
1541 Math::MPFR::Rmpfr_d_div($tmp, 3.88, $rx, $rnd);
1542 Math::MPFR::Rmpfr_add($acc, $acc, $tmp, $rnd);
1543 Math::MPFR::Rmpfr_d_div($tmp, 27.57, $rx*$rx, $rnd);
1544 Math::MPFR::Rmpfr_add($acc, $acc, $tmp, $rnd);
1545 Math::MPFR::Rmpfr_mul($rx, $acc, $sqx/$rx, $rnd);
1546 #Math::MPFR::Rmpfr_mul($rx, $rx, $sqx, $rnd);
1547 } else { # Büthe 2014 7.4
1548 Math::MPFR::Rmpfr_mul($rx, $rx, $sqx, $rnd);
1549 Math::MPFR::Rmpfr_const_pi($sqx, $rnd);
1550 Math::MPFR::Rmpfr_mul_ui($sqx, $sqx, 8, $rnd);
1551 Math::MPFR::Rmpfr_div($rx, $rx, $sqx, $rnd);
1552 }
15211553 Math::MPFR::Rmpfr_sub($rx, $lix, $rx, $rnd);
15221554 my $strval = Math::MPFR::Rmpfr_integer_string($rx, 10, $rnd);
15231555 $result = ($wantbf) ? Math::BigInt->new($strval) : int($strval);
15241556 } else {
15251557 my $lix = LogarithmicIntegral($x);
15261558 my $sqx = sqrt($x);
1527 if (ref($x) eq 'Math::BigFloat') {
1528 my $xdigits = _find_big_acc($x);
1529 $result = $lix - ($fl1*$sqx / (Math::BigFloat->bpi($xdigits)*8));
1559 if ($x < 1e19) {
1560 $result = $lix - ($sqx/$fl1) * (1.94 + 3.88/$fl1 + 27.57/$fl2);
15301561 } else {
1531 $result = $lix - ($fl1*$sqx / PI_TIMES_8);
1532 }
1533 }
1534 } elsif ($x < 599) { # Decent for small numbers
1535 $result = $x / ($fl1 - 0.7);
1536 } elsif ($x < 1332479531) { # Dusart 2010 tweaked
1537 if ($x < 2700) { $a = 0.30; }
1538 elsif ($x < 5500) { $a = 0.90; }
1539 elsif ($x < 19400) { $a = 1.30; }
1540 elsif ($x < 32299) { $a = 1.60; }
1541 elsif ($x < 176000) { $a = 1.80; }
1542 elsif ($x < 315000) { $a = 2.10; }
1543 elsif ($x < 1100000) { $a = 2.20; }
1544 elsif ($x < 4500000) { $a = 2.31; }
1545 elsif ($x < 233000000) { $a = 2.36; }
1546 else { $a = 2.32; }
1547 $result = ($x/$fl1) * ($one + $one/$fl1 + $a/$fl2);
1562 if (ref($x) eq 'Math::BigFloat') {
1563 my $xdigits = _find_big_acc($x);
1564 $result = $lix - ($fl1*$sqx / (Math::BigFloat->bpi($xdigits)*8));
1565 } else {
1566 $result = $lix - ($fl1*$sqx / PI_TIMES_8);
1567 }
1568 }
1569 }
15481570 } else { # Axler 2014 1.4
15491571 if (_MPFR_available()) {
15501572 my $wantbf = (defined $bignum::VERSION || ref($x) =~ /^Math::Big/);
16531675 elsif ($x < 2953652287) { $a = 2.362; }
16541676 else { $a = 2.334; } # Dusart 2010, page 2
16551677 $result = ($x/$fl1) * ($one + $one/$fl1 + $a/$fl2) + $one;
1656 } elsif ($x < 1e14) { # Skewes number lower limit
1657 $a = ($x < 1.1e9) ? 0.031 : ($x < 4.5e9) ? 0.02 : 0.0;
1678 } elsif ($x < 1e19) { # Skewes number lower limit
1679 $a = ($x < 110e7) ? 0.032 : ($x < 1001e7) ? 0.027 : ($x < 10126e7) ? 0.021 : 0.0;
16581680 $result = LogarithmicIntegral($x) - $a * $fl1*sqrt($x)/PI_TIMES_8;
16591681 } elsif ($x < 1.4e25 || Math::Prime::Util::prime_get_config()->{'assume_rh'}) {
16601682 # Schoenfeld / Büthe 2014 Th 7.4
+3
-3
lib/Math/Prime/Util.pm less more
14501450 avoid calling C<prime_count>.
14511451
14521452 These routines use verified tight limits below a range at least C<2^35>.
1453 For larger inputs various methods are used including Kotnick (2008),
1454 Büthe (2014), and Axler (2014).
1453 For larger inputs various methods are used including Dusart (2010),
1454 Büthe (2014,2015), and Axler (2014).
14551455 These bounds do not assume the Riemann Hypothesis.
14561456 If the configuration option C<assume_rh> has been set (it is off by default),
1457 then the Schoenfeld (1976) bounds can be used for large values.
1457 then the Schoenfeld (1976) bounds can be used for very large values.
14581458
14591459
14601460 =head2 prime_count_approx
+4
-4
t/80-pp.t less more
441441 is( prime_count_lower(450), 87, "prime_count_lower(450)" );
442442 is( prime_count_upper(450), 87, "prime_count_upper(450)" );
443443 # Make sure these are about right
444 cmp_closeto( prime_count_lower(1234567), 95327, 10, "prime_count_lower(1234567) in range" );
445 cmp_closeto( prime_count_upper(1234567), 95413, 10, "prime_count_upper(1234567) in range" );
446 cmp_closeto( prime_count_lower(412345678), 21956686, 1000, "prime_count_lower(412345678) in range" );
447 cmp_closeto( prime_count_upper(412345678), 21959328, 1000, "prime_count_upper(412345678) in range" );
444 cmp_closeto( prime_count_lower(1234567), 95360, 60, "prime_count_lower(1234567) in range" );
445 cmp_closeto( prime_count_upper(1234567), 95360, 60, "prime_count_upper(1234567) in range" );
446 cmp_closeto( prime_count_lower(412345678), 21958997, 1500, "prime_count_lower(412345678) in range" );
447 cmp_closeto( prime_count_upper(412345678), 21958997, 1500, "prime_count_upper(412345678) in range" );
448448
449449 ###############################################################################
450450
+6
-4
t/81-bignum.t less more
7575 + 1 # primecount large base small range
7676 + scalar(keys %pseudoprimes)
7777 + 6 # PC lower, upper, approx
78 + 6*2*$extra # more PC tests
78 + 6*3*$extra # more PC tests
7979 + 2*scalar(keys %factors)
8080 + scalar(keys %allfactors)
8181 + 14+3*$extra # moebius, euler_phi, jordan totient, divsum, etc.
234234
235235 ###############################################################################
236236
237 check_pcbounds(31415926535897932384, 716115441142294636, '8e-5', '2e-8');
237 check_pcbounds(31415926535897932384, 716115441142294636, '2e-8', '2e-8');
238238 if ($extra) {
239 check_pcbounds(314159265358979323846, 6803848951392700268, '7e-5', '5e-9');
240 check_pcbounds(31415926535897932384626433, 544551456607147153724423, '4e-5', '3e-11');
239 check_pcbounds(314159265358979323846, 6803848951392700268, '5e-9', '5e-9');
240 check_pcbounds(31415926535897932384626433, 544551456607147153724423, '3e-6', '3e-11');
241 # pi(10^23) = 1925320391606803968923
242 check_pcbounds(10**23, 1925320391607837268776, '5e-10', '5e-10');
241243 }
242244
243245 ###############################################################################
+14
-21
util.c less more
626626 fl1 = logl(n);
627627 fl2 = fl1 * fl1;
628628
629 if (n < UVCONST(1332479531)) { /* Dusart 2010, tweaked */
630 a = (n < 88783) ? 1.89L
629 if (n < UVCONST(52600000)) { /* Dusart 2010, tweaked */
630 a = (n < 88783) ? 1.89L /* n >= 33000 */
631631 : (n < 176000) ? 1.99L
632632 : (n < 315000) ? 2.11L
633633 : (n < 1100000) ? 2.19L
634 : (n < 4500000) ? 2.31L
635 : (n <233000000) ? 2.35L : 2.32L;
634 : (n < 4500000) ? 2.31L : 2.35;
636635 lower = fn/fl1 * (1.0L + 1.0L/fl1 + a/fl2);
637 } else if (BITS_PER_WORD == 32 || n < UVCONST(5589603006000)) {
638 /* Axler 2014 1.4 */
639 long double fl3 = fl2*fl1, fl4 = fl2*fl2, fl5 = fl4*fl1, fl6 = fl4*fl2;
640 lower = fn / (fl1 - 1.0L - 1.0L/fl1 - 2.65L/fl2 - 13.35L/fl3 - 70.3L/fl4 - 455.6275L/fl5 - 3404.4225L/fl6);
641 } else { /* Büthe 2014 7.4 */
636 } else if (fn < 1e19) { /* Büthe 2015 1.9 */
637 lower = _XS_LogarithmicIntegral(fn) - (sqrtl(fn)/fl1) * (1.94L + 3.88L/fl1 + 27.57L/fl2);
638 } else { /* Büthe 2014 v3 7.2 */
642639 lower = _XS_LogarithmicIntegral(fn) - fl1*sqrtl(fn)/25.132741228718345907701147L;
643640 }
644641 return (UV) ceill(lower);
685682 fl1 = logl(n);
686683 fl2 = fl1 * fl1;
687684
688 if (BITS_PER_WORD == 32 || fn <= 821800000.0) {
685 if (BITS_PER_WORD == 32 || fn <= 821800000.0) { /* Dusart 2010, page 2 */
689686 for (i = 0; i < (int)NUPPER_THRESH; i++)
690687 if (n < _upper_thresh[i].thresh)
691688 break;
692 if (i < (int)NUPPER_THRESH) a = _upper_thresh[i].aval;
693 else a = 2.334; /* Dusart 2010, page 2 */
689 a = (i < (int)NUPPER_THRESH) ? _upper_thresh[i].aval : 2.334L;
694690 upper = fn/fl1 * (1.0L + 1.0L/fl1 + a/fl2);
695 } else if (fn < 1e14) { /* Skewes number lower limit */
696 a = (fn < 1100000000) ? 0.031
697 : (fn < 4500000000) ? 0.02 : 0.0;
691 } else if (fn < 1e19) { /* Büthe 2015 1.10 Skewes number lower limit */
692 a = (fn < 1100000000) ? 0.032 /* Empirical */
693 : (fn < 10010000000) ? 0.027 /* Empirical */
694 : (fn < 101260000000) ? 0.021 /* Empirical */
695 : 0.0;
698696 upper = _XS_LogarithmicIntegral(fn) - a * fl1*sqrtl(fn)/25.132741228718345907701147L;
699 #if 0
700 } else if (fn < 5796000000000.0) { /* Axler 2014, theorem 1.3 */
701 long double fl3 = fl2*fl1, fl4 = fl2*fl2, fl5 = fl4*fl1, fl6 = fl4*fl2;
702 upper = fn / (fl1 - 1.0L - 1.0L/fl1 - 3.35L/fl2 - 12.65L/fl3 - 71.7L/fl4 - 466.1275L/fl5 - 3489.8225L/fl6);
703 #endif
704 } else { /* Büthe 2014 7.4 */
697 } else { /* Büthe 2014 7.4 */
705698 upper = _XS_LogarithmicIntegral(fn) + fl1*sqrtl(fn)/25.132741228718345907701147L;
706699 }
707700 return (UV) floorl(upper);