sidef-scripts/Math/count_of_k-powerful_numbers.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 11 February 2020
# https://github.com/trizen

# Fast recursive algorithm for counting the number of k-powerful numbers <= n.
# A positive integer n is considered k-powerful, if for every prime p that divides n, so does p^k.

# Example:
#   2-powerful = a^2 * b^3,             for a,b >= 1
#   3-powerful = a^3 * b^4 * c^5,       for a,b,c >= 1
#   4-powerful = a^4 * b^5 * c^6 * d^7, for a,b,c,d >= 1

# OEIS:
#   https://oeis.org/A001694 -- 2-powerful numbers
#   https://oeis.org/A036966 -- 3-powerful numbers
#   https://oeis.org/A036967 -- 4-powerful numbers
#   https://oeis.org/A069492 -- 5-powerful numbers
#   https://oeis.org/A069493 -- 6-powerful numbers

# See also:
#   https://oeis.org/A118896 -- Number of powerful numbers <= 10^n.

func powerful_count(n, k=2) {

    var count = 0

    func (m,r) {

        var t = iroot(idiv(n,m), r)

        if (r <= k) {
            count += t
            return nil
        }

        t.each_squarefree {|a|
            a.is_coprime(m) || next
            __FUNC__(m * a**r, r-1)
        }

    }(1, 2*k - 1)

    return count
}

for k in (2..10) {
    assert_eq(powerful_count(10**k, k), k.powerful_count(10**k))
    printf("Number of %2d-powerful <= 10^j: {%s}\n", k, (8+k).of {|j| powerful_count(10**j, k) }.join(', '))
}

__END__
Number of  2-powerful <= 10^j: {1, 4, 14, 54, 185, 619, 2027, 6553, 21044, 67231, 214122, 680330, 2158391, 6840384}
Number of  3-powerful <= 10^j: {1, 2, 7, 20, 51, 129, 307, 713, 1645, 3721, 8348, 18589, 41136, 90619, 198767}
Number of  4-powerful <= 10^j: {1, 1, 5, 11, 25, 57, 117, 235, 464, 906, 1741, 3312, 6236, 11654, 21661, 40049}
Number of  5-powerful <= 10^j: {1, 1, 3, 8, 16, 32, 63, 117, 211, 375, 659, 1153, 2000, 3402, 5770, 9713, 16266}
Number of  6-powerful <= 10^j: {1, 1, 2, 6, 12, 21, 38, 70, 121, 206, 335, 551, 900, 1451, 2326, 3706, 5853, 9167}
Number of  7-powerful <= 10^j: {1, 1, 1, 4, 10, 16, 26, 46, 77, 129, 204, 318, 495, 761, 1172, 1799, 2740, 4128, 6200}
Number of  8-powerful <= 10^j: {1, 1, 1, 3, 8, 13, 19, 32, 52, 85, 135, 211, 315, 467, 689, 1016, 1496, 2191, 3214, 4653}
Number of  9-powerful <= 10^j: {1, 1, 1, 2, 6, 11, 16, 24, 38, 59, 94, 145, 217, 317, 453, 644, 919, 1308, 1868, 2651, 3745}
Number of 10-powerful <= 10^j: {1, 1, 1, 1, 5, 9, 14, 21, 28, 43, 68, 104, 155, 227, 322, 447, 621, 858, 1192, 1651, 2279, 3152}