fork download
copy 
  1. // luogu-judger-enable-o2
  2. #include<cstdio>
  3. #include<algorithm>
  4. #include<cmath>
  5. using namespace std;const int N=262144+10;typedef unsigned long long ll;
  6. const int P=65536;const int SF=16;const int msk=65535;ll mod;ll PP;
  7. typedef long double ld;const ld pi=acos(-1.0);
  8. inline ll po(ll a,ll p){ll r=1;for(;p;p>>=1,a=a*a%mod)if(p&1)r=r*a%mod;return r;}
  9. struct cmp
  10. {
  11.     ld r;ld v;
  12.     friend cmp operator +(cmp a,cmp b){return (cmp){a.r+b.r,a.v+b.v};}
  13.     friend cmp operator -(cmp a,cmp b){return (cmp){a.r-b.r,a.v-b.v};}
  14.     friend cmp operator *(cmp a,cmp b){return (cmp){a.r*b.r-a.v*b.v,a.r*b.v+a.v*b.r};}
  15.     void operator /=(const int& len){r/=len;v/=len;}
  16. }rt[2][22][N],tr[N],tr1[N],tr2[N],tr3[N],tr4[N],tr5[N],tr6[N];
  17. int rv[22][N];ll m13[N],m14[N],m23[N],m24[N];
  18. inline void pre()
  19. {
  20.     for(int d=1;d<=18;d++)
  21.         for(int i=1;i<(1<<d);i++)rv[d][i]=(rv[d][i>>1]>>1)|((i&1)<<(d-1));
  22.     for(int d=1,t=1;d<=18;d++,t<<=1)
  23.         for(int i=0;i<(1<<d);i++)rt[0][d][i]=(cmp){cos(pi*i/t),sin(pi*i/t)};
  24.     for(int d=1,t=1;d<=18;d++,t<<=1)
  25.         for(int i=0;i<(1<<d);i++)rt[1][d][i]=(cmp){cos(pi*i/t),-sin(pi*i/t)};
  26. }inline void fft(cmp* a,int len,int d,int o)
  27. {
  28.     for(int i=1;i<len;i++)if(i<rv[d][i])swap(a[i],a[rv[d][i]]);cmp* w;int i;
  29.     for(int k=1,j=1;k<len;k<<=1,j++)
  30.         for(int s=0;s<len;s+=(k<<1))
  31.             for(i=s,w=rt[o][j];i<s+k;i++,++w)
  32.                 {cmp a1=a[i+k]*(*w);a[i+k]=a[i]-a1;a[i]=a[i]+a1;}
  33.     if(o)for(int i=0;i<len;i++)a[i]/=len;
  34. }inline void dbdft(ll* a,int len,int d,cmp* op1,cmp* op2)
  35. {
  36.     for(int i=0;i<len;i++)tr[i]=(cmp){(ld)(a[i]>>SF),(ld)(a[i]&msk)};
  37.     fft(tr,len,d,0);tr[len]=tr[0];
  38.     for(cmp* p1=tr,*p2=tr+len,*p3=op1;p1!=tr+len;++p1,--p2,++p3)
  39.         (*p3)=(cmp){p1->r+p2->r,p1->v-p2->v}*(cmp){0.5,0};
  40.     for(cmp* p1=tr,*p2=tr+len,*p3=op2;p1!=tr+len;++p1,--p2,++p3)
  41.         (*p3)=(cmp){p1->r-p2->r,p1->v+p2->v}*(cmp){0,-0.5};
  42. }inline void dbidft(cmp* tr,int len,int d,ll* a,ll* b)
  43. {
  44.     fft(tr,len,d,1);
  45.     for(int i=0;i<len;i++)a[i]=(ll)(tr[i].r+0.5)%mod;
  46.     for(int i=0;i<len;i++)b[i]=(ll)(tr[i].v+0.5)%mod;
  47. }inline void poly_mul(ll* a,ll* b,ll* c,int len,int d)//以上都是任意模数fft的板子 
  48. {
  49.     dbdft(a,len,d,tr1,tr2);dbdft(b,len,d,tr3,tr4);
  50.     for(int i=0;i<len;i++)tr5[i]=tr1[i]*tr3[i]+(cmp){0,1}*tr2[i]*tr4[i];
  51.     for(int i=0;i<len;i++)tr6[i]=tr2[i]*tr3[i]+(cmp){0,1}*tr1[i]*tr4[i];
  52.     dbidft(tr5,len,d,m13,m24);dbidft(tr6,len,d,m23,m14);
  53.     for(int i=0;i<len;i++)c[i]=m13[i]*PP%mod;
  54.     for(int i=0;i<len;i++)(c[i]+=(m23[i]+m14[i])*P+m24[i])%=mod;
  55. }namespace iter
  56. {
  57.     ll f[N];ll g[N];ll h[N];ll ifac[N];
  58.     inline void ih()
  59.     {
  60.         ifac[0]=ifac[1]=1;
  61.         for(int i=2;i<min((ll)N,mod);i++)ifac[i]=(mod-mod/i)*ifac[mod%i]%mod;
  62.         for(int i=1;i<min((ll)N,mod);i++)(ifac[i]*=ifac[i-1])%=mod;
  63.     }inline void calch(ll del,int cur,ll* ip,ll* op)
  64.     {
  65.         int d=0;int len=1;while(len<=cur+cur+cur)len<<=1,d++;
  66.         for(int i=0;i<=cur;i++)f[i]=ip[i]*ifac[i]%mod*ifac[cur-i]%mod;
  67.         for(int i=cur-1;i>=0;i-=2)f[i]=(mod-f[i])%mod;
  68.         for(int i=0;i<=cur+cur;i++)g[i]=po((del+mod-cur+i)%mod,mod-2); 
  69.         for(int i=cur+1;i<len;i++)f[i]=0;for(int i=cur+cur+1;i<len;i++)g[i]=0;
  70.         poly_mul(f,g,h,len,d);//卷积求出h' 
  71.         ll xs=1;ll p1=del-cur;ll p2=del;
  72.         for(int i=p1;i<=p2;i++)(xs*=i)%=mod;
  73.         for(int i=0;i<=cur;i++,p1++,p2++)//双指针求出系数 
  74.         {
  75.             op[i]=h[i+cur]*xs%mod;
  76.             (xs*=po(p1,mod-2))%=mod,(xs*=(p2+1))%=mod;
  77.         }  
  78.     }
  79. }ll val[N];ll fv1[N];ll fv2[N];
  80. inline void solve(int n)//倍增 
  81. {
  82.     int hb=0;for(int p=n;p;p>>=1)hb++;val[0]=1;
  83.     for(int z=hb,cur=0;z>=0;z--)
  84.     {
  85.         if(cur!=0)//把d乘2 
  86.         {
  87.             iter::calch(cur+1,cur,val,fv1);
  88.             for(int i=0;i<=cur;i++)val[cur+i+1]=fv1[i];val[cur<<1|1]=0;
  89.             iter::calch(cur*po(n,mod-2)%mod,cur<<1,val,fv2);
  90.             cur<<=1;for(int i=0;i<=cur;i++)(val[i]*=fv2[i])%=mod;
  91.         }if((n>>z)&1)//把d加1 
  92.         {
  93.             for(int i=0;i<=cur;i++)(val[i]*=(ll)(n*i)+cur+1)%=mod;cur|=1;val[cur]=1;
  94.             for(int i=1;i<=cur;i++)(val[cur]*=(ll)cur*n+i)%=mod;
  95.         }
  96.     }
  97. }
  98. void solve()
  99. {
  100.     pre();int n;scanf("%d%lld",&n,&mod);iter::ih();
  101.     int bl=sqrt(n);PP=(ll)P*P%mod;solve(bl);ll res=1;
  102.     for(int i=0,id=0;;i+=bl,id++)//分块 
  103.     {
  104.         if((ll)i+bl>n){for(int j=i+1;j<=n;j++)(res*=j)%=mod;break;}
  105.         (res*=val[id])%=mod;
  106.     }printf("%lld\n",res);
  107. }
  108. int main()
  109. {
  110. 	int T;scanf("%d",&T);while(T--)solve();
  111. 	return 0;//拜拜程序~ 
  112. }
Success #stdin #stdout 2.08s 162964KB
comments (?)
stdin
copy 
3
1073741823 2147483629
1073741823 2147481811
2147483628 2147483629
stdout
copy 
0
0
2147483628
https://ideone.com/QluBuZ
hack for FACTMODP
language:
C++14 (gcc 8.3)
created:
3 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!