Luogu6271 [湖北省队互测2014]一个人的数论

这题确实在能力范围之外，现在学反演的瓶颈找到了，就是除了套路之外啥也不敢想，然后应该自己设计的线性筛总是设计不了

（这题想了也没法做……这求多项式的手段真的让我大开眼界了……）

Description

link

给定一个大整数 (n) 的质因数分解形式和一个整数 (d)

求不大于 (n) 且和 (n) 互质的数的 (d) 次方和

(dle 1000)

Solution

上手推式子

[f_d(n)=sum_{p|n} p^dmu(p) sum_{i=1}^{lfloor frac n p floor} i^d ]

然后我们发现这东西并不可做，在 (n) 这么大的情况下开始毫无头绪

以下参考了这篇题解link

设(S_d(n)=sumlimits_{i=1}^n i^d)

那么则有

[f_d(n)=sum_{p|n} p^dmu(p) S_d(lfloor frac n p floor) ]

猜 (S_d(n)) 是一个多项式，那么有如下表示（好像有证明？link)

这一部分等过一段时间填坑吧……（先背个结论）

[S_d(n) =sum_{i=0}^{d+1} a_i n^i ]

我们这里可以用任意插值 (+) 高斯消元得到这个多项式

接着推式子，得到

[f_d(n)=sum_{i=0}^{d+1} a_i imes (sum_{p|n} p^d (frac n p)^i mu (p) ) ]

然后看看后面的部分就是两个函数 (u(x)=x^d mu (x)) 和 (v(x)=x^i) 的狄利克雷卷积

令 (h=u * v)

(u,v) 都是积性函数，所以它们的狄利克雷卷积也是积性函数（可以简单证明一下，最后发现推出来要证 (mu) 是积性函数，所以就证毕了），所以要线性筛

因为它是个积性函数，所以可以把 (n) 分解质因数然后乘起来

(h_i(p^a)) 把 (mu) 的值为 (0) 的部分消掉之后剩下了：

[p^{ai}(1-p^{d-i}) ]

然后是代码实现……（怎么感觉都是知识点门槛很高但是没啥水准）

[inom n m ]

我数据范围开小了卡了 (3h)

Code

#include<bits/stdc++.h>
using namespace std;
#define int long long
namespace yspm{
	inline int read()
	{
		int res=0,f=1; char k;
		while(!isdigit(k=getchar())) if(k=='-') f=-1;
		while(isdigit(k)) res=res*10+k-'0',k=getchar();
		return res*f;
	}
	const int mod=1e9+7,N=1010;
	inline int ksm(int x,int y)
	{
		if(y<0) return ksm(ksm(x,-y),mod-2);
		int res=1; for(;y;y>>=1,(x*=x)%=mod) if(y&1) res=res*x%mod;
		return res;
	}
	int ans,n,d,p[N][2],m[N][N],res[N],sum[N];
	inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
	signed main()
	{
		d=read(); n=read();
		for(int i=1;i<=n;++i) p[i][0]=read(),p[i][1]=read();
		for(int i=0;i<=d+2;++i) 
		{
			int res=0; for(int j=1;j<=i+1;++j) res=add(res,ksm(j,d));
			m[i][d+2]=res; m[i][0]=1;
			for(int j=1;j<=d+1;++j) m[i][j]=m[i][j-1]*(i+1)%mod;
		}
		for(int i=0;i<=d+2;++i) 
		{
			int t=i;
			for(int j=i;j<=d+2;++j) if(m[j][i]){t=j; break;}
			for(int j=i;j<=d+2;++j) swap(m[i][j],m[t][j]);
			for(int j=i+1;j<=d+2;++j)
			{
				int x=m[j][i]*ksm(m[i][i],mod-2)%mod;
				for(int k=i;k<=2+d;++k) m[j][k]=(m[j][k]-m[i][k]*x%mod+mod)%mod; 
			} 
		}
		for(int i=d+2;i>=0;--i) 
		{
			m[i][d+2]=m[i][d+2]*ksm(m[i][i],mod-2)%mod;
			for(int j=i-1;j>=0;--j) m[j][d+2]-=m[i][d+2]*m[j][i]%mod,m[j][d+2]=(m[j][d+2]+mod)%mod;
		}
		for(int i=0;i<=d+1;++i) res[i]=m[i][d+2]; 
		for(int i=0;i<=d+1;++i)
		{
			int tmp=1;
			for(int j=1;j<=n;++j) 
			{
				tmp=tmp*ksm(p[j][0],p[j][1]*i)%mod*((1-ksm(p[j][0],d-i)+mod)%mod)%mod; 
			}
			ans=add(ans,res[i]*tmp%mod);
		} cout<<ans<<endl;
		return 0;
	}
}
signed main(){return yspm::main();}