count_of_smooth_numbers.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 19 May 2020
# https://github.com/trizen

# Count the number of B-smooth numbers <= n.

# See also:
#   https://en.wikipedia.org/wiki/Smooth_number

use 5.036;
use ntheory qw(:all);

sub my_smooth_count ($n, $p) {

    if ($p == 2) {
        return 1 + logint($n, 2);
    }

    my $q = prev_prime($p);

    my $count = 0;
    foreach my $k (0 .. logint($n, $p)) {
        $count += __SUB__->(divint($n, powint($p, $k)), $q);
    }

    return $count;
}

foreach my $p (@{primes(50)}) {
    say "Ψ(10^n, $p) for n <= 10: [", join(', ', map { my_smooth_count(powint(10, $_), $p) } 0 .. 10), "]";
}

__END__
Ψ(10^n, 2) for n <= 10: [1, 4, 7, 10, 14, 17, 20, 24, 27, 30, 34]
Ψ(10^n, 3) for n <= 10: [1, 7, 20, 40, 67, 101, 142, 190, 244, 306, 376]
Ψ(10^n, 5) for n <= 10: [1, 9, 34, 86, 175, 313, 507, 768, 1105, 1530, 2053]
Ψ(10^n, 7) for n <= 10: [1, 10, 46, 141, 338, 694, 1273, 2155, 3427, 5194, 7575]
Ψ(10^n, 11) for n <= 10: [1, 10, 55, 192, 522, 1197, 2432, 4520, 7838, 12867, 20193]
Ψ(10^n, 13) for n <= 10: [1, 10, 62, 242, 733, 1848, 4106, 8289, 15519, 27365, 45914]
Ψ(10^n, 17) for n <= 10: [1, 10, 67, 287, 945, 2579, 6179, 13389, 26809, 50351, 89679]
Ψ(10^n, 19) for n <= 10: [1, 10, 72, 331, 1169, 3419, 8751, 20198, 42950, 85411, 160626]
Ψ(10^n, 23) for n <= 10: [1, 10, 76, 369, 1385, 4298, 11654, 28434, 63768, 133440, 263529]
Ψ(10^n, 29) for n <= 10: [1, 10, 79, 402, 1581, 5158, 14697, 37627, 88415, 193571, 399341]
Ψ(10^n, 31) for n <= 10: [1, 10, 82, 434, 1778, 6070, 18083, 48366, 118599, 270648, 581272]
Ψ(10^n, 37) for n <= 10: [1, 10, 84, 461, 1958, 6952, 21535, 59867, 152482, 361173, 804369]
Ψ(10^n, 41) for n <= 10: [1, 10, 86, 485, 2129, 7833, 25133, 72345, 190767, 467495, 1076462]
Ψ(10^n, 43) for n <= 10: [1, 10, 88, 508, 2300, 8740, 28955, 86086, 234423, 592949, 1408465]
Ψ(10^n, 47) for n <= 10: [1, 10, 90, 529, 2463, 9639, 32876, 100688, 282397, 735425, 1797897]