集训队作业2018: 取名字太难了(FFT)

题意:
大概是求 ∏ i = 1 n ( x + i ) \prod_{i=1}^n(x+i) i=1n(x+i)系数模 p p p意义下的分布。

题解:
分为 ( ∏ i = 1 p − 1 ( x + i ) ) ⌊ n p ⌋ (\prod_{i=1}^{p-1}(x+i))^{\lfloor\frac{n}{p}\rfloor} (i=1p1(x+i))pn,以及 ( ∏ i = 1 n   m o d   p ) ( x + i ) (\prod_{i=1}^{n \bmod p})(x+i) (i=1nmodp)(x+i)

前者等于 ( x p − 1 − 1 ) ⌊ n p ⌋ (x^{p-1}-1)^{\lfloor\frac{n}{p}\rfloor} (xp11)pn,后者最高次数小于 p − 1 p-1 p1(等于把前面的次数加一即可)。

发现两部分系数不影响,可以分开算,后面直接倍增FFT了,前面二项式展开然后Lucas定理合并即可,合并可以用原根+FFT,时间复杂度 O ( p log ⁡ p + p log ⁡ n ) O(p \log p + p \log_n) O(plogp+plogn)

#include <bits/stdc++.h>
using namespace std;
typedef double DB;
typedef long long LL;

const int N=3e5+50, mod=998244353;
LL n,k; int p,G,ind[N],ori[N];
inline int add(int x,int y,int md) {return (x+y>=md) ? (x+y-md) : (x+y);}
inline int dec(int x,int y,int md) {return (x-y<0) ? (x-y+md) : (x-y);} 
inline int mul(int x,int y,int md) {return (long long)x*y%md;}
inline int power(int a,int b,int md,int rs=1) {for(;b;b>>=1,a=mul(a,a,md)) if(b&1) rs=mul(rs,a,md); return rs;}

