【题目链接】
【思路要点】
令 f ( i ) = ( i M i n ( i ) ) k ( i > 1 ) f(i)=(\frac{i}{Min(i)})^k(i>1) f(i)=(Min(i)i)k(i>1) ,即 f ( i ) f(i) f(i) 表示 i i i 次大的因子的 k k k 次方,特别规定 f ( 1 ) = 0 f(1)=0 f(1)=0 。
那么原式即为 ∑ i = 1 N ∑ j = 1 N f ( g c d ( i , j ) ) \sum_{i=1}^{N}\sum_{j=1}^{N}f(gcd(i,j)) ∑i=1N∑j=1Nf(gcd(i,j)) 。
考虑枚举 g c d ( i , j ) gcd(i,j) gcd(i,j) ,那么原式
= ∑ g = 1 N f ( g ) ∑ i = 1 ⌊ N g ⌋ ∑ j = 1 ⌊ N g ⌋ ϵ ( g c d ( i , j ) ) =\sum_{g=1}^{N}f(g)\sum_{i=1}^{\lfloor\frac{N}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{g}\rfloor}\epsilon(gcd(i,j)) =∑g=1Nf(g)∑i=1⌊gN⌋∑j=1⌊gN⌋ϵ(gcd(i,j))
= ∑ g = 1 N f ( g ) ∑ i = 1 ⌊ N g ⌋ ∑ j = 1 ⌊ N g ⌋ ∑ d / i , d / j μ ( d ) =\sum_{g=1}^{N}f(g)\sum_{i=1}^{\lfloor\frac{N}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{g}\rfloor}\sum_{d/i,d/j}\mu(d) =∑g=1Nf(g)∑i=1⌊gN⌋∑j=1⌊gN⌋∑d/i,d/jμ(d)
= ∑ g = 1 N f ( g ) ∑ d = 1 ⌊ N g ⌋ μ ( d ) ⌊ N g d ⌋ 2 =\sum_{g=1}^{N}f(g)\sum_{d=1}^{\lfloor\frac{N}{g}\rfloor}\mu(d)\lfloor\frac{N}{gd}\rfloor^2 =∑g=1Nf(g)∑d=1⌊gN⌋μ(d)⌊gdN⌋2
= ∑ g d = 1 N ( f ∗ μ ) ( g d ) ⌊ N g d ⌋ 2 =\sum_{gd=1}^{N}(f*\mu)(gd)\lfloor\frac{N}{gd}\rfloor^2 =∑gd=1N(f∗μ)(gd)⌊gdN⌋2
对 ⌊ N g d ⌋ \lfloor\frac{N}{gd}\rfloor ⌊gdN⌋ 整数分块,我们剩余的工作是计算函数 ( f ∗ μ ) (f*\mu) (f∗μ) 在每一个 ⌊ N g d ⌋ \lfloor\frac{N}{gd}\rfloor ⌊gdN⌋ 处的前缀和。
令 g ( N ) = ∑ i = 1 N ( f ∗ μ ) ( i ) g(N)=\sum_{i=1}^{N}(f*\mu)(i) g(N)=∑i=1N(f∗μ)(i) ,注意到 ( f ∗ μ ) ∗ 1 = f ∗ ( μ ∗ 1 ) = f (f*\mu)*1=f*(\mu*1)=f (f∗μ)∗1=f∗(μ∗1)=f ,因此
∑ i = 1 N f ( i ) = ∑ i = 1 N ( f ∗ μ ∗ 1 ) ( i ) = ∑ i = 1 N g ( ⌊ N i ⌋ ) \sum_{i=1}^{N}f(i)=\sum_{i=1}^{N}(f*\mu*1)(i)=\sum_{i=1}^{N}g(\lfloor\frac{N}{i}\rfloor) ∑i=1Nf(i)=∑i=1N(f∗μ∗1)(i)=∑i=1Ng(⌊iN⌋)
g ( N ) = ∑ i = 1 N f ( i ) − ∑ i = 2 N g ( ⌊ N i ⌋ ) g(N)=\sum_{i=1}^{N}f(i)-\sum_{i=2}^{N}g(\lfloor\frac{N}{i}\rfloor) g(N)=∑i=1Nf(i)−∑i=2Ng(⌊iN⌋)
上面的式子显然可以杜教筛,剩余的唯一问题在于计算函数 f f f 在每一个 ⌊ N i ⌋ \lfloor\frac{N}{i}\rfloor ⌊iN⌋ 处的前缀和。
分开考虑质数和合数对 ∑ i = 1 N f ( i ) \sum_{i=1}^{N}f(i) ∑i=1Nf(i) 的影响,有
∑ i = 1 N f ( i ) = c n t ( N ) + ∑ p i s a p r i m e , p ≤ N s u m ( ⌊ N p ⌋ , p ) \sum_{i=1}^{N}f(i)=cnt(N)+\sum_{p\ is\ a\ prime,p≤\sqrt{N}}sum(\lfloor\frac{N}{p}\rfloor,p) ∑i=1Nf(i)=cnt(N)+∑p is a prime,p≤Nsum(⌊pN⌋,p)
其中 c n t ( N ) cnt(N) cnt(N) 表示 N N N 以内质数的个数,
s u m ( N , p ) sum(N,p) sum(N,p) 表示 N N N 以内最小质因子大于等于 p p p 的合数的 k k k 次方和。
用 M i n 25 Min25 Min25 筛解决即可。
还有一个小问题,由于模数非质数,求解自然数的 k k k 次方和是需要用第二类斯特林数转化一下。
∑ i = 1 N i k = ∑ i = 1 N ∑ j = 1 k S ( k , j ) ∗ i j ‾ \sum_{i=1}^{N}i^k=\sum_{i=1}^{N}\sum_{j=1}^{k}S(k,j)*i^{\underline{j}} ∑i=1Nik=∑i=1N∑j=1kS(k,j)∗ij
= ∑ j = 1 k S ( k , j ) ∑ i = 1 N i j ‾ =\sum_{j=1}^{k}S(k,j)\sum_{i=1}^{N}i^{\underline{j}} =∑j=1kS(k,j)∑i=1Nij
= ∑ j = 1 k S ( k , j ) ∗ j ! ∑ i = 1 N ( i j ) =\sum_{j=1}^{k}S(k,j)*j!\sum_{i=1}^{N}\binom{i}{j} =∑j=1kS(k,j)∗j!∑i=1N(ji)
= ∑ j = 1 k S ( k , j ) ∗ j ! ( N + 1 j + 1 ) =\sum_{j=1}^{k}S(k,j)*j!\binom{N+1}{j+1} =∑j=1kS(k,j)∗j!(j+1N+1)
= ∑ j = 1 k S ( k , j ) ( N + 1 ) j + 1 ‾ j + 1 =\sum_{j=1}^{k}S(k,j)\frac{(N+1)^{\underline{j+1}}}{j+1} =∑j=1kS(k,j)j+1(N+1)j+1
时间复杂度 O ( N 3 4 + K 2 N ) O(N^{\frac{3}{4}}+K^2\sqrt{N}) O(N43+K2N) 。
【代码】
#include<bits/stdc++.h> using namespace std; const int MAXN = 100005; const int MAXK = 55; template <typename T> void chkmax(T &x, T y) {x = max(x, y); } template <typename T> void chkmin(T &x, T y) {x = min(x, y); } template <typename T> void read(T &x) { x = 0; int f = 1; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == '-') f = -f; for (; isdigit(c); c = getchar()) x = x * 10 + c - '0'; x *= f; } template <typename T> void write(T x) { if (x < 0) x = -x, putchar('-'); if (x > 9) write(x / 10); putchar(x % 10 + '0'); } template <typename T> void writeln(T x) { write(x); puts(""); } unsigned n, k, tot, limit, pre[MAXN], prime[MAXN], f[MAXN], primek[MAXN]; unsigned m, val[MAXN], home1[MAXN], home2[MAXN], cnt[MAXN], sum[MAXN]; unsigned pref[MAXN], s[MAXK][MAXK], ans[MAXN]; unsigned power(unsigned x, unsigned y) { if (y == 0) return 1; unsigned tmp = power(x, y / 2); if (y % 2 == 0) return tmp * tmp; else return tmp * tmp * x; } unsigned calc(unsigned n) { unsigned ans = 0; for (unsigned i = 1; i <= k && i <= n; i++) { unsigned now = 1; for (unsigned j = n + 1; j >= n - i + 1; j--) if (j % (i + 1) == 0) now *= j / (i + 1); else now *= j; ans += now * s[k][i]; } return ans; } void init(unsigned n) { s[0][0] = 1; for (unsigned i = 1; i < MAXK; i++) for (unsigned j = 1; j <= i; j++) s[i][j] = s[i - 1][j - 1] + s[i - 1][j] * j; for (unsigned i = 2; i <= n; i++) { if (f[i] == 0) { prime[++tot] = f[i] = i; primek[tot] = power(i, k); pre[tot] = pre[tot - 1] + primek[tot]; } for (unsigned j = 1; j <= tot && prime[j] <= f[i]; j++) { unsigned tmp = prime[j] * i; if (tmp > n) break; f[tmp] = prime[j]; } } } int main() { read(n), read(k); limit = sqrt(n); init(limit); for (unsigned i = 1, nxt; i <= n; i = nxt) { unsigned tmp = n / i; val[++m] = tmp; cnt[m] = tmp - 1; sum[m] = calc(tmp) - 1; nxt = n / tmp + 1; if (tmp <= limit) home1[tmp] = m; else home2[n / tmp] = m; } for (unsigned i = 1; i <= tot; i++) { for (unsigned j = 1; prime[i] * prime[i] <= val[j]; j++) { unsigned tmp = val[j] / prime[i], pos; if (tmp <= limit) pos = home1[tmp]; else pos = home2[n / tmp]; cnt[j] -= cnt[pos] - (i - 1); pref[j] += sum[pos] - pre[i - 1]; sum[j] -= primek[i] * (sum[pos] - pre[i - 1]); } } for (unsigned i = 1; i <= m; i++) pref[i] += cnt[i]; unsigned finalans = 0; for (unsigned i = m; i >= 1; i--) { ans[i] = pref[i]; for (unsigned j = 2, nxt; j <= val[i]; j = nxt) { unsigned tmp = val[i] / j, pos; if (tmp <= limit) pos = home1[tmp]; else pos = home2[n / tmp]; nxt = val[i] / tmp + 1; ans[i] -= ans[pos] * (nxt - j); } finalans += (ans[i] - ans[i + 1]) * (n / val[i]) * (n / val[i]); } writeln(finalans); return 0; }