题面

前置芝士

（四次\(FFT\)）

题解

有多点求值的做法然而被\(shadowice\)巨巨吊起来打了一顿，所以来学一下倍增

成功同时拿到本题最优解和最劣解……

\(Min_{25}\)牛逼！（据说这是然而我看不懂就是了）

真的快的不要不要的……

和多点求值一样，我们还是设\(s=\sqrt{n}\)，并设多项式

\[g_s(x)=\sum_{i=1}^s(x+i)\]

求出\(g_s(0),g_s(s),g_s(2s),..,g_s((s-1)s)\)的值，最后再乘上\(\prod_{i=s^2+1}^ni\)即可

具体来说我们现在要完成两个操作

已知

\[g_d(0),g_d(s),...g_d(ds)\]

求

\[g_{2d}(0),g_{2d}(s),g_{2d}(2s),...,g_{2d}(2ds)\]

通过这个操作我们可以把\(d\)乘\(2\)

已知

\[g_d(0),g_d(s),...g_d(ds)\]

求

\[g_{d+1}(0),g_{d+1}(s),...g_{d+1}((d+1)s)\]

通过这个操作我们可以把\(d\)加\(1\)

然后就可以用一个类似快速幂的迭代来求出\(g_s\)了

把\(d\)加\(1\)

这个简单一点先说这个好了

对于\(g_{d+1}((d+1)s)\)，直接暴力计算就行了

对于剩下的项，也可以直接暴力计算，因为有

\[g_{d+1}(x)=g_{d}(x)\times (x+d+1)\]

那么一次的复杂度就是\(O(s)\)的

把\(d\)乘\(2\)

这种情况就比较辣手了

具体来说，我们已知

\[g_d(0),g_d(s),..,g_d(ds)\]

需要求

\[g_{2d}(0),g_{2d}(s),g_{2d}(2s),...,g_{2d}(2ds)\]

首先有

\[g_{2d}(x)=g_d(x)g_d(x+d)\]

我们先考虑一下，已知

\[g_d(0),g_d(s),..,g_d(ds)\]

如何求出

\[g_d((d+1)s),g_d((d+2)s),..,g_d((d+d)s)\]

我们可以设\(h(i)=g_d(is)\)，那么问题可以转化为已知

\[h(0),h(1),...,h(d)\]

求

\[h(d+1),h(d+2),...,h(2d)\]

另一个问题就是已知

\[g_d(0),g_d(s),..,g_d(2ds)\]

求

\[g_{2d}(0),g_{2d}(s),...,g_{2d}(2ds)\]

也可以转化为已知

\[h(0),h(s),...,h(2d)\]

求

\[h(d/s),h(d/s+1),h(d/s+2),...,h(d/s+2d)\]

总和起来，问题就是已知

\[h(0),h(1),...,h(d)\]

我们需要求得

\[h(k),k(1+k),...,h(d+k)\]

根据拉格朗日差值公式，有

\[ \begin{aligned} h(m+k) &=\sum_{i=0}^n h(i)\prod_{j\neq i}{m+k-j\over i-j}\\ &=\left(\prod_{j=0}^n (m+k-j)\right)\sum_{i=0}^d {h(i)\over i!(d-i)!(-1)^{d-i}}{1\over m+k-i}\\ \end{aligned} \]

后面明显是个卷积的形式，前面是一段连续区间的乘积，可以用双指针维护

然而万一这里\(m+k-i=0\)怎么办？

因为 \(h(x)\) 本质上代表了一段连续数字的乘积，而根据我们算法的原理， \(h(m+k)\) 和 \(h(i)\) 一定代表了两段不同的数字，所以分母不可能为0

那么我们就可以把\(d\)乘\(2\)了

然而这里有两个细节问题

第一，如果直接求逆元的话是\(O(d\log p)\)的，不过其实有\(O(d)\)求\(n\)个数逆元的方法，感兴趣的可以自己去看

第二，虽然\(m\)很小，但是\(m+k\)很大，我们不能直接暴力卷积。注意到后面那个卷积柿子中当\(i>n\)时可以默认\(h(i)=0\)，所以有用的项其实只有\(m=0,{1\over m+k-d}\)，以及\(m=d,{1\over m+k}\)之间的数，也就是说后面有用的项不超过\(O(d)\)个，所以我们需要把多项式平移一下再来卷积

然后没有然后了

不过可以用威尔逊定理优化一下常数就是了

