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