P5325 【模板】Min_25 筛
lovely_nst
·
·
题解
P5325 【模板】Min_25 筛（PN 筛）
定义 \tt{Powerful\ Number} 为每个质因子项次数至少为 2 的数。
显然的一点，任意一个 \tt{Powerful\ Number} 可以被表示为 a^2b^3。而在 n 以内的 \tt{Powerful\ Number} 的数量即为：
\sum_{a^2b^3\le n} 1\\
=\sum_{a=1}^{\sqrt n}\sqrt[3]{\lfloor\frac n{a^2}\rfloor}\\
\approx \int_1^{\sqrt{n}} \frac {\sqrt[3]{n}}{\sqrt[3]{x}} \text{d} x\\
\approx (3\sqrt[6]{n}-3)\sqrt[3]n\\
=O(\sqrt{n})
这是一个重要的结论，即在 n 以内的 \tt{Powerful\ Number} 的数量是 \sqrt n 级别的。
回到本题，要求积性函数 f(p^k)=p^k(p^k-1)（p 为质数）的前缀和。
设 g(x)=\varphi(x)x。因为 p 为质数，因此有 \varphi(p)=p-1，同时便满足 f(p)=g(p)。
设 f=g*h（狄利克雷卷积），则：
f(x)=\sum_{d|x} g(d)h(\frac xd)\\
f(p)=g(1)h(p)+g(p)h(1)\\
\because f(p)=g(p),g(1)=h(1)=1,\\
\therefore f(p)=h(p)+f(p),h(p)=0.
同时由上也可以得出 h 为积性函数。因此只有当 x 为 \tt{Powerful\ Number} 时才可能满足 h(x)\ne0。继续推式子：
f(p^k)=(g*h)(p^k)\\
=\sum_{d|p^k} g(d)h(\frac {p^k}d)\\
=\sum_{i=0}^k g(p^i)h(p^{k-i})\\
g(p^i)=\varphi(p^i)p^i=(p^i-p^{i-1})p^i=(p-1)p^{2i-1}\\
g(1)=1\\
p^k(p^k-1)=h(p^k)+\sum_{i=1}^k(p-1)p^{2i-1}h(p^{n-i})\\
p^k(p^k-1)=h(p^k)+\sum_{i=0}^{k-1}(p-1)p^{2k-2i-1}h(p^i)\\
h(p^k)=p^k(p^k-1)-\sum_{i=0}^{k-1}(p-1)p^{2k-2i-1}h(p^i)\\
于是我们得到了 h(p^k) 的递推式，直接递推即可算出。
$$
\sum_{ij\le n}h(i)g(j)\\
=\sum_{i=1}^nh(i)\sum_{j=1}^{\lfloor\frac ni\rfloor}g(j)
$$
后半部分可以使用杜教筛算出，而又因为 $\tt{Powerful\ Number}$ 的数量在 $O(\sqrt n)$ 级别，因此只会有 $O(\sqrt n)$ 个 $h(i)$ 有值，直接算即可，时间复杂度为 $O(n^{\frac34})$。
## AC Code
```cpp
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 2e6 + 5;
const int mod = 1e9 + 7 , inv3 = (mod + 1) / 3;
int n , m , k;
int power (int x , int y)
{
int ans = 1;
while (y)
{
if (y & 1) ans = ans * x % mod;
x = x * x % mod;
y >>= 1;
}
return ans;
}
int Sum1 (int n)
{
n %= mod;
return n * (n + 1) / 2 % mod;
}
int Sum2 (int n)
{
n %= mod;
return n * (n + 1) / 2 % mod * (2 * n + 1) % mod * inv3 % mod;
}
namespace PhiSum
{
int sum[N];
unordered_map <int , int> f;
int S (int x)
{
if (x < N) return sum[x];
if (f.count (x)) return f[x];
int ans = 0 , nxt , v;
for (register int i = 2;i <= x;)
{
v = x / i , nxt = x / v;
ans = (ans + S (v) * (Sum1 (nxt) - Sum1 (i - 1))) % mod;
i = nxt + 1;
}
return f[x] = (Sum2 (x) - ans) % mod;
}
bool vis[N];
int pri[N] , ph[N];
void init ()
{
sum[1] = 1;
for (int i = 2;i < N;i ++)
{
if (!vis[i])
pri[++ m] = i , ph[i] = i - 1;
for (int j = 1 , k;j <= m && (k = pri[j] * i) < N;j ++)
{
vis[k] = 1;
if (i % pri[j] == 0)
{
ph[k] = ph[i] * pri[j];
break;
}
ph[k] = ph[i] * ph[pri[j]];
}
sum[i] = sum[i - 1] + ph[i] * i;
sum[i] %= mod;
}
}
}
int a[N] , h[N] , s[N];
signed main ()
{
ios::sync_with_stdio (0);
cin.tie (0) , cout.tie (0);
cin >> n;
PhiSum::init ();
a[++ k] = 1;
h[k] = s[0] = 1;
for (int i = 1;i <= m;i ++)
{
int pr = PhiSum::pri[i] , Mul2 = pr * pr;
if (Mul2 > n) break;
__int128 Mul = Mul2;
for (int j = 2;Mul <= n;j ++)
{
s[j] = Mul % mod * (Mul - 1) % mod;
for (int l = 0;l < j;l ++)
s[j] = (s[j] - (pr - 1) * power (pr , 2 * j - 2 * l - 1) % mod * s[l] % mod);
s[j] %= mod;
Mul *= pr;
}
int K = k;
for (int j = 1;j <= K;j ++)
{
if (a[j] * pr > n) continue;
int mul = a[j] * Mul2;
int cc = 2;
while (mul <= n)
{
a[++ k] = mul;
h[k] = h[j] * s[cc] % mod;
cc ++;
mul *= pr;
}
}
}
int ans = 0;
for (int i = 1;i <= k;i ++)
ans = (ans + h[i] * PhiSum::S (n / a[i])) % mod;
cout << (ans + mod) % mod;
return 0;
}
```