Project Euler 694

题目

Cube-full Divisors

A positive integer \(n\) is considered cube-full, if for every prime \(p\) that divides \(n\), so does \(p^3\). Note that \(1\) is considered cube-full.

Let \(s(n)\) be the function that counts the number of cube-full divisors of \(n\). For example, \(1\), \(8\) and \(16\) are the three cube-full divisors of \(16\). Therefore, \(s(16)=3\).

Let \(S(n)\) represent the summatory function of \(s(n)\), that is \(S(n)=\displaystyle\sum_{i=1}^n s(i)\).

You are given \(S(16) = 19\), \(S(100) = 126\) and \(S(10000) = 13344\).

Find \(S(10^{18})\).

解决方案

令\(N=10^{18},M=\lfloor\sqrt[4]{N}\rfloor.\)计算\(S\)时，我们考虑枚举每一个满立方数\(w\)，那么它们在\(N\)以内出现的次数为\(\left\lfloor\dfrac{N}{w}\right\rfloor\).

通过考察\(n\)分解质因数\(n=\prod_{i=1}^kp_i^{e_i}\)的每个\(e_i\)，可以发现每一个满立方数可以通过\(x^3y^4z^5\)来不重不漏地表出，其中\(x,y\)是两个无平方因子数，并且\(\gcd(x,y)=1.\)

实现时，首先使用筛法将\(M\)以内的所有无平方因子数进行筛选出来。前两层循环分别枚举两个无平方因子数\(z\)和\(y\)。注意判断\(\gcd(z,y)=1\)时，只需要判断\(zy\)是否为无平方因子数即可，接下来直接枚举\(x\)，从而构造出一个满立方因子数\(x^3y^4z^5\)。

代码

from itertools import count
from tools import int_sqrt

N = 10 ** 18
M = int(N ** (1 / 4))
square_free = [1] * (M + 1)
for i in range(2, int_sqrt(M) + 1):
    k = i * i
    for j in range(k, M + 1, k):
        square_free[j] = 0
ans = 0
for z in count(1, 1):
    z5 = z ** 5
    if z ** 5 > N:
        break
    if square_free[z]:
        for y in count(1, 1):
            y4 = y ** 4
            mul = y4 * z5
            if mul > N:
                break
            if square_free[y] and square_free[z * y]:
                for x in count(1, 1):
                    x3 = x ** 3
                    if mul * x3 > N:
                        break
                    ans += N // (mul * x3)
print(ans)