1. 程式人生 > >青春野狼不做理性小魔女的夢 - 莫比烏斯反演 - 拉格朗日插值 - 杜教篩

青春野狼不做理性小魔女的夢 - 莫比烏斯反演 - 拉格朗日插值 - 杜教篩

題目大意:
給定 { A k } \{A_k\} (有些位置是-1表示一會要由你決定)。
對每個 m

[ 1 , n ] m\in[1,n] ,求有多少方案(將 1
-1
變為 [ 0 , m ) [0,m) 之間的整數)使得
i = 1 k A i x i = 1 (   m o d     m ) \sum_{i=1}^kA_ix_i=1(\bmod\ m)
存在解。
題解:
顯然有解當且僅當所有的A和m的gcd是1。
容斥反演一波分類討論一下發現要寫個拉格朗日插值和杜教篩。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define lint long long
#define mod 1000000007
#define ull unsigned lint
#define db long double
#define pb push_back
#define mp make_pair
#define fir first
#define sec second
#define gc getchar()
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
typedef pair<int,int> pii;
typedef set<int>::iterator sit;
inline int inn() { int x;scanf("%d",&x);return x; }
const int MXK=60,N=6000000;
int gcd(int a,int b) { return a?gcd(b%a,a):b; }
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans; }
int a[MXK];
namespace FS_space{
	int k,t,y[MXK],ys[MXK],fac[MXK],facinv[MXK],pre[MXK],suf[MXK];
	inline int prelude(int _k)
	{
		k=_k,t=k+2;
		rep(i,1,t) y[i]=fast_pow(i,k)+y[i-1],(y[i]>=mod?y[i]-=mod:0);
		rep(i,fac[0]=1,t) fac[i]=(lint)fac[i-1]*i%mod;
		facinv[t]=fast_pow(fac[t],mod-2);
		for(int i=t-1;i>=0;i--) facinv[i]=(i+1ll)*facinv[i+1]%mod;
		rep(i,1,t)
		{
			ys[i]=y[i],y[i]=(lint)y[i]*facinv[i-1]%mod*facinv[t-i]%mod;
			if((t-i)&1) (y[i]?y[i]=mod-y[i]:0);
		}
		return 0;
	}
	inline int Fs(int n)
	{
//		int Ans=0;rep(i,1,n) (Ans+=fast_pow(i,k))%=mod;return Ans;
		if(n<=t) return ys[n];pre[0]=suf[t+1]=1;int ans=0;
		rep(i,1,t) pre[i]=(lint)pre[i-1]*(n-i)%mod;
		for(int i=t;i;i--) suf[i]=(lint)suf[i+1]*(n-i)%mod;
		rep(i,1,t) ans=(ans+(lint)y[i]*pre[i-1]%mod*suf[i+1])%mod;
		return ans;
	}
}using FS_space::Fs;
namespace MS_space{
	int mxn,mu[N],ms[N],p[N],np[N];
	unordered_map<int,int> sav;
	inline int prelude(int n)
	{
		mu[1]=ms[1]=1,mxn=n;
		for(int i=2,c=0;i<=n;i++)
		{
			if(!np[i]) p[++c]=i,mu[i]=-1;
			ms[i]=ms[i-1]+mu[i];
			for(int j=1,x;j<=c&&p[j]<=n/i;j++)
			{
				np[x=p[j]*i]=1,mu[x]=-mu[i];
				if(i%p[j]==0) { mu[x]=0;break; }
			}
		}
		return 0;
	}
	int Ms(int n)
	{
		if(n<=mxn) return ms[n];int ans=0;
		if(sav.count(n)) return sav[n];
		for(int s=2,t;s<=n;s=t+1)
			t=n/(n/s),ans+=(t-s+1)*Ms(n/s);
		return sav[n]=1-ans;
	}
	inline int Ms(int s,int t)
	{
		int ans=Ms(t)-Ms(s-1);
		ans%=mod,(ans<0?ans+=mod:0);
		return ans;
	}
}using MS_space::Ms;
int p[100010];
int main()
{
	int k=inn(),n=inn(),g=0,cnt=0;
	rep(i,1,k) a[i]=inn();
	rep(i,1,k) if(a[i]>=0) g=gcd(g,a[i]);else cnt++;
	FS_space::prelude(cnt),MS_space::prelude(N-1);
	if(!g)//\sum_{d=1}^n \mu(d)F_cnt(n/d)
	{
		int ans=0;
		for(int s=1,t;s<=n;s=t+1)
			t=n/(n/s),ans=(ans+(lint)Ms(s,t)*Fs(n/s))%mod;
		return !printf("%d\n",ans);
	}
	else{
		int c=0,ans=0;
		rep(i,2,g/i) if(g%i==0)
			for(p[c++]=i,g/=i;g%i==0;g/=i);
		if(g>1) p[c++]=g;
		int all=(1<<c)-1;
		rep(i,0,all)
		{
			int d=1,s=1;
			rep(j,0,c-1) if((i>>j)&1) d*=p[j],s*=-1;
//			debug(d)sp,debug(s)sp,debug(n/d)sp,debug(cnt)sp,debug(Fs(n/d))ln;
			if(d>n) continue;
			if(s>0) ans+=Fs(n/d),(ans>=mod?ans-=mod:0);
			else ans-=Fs(n/d),(ans<0?ans+=mod:0);
		}
		return !printf("%d\n",ans);
	}
	return 0;
}