考虑分成
B
=
n
B=\sqrt n
B=n块分别求
即设
f
(
d
,
x
)
=
∏
i
=
1
d
(
x
+
i
)
f(d,x)=\prod_{i=1}^d(x+i)
f(d,x)=∏i=1d(x+i)
则考虑求出
f
(
B
,
0
)
,
f
(
B
,
2
B
)
.
.
.
.
f
(
B
,
B
∗
B
)
f(B,0),f(B,2B)....f(B,B*B)
f(B,0),f(B,2B)....f(B,B∗B)
剩下的一点直接乘即可
对于这个考虑倍增
即考虑已知
f
(
d
,
0
)
,
f
(
d
,
B
)
.
.
.
f
(
d
,
d
B
)
f(d,0),f(d,B)...f(d,dB)
f(d,0),f(d,B)...f(d,dB)
求
f
(
2
d
,
0
)
,
f
(
2
d
,
B
)
.
.
.
.
f
(
2
d
,
2
d
B
)
f(2d,0),f(2d,B)....f(2d,2dB)
f(2d,0),f(2d,B)....f(2d,2dB)
由于有
f
(
2
d
,
x
)
=
f
(
d
,
x
)
f
(
d
,
d
+
x
)
f(2d,x)=f(d,x)f(d,d+x)
f(2d,x)=f(d,x)f(d,d+x)
那么只需要知道
f
(
d
,
0
)
.
.
.
f
(
d
,
2
d
B
)
f(d,0)...f(d,2dB)
f(d,0)...f(d,2dB)
和
f
(
d
,
d
)
,
f
(
d
,
2
d
)
.
.
.
f
(
d
,
2
d
B
+
d
)
f(d,d),f(d,2d)...f(d,2dB+d)
f(d,d),f(d,2d)...f(d,2dB+d)
设
h
(
x
)
=
f
(
d
,
B
x
)
h(x)=f(d,Bx)
h(x)=f(d,Bx)
那么就是已知
h
(
0
)
,
h
(
1
)
.
.
.
h
(
d
)
h(0),h(1)...h(d)
h(0),h(1)...h(d)
求
h
(
d
+
1
)
,
h
(
d
+
2
)
.
.
h
(
2
d
+
1
)
h(d+1),h(d+2)..h(2d+1)
h(d+1),h(d+2)..h(2d+1)和
h
(
d
B
)
,
h
(
d
B
+
1
)
.
.
.
h
(
d
B
+
2
d
)
h(\frac d B),h(\frac d B+1)...h(\frac d B+2d)
h(Bd),h(Bd+1)...h(Bd+2d)
考虑
已知
h
(
0
)
,
h
(
1
)
.
.
.
h
(
k
)
h(0),h(1)...h(k)
h(0),h(1)...h(k)
求
h
(
p
)
,
h
(
p
+
1
)
.
.
.
h
(
p
+
k
)
h(p),h(p+1)...h(p+k)
h(p),h(p+1)...h(p+k)
考虑利用拉格朗日插值
h
(
p
+
n
)
=
∑
i
=
0
k
h
(
i
)
∏
j
=
0
,
j
≠
i
k
p
+
n
−
j
i
−
j
=
∏
j
=
0
k
(
p
+
n
−
j
)
(
∑
i
=
0
k
h
(
i
)
p
+
n
−
i
∏
j
=
0
,
j
≠
i
k
1
i
−
j
)
=
∏
j
=
0
k
(
p
+
n
−
j
)
(
∑
i
=
0
k
1
p
+
n
−
i
∗
h
(
i
)
i
!
(
k
−
i
)
!
(
−
1
)
k
−
i
)
h(p+n)=\sum_{i=0}^kh(i)\prod_{j=0,j\not =i}^k\frac{p+n-j}{i-j}\\ =\prod_{j=0}^{k}(p+n-j)(\sum_{i=0}^k\frac{h(i)}{p+n-i}\prod_{j=0,j\not=i}^k\frac 1{i-j})\\ =\prod_{j=0}^{k}(p+n-j)(\sum_{i=0}^k\frac{1}{p+n-i}*\frac{h(i)}{i!(k-i)!(-1)^{k-i}})
h(p+n)=i=0∑kh(i)j=0,j=i∏ki−jp+n−j=j=0∏k(p+n−j)(i=0∑kp+n−ih(i)j=0,j=i∏ki−j1)=j=0∏k(p+n−j)(i=0∑kp+n−i1∗i!(k−i)!(−1)k−ih(i))
后面可以卷积求
即设
g
(
x
)
=
∑
i
1
p
−
k
+
i
x
i
,
f
(
x
)
=
∑
i
h
(
i
)
i
!
(
k
−
i
)
!
(
−
1
)
k
−
i
,
H
(
x
)
=
∑
i
h
(
i
+
p
−
k
)
∏
j
=
0
k
(
p
+
n
−
j
)
g(x)=\sum_i\frac{1}{p-k+i}x^i,f(x)=\sum_{i}\frac{h(i)}{i!(k-i)!(-1)^{k-i}},H(x)=\sum_i\frac{h(i+p-k)}{\prod_{j=0}^k(p+n-j)}
g(x)=∑ip−k+i1xi,f(x)=∑ii!(k−i)!(−1)k−ih(i),H(x)=∑i∏j=0k(p+n−j)h(i+p−k)
那么
H
=
f
∗
g
H=f*g
H=f∗g
这里面所有求逆都可以离线下来做到线性
这样能减少很多常数
(
10
s
→
3
s
)
(10s\rightarrow 3s)
(10s→3s)
这样可以做到 d → d ∗ 2 d\rightarrow d*2 d→d∗2
再考虑
d
→
d
+
1
d\rightarrow d+1
d→d+1
即
f
(
d
,
x
)
∗
(
x
+
d
+
1
)
→
f
(
d
+
1
,
x
)
f(d,x)*(x+d+1)\rightarrow f(d+1,x)
f(d,x)∗(x+d+1)→f(d+1,x)
另外多一个
f
(
d
,
d
B
+
B
)
f(d,dB+B)
f(d,dB+B)，直接算即可
复杂度为
O
(
n
l
o
g
n
)
O(\sqrt nlogn)
O(nlogn)
另外可以用威尔逊定理优化常数
由于模数要用 M T T MTT MTT
#include<bits/stdc++.h>
using namespace std;
#define cs const
#define re register
#define pb push_back
#define pii pair<int,int>
#define ll long long
#define fi first
#define se second
#define bg begin
cs int RLEN=1<<20|1;
inline char gc(){
static char ibuf[RLEN],*ib,*ob;
(ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
return (ib==ob)?EOF:*ib++;
}
inline int read(){
char ch=gc();
int res=0;bool f=1;
while(!isdigit(ch))f^=ch=='-',ch=gc();
while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
return f?res:-res;
}
inline ll readll(){
char ch=gc();
ll res=0;bool f=1;
while(!isdigit(ch))f^=ch=='-',ch=gc();
while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
return f?res:-res;
}
inline int readstring(char *s){
int top=0;char ch=gc();
while(isspace(ch))ch=gc();
while(!isspace(ch)&&ch!=EOF)s[++top]=ch,ch=gc();
return top;
}
template<typename tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<typename tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
int mod;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline int dec(int a,int b){return a-b<0?a-b+mod:a-b;}
inline int mul(int a,int b){static ll r;r=(ll)a*b;return (r>=mod)?(r%mod):r;}
inline void Add(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){static ll r;r=(ll)a*b;a=(r>=mod)?(r%mod):r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
#define poly vector<int>
namespace Poly{
cs int C=17,M=(1<<C)|5,T=(1<<16)-1;
int PP,P;
struct plx{
double x,y;
plx(double _x=0,double _y=0):x(_x),y(_y){}
friend inline plx operator +(cs plx &a,cs plx &b){
return plx(a.x+b.x,a.y+b.y);
}
friend inline plx operator -(cs plx &a,cs plx &b){
return plx(a.x-b.x,a.y-b.y);
}
friend inline plx operator *(cs plx &a,cs plx &b){
return plx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
friend inline plx operator /(cs plx &a,cs int &b){
return plx(a.x/b,a.y/b);
}
inline plx conj()cs{return plx(x,-y);}
};
plx *w[C+1];
cs double pi=acos(-1);
int rev[M];
inline void init_rev(int lim){
for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void init_w(){
for(int i=1;i<=C;i++)w[i]=new plx[(1<<(i-1))|1];
for(int i=0,l=1<<(C-1);i<l;i++)
w[C][i]=plx(cos(i*pi/l),sin(i*pi/l));
for(int i=C-1;i;i--)
for(int j=0,l=1<<(i-1);j<l;j++)w[i][j]=w[i+1][j<<1];
}
inline void fft(plx *f,int lim,int kd){
for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
plx a0,a1;
for(int mid=1,l=1;mid<lim;mid<<=1,l++)
for(int i=0;i<lim;i+=mid<<1)
for(int j=0;j<mid;j++)
a0=f[i+j],a1=w[l][j]*f[i+j+mid],f[i+j]=a0+a1,f[i+j+mid]=a0-a1;
if(kd==-1){
reverse(f+1,f+lim);
for(int i=0;i<lim;i++)f[i]=f[i]/lim;
}
}
inline void Mul(int *ret,int *A,int *B,int lim){
static plx a[M],b[M],c[M],d[M],da,db,dc,dd;
init_rev(lim);
for(int i=0;i<lim;i++)a[i]=plx(A[i]&T,A[i]>>16);
for(int i=0;i<lim;i++)b[i]=plx(B[i]&T,B[i]>>16);
fft(a,lim,1),fft(b,lim,1);
for(int i=0;i<lim;i++){
int j=(lim-i)&(lim-1);
da=(a[i]+a[j].conj())*plx(0.5,0);
db=(a[j].conj()-a[i])*plx(0,0.5);
dc=(b[i]+b[j].conj())*plx(0.5,0);
dd=(b[j].conj()-b[i])*plx(0,0.5);
c[i]=(da*dc)+((da*dd)*plx(0,1));
d[i]=(db*dd)+((db*dc)*plx(0,1));
}
fft(c,lim,-1),fft(d,lim,-1);
for(int i=0;i<lim;i++){
ll da=(ll)(d[i].x+0.5)%mod,db=(ll)(d[i].y+0.5)%mod,dc=(ll)(c[i].y+0.5)%mod,dd=(ll)(c[i].x+0.5)%mod;
ret[i]=((da*PP%mod)+((db+dc)*P%mod)+dd)%mod;
}
}
}
using namespace Poly;
cs int N=(1<<18)|8;
int fac[N],ifac[N],p[N],ip[N],t[N];
inline void init_inv(int n){
fac[0]=ifac[0]=1;
for(int i=1;i<=n;i++)fac[i]=mul(fac[i-1],i);
ifac[n]=Inv(fac[n]);
for(int i=n-1;i;i--)ifac[i]=mul(ifac[i+1],i+1);
}
inline void calc(int del,int n,int *f,int *g){
static int a[N],b[N];
int len=1;while(len<=n+n)len<<=1;
for(int i=0;i<=n;i++){
a[i]=mul(f[i],mul(ifac[i],ifac[n-i]));
if((n-i)&1)a[i]=dec(0,a[i]);
}
for(int i=0;i<=n+n;i++)t[i]=dec(add(del,i),n);
p[0]=t[0];
for(int i=1;i<=n+n;i++)p[i]=mul(t[i],p[i-1]);
ip[n+n]=Inv(p[n+n]);
for(int i=n+n;i;i--)ip[i-1]=mul(ip[i],t[i]);
b[0]=ip[0];
for(int i=1;i<=n+n;i++)b[i]=mul(ip[i],p[i-1]);
for(int i=n+1;i<len;i++)a[i]=0;
for(int i=n+n+1;i<len;i++)b[i]=0;
Mul(a,a,b,len);
int coef=1,l=dec(del,n),r=del;
for(int i=l;i<=r;i++)Mul(coef,i);
for(int i=0;i<=n;i++)t[i]=i+l;
p[0]=t[0];
for(int i=1;i<=n;i++)p[i]=mul(t[i],p[i-1]);
ip[n]=Inv(p[n]);
for(int i=n;i;i--)ip[i-1]=mul(ip[i],t[i]);
t[0]=ip[0];
for(int i=1;i<=n;i++)t[i]=mul(ip[i],p[i-1]);
for(int i=0;i<=n;i++,l++,r++){
g[i]=mul(a[i+n],coef);
Mul(coef,mul(t[i],r+1));
}
}
int f[N],g[N];
inline void getfac(int b){
int ti=0;f[0]=1;
for(int i=b;i;i>>=1)ti++;
for(int bs=0;ti>=0;ti--){
if(bs){
calc(bs+1,bs,f,g);
for(int i=0;i<=bs;i++)f[bs+i+1]=g[i],g[i]=0;f[bs+bs+1]=0;
calc(mul(bs,Inv(b)),bs<<1,f,g),bs<<=1;
for(int i=0;i<=bs;i++)Mul(f[i],g[i]),g[i]=0;
}
if((b>>ti)&1){
for(int i=0;i<=bs;i++)
Mul(f[i],(1ll*i*b+bs+1)%mod);
bs++,f[bs]=1;
for(int i=1,v=mul(bs,b);i<=bs;i++)Mul(f[bs],add(v,i));
}
}
}
inline int solve(int d){
int b=sqrt(d);
init_inv(b),getfac(b);
int res=1;
for(int i=0,p=0;;i+=b,p++){
if(i+b>=d){
for(int j=i+1;j<=d;j++)Mul(res,j);
return res;
}
Mul(res,f[p]);
}
}
inline int calc(int n){
mod=read();
P=ksm(2,16),PP=mul(P,P);
if(n>mod-1-n){
int res=Inv(solve(mod-1-n));
return ((mod-n-1)&1)?res:mod-res;
}
return solve(n);
}
signed main(){
#ifdef Stargazer
freopen("lx.in","r",stdin);
#endif
init_w();
int T=read();
while(T--)
cout<<calc(read())<<'\n';
}
扫一扫
596
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
