作为 min-25 筛的模板题，本题可以使用 min-25 筛 / 洲阁筛 等亚线性积性函数求和的方法在 O\left(\frac{n^{\frac{3}{4}}}{\log n}\right) 的时间复杂度内解决。树状数组优化后的 min-25 筛可以做到 O\left(\frac{n^{\frac{3}{4}}}{\log n}\right)

算法一

这种解法的思路参考自 zzq's blog .

定义：Powerful Number

若正整数 x 满足它所有质因子的幂次都 \geq 2 ，则称 x 是一个 Powerful Number .

性质

[1,n] 间的 Powerful Number 有 O(\sqrt n) 个。

证明：

由于每个质因子间是独立的，显然所有 Powerful Number 都可以表示成 a^2b^3 的形式。考虑枚举 a ，则 n 以内这种数的数量为

O\left(\sum_{a=1}^{\sqrt n}\left(\frac n{a^2}\right)^{\frac 1 3}\right)

=O\left(\int_1^{\sqrt n}\left(\frac n{x^2}\right)^{\frac 1 3}dx\right)

=O(\sqrt n)

这个性质保证了该算法的复杂度。

本题思路

设 f(x) 为要求的积性函数。对于质数 p ，有 f(p)=p(p-1) . 设数论函数 g(x)=x\varphi(x) ，显然 g(x) 是积性函数，且满足 g(p)=f(p) . 设数论函数 h(x) 满足 f=g*h ，其中 * 表示狄利克雷卷积。则有

f(p)=g(1)h(p)+g(p)h(1)

因为 f,g 均为积性函数，所以 h 也为积性函数。所以 g(1)=h(1)=1 ，又因为 g(p)=f(p) ，可以得出 h(p)=0 . 由于 h 是积性函数，所以仅当 x 是 Powerful Number 时，h(x) 才可能不为 0 .

考虑答案。答案为

\sum_{i=1}^{n}f(i)=\sum_{ij\leq n}g(i)h(j)=\sum_{i=1}^{n}h(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}g(j)

g 本身是简单的积性函数，可以使用杜教筛在 O(n^{\frac{2}{3}}) 的时间内筛出所有 \lfloor\frac{n}{i}\rfloor 位置上的前缀和。接下来只需考虑 i 为 Powerful Number 时，h(i) 如何计算。因为 h 是积性函数，所以我们只需要求出 h(p^k) . 由狄利克雷卷积，有

f(p^k)=\sum_{i=0}^{k}g(p^i)h(p^{k-i})

根据欧拉函数的有关性质，g(p^k)=p^k\varphi(p^k)=p^k(p-1)p^{k-1}=(p-1)p^{2k-1} . 将 f,g 展开，得

p^k(p^k-1)=\sum_{i=0}^{k}(p-1)p^{2i-1}h(p^{k-i})

归纳可证 h(p^k)=(k-1)(p^{k+1}-p^k) .

搜索 [1,\sqrt n] 内质数，并从 2 开始枚举它们的幂次，即可得到 [1,n] 内的所有 Powerful Number . 在这个过程中，我们可以同时得到这些数的 h 函数值。事实上，如果 h(p^k) 不像本题一样可以快速求出，因为 \sqrt n 以内的质数只有 O\left(\frac{\sqrt n}{\log n}\right) 个且 p 的幂次是 O(\log n) 的，也可以暴力地计算出它们。

这个做法的时间复杂度是 O(n^{\frac{2}{3}}) 的，其瓶颈在于杜教筛 g(x) 的前缀和。渐进意义上这个复杂度低于 min-25 筛的复杂度 O\left(\frac{n^{\frac{3}{4}}}{\log n}\right)，但在本题的数据范围下表现不如后者。在洛谷上不开启 O2 优化的情况下，最慢的测试点用了 1.5s 的时间，明显是不如实现正常的 min-25 筛的。但是这个算法也有它的优点：实现简单；当 g 的前缀和可以快速求出时复杂度可以降低到 O(\sqrt n) ，有明显优势。

代码

#include <cstdio>
#include <cmath>
#include <cstring>
typedef long long ll;
const int MAXN=5500000;
const int MOD=1E9+7;
const int INV2=500000004, INV6=166666668;
ll n; int lim, cp, ans;
int pri[MAXN], phi[MAXN];
bool flag[MAXN];
int g1[MAXN], g2[MAXN];
void prework()
{
    memset(flag, 0, sizeof flag);
    phi[1]=1;
    for (int i=2; i<=lim; i++)
    {
        if (!flag[i]) pri[++cp]=i, phi[i]=i-1;
        for (int j=1; j<=cp&&i*pri[j]<=lim; j++)
        {
            flag[i*pri[j]]=1;
            if (i%pri[j]==0)
            {
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
    for (int i=1; i<=lim; i++) g1[i]=(g1[i-1]+1ll*i*phi[i])%MOD;
}
int djs(ll m)
{
    if (m<=lim) return g1[m];
    if (g2[n/m]) return g2[n/m];
    int res=(m%MOD)*((m+1)%MOD)%MOD*((2*m+1)%MOD)%MOD*INV6%MOD;
    for (ll i=2, j; i<=m; i=j+1)
    {
        j=m/(m/i);
        res=(res-((i+j)%MOD)*((j-i+1)%MOD)%MOD*INV2%MOD*djs(m/i)%MOD+MOD)%MOD;
    }
    return g2[n/m]=res;
}
void dfs(int k, ll m, int h)
{
    if (k>cp||m*pri[k]>n)
    {
        if (n/m<=lim) ans=(ans+1ll*h*g1[n/m])%MOD;
        else ans=(ans+1ll*h*g2[n/(n/m)])%MOD;
        return;
    }
    ll p=1ll*pri[k]*pri[k];
    dfs(k+1, m, h);
    for (int e=2; m*p<=n; p*=pri[k], e++)
        dfs(k+1, m*p, p%MOD*(pri[k]-1)%MOD*(e-1)%MOD*h%MOD);
}
int main()
{
//  freopen("P5325.in", "r", stdin);
//  freopen("P5325.out", "w", stdout);
    scanf("%lld", &n);
    lim=pow(n, 2.0/3);
    prework(); djs(n);
    dfs(1, 1, 1);
    printf("%d\n", ans);
    return 0;
}

算法二