「51nod 1847」奇怪的数学题

51nod 1847

题意

给出 \(n,k\)

\[\sum_{i=1}^n\sum_{j=1}^n sgcd(i,j)^k\]

其中 \(sgcd(i,j)\) 表示 \(i,j\) 的次大公约数,特殊地,\(sgcd(1,1)=0\)

\(2^{32}\) 取模

\(n\le 10^9, k\le 50\)


做法

\(f(x)\) 表示 \(x\) 的次大因数,特殊地,\(f(1)=0\)

于是

\[ \begin{align} & \sum_{i=1}^n\sum_{j=1}^n sgcd(i,j)^k \\ = & \sum_{i=1}^n\sum_{j=1}^n f(\gcd(i,j))^k \\ = & \sum_{g=1}^n f(g)^k \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{g} \rfloor} [\gcd(i,j)=1] \\ = & \sum_{g=1}^n f(g)^k \left(2 \left( \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} \varphi(i) \right) - 1 \right) \end{align} \]

显然右边的部分可以Min_25筛/杜教筛处理

Min_25筛参考这里

我们需要对于每个 \(\lfloor \frac{n}{x} \rfloor\) 求出 \(\sum_{g=1}^{\lfloor \frac{n}{x} \rfloor} f(g)^k\)

显然对于 \(g\ne 1\)\(f(g)=\frac{g}{g\text{ 的最小质因子}}\),并且对于合数 \(g\)\(g\text{ 的最小质因子}\le \sqrt{g}\)

每个质数的贡献是 \(1\),这可以在Min_25筛的第一部分预处理质数时计算

事实上这里不需要第二部分,考虑对于一个 \(\lfloor \frac{n}{x} \rfloor\)\(\lfloor \frac{n}{x} \rfloor\) 以内的所有合数都恰好被筛到一次,直接在第一部分内统计答案就好了

还有一点是这里模数不是质数,在Min_25筛预处理时需要求 \(\sum_{i=1}^n i^k\),不能用插值之类的做

因为

\[ \begin{align} x^k & = \sum_{i=0}^k \begin{Bmatrix} k \\ i \end{Bmatrix} x^{\underline{i}} \\ & =\sum_{i=0}^k \begin{Bmatrix} k \\ i \end{Bmatrix} \binom{x}{i}i! \end{align} \]

其中 \(\begin{Bmatrix} n \\ k \end{Bmatrix}\) 表示第二类斯特林数,上式可以用组合意义理解(膜yx

所以

\[ \begin{align} \sum_{x=1}^n x^k & = \sum_{x=1}^n \sum_{i=0}^k \begin{Bmatrix} k \\ i \end{Bmatrix} \binom{x}{i}i! \\ & = \sum_{i=0}^k i! \begin{Bmatrix} k \\ i \end{Bmatrix} \sum_{x=1}^n \binom{x}{i} \\ & = \sum_{i=0}^k i! \begin{Bmatrix} k \\ i \end{Bmatrix} \binom{n+1}{i+1} \\ & = \sum_{i=0}^k \begin{Bmatrix} k \\ i \end{Bmatrix} \frac{(n+1)^{\underline{i+1}}}{i+1} \end{align} \]

显然最后的下降幂中必有一项被 \(i+1\) 整除,所以可以 \(\mathcal O(k^2)\) 地计算一个 \(k\) 次幂前缀和

总复杂度 \(\mathcal O(\frac{n^{\frac{3}{4}}}{\log n} + \sqrt{n}k^2)\)\(\mathcal O(\frac{n^{\frac{3}{4}}}{\log n} + \sqrt{n}k^2 + n^{\frac{2}{3}})\)


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include<cstdio>
#include<algorithm>
#include<ctype.h>
#include<string.h>
#include<math.h>

using namespace std;
#define ll long long

const int N = 31627, K = 52;
int n, k, id, cnt, sn, prime[N], pk[N], a[N<<1];
unsigned Ans, s[K][K], g[N<<1], f[N<<1], h[N<<1], F[N<<1], ans[N<<1];
inline int Id(int x){ return x<=sn?x:id-n/x+1;}
int main() {
scanf("%d%d", &n, &k), sn=sqrt(n);
s[1][1]=1;
for(int i=2; i<=k; ++i) for(int j=1; j<=i; ++j) s[i][j]=s[i-1][j-1]+j*s[i-1][j];
for(int i=1; i<=n; i=a[id]+1){
a[++id]=n/(n/i);
for(int j=0; j<=k; ++j){
unsigned x=s[k][j];
for(int t=0; t<=j; ++t)
if((a[id]+1-t)%(j+1)) x*=a[id]+1-t;
else x*=(a[id]+1-t)/(j+1);
g[id]+=x;
}
--g[id], f[id]=(ll)a[id]*(a[id]+1)/2-1, h[id]=a[id]-1;
}
for(int i=2; i<=sn; ++i) if(h[i]!=h[i-1]){
prime[++cnt]=i;
for(int j=pk[cnt]=1; j<=k; ++j) pk[cnt]*=i;
for(int j=id; i*i<=a[j]; --j){
int p=Id(a[j]/i);
g[j]-=pk[cnt]*(g[p]-g[i-1]);
ans[j]+=g[p]-g[i-1];
f[j]-=i*(f[p]-f[i-1]);
h[j]-=h[p]-(cnt-1);
}
}
for(int i=1; i<=id; ++i) F[i]=(f[i]-=h[i]);
for(int i=cnt; i; --i)
for(int j=id; prime[i]*prime[i]<=a[j]; --j)
for(unsigned x=prime[i], y=x-1; (ll)x*prime[i]<=a[j]; x*=prime[i], y*=prime[i])
F[j]+=y*(F[Id(a[j]/x)]-f[prime[i]]+prime[i]);
for(int i=2; i<=id; ++i) Ans+=(ans[i]-ans[i-1]+h[i]-h[i-1])*(2*F[id-i+1]+1);
return printf("%u\n", Ans), 0;
}