多項式乘法
FFT
因為有浮點數參與運算，所以可能會出現精度的問題
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const double pi=acos(-1);
const int maxn=4e6+5;
struct fs{
double x,y;
friend fs operator + (rg const fs& A,rg const fs& B){
return (fs){A.x+B.x,A.y+B.y};
}
friend fs operator - (rg const fs& A,rg const fs& B){
return (fs){A.x-B.x,A.y-B.y};
}
friend fs operator * (rg const fs& A,rg const fs& B){
return (fs){A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
}
friend fs operator / (rg const fs& A,rg const double& B){
return (fs){A.x/B,A.y/B};
}
}a[maxn],b[maxn],w[25][maxn];
int wz[maxn],n,m;
void fft(rg fs A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg fs x=A[j+k],y=A[j+k+len]*w[t0][k];
A[j+k]=x+y;
A[j+k+len]=x-y;
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
for(rg int i=0;i<lim;i++) A[i]=A[i]/(double)lim;
}
}
int main(){
n=read(),m=read();
for(rg int i=0;i<=n;i++) a[i].x=read();
for(rg int i=0;i<=m;i++) b[i].x=read();
rg int lim=1,bit=0;
for(;lim<=n+m;lim<<=1) bit++;
for(rg int i=0;i<lim;i++){
wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
}
for(rg int i=0,len=1;len<lim;len<<=1,i++){
for(rg int j=0;j<len;j++){
w[i][j]=(fs){cos(pi*j/len),sin(pi*j/len)};
}
}
fft(a,lim,1),fft(b,lim,1);
for(rg int i=0;i<lim;i++) a[i]=a[i]*b[i];
fft(a,lim,-1);
for(rg int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x+0.5));
printf("\n");
return 0;
}
NTT
用原根代替複數
但必須保證模數 \(p = a \cdot 2^k + 1\)
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e6+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,m,w[23][maxn];
int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int main(){
n=read(),m=read();
for(rg int i=0;i<=n;i++) a[i]=read();
for(rg int i=0;i<=m;i++) b[i]=read();
rg int lim=1,bit=0;
for(;lim<=n+m;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0,k=1;k<lim;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
ntt(a,lim,1),ntt(b,lim,1);
for(rg int i=0;i<lim;i++) a[i]=mulmod(a[i],b[i]);
ntt(a,lim,-1);
for(rg int i=0;i<=n+m;i++) printf("%d ",a[i]);
printf("\n");
return 0;
}
\(FFT\) 和 \(NTT\) 預處理單位根要快不少
預處理單位根有兩種寫法，一種是 \(O(nlogn)\) 的，另一種是 \(O(n)\) 的
第一種因為訪問記憶體連續實測更快
任意模數多項式乘法（MTT）
當多項式乘法進行取模的數並不是我們喜歡的模數時，就要用到 \(MTT\)
實現的方法有很多，這裡介紹一種只需要 \(2\) 次 \(DFT\) 和 \(2\) 次 \(IDFT\) 的優秀做法
首先，我們對原來的多項式拆係數
令\(k=sqrt(mod)\)，則
\[\begin{aligned} A(x)=kA_0(x)+A_1(x)\\ B(x)=kB_0(x)+B_1(x)\end{aligned}
\]
即拆成了 \(A(x)/k\) 和 \(A(x)\%k\) 兩部分
那麼
\[\begin{aligned} A(x)B(x)&=(kA_0(x)+A_1(x))(kB_0(x)+B_1(x))\\ &=k^2A_0(x)B_0(x)+A_1(x)B_1(x)+kA_0(x)B_1(x)+kA_1(x)B_0(x) \end{aligned}
\]
這樣就有一個\(4\) 次 \(DFT\) 和 \(4\) 次 \(IDFT\) 的暴力做法
考慮優化，設
\[\begin{aligned} P(x)=A_0(x)+iA_1(x)\\ Q(x)=A_0(x)-iA_1(x)\end{aligned}
\]
由於 \(A, B\) 的虛部都為 \(0\)
\(P,Q\) 的每一項係數都互為共軛，同樣每一個點值也互為共軛
那麼只需對 \(P\) 做一次 \(DFT\) ，就可以通過共軛 \(O(n)\) 求出 \(Q\) 的點值表示
具體地說，在點值表示中 \(P\) 的第 \(k\) 項與 \(Q\) 的第 \(n-k\) 項是共軛
要特判下標為 \(0\) 的情況
求出了 \(P\) 和 \(Q\) ，我們就能用二元一次方程的知識推出 \(A_0\) 和 \(A_1\)
對於 \(B\) 也是同理
此時，我們進行了兩次 \(DFT\) 得到了所有多項式的點值表示
下面就是怎麼把它們變成係數表示
仍然借鑑上面的思路，構造兩個多項式
\(P(x)=A_0(x)B_0(x)+iA_1(x)B_0(x)\)
\(Q(x)=A_0(x)B_1(x)+iA_1(x)B_1(x)\)
通過已知的點值求出此時 \(P, Q\) 的點值
然後分別對 \(P, Q\) 做 \(IDFT\)
由於 \(A_0 B_0,A_0 B_1,A_1 B_0,A_1 B_1\)這四個多項式捲起來後的係數表示中虛部一定為 \(0\)
那麼此時 \(P\) 的實部和虛部就分別為 \(A_0(x) B_0(x)\)和 \(A_1(x) B_0(x)\)
同樣 \(Q\) 的實部和虛部就分別為 \(A_0(x) B_1(x)\) 和 \(A_1(x) B_1(x)\)。
這樣，我們又進行了兩次 \(IDFT\) 得到了所有多項式的係數表示
程式碼
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5;
const double pi=acos(-1);
struct fs{
double x,y;
friend fs operator + (rg const fs& A,rg const fs& B){
return (fs){A.x+B.x,A.y+B.y};
}
friend fs operator - (rg const fs& A,rg const fs& B){
return (fs){A.x-B.x,A.y-B.y};
}
friend fs operator * (rg const fs& A,rg const fs& B){
return (fs){A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
}
friend fs operator / (rg const fs& A,rg const double& B){
return (fs){A.x/B,A.y/B};
}
}a0[maxn],b0[maxn],a1[maxn],b1[maxn],w[25][maxn],p[maxn],q[maxn];
const fs I=(fs){0,1};
int wz[maxn],n,m,mod,blo;
void fft(rg fs A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg fs x=A[j+k],y=A[j+k+len]*w[t0][k];
A[j+k]=x+y;
A[j+k+len]=x-y;
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
for(rg int i=0;i<lim;i++) A[i]=A[i]/(double)lim;
}
}
void fftfft(rg fs A[],rg fs B[],rg int lim){
for(rg int i=0;i<lim;i++){
A[i]=A[i]+I*B[i];
}
fft(A,lim,1);
B[0]=(fs){A[0].x,-A[0].y};
for(rg int i=1;i<lim;i++) B[i]=(fs){A[lim-i].x,-A[lim-i].y};
rg fs t1,t2;
for(rg int i=0;i<lim;i++){
t1=A[i],t2=B[i];
A[i]=(t1+t2)/2.0;
B[i]=(t2-t1)*I/2.0;
}
}
long long get(rg double now){
if(now>=0) return (long long)(now+0.5)%mod;
else return (long long)(now-0.5)%mod;
}
int main(){
n=read(),m=read(),mod=read();
blo=sqrt(mod)+1;
rg int aa;
for(rg int i=0;i<=n;i++){
aa=read();
aa%=mod;
a0[i].x=aa/blo;
a1[i].x=aa%blo;
}
for(rg int i=0;i<=m;i++){
aa=read();
aa%=mod;
b0[i].x=aa/blo;
b1[i].x=aa%blo;
}
rg int bit=0,lim=1;
for(;lim<=n+m;lim<<=1) bit++;
for(rg int i=0;i<lim;i++){
wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
}
for(rg int i=0,len=1;len<lim;len<<=1,i++){
for(rg int j=0;j<len;j++){
w[i][j]=(fs){cos(pi*j/len),sin(pi*j/len)};
}
}
fftfft(a0,a1,lim),fftfft(b0,b1,lim);
for(rg int i=0;i<lim;i++){
p[i]=a0[i]*b0[i]+I*a1[i]*b0[i];
q[i]=a0[i]*b1[i]+I*a1[i]*b1[i];
}
fft(p,lim,-1),fft(q,lim,-1);
for(rg int i=0;i<=n+m;i++){
long long ans=1LL*blo*blo%mod*get(p[i].x)%mod+1LL*blo*get(q[i].x+p[i].y)%mod+get(q[i].y)%mod;
printf("%lld ",ans%mod);
}
return 0;
}
多項式乘法逆
如果 \(F(x)\) 只剩一項
顯然 \(G_0\) 就是 \(F_0\) 的逆元
如果有多項，可以遞迴求解
已知 \(F(x)H(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
又因為 \(F(x)G(x) \equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
那麼 \(F(x)(G(x)-H(x))\equiv 1 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
即 \(G(x)-H(x)\equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
兩邊同時平方
\(G(x)^2+H(x)^2-2G(x)H(x)\equiv 0 \pmod{x^n}\)
同時乘上 \(F(x)\)
\(G(x)+F(x)H(x)^2-2H(x)\equiv 0 \pmod{x^n}\)
所以
\(G(x)=2H(x)-F(x)H(x)^2(mod{x^n})\)
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getinv(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("\n");
return 0;
}
多項式開根
如果 \(F(x)\) 只剩一項
顯然如果 \(F_0=1\), \(G_0\) 就是 \(1\)
如果有多項，可以遞迴求解
已知 \(H(x)^2=F(x) \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
又因為 \(G(x)^2=F(x) \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
那麼 \(G(x)=H(x)\pmod{x^{\lceil \frac{n}{2} \rceil}}\)
即 \(G(x)-H(x)\equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
兩邊同時平方
\(G(x)^2+H(x)^2-2G(x)H(x)\equiv 0 \pmod{x^n}\)
\(F(x)+H(x)^2-2G(x)H(x)\equiv 0 \pmod{x^n}\)
所以
\(G(x)=\frac{F(x)+H(x)^2}{H(x)}(mod{x^n})\)
套一個多項式求逆即可
程式碼
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fsqrt1[maxn],fsqrt2[maxn];
void getsqrt(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=1;
return;
}
getsqrt(A,B,(deg+1)>>1);
getinv(B,fsqrt1,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) fsqrt2[i]=A[i];
for(rg int i=deg;i<lim;i++) fsqrt2[i]=0;
ntt(fsqrt1,lim,1),ntt(fsqrt2,lim,1);
for(rg int i=0;i<lim;i++) fsqrt1[i]=mulmod(fsqrt1[i],fsqrt2[i]);
ntt(fsqrt1,lim,-1);
rg int ny=ksm(2,mod-2);
for(rg int i=0;i<lim;i++) B[i]=mulmod(addmod(fsqrt1[i],B[i]),ny);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getsqrt(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("\n");
return 0;
}
多項式求導
\[F(x)=\sum^{n}_{i=0}a_ix^i
\]
\[F(x+dx)=\sum^{n}_{i=0}a_i(x+dx)^i
\]
\[F(x+dx)=\sum^{n}_{i=0}a_i(\binom{i}{0}x^i+\binom{i}{1}x^{i-1}dx+\binom{i}{2}x^{i-2}dx^2+...+\binom{i}{i-1}xdx^{i-1}+\binom{i}{i}dx^i)
\]
\[F(x+dx)-F(x)=\sum^{n}_{i=1}a_i(\binom{i}{1}x^{i-1}dx+\binom{i}{2}x^{i-2}dx^2+...+\binom{i}{i-1}xdx^{i-1}+\binom{i}{i}dx^i)
\]
\(F'(x)=\frac{F(x+dx)-F(x)}{dx}=\sum^{n}_{i=1}a_i\binom{i}{1}x^{i-1}\)
因為 \(dx\) 趨進於無窮小，所以後面的幾項都變成了 \(0\)
\(F'(x)=\sum^{n}_{i=0}a_{i + 1}(i+1)x^i\)
程式碼
void getdao(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i-1]=mulmod(A[i],i);
}
B[deg-1]=0;
}
多項式積分
求導的逆運算
程式碼
void getfen(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i]=mulmod(A[i-1],ksm(i,mod-2));
}
B[0]=0;
}
多項式對數函式（多項式 ln）
對兩邊同時求導
\(G(x)=F(A(x)),F(x)=ln(x)\)
\(G'(x)=F'(A(x))A'(x)=\frac{A'(x)}{A(x)}\)
所以對 \(A\) 求導求個逆在積一下分就算出了 \(G\)
程式碼
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fsqrt1[maxn],fsqrt2[maxn];
void getsqrt(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=1;
return;
}
getsqrt(A,B,(deg+1)>>1);
getinv(B,fsqrt1,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) fsqrt2[i]=A[i];
for(rg int i=deg;i<lim;i++) fsqrt2[i]=0;
ntt(fsqrt1,lim,1),ntt(fsqrt2,lim,1);
for(rg int i=0;i<lim;i++) fsqrt1[i]=mulmod(fsqrt1[i],fsqrt2[i]);
ntt(fsqrt1,lim,-1);
rg int ny=ksm(2,mod-2);
for(rg int i=0;i<lim;i++) B[i]=mulmod(addmod(fsqrt1[i],B[i]),ny);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i-1]=mulmod(A[i],i);
}
B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i]=mulmod(A[i-1],ksm(i,mod-2));
}
B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
memset(fln1,0,maxn<<2);
memset(fln2,0,maxn<<2);
getdao(A,fln1,deg);
getinv(A,fln2,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
ntt(fln1,lim,1),ntt(fln2,lim,1);
for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
ntt(fln1,lim,-1);
getfen(fln1,B,deg);
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getln(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("\n");
return 0;
}
多項式指數函式（多項式 exp）
\(F(x)\equiv F_0(x)(1-\ln F_0(x)+A(x))\pmod{x^n}\)
\(F_0(x)\) 指的是在模 \(x^{n\over 2}\) 意義下的答案
程式碼
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fsqrt1[maxn],fsqrt2[maxn];
void getsqrt(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=1;
return;
}
getsqrt(A,B,(deg+1)>>1);
getinv(B,fsqrt1,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) fsqrt2[i]=A[i];
for(rg int i=deg;i<lim;i++) fsqrt2[i]=0;
ntt(fsqrt1,lim,1),ntt(fsqrt2,lim,1);
for(rg int i=0;i<lim;i++) fsqrt1[i]=mulmod(fsqrt1[i],fsqrt2[i]);
ntt(fsqrt1,lim,-1);
rg int ny=ksm(2,mod-2);
for(rg int i=0;i<lim;i++) B[i]=mulmod(addmod(fsqrt1[i],B[i]),ny);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i-1]=mulmod(A[i],i);
}
B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i]=mulmod(A[i-1],ksm(i,mod-2));
}
B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
memset(fln1,0,maxn<<2);
memset(fln2,0,maxn<<2);
getdao(A,fln1,deg);
getinv(A,fln2,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
ntt(fln1,lim,1),ntt(fln2,lim,1);
for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
ntt(fln1,lim,-1);
getfen(fln1,B,deg);
}
int fexp1[maxn],fexp2[maxn];
void getexp(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=1;
return;
}
getexp(A,B,(deg+1)>>1);
getln(B,fexp1,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) fexp2[i]=A[i];
ntt(fexp1,lim,1),ntt(fexp2,lim,1),ntt(B,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],addmod(delmod(1,fexp1[i]),fexp2[i]));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getexp(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("\n");
return 0;
}