k-powerful_numbers_in_range.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 28 February 2021
# https://github.com/trizen

# Fast recursive algorithm for generating all the k-powerful numbers in a given range [a,b].
# 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

func powerful_numbers(A, B, k=2) {

    var powerful = []

    func (m,r) {

        var from = 1
        var upto = iroot(idiv(B, m), r)

        if (r <= k) {

            if (A > m) {

                # Optimization by Dana Jacobsen (from Math::Prime::Util::PP)
                with (idiv_ceil(A,m)) {|l|
                    if ((l >> r) == 0) {
                        from = 2
                    }
                    else {
                        from = l.iroot(r)
                        from++ if (ipow(from, r) != l)
                    }
                }
            }

            return nil if (from > upto)

            for j in (from .. upto) {
                powerful << (m * ipow(j, r))
            }

            return nil
        }

        for j in (from .. upto) {

            j.is_coprime(m) || next
            j.is_squarefree || next

            __FUNC__(m * ipow(j, r), r-1)
        }
    }(1, 2*k - 1)

    powerful.sort
}

var a = 1e5.irand
var b = 1e7.irand

for k in (2..5) {
    say "Testing: k = #{k}"
    assert_eq(
        powerful_numbers(a, b, k),
        k.powerful(a, b)
    )
}

say powerful_numbers(1e6 - 1e4, 1e6, 2)     #=> [990025, 990125, 990584, 991232, 992016, 994009, 995328, 996004, 996872, 998001, 998784, 1000000]