A371185 Indices of the cubefull numbers in the sequence of powerful numbers.

1, 3, 5, 7, 8, 11, 13, 17, 18, 23, 25, 26, 30, 34, 38, 41, 42, 45, 49, 54, 55, 61, 63, 72, 77, 78, 80, 82, 83, 87, 89, 93, 99, 105, 106, 113, 115, 116, 127, 128, 130, 137, 140, 148, 151, 153, 161, 164, 166, 179, 185, 186, 188, 192, 196, 201, 206, 221, 227, 234

Offset: 1

Author

Amiram Eldar, Mar 14 2024

Keywords

• nonn
• easy

Examples

			The first 5 powerful numbers are 1, 4, 8, 9, and 16. The 1st, 3rd, and 5th, 1, 8, and 16, are also cubefull numbers. Therefore, the first 3 terms of this sequence are 1, 3, and 5.
		

Links

• Amiram Eldar, Table of n, a(n) for n = 1..10000
• Index entries for sequences related to powerful numbers.

Crossrefs

Cf. A001694, A036966.

Cf. A090699, A362974.

Similar sequences: A361936, A371186.

Programs

Mathematica

powQ[n_, emin_] := n == 1 || AllTrue[FactorInteger[n], Last[#] >= emin &]; Position[Select[Range[20000], powQ[#, 2] &], _?(powQ[#1, 3] &), Heads -> False] // Flatten

PARI

ispow(n, emin) = n == 1 || vecmin(factor(n)[, 2]) >= emin;
lista(kmax) = {my(f, c = 0); for(k = 1, kmax, if(ispow(k, 2), c++; if(ispow(k, 3), print1(c, ", "))));}

Python

from math import isqrt, gcd
from sympy import mobius, integer_nthroot, factorint
def A371185(n):
    def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
    def bisection(f,kmin=0,kmax=1):
        while f(kmax) > kmax: kmax <<= 1
        while kmax-kmin > 1:
            kmid = kmax+kmin>>1
            if f(kmid) <= kmid:
                kmax = kmid
            else:
                kmin = kmid
        return kmax
    def f(x):
        c = n+x
        for w in range(1,integer_nthroot(x,5)[0]+1):
            if all(d<=1 for d in factorint(w).values()):
                for y in range(1,integer_nthroot(z:=x//w**5,4)[0]+1):
                    if gcd(w,y)==1 and all(d<=1 for d in factorint(y).values()):
                        c -= integer_nthroot(z//y**4,3)[0]
        return c
    c, l, m = 0, 0, bisection(f,n,n)
    j = isqrt(m)
    while j>1:
        k2 = integer_nthroot(m//j**2,3)[0]+1
        w = squarefreepi(k2-1)
        c += j*(w-l)
        l, j = w, isqrt(m//k2**3)
    c += squarefreepi(integer_nthroot(m,3)[0])-l
    return c # Chai Wah Wu, Sep 12 2024

Formula

A001694(a(n)) = A036966(n).

a(n) ~ c * n^(3/2), where c = A090699 / A362974^(3/2) = 0.216089803749...