题面
解析
定义$f(i)$为$i$的次大因数,特别地$f(1)=0$
那么我们就是去求$$\sum_{i=1}^{n}\sum_{j=1}^{n}f^{m}(gcd(i, j))$$
这种东西的套路就是枚举$gcd$然后用欧拉函数处理, \begin{align*}\sum_{i=1}^{n}\sum_{j=1}^{n}f^{m}(gcd(i, j)) &= \sum_{i=1}^{n}\sum_{j=1}^{\left \lfloor \frac{n}{i}\right \rfloor} \sum_{k=1}^{\left \lfloor \frac{n}{i}\right \rfloor}[gcd(j, k)==1]f(i)\\&=\sum_{i=1}^{n}f^{m}(i)((2\sum_{j=1}^{\left \lfloor \frac{n}{i}\right \rfloor}\varphi(j)) - 1)\end{align*}
求$\varphi(i)$的前缀和用杜教筛即可,问题在于如何求$f(i)$的前缀和
考虑$min25$筛的过程,它把所有数对前缀和的贡献分为了质数的贡献与合数的贡献,这道题中质数的贡献显然是1,因此,可以构造一个完全积性函数$c(i)=1$处理质数的贡献。接下来处理合数的贡献,考虑在$min25$的第一步处理过程中,每一个合数都会被它最小的质因子$minp$筛去。因此在枚举第$i$个质数$pri[i]$时就会有以$pri[i]$为最小质因子的所有合数的贡献。定义完全积性函数$b(i)=i^{m}$,$g$数组为筛法第一步中处理的函数,$h(i)$表示$1$到$i$中合数对答案的贡献,用公式写来,就是这样:$$g(n, i) = g(n, i)-b(pri[i])*(g(\left \lfloor \frac{n}{pri[i]}\right \rfloor, i-1) - \sum_{j=1}^{i-1}b(pri[j]))\\h(n)=h(n)+g(\left \lfloor \frac{n}{pri[i]}\right \rfloor, i-1) - \sum_{j=1}^{i-1}b(pri[j])$$
$b(i)$的前缀和在线性筛时预处理即可,问题在于如何预处理$g(n, 0)=\sum_{i=1}^{n}i^{m}$,考虑用第二类斯特林数,原式变为:$$\sum_{i=1}^{n}\sum_{j=0}^{m}\begin{Bmatrix}m\\j\end{Bmatrix}\cdot j!\cdot \binom{i}{j}$$
上式$j$可以直接枚举到$m$, 若$j$大于$m$,$\begin{Bmatrix}m\\j\end{Bmatrix}=0$;若$j$大于$i$,$C_{i}^{j}=0$,两种情况均不会对求和造成影响。
套路交换求和号,变为:$$\sum_{j=0}^{m}\begin{Bmatrix}m\\j\end{Bmatrix}\cdot j!\sum_{i=0}^{n}\binom{i}{j}$$
注意到$\sum_{i=0}^{n}\binom{i}{j}=\binom{n+1}{j+1}$,代入得:$$\sum_{j=0}^{m}\begin{Bmatrix}m\\j\end{Bmatrix}\cdot j! \cdot \binom{n+1}{j+1}$$
展开$C_{i}^{j}=0$, 将$j!$代入,最终变为$$\sum_{j=0}^{m}\begin{Bmatrix}m\\j\end{Bmatrix}\frac{\prod_{i=n-j+1}^{n+1}i}{j+1}$$
注意$j+1$可能没有逆元,需要先把$j+1$除了才能取模
最后做一次数论分块统计答案
在做这道题时才意识到的细节,在数论分块时,存下来的值恰好为每一块区间的右端点...
代码:
#include<bits/stdc++.h> using namespace std; typedef long long ll; typedef unsigned int ui; const int maxn = 1000005; int n, m, mp[maxn], Mp[maxn], a[maxn]; ui s2[60][60], b[maxn], c[maxn]; ll pri[maxn]; int phi[maxn], cnt; ui f[maxn], g[maxn], s[maxn]; bool notp[maxn]; map<int, ui> F; ui qpow(ui x, int y) { ui ret = 1; while(y) { if(y&1) ret = ret * x; x = x * x; y >>= 1; } return ret; } void Euler() { phi[1] = 1; for(int i = 2; i <= 1000000; ++i) { if(!notp[i]) { pri[++cnt] = i; phi[i] = i - 1; s[cnt] = s[cnt-1] + qpow(i, m); } for(int j = 1; j <= cnt; ++j) { if(pri[j] * i > 1000000) break; notp[pri[j]*i] = 1; if(i % pri[j] == 0) { phi[pri[j]*i] = pri[j] * phi[i]; break; } phi[pri[j]*i] = (pri[j] - 1) * phi[i]; } } for(int i = 1; i <= 1000000; ++i) f[i] = f[i-1] + phi[i]; } void init() { s2[0][0] = 1; for(int i = 1; i <= m; ++i) for(int j = 1; j <= i; ++j) s2[i][j] = s2[i-1][j-1] + j * s2[i-1][j]; } int get_id(int x) { return x <= 1000000? mp[x]: Mp[n/x]; } int gcd(int a, int b) { return b == 0? a: gcd(b, a % b); } int stak[60]; ui get_mi(int x) { ui ret = 0, mul; int sj = min(x, m), t, gc; for(int i = 1; i <= sj; ++i) { mul = s2[m][i]; for(int j = 1; j <= i + 1; ++j) stak[j] = x + 2 - j; t = i + 1; for(int j = 1; j <= i + 1; ++j) { gc = gcd(stak[j], t); stak[j] /= gc; t /= gc; if(t == 1) break; } for(int j = 1; j <= i + 1; ++j) mul = mul * stak[j]; ret += mul; } return ret; } ui dfs(int x) { if(x <= 1000000) return f[x]; if(F.find(x) != F.end()) return F[x]; ll tmp = (1LL * x * (x + 1)) >> 1; ui ret = (ui)tmp; for(int l = 2, r; l <= x; l = r + 1) { r = x / (x / l); ret -= (r - l + 1) * dfs(x / l); } return F[x] = ret; } int main() { scanf("%d%d", &n, &m); Euler(); init(); int t, tot = 0, id, now; for(int l = 1, r; l <= n; l = r + 1) { t = n / l; r = n / t; a[++tot] = t; if(t <= 1000000) mp[t] = tot; else Mp[n/t] = tot; b[tot] = get_mi(t) - 1; c[tot] = t - 1; } for(int i = 1; i <= cnt && pri[i] * pri[i] <= n; ++i) for(int j = 1; j <= tot && pri[i] * pri[i] <= a[j]; ++j) { id = get_id(a[j] / pri[i]); b[j] -= (b[id] - s[i-1]) * (s[i] - s[i-1]); g[j] += b[id] - s[i-1]; c[j] -= c[id] - i + 1; } ui ans = 0; for(int l = 1, r; l <= n; l = r + 1) { r = n / (n / l); now = get_id(r); id = get_id(l - 1); ans += ((g[now] + c[now]) - (g[id] + c[id])) * (2 * dfs(n / l) - 1); } printf("%u", ans); return 0; }
原文地址:https://www.cnblogs.com/Joker-Yza/p/11969074.html