三个模数NTT,中国剩余定理转换答案即可
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e6+10;
const ll mod1=998244353,mod2=1004535809,mod3=469762049,G=3;
ll mul(ll a,ll k,ll mod)
{ll ans=1;a%=mod;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll mod1_2=mod1*mod2;
ll inv1=mul(mod1,mod2-2,mod2),inv2=mul(mod1_2,mod3-2,mod3);
int n,m;ll p;
struct Int{ll a[4];Int(){}Int(ll x){a[1]=a[2]=a[3]=x;}Int(ll x,ll y,ll z){a[1]=x,a[2]=y,a[3]=z;}inline Int operator*(const Int x)const{return Int(a[1]*x.a[1]%mod1,a[2]*x.a[2]%mod2,a[3]*x.a[3]%mod3);}inline Int operator+(const Int x)const{return Int((a[1]+x.a[1])%mod1,(a[2]+x.a[2])%mod2,(a[3]+x.a[3])%mod3);}inline Int operator-(const Int x)const{return Int((a[1]-x.a[1]+mod1)%mod1,(a[2]-x.a[2]+mod2)%mod2,(a[3]-x.a[3]+mod3)%mod3);}ll get(){ll ans=(a[2]-a[1]+mod2)*inv1%mod2*mod1+a[1];ans=(a[3]-ans%mod3+mod3)%mod3*inv2%mod3*(mod1_2%p)%p+ans;return ans%p;}
};
int lim,rev[N];
Int wn[N];
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);wn[0]=Int(1);Int inv=Int(mul(G,(mod1-1)/lim,mod1),mul(G,(mod2-1)/lim,mod2),mul(G,(mod3-1)/lim,mod3));for(int i=1;i<=lim;i++)wn[i]=wn[i-1]*inv;
}
void NTT(Int *A,int flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){int t=lim/mid>>1;for(int i=0;i<lim;i+=mid<<1)for(int j=0;j<mid;j++){Int w=wn[t*j];Int X=A[i+j],Y=A[i+j+mid]*w;A[i+j]=X+Y,A[i+j+mid]=X-Y;}}if(flag)return ;Int inv=Int(mul(lim,mod1-2,mod1),mul(lim,mod2-2,mod2),mul(lim,mod3-2,mod3));reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=A[i]*inv;
}
ll fw;
Int A[N],B[N];
int main()
{scanf("%d%d%lld",&n,&m,&p);n++;m++;for(int i=0;i<n;i++)scanf("%lld",&fw),A[i]=Int(fw);for(int i=0;i<m;i++)scanf("%lld",&fw),B[i]=Int(fw);init(n+m);NTT(A,1),NTT(B,1);for(int i=0;i<lim;i++)A[i]=A[i]*B[i];NTT(A,0);for(int i=0;i<n+m-1;i++){printf("%lld ",A[i].get());}
}
CDQ分治的思路
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e5+10;
const int mod=998244353,G=3;
ll quick_mul(ll a,ll k)
{ll ans=1;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll inv(ll a){return quick_mul(a,mod-2);}
int lim,rev[N];
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*w%mod){ll X=A[i+j],Y=1ll*A[i+j+mid]*fw%mod;A[i+j]=(X+Y)%mod,A[i+j+mid]=(X-Y+mod)%mod;}}}if(flag)return ;ll inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=A[i]*inv_%mod;
}
int F[N],H[N];
void solve(int *A,int *B,int l,int r)
{if(l==r){if(!l)B[l]=1;return ;}int mid=(l+r)>>1;solve(A,B,l,mid);init(r-l+1);for(int i=0;i<lim;i++)F[i]=H[i]=0;for(int i=l;i<=mid;i++)F[i-l]=B[i];for(int i=1;i<=r-l;i++)H[i]=A[i];NTT(F,1);NTT(H,1);for(int i=0;i<lim;i++)F[i]=1ll*F[i]*H[i]%mod;NTT(F,0);for(int i=mid+1;i<=r;i++)B[i]=(B[i]+F[i-l])%mod;solve(A,B,mid+1,r);
}
int n,a[N],b[N];
int main()
{scanf("%d",&n);for(int i=1;i<n;i++)scanf("%d",&a[i]);solve(a,b,0,n-1);for(int i=0;i<n;i++)printf("%d ",b[i]);
}
牛顿迭代老好用了
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e5+10,mod=998244353,G=3;
ll quick_mul(ll a,ll k)
{ll ans=1;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll inv(ll a){return quick_mul(a,mod-2);}
int lim,rev[N];
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*w%mod){ll x=A[i+j],y=A[i+j+mid]*fw%mod;A[i+j]=(x+y)%mod;A[i+j+mid]=(x-y+mod)%mod; }}}if(flag)return ;ll inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=A[i]*inv_%mod;
}
int F[N];
void poly_inv(int *A,int *B,int len)
{if(len==1){B[0]=inv(1ll*A[0]);return ;}poly_inv(A,B,(len+1)>>1);init(len<<1);for(int i=0;i<len;i++)F[i]=A[i];for(int i=len;i<lim;i++)F[i]=0;NTT(F,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=(2ll-1ll*F[i]*B[i]%mod+mod)*B[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;
}
int n;
int a[N],b[N];
int main()
{scanf("%d",&n);for(int i=0;i<n;i++)scanf("%d",&a[i]);poly_inv(a,b,n);for(int i=0;i<n;i++)printf("%d ",b[i]);
}
求导积分记一下
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e5+10,mod=998244353,G=3;
ll quick_mul(ll a,ll k)
{ll ans=1;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll inv(ll a){return quick_mul(a,mod-2);}
int lim,rev[N];
void pr(int *A){for(int i=0;i<lim;i++)printf("%d ",A[i]);printf("n");}
void Direv(int *A,int *B,int n){for(int i=1;i<n;i++)B[i-1]=1ll*A[i]*i%mod;B[n-1]=0;}
void Inter(int *A,int *B,int n){for(int i=1;i<n;i++)B[i]=A[i-1]*inv(i)%mod;B[0]=0;}
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*w%mod){ll x=A[i+j],y=A[i+j+mid]*fw%mod;A[i+j]=(x+y)%mod,A[i+j+mid]=(x-y+mod)%mod;}}}if(flag)return ;ll inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=A[i]*inv_%mod;
}
int F[N];
void poly_inv(int *A,int *B,int len)
{if(len==1){B[0]=inv(1ll*A[0]);return ;}poly_inv(A,B,(len+1)>>1);init(len<<1);for(int i=0;i<len;i++)F[i]=A[i];for(int i=len;i<lim;i++)F[i]=0;NTT(F,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=(2ll-1ll*F[i]*B[i]%mod+mod)%mod*B[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;
}
int H1[N],H2[N];
void poly_ln(int *A,int *B,int len)
{Direv(A,H1,len);poly_inv(A,H2,len);init(len<<1);NTT(H1,1);NTT(H2,1);for(int i=0;i<lim;i++)H1[i]=1ll*H1[i]*H2[i]%mod;NTT(H1,0);Inter(H1,B,len);
}
int n,a[N],b[N];
int main()
{scanf("%d",&n);for(int i=0;i<n;i++)scanf("%d",&a[i]);poly_ln(a,b,n);for(int i=0;i<n;i++)printf("%d ",b[i]);
}
临时数组千万要清零
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e5+10,mod=998244353,G=3;
ll quick_mul(ll a,ll k)
{ll ans=1;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll inv(ll a){return quick_mul(a,mod-2);}
int lim,rev[N];
void pr(int *A){for(int i=0;i<lim;i++)printf("%d ",A[i]);printf("n");}
void Direv(int *A,int *B,int len){for(int i=1;i<len;i++)B[i-1]=1ll*A[i]*i%mod;B[len-1]=0;}
void Inter(int *A,int *B,int len){for(int i=1;i<len;i++)B[i]=inv(i)*A[i-1]%mod;B[0]=0;}
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*w%mod){ll x=A[i+j],y=fw*A[i+j+mid]%mod;A[i+j]=(x+y)%mod,A[i+j+mid]=(x-y+mod)%mod;}}}if(flag)return ;ll inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=inv_*A[i]%mod;
}
int F[N];
void poly_inv(int *A,int *B,int len)
{if(len==1){B[0]=inv(A[0]);return ;}poly_inv(A,B,(len+1)>>1);init(len<<1);for(int i=0;i<len;i++)F[i]=A[i];for(int i=len;i<lim;i++)F[i]=0;NTT(F,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=(2ll-1ll*F[i]*B[i]%mod+mod)%mod*B[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F[i]=0;
}
int H1[N],H2[N];
void poly_ln(int *A,int *B,int len)
{Direv(A,H1,len);poly_inv(A,H2,len);init(len<<1);NTT(H1,1);NTT(H2,1);for(int i=0;i<lim;i++)H1[i]=1ll*H1[i]*H2[i]%mod;NTT(H1,0);Inter(H1,B,len);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)H1[i]=H2[i]=0;
}
int F1[N],F2[N];
void poly_exp(int *A,int *B,int len)
{if(len==1){B[0]=1ll;return ;}poly_exp(A,B,(len+1)>>1);init(len<<1);poly_ln(B,F1,len);for(int i=0;i<len;i++)F1[i]=(1ll*A[i]-F1[i]+mod)%mod;for(int i=len;i<lim;i++)F1[i]=0;F1[0]=(F1[0]+1)%mod;NTT(F1,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=1ll*F1[i]*B[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F1[i]=0;
}
int n,a[N],b[N];
int main()
{scanf("%d",&n);for(int i=0;i<n;i++)scanf("%d",&a[i]);poly_exp(a,b,n);for(int i=0;i<n;i++)printf("%d ",b[i]);
}
e x p ( k ∗ l n ( A ( x ) ) ) exp(k*ln(A(x))) exp(k∗ln(A(x)))
注意此时k是作为系数来乘,应%mod而不是%mod-1
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e5+10,mod=998244353,G=3;
int read()
{int ans=0,f=1;char s=getchar();while(!isdigit(s)){if(s=='-')f=-1;s=getchar();}while(isdigit(s)){ans=(1ll*ans*10%mod+s-'0'+mod)%mod;s=getchar();}return ans*f;
}
ll quick_mul(ll a,ll k)
{ll ans=1;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll inv(ll a){return quick_mul(a,mod-2);}
int lim,rev[N];
void Direv(int *A,int *B,int n){for(int i=1;i<n;i++)B[i-1]=1ll*A[i]*i%mod;B[n-1]=0;}
void Inter(int *A,int *B,int n){for(int i=1;i<n;i++)B[i]=inv(i)*A[i-1]%mod;B[0]=0;}
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*w%mod){ll x=A[i+j],y=fw*A[i+j+mid]%mod;A[i+j]=(x+y)%mod,A[i+j+mid]=(x-y+mod)%mod;}}}if(flag)return ;ll inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=inv_*A[i]%mod;
}
int F1[N];
void poly_inv(int *A,int *B,int len)
{if(len==1){B[0]=inv(A[0]);return ;}poly_inv(A,B,(len+1)>>1);init(len<<1);for(int i=0;i<len;i++)F1[i]=A[i];for(int i=len;i<lim;i++)F1[i]=0;NTT(F1,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=(2ll-1ll*F1[i]*B[i]%mod+mod)%mod*B[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F1[i]=0;
}
int F2[N],F3[N];
void poly_ln(int *A,int *B,int len)
{Direv(A,F2,len);poly_inv(A,F3,len);init(len<<1);NTT(F2,1);NTT(F3,1);for(int i=0;i<lim;i++)F2[i]=1ll*F2[i]*F3[i]%mod;NTT(F2,0);Inter(F2,B,len);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F2[i]=F3[i]=0;
}
int F4[N];
void poly_exp(int *A,int *B,int len)
{if(len==1){B[0]=1;return ;}poly_exp(A,B,(len+1)>>1);poly_ln(B,F4,len);init(len<<1);for(int i=0;i<len;i++)F4[i]=(A[i]-F4[i]+mod)%mod;F4[0]=(F4[0]+1)%mod;NTT(B,1);NTT(F4,1);for(int i=0;i<lim;i++)B[i]=1ll*B[i]*F4[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F4[i]=0;
}
int F5[N];
void poly_mul(int *A,int *B,int len,int k)
{poly_ln(A,F5,len);for(int i=0;i<len;i++)F5[i]=1ll*F5[i]*k%mod;poly_exp(F5,B,len);
}
int n,k,a[N],b[N];
int main()
{n=read();k=read();for(int i=0;i<n;i++)a[i]=read();poly_mul(a,b,n,k);for(int i=0;i<n;i++)printf("%d ",b[i]);
}
不保证 a 0 = 1 a_0=1 a0=1,就全部除以 a 0 a_0 a0在快速幂,最后乘上 a 0 k a_0^k a0k,这个k就应该%mod-1了
要是 a 0 = 0 a_0=0 a0=0需要将多项式左移,最后再移回去
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=4e5+10,mod=998244353,G=3;
int num(char *s,int p)
{int ans=0;for(int i=0;s[i];i++)ans=(1ll*ans*10%p+s[i]-'0'+p)%p;return ans;
}
ll quick_mul(ll a,ll k)
{ll ans=1;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll inv(ll a){return quick_mul(a,mod-2);}
int lim,rev[N];
ll inv_num[N];
void pre_work()
{for(int i=1;i<N-10;i++)inv_num[i]=inv(i);
}
void Direv(int *A,int *B,int n){for(int i=1;i<n;i++)B[i-1]=1ll*A[i]*i%mod;B[n-1]=0;}
void Inter(int *A,int *B,int n){for(int i=1;i<n;i++)B[i]=inv_num[i]*A[i-1]%mod;B[0]=0;}
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*w%mod){ll x=A[i+j],y=fw*A[i+j+mid]%mod;A[i+j]=(x+y)%mod;A[i+j+mid]=(x-y+mod)%mod;}}}if(flag)return ;ll inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=inv_*A[i]%mod;
}
int F1[N];
void poly_inv(int *A,int *B,int len)
{if(len==1){B[0]=inv(A[0]);return ;}poly_inv(A,B,(len+1)>>1);init(len<<1);for(int i=0;i<len;i++)F1[i]=A[i];for(int i=len;i<lim;i++)F1[i]=0;NTT(F1,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=(2ll-1ll*F1[i]*B[i]%mod+mod)%mod*B[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F1[i]=0;
}
int F2[N],F3[N];
void poly_ln(int *A,int *B,int len)
{Direv(A,F2,len);poly_inv(A,F3,len);init(len<<1);NTT(F2,1);NTT(F3,1);for(int i=0;i<lim;i++)F2[i]=1ll*F2[i]*F3[i]%mod;NTT(F2,0);Inter(F2,B,len);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F2[i]=F3[i]=0;
}
int F4[N];
void poly_exp(int *A,int *B,int len)
{if(len==1){B[0]=1;return ;}poly_exp(A,B,(len+1)>>1);poly_ln(B,F4,len);init(len<<1);for(int i=0;i<len;i++)F4[i]=(A[i]-F4[i]+mod)%mod;for(int i=len;i<lim;i++)F4[i]=0;F4[0]=(F4[0]+1)%mod;NTT(F4,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=1ll*B[i]*F4[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F4[i]=0;
}
int F5[N],F6[N];
void poly_ex_mul(int *A,int *B,int len,char *s)
{bool f=1;int g=0;for(int i=0;i<len;i++)if(A[i]){g=i;f=0;break;}if(f){if(num(s,mod-1)==0)B[0]=1;return ;}ll kt=0;for(int i=0;s[i];i++)if((ll)(kt*10+s[i]-'0')<=mod)kt=kt*10+s[i]-'0';if(1ll*g*kt>=1ll*len)return ;int n=len-g,k1=num(s,mod-1),k2=num(s,mod);ll inv_=inv(A[g]),x=quick_mul(A[g],k1);for(int i=0;i<n;i++)F5[i]=inv_*A[i+g]%mod;poly_ln(F5,F6,n);for(int i=0;i<n;i++)F6[i]=1ll*k2*F6[i]%mod;poly_exp(F6,B,n);for(int i=0;i<n;i++)B[i]=x*B[i]%mod;ll num=min(1ll*g*kt,1ll*len);for(int i=len-1;i>=num;i--)B[i]=B[i-num];for(int i=num-1;i>=0;i--)B[i]=0;
}
int n,a[N],b[N];
char s[N];
int main()
{pre_work();scanf("%d%s",&n,s);for(int i=0;i<n;i++)scanf("%d",&a[i]);poly_ex_mul(a,b,n,s);for(int i=0;i<n;i++)printf("%d ",b[i]);
}
远古代码了。。。
应该可以快速幂?
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=8e5+10;
const ll G=3,mod=998244353;
void pr(ll *a){for(int i=0;i<=10;i++)printf("%lld ",a[i]);printf("n");}
ll inv2;
ll quick_mul(ll a,ll k)
{ll as=1;while(k){if(k&1)as=as*a%mod;a=a*a%mod;k>>=1;}return as;
}
ll powM(ll a){return quick_mul(a,mod-2);}
void dao(ll *A,int n)
{for(int i=1;i<n;i++)A[i]=A[i-1]*i%mod;A[n-1]=0;
}//求导
void jifen(ll *A,int n)
{for(int i=n;i;i--)A[i]=A[i-1]*powM(i)%mod;A[0]=0;
}//积分
int lim;
int rev[N];
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
}
void NTT(ll *A,int flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll wn=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*wn%mod){ll x=A[i+j],y=A[i+j+mid]*fw%mod;A[i+j]=(x+y)%mod,A[i+j+mid]=(x-y+mod)%mod;}}}if(flag==1)return ;ll inv=powM(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=A[i]*inv%mod;
}
ll ls[N];
void inv_poly(ll *a,ll *b,int len)
{if(len==1){b[0]=powM(a[0]);return ;}inv_poly(a,b,(len+1)>>1);init(len<<1);for(int i=0;i<len;i++) ls[i]=a[i];for(int i=len;i<lim;i++)ls[i]=0;NTT(ls,1);NTT(b,1);for(int i=0;i<lim;i++)b[i]=(2ll-ls[i]*b[i]%mod+mod)%mod*b[i]%mod;NTT(b,-1);for(int i=len;i<lim;i++)b[i]=0;
}
ll lss[N];
void sqr_poly(ll *a,ll *b,int len)
{if(len==1){b[0]=1;return ;}sqr_poly(a,b,(len+1)>>1);for(int i=0;i<len<<1;i++)lss[i]=0;//lss需要清零inv_poly(b,lss,len);init(len<<1);for(int i=0;i<len;i++)ls[i]=a[i];for(int i=len;i<lim;i++)ls[i]=0;NTT(ls,1);NTT(lss,1);NTT(b,1);for(int i=0;i<lim;i++)b[i]=(b[i]+ls[i]*lss[i]%mod)%mod*inv2%mod;NTT(b,-1);for(int i=len;i<lim;i++)b[i]=0;
}
int n;
ll a[N],b[N];
int main()
{scanf("%d",&n);inv2=powM(2ll);for(int i=0;i<n;i++)scanf("%lld",&a[i]);sqr_poly(a,b,n);for(int i=0;i<n;i++)printf("%lld ",b[i]);
}
学会背 了二次剩余
#include<bits/stdc++.h>
using namespace std;
const int N=8e5+10;
const int G=3,mod=998244353;
void pr(int *a){for(int i=0;i<=10;i++)printf("%d ",a[i]);printf("n");}
int inv2;
int add(int a,int b) {return a+b>=mod?a+b-mod:a+b;}
int del(int a,int b) {return a<b?a+mod-b:a-b;}
int mul(int a,int b) {return 1ll*a*b%mod;}
int quick_mul(int a,int k,int ans=1){for(;k;k>>=1,a=mul(a,a))if(k&1)ans=mul(ans,a);return ans;}
int w;
struct nd{int x,y;nd(){}nd(int a,int b){x=a;y=b;}nd operator*(const nd &a)const{return nd(add(mul(x,a.x),mul(mul(y,a.y),w)),add(mul(x,a.y),mul(y,a.x)));}
};
nd quick_mul(nd a,int k,nd ans=nd(1,0)){for(;k;k>>=1,a=a*a)if(k&1)ans=ans*a;return ans;}
bool check_sqrt(int a){return quick_mul(a,(mod-1)/2)==1;}
int get_w(int n)
{for(int a=rand();;a=rand())if(!check_sqrt(del(mul(a,a),n)))return w=del(mul(a,a),n),a;
}
int get_sqrt(int n)
{int a=get_w(n);int x1=(quick_mul(nd(a,1),(mod+1)/2)).x;int x2=mod-x1;return x1<x2?x1:x2;
}
int inv(int a){return quick_mul(a,mod-2);}
int lim,rev[N];
void Direv(int *A,int *B,int n){for(int i=1;i<n;i++)B[i-1]=mul(A[i],i);B[n-1]=0;}
void Inter(int *A,int *B,int n){for(int i=1;i<n;i++)B[i]=mul(inv(i),A[i-1]);B[0]=0;}
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){int w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){int fw=1;for(int j=0;j<mid;j++,fw=mul(fw,w)){int x=A[i+j],y=mul(fw,A[i+j+mid]);A[i+j]=add(x,y),A[i+j+mid]=del(x,y);}}}if(flag)return ;int inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=mul(inv_,A[i]);
}
void poly_inv(int *A,int *B,int len)
{if(len==1){B[0]=inv(A[0]);return ;}static int F[N];poly_inv(A,B,(len+1)>>1);for(int i=0;i<len;i++)F[i]=A[i];init(len<<1);NTT(F,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=mul(del(2,mul(F[i],B[i])),B[i]);NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F[i]=0;
}
void poly_sqrt(int *A,int *B,int len)
{if(len==1){B[0]=get_sqrt(A[0]);return ;}static int F[N],G[N];poly_sqrt(A,B,(len+1)>>1);poly_inv(B,F,len);for(int i=0;i<len;i++)G[i]=A[i];init(len<<1);NTT(G,1);NTT(B,1);NTT(F,1);for(int i=0;i<lim;i++)B[i]=mul((add(mul(G[i],F[i]),B[i])),inv2);NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F[i]=G[i]=0;
}
int n;
int a[N],b[N];
int main()
{scanf("%d",&n);inv2=inv(2ll);for(int i=0;i<n;i++)scanf("%d",&a[i]);poly_sqrt(a,b,n);for(int i=0;i<n;i++)printf("%d ",b[i]);
}
翻转小心点吧
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e5+10,mod=998244353,G=3;
ll quick_mul(ll a,ll k)
{ll ans=1;while(k){if(k&1)ans=ans*a%mod;a=a*a%mod;k>>=1;}return ans;
}
ll inv(ll a){return quick_mul(a,mod-2);}int lim,rev[N];
void Direv(int *A,int *B,int n){for(int i=1;i<n;i++)B[i-1]=1ll*A[i]*i%mod;B[n-1]=0;}
void Inter(int *A,int *B,int n){for(int i=1;i<n;i++)B[i]=inv(i)*A[i-1]%mod;B[0]=0;}
void init(int n)
{for(lim=1;lim<n;lim<<=1);for(int i=0;i<lim;i++)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
}
void NTT(int *A,bool flag)
{for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(int mid=1;mid<lim;mid<<=1){ll w=quick_mul(G,(mod-1)/(mid<<1));for(int i=0;i<lim;i+=mid<<1){ll fw=1;for(int j=0;j<mid;j++,fw=fw*w%mod){ll x=A[i+j],y=fw*A[i+j+mid]%mod;A[i+j]=(x+y)%mod,A[i+j+mid]=(x-y+mod)%mod;}}}if(flag)return ;ll inv_=inv(lim);reverse(A+1,A+lim);for(int i=0;i<lim;i++)A[i]=inv_*A[i]%mod;
}
int F1[N];
void poly_inv(int *A,int *B,int len)
{if(len==1){B[0]=inv(A[0]);return ;}poly_inv(A,B,(len+1)>>1);init(len<<1);for(int i=0;i<len;i++)F1[i]=A[i];for(int i=len;i<lim;i++)F1[i]=0;NTT(F1,1);NTT(B,1);for(int i=0;i<lim;i++)B[i]=(2ll-1ll*F1[i]*B[i]%mod+mod)%mod*B[i]%mod;NTT(B,0);for(int i=len;i<lim;i++)B[i]=0;for(int i=0;i<lim;i++)F1[i]=0;
}
int n,m,a[N],b[N],ar[N],br[N],bt[N];
int main()
{scanf("%d%d",&n,&m);for(int i=0;i<=n;i++)scanf("%d",&a[i]),ar[n-i]=a[i];for(int i=0;i<=m;i++)scanf("%d",&b[i]),br[m-i]=b[i];for(int i=n-m+2;i<=m;i++)br[i]=0;poly_inv(br,bt,n-m+1);init((n*2-m+1));NTT(bt,1);NTT(ar,1);for(int i=0;i<lim;i++)bt[i]=1ll*bt[i]*ar[i]%mod;NTT(bt,0);reverse(bt,bt+n-m+1);for(int i=n-m+1;i<lim;i++)bt[i]=0;for(int i=0;i<=n-m;i++)printf("%d ",bt[i]);printf("n");init((n<<1));NTT(bt,1);NTT(b,1);for(int i=0;i<lim;i++)bt[i]=1ll*bt[i]*b[i]%mod;NTT(bt,0);for(int i=0;i<m;i++)printf("%d ",(a[i]-bt[i]+mod)%mod);
}
本文发布于:2024-02-02 10:26:27,感谢您对本站的认可!
本文链接:https://www.4u4v.net/it/170684078643187.html
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。
留言与评论(共有 0 条评论) |