Skip to content

Commit 4534977

Browse files
trizentrizen
committed
new file: Math/prime_counting_liouville_formula.pl
new file: Math/prime_counting_mertens_formula.pl
1 parent a35048a commit 4534977

6 files changed

+181
-28
lines changed

‎Math/count_of_k-almost_primes.pl‎

Lines changed: 5 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -19,9 +19,8 @@
1919
# https://oeis.org/A082996 -- count of 4-almost primes
2020
# https://oeis.org/A126280 -- Triangle read by rows: T(k,n) is number of numbers <= 10^n that are products of k primes.
2121

22-
use 5.020;
22+
use 5.036;
2323
use ntheory qw(:all);
24-
use experimental qw(signatures);
2524

2625
sub k_prime_count ($n, $k) {
2726

@@ -41,14 +40,14 @@ ($n, $k)
4140

4241
if ($k == 2) {
4342

44-
foreach my $q (@{primes($p, $s)}) {
45-
$count += prime_count(divint($n, mulint($m, $q))) - $j++;
46-
}
43+
forprimes {
44+
$count += prime_count(divint($n, mulint($m, $_))) - $j++;
45+
} $p, $s;
4746

4847
return;
4948
}
5049

51-
for (my $q = $p ; $q <= $s ; $q = next_prime($q)) {
50+
foreach my $q (@{primes($p, $s)}) {
5251
__SUB__->($m * $q, $q, $k - 1, $j++);
5352
}
5453
}->(1, 2, $k);

‎Math/count_of_squarefree_k-almost_primes.pl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -39,9 +39,9 @@ ($n, $k)
3939

4040
if ($k == 2) {
4141

42-
for (; $p <= $s ; $p = next_prime($p)) {
43-
$count += prime_count(divint($n, mulint($m, $p))) - $j++;
44-
}
42+
forprimes {
43+
$count += prime_count(divint($n, mulint($m, $_))) - $j++;
44+
} $p, $s;
4545

4646
return;
4747
}
@@ -63,7 +63,7 @@ ($n, $k)
6363
my $upto = pn_primorial($k) + int(rand(1e5));
6464

6565
my $x = squarefree_almost_prime_count($upto, $k);
66-
my $y = scalar grep { is_square_free($_) } @{almost_primes($k, $upto)};
66+
my $y = scalar grep { is_square_free($_) } @{almost_primes($k, 1, $upto)};
6767

6868
say "Testing: $k with n = $upto -> $x";
6969

‎Math/partial_sums_of_exponential_prime_omega_functions.pl‎

Lines changed: 9 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,8 @@
99
# S2(n) = Sum_{k=1..n} v^omega(k)
1010
# S3(n) = Sum_{k=1..n} v^omega(k) * mu(k)^2
1111

12-
use 5.020;
13-
use warnings;
14-
12+
use 5.036;
1513
use ntheory qw(:all);
16-
use experimental qw(signatures);
1714

1815
sub squarefree_almost_prime_count ($k, $n) {
1916

@@ -27,29 +24,23 @@ ($k, $n)
2724

2825
my $count = 0;
2926

30-
sub ($m, $p, $k, $j = 0) {
27+
sub ($m, $p, $k, $j = 1) {
3128

3229
my $s = rootint(divint($n, $m), $k);
3330

3431
if ($k == 2) {
3532

36-
foreach my $q (@{primes($p, $s)}) {
37-
38-
++$j;
39-
40-
if (modint($m, $q) != 0) {
41-
$count += prime_count(divint($n, mulint($m, $q))) - $j;
42-
}
43-
}
33+
forprimes {
34+
$count += prime_count(divint($n, mulint($m, $_))) - $j++;
35+
} $p, $s;
4436

4537
return;
4638
}
4739

48-
foreach my $p (@{primes($p, $s)}) {
49-
if (modint($m, $p) != 0) {
50-
__SUB__->(mulint($m, $p), $p, $k - 1, $j);
51-
}
52-
++$j;
40+
for (; $p <= $s ; ++$j) {
41+
my $r = next_prime($p);
42+
__SUB__->(mulint($m, $p), $r, $k - 1, $j + 1);
43+
$p = $r;
5344
}
5445
}->(1, 2, $k);
5546

‎Math/prime_counting_liouville_formula.pl‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
#!/usr/bin/perl
2+
3+
# Author: Trizen
4+
# Date: 17 July 2025
5+
# https://github.com/trizen
6+
7+
# A sublinear algorithm for computing the Prime Counting function `pi(n)`,
8+
# based on the Liouville function and the number of k-almost primes <= n, for `k >= 2`.
9+
10+
# See also:
11+
# https://en.wikipedia.org/wiki/Mertens_function
12+
# https://en.wikipedia.org/wiki/M%C3%B6bius_function
13+
14+
use 5.036;
15+
use ntheory qw(:all);
16+
17+
sub k_prime_count ($k, $n) {
18+
19+
if ($k == 1) {
20+
return my_prime_count($n);
21+
}
22+
23+
my $count = 0;
24+
25+
sub ($m, $p, $k, $j = 0) {
26+
27+
my $s = rootint(divint($n, $m), $k);
28+
29+
if ($k == 2) {
30+
31+
forprimes {
32+
$count += my_prime_count(divint($n, mulint($m, $_))) - $j++;
33+
} $p, $s;
34+
35+
return;
36+
}
37+
38+
foreach my $q (@{primes($p, $s)}) {
39+
__SUB__->($m * $q, $q, $k - 1, $j++);
40+
}
41+
}->(1, 2, $k);
42+
43+
return $count;
44+
}
45+
46+
sub my_prime_count ($n) {
47+
48+
state $pi_table = [0, 0, 1, 2, 2]; # a larger lookup table helps a lot!
49+
50+
if ($n < 0) {
51+
return 0;
52+
}
53+
54+
if (defined($pi_table->[$n])) {
55+
return $pi_table->[$n];
56+
}
57+
58+
my $M = sumliouville($n);
59+
60+
foreach my $k (2 .. logint($n, 2)) {
61+
$M -= (-1)**$k * k_prime_count($k, $n);
62+
}
63+
64+
return ($pi_table->[$n] //= 1 - $M);
65+
}
66+
67+
foreach my $n (1..7) { # takes ~3 seconds
68+
say "pi(10^$n) = ", my_prime_count(10**$n);
69+
}
70+
71+
__END__
72+
pi(10^1) = 4
73+
pi(10^2) = 25
74+
pi(10^3) = 168
75+
pi(10^4) = 1229
76+
pi(10^5) = 9592
77+
pi(10^6) = 78498
78+
pi(10^7) = 664579

‎Math/prime_counting_mertens_formula.pl‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
#!/usr/bin/perl
2+
3+
# Author: Trizen
4+
# Date: 17 July 2025
5+
# https://github.com/trizen
6+
7+
# A sublinear algorithm for computing the Prime Counting function `pi(n)`, based on the
8+
# Mertens function and the number of squarefree k-almost primes <= n, for `k >= 2`.
9+
10+
# See also:
11+
# https://en.wikipedia.org/wiki/Mertens_function
12+
# https://en.wikipedia.org/wiki/M%C3%B6bius_function
13+
14+
use 5.036;
15+
use ntheory qw(:all);
16+
17+
sub squarefree_almost_prime_count ($k, $n) {
18+
19+
if ($k == 0) {
20+
return (($n <= 0) ? 0 : 1);
21+
}
22+
23+
if ($k == 1) {
24+
return my_prime_count($n);
25+
}
26+
27+
my $count = 0;
28+
29+
sub ($m, $p, $k, $j = 1) {
30+
31+
my $s = rootint(divint($n, $m), $k);
32+
33+
if ($k == 2) {
34+
35+
forprimes {
36+
$count += my_prime_count(divint($n, mulint($m, $_))) - $j++;
37+
} $p, $s;
38+
39+
return;
40+
}
41+
42+
foreach my $q (@{primes($p, $s)}) {
43+
__SUB__->(mulint($m, $q), $q + 1, $k - 1, ++$j);
44+
}
45+
}
46+
->(1, 2, $k);
47+
48+
return $count;
49+
}
50+
51+
sub my_prime_count ($n) {
52+
53+
state $pi_table = [0, 0, 1, 2, 2]; # a larger lookup table helps a lot!
54+
55+
if ($n < 0) {
56+
return 0;
57+
}
58+
59+
if (defined($pi_table->[$n])) {
60+
return $pi_table->[$n];
61+
}
62+
63+
my $M = mertens($n);
64+
65+
foreach my $k (2 .. logint($n, 2)) {
66+
$M -= (-1)**$k * squarefree_almost_prime_count($k, $n);
67+
}
68+
69+
return ($pi_table->[$n] //= 1 - $M);
70+
}
71+
72+
foreach my $n (1 .. 7) { # takes ~1 second
73+
say "pi(10^$n) = ", my_prime_count(10**$n);
74+
}
75+
76+
__END__
77+
pi(10^1) = 4
78+
pi(10^2) = 25
79+
pi(10^3) = 168
80+
pi(10^4) = 1229
81+
pi(10^5) = 9592
82+
pi(10^6) = 78498
83+
pi(10^7) = 664579

‎README.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1065,6 +1065,8 @@ A nice collection of day-to-day Perl scripts.
10651065
* [Prime 41](./Math/prime_41.pl)
10661066
* [Prime abundant sequences](./Math/prime_abundant_sequences.pl)
10671067
* [Prime count smooth sum](./Math/prime_count_smooth_sum.pl)
1068+
* [Prime counting liouville formula](./Math/prime_counting_liouville_formula.pl)
1069+
* [Prime counting mertens formula](./Math/prime_counting_mertens_formula.pl)
10681070
* [Prime factorization concept](./Math/prime_factorization_concept.pl)
10691071
* [Prime factors of binomial coefficients](./Math/prime_factors_of_binomial_coefficients.pl)
10701072
* [Prime factors of binomial product](./Math/prime_factors_of_binomial_product.pl)

0 commit comments

Comments
 (0)