[51nod 1847] 奇怪的数学题

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

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

一、题目

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

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

(1leq nleq 10^9,1leq kleq 50)

二、解法

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

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

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

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

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

(f(i)=frac{i}{minp}) ，其中 (minp) 表示 (i) 的最小质因子，这个东西似乎不是积性函数，直接套 ( t 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) 以内质数的个数，所以魔改一下 ( t 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}^nsum_{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) 来判断。这道题理应写自然溢出，但是由于我一写就炸所以还是用了 ( t#definespace intspace longspace 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
",(ans+MOD)%MOD);
}