Therefore we can solve the problem if we can find all the cube-full divisors less than 10^18.

First I noticed that (2*3*5*7*11*13*17*19)^3 > 10^18 so a cube-full divisor will have a maximum of 8 prime factors.

I brute force searched for all of them by first finding all prime powers for primes less than (10^18)^(1/3) = 10^6 as this is the largest prime factor a cube-full divisor can have, that is d = [2^3, 2^4, ..., 3^3, 3^4, ...].

Now by multiplying every element in d with every element in d we get all cube-full divisors with 1 or 2 prime factors, store this in a list called d1, then by multiplying this list with itself again we get all cube-full divisors with 1,2,3, and 4 divisors, and finally doing this one more time we get all cube-full divisors with 1,2,3,4,5,6,7, and 8 prime factors, which is what we need.

My code takes around ~130s to run, the forum is extremely informative on how to get fast recursive solutions.

EDIT: After solving Problem 641 I came back to this one to solve it with a recursive solution, now it takes 8 seconds to solve the problem

One interesting mathematical point is that we can notice that S(n)/n is clearly converging to some number:

S(10^6) = 1,339,485, S(10^7) = 13397119, ..., S(10^14) = 133978415161307, why is this?