[51nod 1847] 奇怪的数学题

爽到了爽到了，真的给爷爽到了！！！！！

时限：\(\tt 1500ms\) ，我的代码：\(\tt 1484ms\)

一、题目

\(\sum_{i=1}^n\sum_{j=1}^nsgcd(i,j)^k\)

其中 \(sgcd(i,j)\) 表示 \((i,j)\) 的所有公约数中第二大的数，输出答案模 \(2^{32}\) 的结果。

\(1\leq n\leq 10^9,1\leq k\leq 50\)

二、解法

首先可以推式子，设 \(f(i)\) 表示 \(i\) 的次大公约数：

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

然后就是和 这道题 一模一样的反演方法：

\[\sum_{T=1}^n(\frac{n}{T})^2\sum_{d|T}f(d)^k\mu(\frac{T}{d}) \]

外面还是用整除分块，里面还是用杜教筛，令辅助函数 \(g\) 为 \(I\) ，现在唯一的问题是解决 \(\sum_{i=1}^n f(i)^k\) ，其他的地方都跟那道题差不了多少。

\(f(i)=\frac{i}{minp}\) ，其中 \(minp\) 表示 \(i\) 的最小质因子，这个东西似乎不是积性函数，直接套 \(\tt Min25\) 好像不行。仔细回忆一下我们算法的全过程，还有一个地方是用到了最小质因子的，就是筛 \(g\) 的那个部分，记得这个柿子吗？

\[g(n,j)=g(n,j-1)-f(p_j)(g(n/p_j,j-1)-g(p_{j-1},j-1)) \]

后面是算所谓合数，\(p_j\) 就是这些合数的最小质因子了，我们可以利用这个部分算 \(f\) ，定义 \(G(n)\) 为 \(n\) 以内合数的 \(f\) 函数值，最小质因子是都被这个部分枚举到了的，我们不把最小质因子的函数值算进去就可以了：

\[G(n)=\sum g(n/p_j,j-1)-g(p_{j-1},j-1) \]

这个 \(g\) 表示的函数值的 \(x^k\) ，最后 \(f\) 的求和可以表示成 \(pcnt[n]+G[n]\) ，\(pcnt[n]\) 表示 \(n\) 以内质数的个数，所以魔改一下 \(\tt Min25\) 的筛法就是这个样子的：

for(int i=1;i<=cnt;i++)
	for(int j=1;j<=tot && p[i]*p[i]<=w[j];j++)
	{
		int k=id(w[j]/p[i]);
		g1[j]=(g1[j]-g1[k]+i-1)%MOD;//这个筛的是质数个数
		g2[j]=(g2[j]-pk[i]*(g2[k]-g2[id(p[i-1])]))%MOD;//表示函数值为x^k,pk[i]是质数pi的k次方
		G[j]=(G[j]+g2[k]-g2[id(p[i-1])])%MOD;//也就是对于所有最小质因子的答案求和
	}

处理出一开始的 \(g\) 需要快速算：\(\sum_{i=1}^ni^k\) ，这个经典问题用第二类斯特林数就行了：

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

\[=\sum_{j=1}^kS(k,j)j!\sum_{i=1}^nC(i,j) \]

\[=\sum_{j=1}^kS(k,j)j!C(n+1,j+1) \]

\[=\sum_{j=1}^kS(k,j)\frac{(n+1)!}{(n-j)!j+1} \]

因为模数没有逆元所以只能暴力把他们乘起来，遇到能除 \(j+1\) 的因子时再除，可以通过模 \(j+1\) 余数为 \(0\) 来判断。这道题理应写自然溢出，但是由于我一写就炸所以还是用了 \(\tt\#define\space int\space long\space long\)

#include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;
const int N = 200000;
const int M = 200005;
#define int long long
const int MOD = 4294967296ll;
int read()
{
	int x=0,f=1;char c;
	while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
	while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
	return x*f;
}
int n,k,cnt,ans,p[M],pk[M],vis[M],s[60][60],pd[M];
int sqr,tot,w[M],g1[M],g2[M],G[M],id1[M],id2[M],sum[M];
int qkpow(int a,int b)
{
	int r=1;
	while(b>0)
	{
		if(b&1) r=r*a%MOD;
		a=a*a%MOD;
		b>>=1;
	}
	return r;
}
void init(int n)
{
	for(int i=2;i<=n;i++)
	{
		if(!vis[i])
		{
			p[++cnt]=i;
			pk[cnt]=qkpow(i,k); 
		}
		for(int j=1;j<=cnt && i*p[j]<=n;j++)
		{
			vis[i*p[j]]=1;
			if(i%p[j]==0) break;
		}
	}
	s[0][0]=1;
	for(int i=1;i<=k;i++)
		for(int j=1;j<=i;j++)
			s[i][j]=(s[i-1][j-1]+s[i-1][j]*j)%MOD;
}
int cal(int n)
{
	int ret=0,tmp=0,c=0;
	for(int i=1;i<=k;i++)
	{
		c=1;tmp=(n-i+1)%(i+1);
		for(int j=n-i+1;j<=n+1;j++,tmp++)
		{
			if(tmp>=i+1) tmp-=i+1;
			if(!tmp) c=c*j/(i+1)%MOD;
			else c=c*j%MOD;
		}
		ret=(ret+c*s[k][i])%MOD;
	}
	return ret;
}
int id(int x)
{
	if(x<=sqr) return id1[x];
	return id2[n/x];
}
int get(int n)
{
	if(n<=1) return 0;
	if(pd[id(n)]) return sum[id(n)];
	int ans=g1[id(n)]+G[id(n)];
	for(int l=2,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans=(ans-(r-l+1)*get(n/l))%MOD;
	}
	sum[id(n)]=ans;pd[id(n)]=1;
	return ans;
}
signed main()
{
	n=read();k=read();
	init(N);sqr=sqrt(n);
	for(int l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		w[++tot]=n/l;
		g1[tot]=w[tot]-1;
		g2[tot]=cal(w[tot])-1;
		if(n/l<=sqr) id1[n/l]=tot;
		else id2[n/(n/l)]=tot;
	}
	for(int i=1;i<=cnt;i++)
		for(int j=1;j<=tot && p[i]*p[i]<=w[j];j++)
		{
			int k=id(w[j]/p[i]);
			g1[j]=(g1[j]-g1[k]+i-1)%MOD;
			g2[j]=(g2[j]-pk[i]*(g2[k]-g2[id(p[i-1])]))%MOD;
			G[j]=(G[j]+g2[k]-g2[id(p[i-1])])%MOD;
		}
	for(int l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans=(ans+(n/l)*(n/l)%MOD*(get(r)-get(l-1)))%MOD;
	}
	printf("%lld\n",(ans+MOD)%MOD);
}