「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}})$

---

代码

#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;
}