//minamoto#include
    
     #define R register#define ll long long#define fp(i,a,b) for(R int i=(a),I=(b)+1;i
     
      I;--i)#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)using namespace std;const int N=(1<<17)+5;int P;inline int add(R int x,R int y){return 0ll+x+y>=P?0ll+x+y-P:x+y;}inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}int ksm(R int x,R int y){    R int res=1;    for(;y;y>>=1,x=mul(x,x))(y&1)?res=mul(res,x):0;    return res;}const double Pi=acos(-1.0);struct cp{    double x,y;    inline cp(){}    inline cp(R double xx,R double yy):x(xx),y(yy){}    inline cp operator +(const cp &b)const{return cp(x+b.x,y+b.y);}    inline cp operator -(const cp &b)const{return cp(x-b.x,y-b.y);}    inline cp operator *(const cp &b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}    inline cp operator *(const double &b)const{return cp(x*b,y*b);}    inline cp operator ~()const{return cp(x,-y);}}w[2][N];int r[21][N],ifac[N],lg[N],inv[N];double iv[21];void Pre(){    iv[0]=1;    fp(d,1,17){        fp(i,0,(1<
      
       >1]>>1)|((i&1)<<(d-1));        lg[1<
       
        >16,a[i]&65535),g[i]=cp(b[i]>>16,b[i]&65535);    fp(i,len,lim-1)f[i]=g[i]=cp(0,0);    FFT(f,1),FFT(g,1);    fp(i,0,lim-1){        cp t,f0,f1,g0,g1;        t=~f[i?lim-i:0],f0=(f[i]-t)*cp(0,-0.5),f1=(f[i]+t)*0.5;        t=~g[i?lim-i:0],g0=(g[i]-t)*cp(0,-0.5),g1=(g[i]+t)*0.5;        p[i]=f1*g1,q[i]=f1*g0+f0*g1+f0*g0*cp(0,1);    }    FFT(p,0),FFT(q,0);    fp(i,0,lim-1)c[i]=((((ll)(p[i].x+0.5)%P<<16)%P<<16)+((ll)(q[i].x+0.5)<<16)+((ll)(q[i].y+0.5)))%P;}void calc(int *a,int *b,int n,int k){    static int f[N],g[N],h[N],sum[N],isum[N];    int len=1;while(len<=n+n)len<<=1;    fp(i,0,n)f[i]=mul(a[i],mul(ifac[i],ifac[n-i]));    for(R int i=n-1;i>=0;i-=2)f[i]=P-f[i];    int t=dec(k,n);    fp(i,0,n+n)g[i]=add(i,t);    sum[0]=g[0];fp(i,1,n+n)sum[i]=mul(sum[i-1],g[i]);    isum[n+n]=ksm(sum[n+n],P-2);    fd(i,n+n,1)isum[i-1]=mul(isum[i],g[i]);    fp(i,1,n+n)g[i]=mul(isum[i],sum[i-1]);g[0]=isum[0];    fp(i,n+1,len-1)f[i]=0;fp(i,n+n+1,len-1)g[i]=0;        MTT(f,g,len,h);    int res=1,p1=k-n,p2=k;    fp(i,p1,p2)res=1ll*res*i%P;    res=dec(res,0);        fp(i,0,n)g[i]=(0ll+P+p1+i)%P;    sum[0]=g[0];fp(i,1,n)sum[i]=mul(sum[i-1],g[i]);    isum[n]=ksm(sum[n],P-2);    fd(i,n,1)isum[i-1]=mul(isum[i],g[i]);    fp(i,1,n)g[i]=mul(isum[i],sum[i-1]);g[0]=isum[0];        for(R int i=0;i<=n;p2=add(p2,1),++i)        b[i]=mul(h[i+n],res),res=mul(res,mul(g[i],p2+1));}int solve(int bl){    static int a[N],b[N],c[N];    int s=0;for(int p=bl;p;p>>=1)++s;a[0]=1,--s;    int qwq=ksm(bl,P-2);    for(int p=0;s>=0;--s){        if(p){            calc(a,b,p,p+1);            fp(i,0,p)a[p+i+1]=b[i];a[p<<1|1]=0;            calc(a,b,p<<1,mul(p,qwq));            p<<=1;fp(i,0,p)a[i]=mul(a[i],b[i]);        }        if(bl>>s&1){            fp(i,0,p)a[i]=mul(a[i],(1ll*bl*i+p+1)%P);            p|=1,a[p]=1;            fp(i,1,p)a[p]=mul(a[p],(1ll*bl*p+i)%P);        }    }    int res=1;    fp(i,0,bl-1)res=mul(res,a[i]);    return res;}int GetFac(int n){    int s=sqrt(n),res=solve(s);    fp(i,s*s+1,n)res=mul(res,i);    return res;}int Fac(int n){    if(n>P-1-n){        int res=ksm(GetFac(P-1-n),P-2);        return (P-1-n)&1?res:P-res;    }    return GetFac(n);}int n;int main(){//  freopen("testdata.in","r",stdin);    scanf("%d%d",&n,&P),Pre();    printf("%d\n",Fac(n));    return 0;}