new file: Math/prime_counting_from_almost_primes.pl

new file:   Math/prime_counting_from_squarefree_almost_primes.pl

Author: trizen

Math/prime_counting_from_almost_primes.pl

@@ -0,0 +1,77 @@
+#!/usr/bin/perl
+
+# Author: Daniel "Trizen" Șuteu
+# Date: 27 August 2025
+# https://github.com/trizen
+
+# A sublinear algorithm for computing the Prime Counting function `pi(n)`,
+# based on the number of k-almost primes <= n, for `k >= 2`, which can be computed in sublinear time.
+
+# See also:
+#   https://mathworld.wolfram.com/AlmostPrime.html
+
+use 5.036;
+use ntheory qw(:all);
+
+sub k_prime_count ($k, $n) {
+
+    if ($k == 1) {
+        return my_prime_count($n);
+    }
+
+    my $count = 0;
+
+    sub ($m, $p, $k, $j = 0) {
+
+        my $s = rootint(divint($n, $m), $k);
+
+        if ($k == 2) {
+
+            forprimes {
+                $count += my_prime_count(divint($n, mulint($m, $_))) - $j++;
+            } $p, $s;
+
+            return;
+        }
+
+        foreach my $q (@{primes($p, $s)}) {
+            __SUB__->($m * $q, $q, $k - 1, $j++);
+        }
+    }->(1, 2, $k);
+
+    return $count;
+}
+
+sub my_prime_count ($n) {
+
+    state $pi_table = [0, 0, 1, 2, 2];      # a larger lookup table helps a lot!
+
+    if ($n < 0) {
+        return 0;
+    }
+
+    if (defined($pi_table->[$n])) {
+        return $pi_table->[$n];
+    }
+
+    my $M = $n - 1;
+
+    foreach my $k (2 .. logint($n, 2)) {
+        $M -= k_prime_count($k, $n);
+    }
+
+    return ($pi_table->[$n] //= $M);
+}
+
+foreach my $n (1..7) {    # takes ~3 seconds
+    say "pi(10^$n) = ", my_prime_count(10**$n);
+}
+
+__END__
+pi(10^1) = 4
+pi(10^2) = 25
+pi(10^3) = 168
+pi(10^4) = 1229
+pi(10^5) = 9592
+pi(10^6) = 78498
+pi(10^7) = 664579

Math/prime_counting_from_squarefree_almost_primes.pl

@@ -0,0 +1,82 @@
+#!/usr/bin/perl
+
+# Author: Daniel "Trizen" Șuteu
+# Date: 27 August 2025
+# https://github.com/trizen
+
+# A sublinear algorithm for computing the Prime Counting function `pi(n)`,
+# based on the number of squarefree k-almost primes <= n, for `k >= 2`, which can be computed in sublinear time.
+
+# See also:
+#   https://mathworld.wolfram.com/AlmostPrime.html
+
+use 5.036;
+use ntheory qw(:all);
+
+sub squarefree_almost_prime_count ($k, $n) {
+
+    if ($k == 0) {
+        return (($n <= 0) ? 0 : 1);
+    }
+
+    if ($k == 1) {
+        return my_prime_count($n);
+    }
+
+    my $count = 0;
+
+    sub ($m, $p, $k, $j = 1) {
+
+        my $s = rootint(divint($n, $m), $k);
+
+        if ($k == 2) {
+
+            forprimes {
+                $count += my_prime_count(divint($n, mulint($m, $_))) - $j++;
+            } $p, $s;
+
+            return;
+        }
+
+        foreach my $q (@{primes($p, $s)}) {
+            __SUB__->(mulint($m, $q), $q + 1, $k - 1, ++$j);
+        }
+      }
+      ->(1, 2, $k);
+
+    return $count;
+}
+
+sub my_prime_count ($n) {
+
+    state $pi_table = [0, 0, 1, 2, 2];      # a larger lookup table helps a lot!
+
+    if ($n < 0) {
+        return 0;
+    }
+
+    if (defined($pi_table->[$n])) {
+        return $pi_table->[$n];
+    }
+
+    my $M = powerfree_count($n, 2) - 1;
+
+    foreach my $k (2 .. exp(LambertW(log($n))) + 1) {
+        $M -= squarefree_almost_prime_count($k, $n);
+    }
+
+    return ($pi_table->[$n] //= $M);
+}
+
+foreach my $n (1..7) {    # takes ~1 second
+    say "pi(10^$n) = ", my_prime_count(10**$n);
+}
+
+__END__
+pi(10^1) = 4
+pi(10^2) = 25
+pi(10^3) = 168
+pi(10^4) = 1229
+pi(10^5) = 9592
+pi(10^6) = 78498
+pi(10^7) = 664579

Math/prime_counting_liouville_formula.pl

@@ -8,8 +8,7 @@
 # based on the Liouville function and the number of k-almost primes <= n, for `k >= 2`.
 
 # See also:
-#   https://en.wikipedia.org/wiki/Mertens_function
-#   https://en.wikipedia.org/wiki/M%C3%B6bius_function
+#   https://mathworld.wolfram.com/AlmostPrime.html
 
 use 5.036;
 use ntheory qw(:all);

Math/prime_counting_mertens_formula.pl

@@ -62,7 +62,7 @@ ($n)
 
     my $M = mertens($n);
 
-    foreach my $k (2 .. logint($n, 2)) {
+    foreach my $k (2 .. exp(LambertW(log($n))) + 1) {
         $M -= (-1)**$k * squarefree_almost_prime_count($k, $n);
     }
 

README.md

@@ -1067,6 +1067,8 @@ A nice collection of day-to-day Perl scripts.
     * [Prime 41](./Math/prime_41.pl)
     * [Prime abundant sequences](./Math/prime_abundant_sequences.pl)
     * [Prime count smooth sum](./Math/prime_count_smooth_sum.pl)
+    * [Prime counting from almost primes](./Math/prime_counting_from_almost_primes.pl)
+    * [Prime counting from squarefree almost primes](./Math/prime_counting_from_squarefree_almost_primes.pl)
     * [Prime counting liouville formula](./Math/prime_counting_liouville_formula.pl)
     * [Prime counting mertens formula](./Math/prime_counting_mertens_formula.pl)
     * [Prime factorization concept](./Math/prime_factorization_concept.pl)