namespace FFT {
	const int N=3e6+50, M=sqrt(mod);
	const DB PI2=acos(-1.0)*2;
	int k,pos[N];
	struct CP {
		DB r,i;
		CP(DB r=0,DB i=0) : r(r),i(i) {}
		friend inline CP operator +(const CP &a,const CP &b) {return CP(a.r+b.r,a.i+b.i);}
		friend inline CP operator -(const CP &a,const CP &b) {return CP(a.r-b.r,a.i-b.i);}
		friend inline CP operator *(const CP &a,const CP &b) {return CP(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
	} A[N],B[N],C[N],w[N];
	inline void clear(CP *a) {for(int i=0;i<k;i++) a[i]=CP(0,0);}
	inline void init(int nn) {
		for(k=1;k<=nn;k<<=1);
		for(int i=0;i<k;i++) pos[i]=(i&1) ? (pos[i>>1]>>1)^(k>>1) : (pos[i>>1]>>1);
	}
	inline void mul(CP *a,CP *b,CP *c) {
		for(int i=0;i<k;i++) c[i]=a[i]*b[i];
	}
	inline void dft(CP *a,int o=1) {
		for(int i=1;i<k;i++) 
			if(pos[i]>i) swap(a[pos[i]],a[i]);
		for(int bl=1;bl<k;bl<<=1) {
			int tl=bl<<1;
			for(int i=0;i<bl;i++) w[i]=CP(cos(PI2/tl*i),sin(PI2/tl*i*o));
			for(int bg=0;bg<k;bg+=tl)
				for(int j=0;j<bl;j++) {
					CP &t1=a[bg+j], &t2=a[bg+j+bl], t=t2*w[j];
					t2=t1-t; t1=t1+t;
				}
		}
		if(!~o) for(int i=0;i<k;i++) a[i].r/=k, a[i].i/=k;
	}
} 

using FFT::CP; 
using FFT::A;
using FFT::B;
using FFT::C;
using FFT::M;

struct poly {
	vector <int> a;
	poly(int d=0,int t=0) {a.resize(d+1); a[d]=t;}
	inline int deg() const {return a.size()-1;}
	inline int& operator [](const int &i) {return a[i];}
	inline const int& operator[](const int &i) const {return a[i];}
	friend inline poly mul(const poly &a,const poly &b,const int &md) {
		poly c(a.deg()+b.deg()); 
		FFT::init(c.deg());
		FFT::clear(A); FFT::clear(B);
		for(int i=0;i<=a.deg();i++) A[i]=CP(a[i]/M,a[i]%M);
		FFT::dft(A);
		for(int i=0;i<=b.deg();i++) B[i]=CP(b[i]/M,0);
		FFT::dft(B);
		FFT::mul(A,B,C);
		FFT::dft(C,-1);
		for(int i=0;i<=c.deg();i++) c[i]=((LL)(C[i].r+0.5)%md*M%md*M%md+(LL)(C[i].i+0.5)%md*M%md)%md;
		FFT::clear(B);
		for(int i=0;i<=b.deg();i++) B[i]=CP(b[i]%M,0);
		FFT::dft(B);
		FFT::mul(A,B,C);
		FFT::dft(C,-1);
		for(int i=0;i<=c.deg();i++) 
			c[i]=add(c[i],((LL)(C[i].r+0.5)%md*M%md+(LL)(C[i].i+0.5))%md,md);
		return c;
	}
	inline void pt() {
		for(int i=0;i<=deg();i++)
			cout<<a[i]<<' ';
		cout<<'\n';
	}
};

struct combin {
	int fac[N],ifac[N];
	inline void init() {
		fac[0]=1;
		for(int i=1;i<p;i++) fac[i]=mul(fac[i-1],i,p);
		ifac[p-1]=power(fac[p-1],p-2,p);
		for(int i=p-2;~i;i--) ifac[i]=mul(ifac[i+1],i+1,p);
	}
	inline int C(int a,int b) {return mul(fac[a],mul(ifac[b],ifac[a-b],p),p);}
} cb;

namespace ORD {
	vector <int> ff;
	inline bool check(int v) {
		for(auto u:ff)
			if(power(v,(p-1)/u,p)==1) return false;
		return true;
	}
	inline void find_G() {
		int nn=p-1;
		for(int i=2;i*i<=nn;i++)
			if(!(nn%i)) {
				while(!(nn%i)) nn/=i;
				ff.push_back(i);
			}
		if(nn!=1) ff.push_back(nn);
		for(G=1;;G++)
			if(check(G)) break;
	}
}

inline poly calc_fac(int nn) {
	if(!nn) return poly(0,1);
	if(nn&1) {
		poly t=poly(1,1); t[0]=nn;
		return mul(calc_fac(nn-1),t,p);
	} else {
		int l=nn/2;
		poly f=calc_fac(l);
	//	f.pt();
		poly c(l,0),d(l,0);
		for(int i=0,pw=1;i<=l;i++) c[i]=mul(f[i],mul(pw,cb.fac[i],p),p), pw=mul(pw,l,p);
		for(int i=0;i<=l;i++) d[l-i]=cb.ifac[i];
		poly g=mul(c,d,p);
		poly h(l,0);
		const int inv=power(l,p-2,p);
		for(int i=0,iv=1;i<=l;i++) h[i]=mul(iv,mul(g[i+l],cb.ifac[i],p),p), iv=mul(iv,inv,p);
		return mul(f,h,p);
	}
}

inline poly calc_pow(LL k) {
	poly f(0,1);
	long long base=1;
	while(k) {
		int z=k%p;
		k/=p;
		poly h(p,0);
		for(int i=0;i<=z;i++) {
			int tag=(base*i%2) ? (p-cb.C(z,i)) : cb.C(z,i);
			h[ind[tag]]++;
		}
		f=mul(f,h,mod);
		if(f.deg()>=p-1) {
			for(int i=p-1;i<=f.deg();i++)
				f[i%(p-1)]=add(f[i%(p-1)],f[i],mod);
			f.a.resize(p-1);
		} 
		while(f.a.size()>1 && !f.a.back()) f.a.pop_back();
		base=base*p;
	}
	return f;
}

int main() {
	cin>>n>>p;
	cb.init();
	ORD::find_G();
	
	for(int i=0,pw=1;i<p-1;i++) ind[pw]=i, ori[i]=pw, pw=mul(pw,G,p);
	k=((n%p)==p-1) ? n/p+1 : n/p;
	poly f1=calc_fac(max(n-k*p,0ll));
	poly f2=calc_pow(k);
	
	poly f3(p-2,0);
	for(int i=0;i<=f1.deg();i++) if(f1[i]) ++f3[ind[f1[i]]];
	
	poly f=mul(f3,f2,mod);
	if(f.deg()>=p-1) {
		for(int i=p-1;i<=f.deg();i++)
			f[i%(p-1)]=add(f[i%(p-1)],f[i],mod);
	} f.a.resize(p-1);
	int total=(n+1)%mod;
	for(int i=0;i<p-1;i++) total=dec(total,f[i],mod);
	cout<<total<<'\n';
	for(int i=1;i<p;i++) cout<<f[ind[i]]<<'\n';
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值