這題確實在能力範圍之外，現在學反演的瓶頸找到了，就是除了套路之外啥也不敢想，然後應該自己設計的線性篩總是設計不了

（這題想了也沒法做……這求多項式的手段真的讓我大開眼界了……）

Description

link

給定一個大整數 \(n\) 的質因數分解形式和一個整數 \(d\)

求不大於 \(n\) 且和 \(n\) 互質的數的 \(d\) 次方和

\(d\le 1000\)

Solution

上手推式子

\(f_d(n)=\sum_{p|n} p^d\mu(p) \sum_{i=1}^{\lfloor \frac n p\rfloor} i^d\)

然後我們發現這東西並不可做，在 \(n\) 這麼大的情況下開始毫無頭緒

以下參考了這篇題解link

設$S_d(n)=\sum\limits_^{n i}d$

那麼則有

\(f_d(n)=\sum_{p|n} p^d\mu(p) S_d(\lfloor \frac n p\rfloor)\)

猜 \(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 \times (\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})\)

然後是程式碼實現……（怎麼感覺都是知識點門檻很高但是沒啥水準）

\(\binom 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();}