new file: Math/count_of_rough_numbers_recursive.pl

modified:   Math/meissel_lehmer_prime_count.pl
modified:   Math/prime_counting_from_squarefree_almost_primes.pl

Author: trizen

Math/count_of_rough_numbers.pl

@@ -10,11 +10,13 @@
 #   https://en.wikipedia.org/wiki/Rough_number
 
 use 5.020;
-use integer;
 use ntheory qw(:all);
 use experimental qw(signatures);
 
-sub rough_count ($n, $p) {
+sub my_rough_count ($n, $p) {
+
+    my %cache;
+
     sub ($n, $p) {
 
         if ($p > sqrtint($n)) {
@@ -26,12 +28,17 @@ ($n, $p)
         }
 
         if ($p == 3) {
-            my $t = $n / 3;
+            my $t = divint($n, 3);
             return ($t - ($t >> 1));
         }
 
+        my $key = "$n,$p";
+
+        return $cache{$key}
+            if exists $cache{$key};
+
         my $u = 0;
-        my $t = $n / $p;
+        my $t = divint($n, $p);
 
         for (my $q = 2 ; $q < $p ; $q = next_prime($q)) {
 
@@ -46,12 +53,12 @@ ($n, $p)
             }
         }
 
-        $t - $u;
+        $cache{$key} = $t - $u;
     }->($n * $p, $p);
 }
 
 foreach my $p (@{primes(30)}) {
-    say "Φ(10^n, $p) for n <= 10: [", join(', ', map { rough_count(powint(10, $_), $p) } 0 .. 10), "]";
+    say "Φ(10^n, $p) for n <= 10: [", join(', ', map { my_rough_count(powint(10, $_), $p) } 0 .. 10), "]";
 }
 
 __END__

Math/count_of_rough_numbers_recursive.pl

@@ -0,0 +1,58 @@
+#!/usr/bin/perl
+
+# Daniel "Trizen" Șuteu
+# Date: 05 September 2025
+# https://github.com/trizen
+
+# Count the number of B-rough numbers <= n.
+
+# See also:
+#   https://en.wikipedia.org/wiki/Rough_number
+
+use 5.020;
+use ntheory qw(:all);
+use experimental qw(signatures);
+
+sub my_rough_count($n, $k) {
+
+    my @P = @{primes($k - 1)};
+
+    return $n if (@P == 0);
+
+    my %cache;
+
+    sub ($n, $a) {
+
+        my $key = "$n,$a";
+
+        return $cache{$key}
+            if exists $cache{$key};
+
+        # Initial count: odd numbers ≤ n
+        my $count = $n - ($n >> 1);
+
+        # Inclusion-Exclusion principle
+        for my $j (1 .. $a - 1) {
+            last if ($P[$j] > $n);
+            $count -= __SUB__->(divint($n, $P[$j]), $j);
+        }
+
+        $cache{$key} = $count;
+    }->($n, scalar @P);
+}
+
+foreach my $p (@{primes(30)}) {
+    say "Φ(10^n, $p) for n <= 10: [", join(', ', map { my_rough_count(powint(10, $_), $p) } 0 .. 10), "]";
+}
+
+__END__
+Φ(10^n,  2) for n <= 10: [1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000, 10000000000]
+Φ(10^n,  3) for n <= 10: [1, 5, 50, 500, 5000, 50000, 500000, 5000000, 50000000, 500000000, 5000000000]
+Φ(10^n,  5) for n <= 10: [1, 3, 33, 333, 3333, 33333, 333333, 3333333, 33333333, 333333333, 3333333333]
+Φ(10^n,  7) for n <= 10: [1, 2, 26, 266, 2666, 26666, 266666, 2666666, 26666666, 266666666, 2666666666]
+Φ(10^n, 11) for n <= 10: [1, 1, 22, 228, 2285, 22857, 228571, 2285713, 22857142, 228571428, 2285714285]
+Φ(10^n, 13) for n <= 10: [1, 1, 21, 207, 2077, 20779, 207792, 2077921, 20779221, 207792207, 2077922077]
+Φ(10^n, 17) for n <= 10: [1, 1, 20, 190, 1917, 19181, 191808, 1918081, 19180820, 191808190, 1918081917]
+Φ(10^n, 19) for n <= 10: [1, 1, 19, 179, 1806, 18053, 180524, 1805251, 18052535, 180525355, 1805253568]
+Φ(10^n, 23) for n <= 10: [1, 1, 18, 170, 1711, 17103, 171021, 1710234, 17102401, 171024023, 1710240224]
+Φ(10^n, 29) for n <= 10: [1, 1, 17, 163, 1634, 16361, 163586, 1635877, 16358819, 163588196, 1635881952]

