count_of_rough_numbers.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 20 July 2020
# 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 integer;
use ntheory qw(:all);
use experimental qw(signatures);

sub rough_count ($n, $p) {
    sub ($n, $p) {

        if ($p > sqrtint($n)) {
            return 1;
        }

        if ($p == 2) {
            return ($n >> 1);
        }

        if ($p == 3) {
            my $t = $n / 3;
            return ($t - ($t >> 1));
        }

        my $u = 0;
        my $t = $n / $p;

        for (my $q = 2 ; $q < $p ; $q = next_prime($q)) {

            my $v = __SUB__->($t - ($t % $q), $q);

            if ($v == 1) {
                $u += prime_count($q, $p - 1);
                last;
            }
            else {
                $u += $v;
            }
        }

        $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), "]";
}

__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]