k-odd-powerful_numbers_in_range.sf

#!/usr/bin/ruby

# Author: Trizen
# Date: 28 February 2021
# Edit: 23 February 2024
# https://github.com/trizen

# Fast recursive algorithm for generating all the odd 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

# See also:
#   https://oeis.org/A062739

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

    var odd_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) {
                odd_powerful << (m * ipow(j, r)) if j.is_odd
            }

            return nil
        }

        for j in (from .. upto) {

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

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

    odd_powerful.sort
}

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

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

say odd_powerful_numbers(1e6 - 2e4, 1e6, 2)     #=> [982081, 984987, 985527, 986049, 990025, 990125, 994009, 998001]