Math/meissel_lehmer_prime_count.pl

@@ -19,9 +19,7 @@
 my %pi_cache;
 
 # Recursive φ(n, a): numbers <= n not divisible by first a primes
-sub recursive_rough_count ($n, $k) {
-
-    my @P = @{primes($k)};
+sub recursive_rough_count ($n, $P) {
 
     sub ($n, $a) {
 
@@ -30,13 +28,17 @@ ($n, $k)
         return $phi_cache{$key}
           if exists $phi_cache{$key};
 
-        if ($a == 0) {
-            return $phi_cache{$key} = $n;
+        my $count = $n - ($n >> 1);
+
+        foreach my $j (1 .. $a - 1) {
+            my $np = divint($n, $P->[$j]);
+            last if ($np == 0);
+            $count -= __SUB__->($np, $j);
         }
 
-        $phi_cache{$key} = __SUB__->($n, $a - 1) - __SUB__->(divint($n, $P[$a - 1]), $a - 1);
+        $phi_cache{$key} = $count;
       }
-      ->($n, scalar @P);
+      ->($n, scalar @$P);
 }
 
 # P2 correction term
@@ -69,7 +71,7 @@ ($n)
     my $a    = scalar @P;
     my $p_a  = $P[-1];
 
-    my $phi = recursive_rough_count($n, $p_a + 1);
+    my $phi = recursive_rough_count($n, \@P);
     my $p2  = P2($n, $a, $p_a);
 
     my $result = $phi + $a - 1 - $p2;

Math/prime_counting_from_squarefree_almost_primes.pl

@@ -33,7 +33,8 @@ ($k, $n)
 
             forprimes {
                 $count += my_prime_count(divint($n, mulint($m, $_))) - $j++;
-            } $p, $s;
+            }
+            $p, $s;
 
             return;
         }
@@ -49,14 +50,20 @@ ($k, $n)
 
 sub my_prime_count ($n) {
 
-    state $pi_table = [0, 0, 1, 2, 2];      # a larger lookup table helps a lot!
+    state %cache = (    # a larger lookup table helps a lot!
+                     0 => 0,
+                     1 => 0,
+                     2 => 1,
+                     3 => 2,
+                     4 => 2,
+                   );
 
     if ($n < 0) {
         return 0;
     }
 
-    if (defined($pi_table->[$n])) {
-        return $pi_table->[$n];
+    if (exists $cache{$n}) {
+        return $cache{$n};
     }
 
     my $M = powerfree_count($n, 2) - 1;
@@ -65,10 +72,10 @@ ($n)
         $M -= squarefree_almost_prime_count($k, $n);
     }
 
-    return ($pi_table->[$n] //= $M);
+    $cache{$n} //= $M;
 }
 
-foreach my $n (1..7) {    # takes ~1 second
+foreach my $n (1 .. 7) {    # takes ~1 second
     say "pi(10^$n) = ", my_prime_count(10**$n);
 }
 

README.md

@@ -743,6 +743,7 @@ A nice collection of day-to-day Perl scripts.
     * [Count of perfect powers](./Math/count_of_perfect_powers.pl)
     * [Count of prime power](./Math/count_of_prime_power.pl)
     * [Count of rough numbers](./Math/count_of_rough_numbers.pl)
+    * [Count of rough numbers recursive](./Math/count_of_rough_numbers_recursive.pl)
     * [Count of smooth numbers](./Math/count_of_smooth_numbers.pl)
     * [Count of smooth numbers mpz](./Math/count_of_smooth_numbers_mpz.pl)
     * [Count of smooth numbers mpz 2](./Math/count_of_smooth_numbers_mpz_2.pl)