前置芝士:麦克拉伦展开,泰勒展开,牛顿迭代,快速傅里叶变换,标记($mod x^n$)。
正片:
1.多项式加法。
2.多项式减法。
3.多项式乘法。
4.多项式求逆。
5.多项式带余除法。
6.多项式开根。
7.多项式对数函数。
8.多项式指数函数。
9.多项式的k次幂/k次方根
10.多项式插值
正片开始:
1.多项式加法:(略)
2.多项式减法:(略)
3.多项式乘法:(FFT即可)
4.多项式求逆:
设原多项式为$H(x)$,所求为$F(x)(mod ; x^t)$
假设对于所求多项式$G(x)equiv F(x)(mod ; x^n)$的$G(x)$已经得到,那么:
$G(x)H(x)equiv 1(mod ; x^n)$
$G(x)H(x)-1equiv 0(mod ; x^n)$
$[G(x)H(x)-1]^2equiv 0(mod ; x^{2n})$
$G(x)^2H(x)^2-2G(x)H(x)+1equiv 0(mod ; x^{2n})$
$G(x)H(x)[2-G(x)H(x)]equiv 1(mod ; x^{2n})$
$frac{1}{H(x)}equiv G(x)[2-G(x)H(x)](mod ; x^{2n})$
$F(x)equiv G(x)[2-G(x)H(x)](mod ; x^{2n})$
这样就可以倍增地得到F(x),边界条件为常数项为原多项式常数项的逆元。
代码:
1 #include<cstdio>
2 #include<cstring>
3 #include<algorithm>
4 typedef long long lnt;
5 const lnt mod=998244353;
6 lnt A[1000000];
7 lnt B[1000000];
8 lnt C[1000000];
9 int pos[1000000];
10 int len;
11 int n;
12 lnt ksm(lnt a,lnt b)
13 {
14 lnt ans=1;
15 while(b)
16 {
17 if(b&1)ans=ans*a%mod;
18 a=a*a%mod;
19 b=b/2;
20 }
21 return ans;
22 }
23 void NTT(lnt *a,int flag)
24 {
25 for(int i=0;i<len;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]);
26 for(int i=2;i<=len;i<<=1)
27 {
28 lnt wn=ksm(3,(mod-1)/i);
29 for(int j=0;j<len;j+=i)
30 {
31 lnt w=1,t;
32 for(int k=0;k<(i>>1);k++,w=w*wn%mod)
33 {
34 t=a[j+k+(i>>1)]*w%mod;
35 a[j+k+(i>>1)]=(a[j+k]-t)%mod;
36 a[j+k]=(a[j+k]+t)%mod;
37 }
38 }
39 }
40 if(flag==-1)std::reverse(a+1,a+len);
41 return ;
42 }
43 void INV(lnt *a,lnt *b,int lim)
44 {
45 if(lim==1)
46 {
47 b[0]=ksm(a[0],mod-2);
48 return ;
49 }
50 INV(a,b,lim-1);
51 len=1<<lim;
52 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1));
53 memcpy(C,a,sizeof(lnt)*(len>>1));
54 NTT(C,1),NTT(b,1);
55 for(int i=0;i<len;i++)C[i]=b[i]*(2-C[i]*b[i]%mod+mod)%mod;
56 NTT(C,-1);
57 lnt inv=ksm(len,mod-2);
58 for(int i=0;i<len>>1;i++)b[i]=C[i]*inv%mod;
59 memset(b+(len>>1),0,sizeof(lnt)*(len>>1));
60 return ;
61 }
62 int main()
63 {
64 scanf("%d",&n);
65 int l=0;
66 while((1<<l)<(n<<1))l++;
67 for(int i=0;i<n;i++)scanf("%lld",&A[i]);
68 INV(A,B,l);
69 for(int i=0;i<n;i++)printf("%lld ",(B[i]%mod+mod)%mod);
70 puts("");
71 return 0;
72 }
5.多项式带余除法:
余项不容易搞定,那么就想办法将余项消掉。
构造一种巧妙的方法。
定义$A^R(x)$为$A(x)$的系数反转的多项式。
例如:$A(x)=x+2$,则$A^R(x)=2x+1$,
这样就可以将$R(x)$的系数转上去被余掉,考虑如下方法实现:
原多项式关系为$A(x)=B(x)C(x)+R(x)$,$A(x)$的最高次为$n$,$B(x)为$m$,$C(x)$为$n-m$,$R(x)$为$m-1$。
$A(x)=x^nA^R(frac{1} {x})$
$x^nA^R(frac{1} {x})=x^mB^R(frac{1}{x})x^{n-m}C^R(frac{1}{x})+x^{m-1}R^R(frac{1}{x})x^{n-m+1}$
这时发现$C(x)$的次数为$n-m$所以整个式子在$mod ;x^{n-m+1}$下求出的$C(x)$是无损的。
而且正好将$R(x)$忽略了。就只需要多项式求逆解出$C^R(x)$反转系数得到$C(x)$,最后回代得到$R(x)$就好了。
代码:
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 typedef long long lnt; 5 const lnt mod=998244353; 6 lnt a[1000000]; 7 lnt b[1000000]; 8 lnt d[1000000]; 9 lnt r[1000000]; 10 lnt ar[1000000]; 11 lnt br[1000000]; 12 lnt ibr[1000000]; 13 lnt A[1000000]; 14 lnt B[1000000]; 15 int n,m; 16 int len; 17 int pos[1000000]; 18 lnt ksm(lnt a,lnt b) 19 { 20 lnt ans=1; 21 while(b) 22 { 23 if(b&1)ans=ans*a%mod; 24 a=a*a%mod; 25 b=b/2; 26 } 27 return ans; 28 } 29 void NTT(lnt *a,int flag) 30 { 31 for(int i=0;i<len;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]); 32 for(int i=2;i<=len;i<<=1) 33 { 34 lnt wn=ksm(3,(mod-1)/i); 35 for(int j=0;j<len;j+=i) 36 { 37 lnt w=1,t; 38 for(int k=0;k<(i>>1);k++,w=w*wn%mod) 39 { 40 t=a[j+k+(i>>1)]*w%mod; 41 a[j+k+(i>>1)]=(a[j+k]-t)%mod; 42 a[j+k]=(a[j+k]+t)%mod; 43 } 44 } 45 } 46 if(flag==-1)std::reverse(a+1,a+len); 47 return ; 48 } 49 void Get_inv(lnt *a_,lnt *b_,int lim) 50 { 51 if(lim==1) 52 { 53 b_[0]=ksm(a_[0],mod-2); 54 return ; 55 } 56 Get_inv(a_,b_,lim-1); 57 len=1<<lim; 58 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 59 memcpy(A,a_,sizeof(lnt)*(len>>1)); 60 NTT(A,1),NTT(b_,1); 61 for(int i=0;i<len;i++)A[i]=b_[i]*(2ll-A[i]*b_[i]%mod+mod)%mod; 62 NTT(A,-1); 63 lnt INV=ksm(len,mod-2); 64 for(int i=0;i<len>>1;i++)b_[i]=A[i]*INV%mod; 65 memset(b_+(len>>1),0,sizeof(lnt)*(len>>1)); 66 return ; 67 } 68 int main() 69 { 70 scanf("%d%d",&n,&m); 71 for(int i=0;i<=n;i++)scanf("%lld",&a[i]); 72 for(int i=0;i<=m;i++)scanf("%lld",&b[i]); 73 std::reverse(a,a+n+1);memcpy(ar,a,sizeof(a)); 74 std::reverse(b,b+m+1);memcpy(br,b,sizeof(b)); 75 std::reverse(a,a+n+1);std::reverse(b,b+m+1); 76 memset(br+n-m+1,0,sizeof(lnt)*m); 77 memset(ar+n-m+1,0,sizeof(lnt)*m); 78 int lim=0; 79 while((1<<lim)<((n-m+1)<<1))lim++; 80 Get_inv(br,ibr,lim); 81 len=1<<lim; 82 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 83 NTT(ar,1);NTT(ibr,1); 84 for(int i=0;i<len;i++)d[i]=ar[i]*ibr[i]%mod; 85 NTT(d,-1); 86 lnt INV=ksm(len,mod-2); 87 for(int i=0;i<len;i++)d[i]=d[i]*INV%mod; 88 std::reverse(d,d+n-m+1); 89 for(int i=0;i<=n-m;i++)printf("%lld ",(d[i]%mod+mod)%mod); 90 puts(""); 91 memset(d+n-m+1,0,sizeof(lnt)*n); 92 lim=0; 93 while((1<<lim)<(n<<1))lim++; 94 len=1<<lim; 95 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 96 NTT(d,1),NTT(a,1),NTT(b,1); 97 for(int i=0;i<len;i++)r[i]=(a[i]-b[i]*d[i]%mod+mod)%mod; 98 NTT(r,-1);INV=ksm(len,mod-2); 99 for(int i=0;i<len;i++)r[i]=r[i]*INV%mod; 100 for(int i=0;i<m;i++)printf("%lld ",(r[i]%mod+mod)%mod); 101 puts(""); 102 return 0; 103 }
6.多项式开根
这里先不介绍那个指对数混合的方法。
继续推式子:
假设已经知道了$F_0(x)^2equiv G(x)(mod;x^n)$
那么:$[F_0(x)^2-G(x)]^2equiv0(mod;x^{2n})$
$[F_0(x)^2+G(x)]^2equiv 4G(x)F_0(x)^2(mod;x^{2n})$
$[frac{F_0(x)^2+G(x)}{2F_0(x)}]^2equiv G(x)(mod;x^{2n})$
所以:$frac{F_0(x)^2+G(x)}{2F_0(x)}equiv F(x)(mod;x^{2n})$
代码:
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 typedef long long lnt; 5 const lnt mod=998244353; 6 lnt A[1000000]; 7 lnt B[1000000]; 8 lnt C[1000000]; 9 lnt f[1000000]; 10 lnt g[1000000]; 11 int n; 12 int len; 13 int pos[1000000]; 14 lnt ksm(lnt a,lnt b) 15 { 16 lnt ans=1; 17 while(b) 18 { 19 if(b&1)ans=ans*a%mod; 20 a=a*a%mod; 21 b=b/2; 22 } 23 return ans; 24 } 25 void NTT(lnt *a,int flag) 26 { 27 for(int i=0;i<len;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]); 28 for(int i=2;i<=len;i<<=1) 29 { 30 lnt wn=ksm(3,(mod-1)/i); 31 for(int j=0;j<len;j+=i) 32 { 33 lnt w=1,t; 34 for(int k=0;k<(i>>1);k++,w=w*wn%mod) 35 { 36 t=a[j+k+(i>>1)]*w%mod; 37 a[j+k+(i>>1)]=(a[j+k]-t)%mod; 38 a[j+k]=(a[j+k]+t)%mod; 39 } 40 } 41 } 42 if(flag==-1)std::reverse(a+1,a+len); 43 return ; 44 } 45 void Get_inv(lnt *a,lnt *b,int lim) 46 { 47 if(lim==1) 48 { 49 b[0]=1; 50 return ; 51 } 52 Get_inv(a,b,lim-1); 53 len=1<<lim; 54 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 55 memcpy(A,a,sizeof(lnt)*(len>>1)); 56 NTT(A,1),NTT(b,1); 57 for(int i=0;i<len;i++)A[i]=b[i]*(2ll-A[i]*b[i]%mod+mod)%mod; 58 NTT(A,-1); 59 lnt INV=ksm(len,mod-2); 60 for(int i=0;i<(len>>1);i++)b[i]=A[i]*INV%mod; 61 memset(b+(len>>1),0,sizeof(lnt)*(len>>1)); 62 return ; 63 } 64 void Get_sqrt(lnt *a,lnt *b,int lim) 65 { 66 if(lim==1) 67 { 68 b[0]=1; 69 return ; 70 } 71 Get_sqrt(a,b,lim-1); 72 memset(C,0,sizeof(lnt)*(1<<lim)); 73 memset(B,0,sizeof(lnt)*(1<<lim)); 74 memset(A,0,sizeof(lnt)*(1<<lim)); 75 Get_inv(b,C,lim); 76 len=1<<lim; 77 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 78 memcpy(B,a,sizeof(lnt)*(len>>1)); 79 NTT(C,1);NTT(B,1);NTT(b,1); 80 lnt inv2=ksm(2,mod-2); 81 for(int i=0;i<len;i++)B[i]=inv2*((b[i]+B[i]*C[i]%mod)%mod)%mod; 82 NTT(B,-1); 83 lnt INV=ksm(len,mod-2); 84 for(int i=0;i<(len>>1);i++)b[i]=B[i]*INV%mod; 85 memset(b+(len>>1),0,sizeof(lnt)*(len>>1)); 86 return ; 87 } 88 int main() 89 { 90 scanf("%d",&n); 91 for(int i=0;i<n;i++)scanf("%lld",&f[i]),f[i]%=mod; 92 int lim=0; 93 while((1<<lim)<(n<<1))lim++; 94 Get_sqrt(f,g,lim); 95 for(int i=0;i<n;i++)printf("%lld ",(g[i]%mod+mod)%mod); 96 puts(""); 97 return 0; 98 }
7.多项式对数函数(ln)
考虑$f(x)=ln(g(x))$
则$f'(x)=frac{1}{g(x)}g'(x)$
多项式求逆即可。
代码:
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 typedef long long lnt; 5 const lnt mod=998244353; 6 lnt A[1000010]; 7 lnt B[1000010]; 8 lnt C[1000010]; 9 lnt f[1000010]; 10 int pos[1000000]; 11 int len; 12 int n; 13 lnt ksm(lnt a,lnt b) 14 { 15 lnt ans=1; 16 while(b) 17 { 18 if(b&1)ans=ans*a%mod; 19 a=a*a%mod; 20 b=b/2; 21 } 22 return ans; 23 } 24 void NTT(lnt *a,int flag) 25 { 26 for(int i=0;i<len;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]); 27 for(int i=2;i<=len;i<<=1) 28 { 29 lnt wn=ksm(3,(mod-1)/i); 30 for(int j=0;j<len;j+=i) 31 { 32 lnt w=1,t; 33 for(int k=0;k<(i>>1);k++,w=w*wn%mod) 34 { 35 t=a[j+k+(i>>1)]*w%mod; 36 a[j+k+(i>>1)]=(a[j+k]-t)%mod; 37 a[j+k]=(a[j+k]+t)%mod; 38 } 39 } 40 } 41 if(flag==-1)std::reverse(a+1,a+len); 42 return ; 43 } 44 void Get_inv(lnt *a,lnt *b,int lim) 45 { 46 if(lim==1) 47 { 48 b[0]=ksm(a[0],mod-2); 49 return ; 50 } 51 Get_inv(a,b,lim-1); 52 len=1<<lim; 53 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 54 memcpy(A,a,sizeof(lnt)*(len>>1)); 55 NTT(A,1),NTT(b,1); 56 for(int i=0;i<len;i++)A[i]=b[i]*(2ll-A[i]*b[i]%mod+mod)%mod; 57 NTT(A,-1); 58 lnt INV=ksm(len,mod-2); 59 for(int i=0;i<(len>>1);i++)b[i]=A[i]*INV%mod; 60 memset(b+(len>>1),0,sizeof(lnt)*(len>>1)); 61 return ; 62 } 63 void Derivation(lnt *a,lnt *b) 64 { 65 b[len-1]=0; 66 for(int i=0;i<len-1;i++)b[i]=a[i+1]*(i+1)%mod; 67 return ; 68 } 69 void Integral(lnt *a,lnt *b) 70 { 71 b[0]=0; 72 for(int i=0;i<len-1;i++)b[i+1]=a[i]*ksm(i+1,mod-2)%mod; 73 return ; 74 } 75 void ln(lnt *a) 76 { 77 int lim=0; 78 while((1<<lim)<(n<<1))lim++; 79 len=1<<lim; 80 Derivation(a,C); 81 Get_inv(a,B,lim); 82 len=1<<lim; 83 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 84 NTT(B,1),NTT(C,1); 85 for(int i=0;i<len;i++)C[i]=C[i]*B[i]%mod; 86 NTT(C,-1); 87 lnt INV=ksm(len,mod-2); 88 for(int i=0;i<len;i++)C[i]=C[i]*INV%mod; 89 Integral(C,a); 90 return ; 91 } 92 int main() 93 { 94 scanf("%d",&n); 95 for(int i=0;i<n;i++)scanf("%lld",&f[i]); 96 ln(f); 97 for(int i=0;i<n;i++)printf("%lld ",(f[i]%mod+mod)%mod); 98 puts(""); 99 return 0; 100 }
8.多项式指数函数(exp)
这个就需要牛顿迭代了。
假设我们原多项式为$H(x)$,所求多项式为$F(x)$
假设我们有一个函数$G(F(x))equiv 0 (mod; x^n)$
我们已经得到$G(F_0(x))equiv 0(mod;x^n)$
那么写成泰勒展开,就得到了:
$G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))equiv 0 (mod; x^{2n})$
这里$G(x)=lnF(x)-H(x)$
代入上式得到:
$F(x)equiv F_0(x)(1-ln(F_0(x))+H(x))(mod;x^{2n})$
代码:
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 typedef long long lnt; 5 const lnt mod=998244353; 6 lnt A[1000000]; 7 lnt B[1000000]; 8 lnt C[1000000]; 9 lnt D[1000000]; 10 lnt E[1000000]; 11 lnt f[1000000]; 12 lnt g[1000000]; 13 int n; 14 int len; 15 int pos[1000000]; 16 lnt ksm(lnt a,lnt b) 17 { 18 lnt ans=1; 19 while(b) 20 { 21 if(b&1)ans=ans*a%mod; 22 a=a*a%mod; 23 b=b/2; 24 } 25 return ans; 26 } 27 void NTT(lnt *a,int flag) 28 { 29 for(int i=0;i<len;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]); 30 for(int i=2;i<=len;i<<=1) 31 { 32 lnt wn=ksm(3,(mod-1)/i); 33 for(int j=0;j<len;j+=i) 34 { 35 lnt w=1,t; 36 for(int k=0;k<(i>>1);k++,w=w*wn%mod) 37 { 38 t=a[j+k+(i>>1)]*w%mod; 39 a[j+k+(i>>1)]=(a[j+k]-t)%mod; 40 a[j+k]=(a[j+k]+t)%mod; 41 } 42 } 43 } 44 if(flag==-1)std::reverse(a+1,a+len); 45 return ; 46 } 47 void Get_inv(lnt *a,lnt *b,int lim) 48 { 49 if(lim==1) 50 { 51 b[0]=ksm(a[0],mod-2); 52 return ; 53 } 54 Get_inv(a,b,lim-1); 55 len=1<<lim; 56 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 57 memcpy(A,a,sizeof(lnt)*(len>>1)); 58 NTT(A,1),NTT(b,1); 59 for(int i=0;i<len;i++)A[i]=b[i]*(2ll-A[i]*b[i]%mod+mod)%mod; 60 NTT(A,-1); 61 lnt INV=ksm(len,mod-2); 62 for(int i=0;i<len>>1;i++)b[i]=A[i]*INV%mod; 63 memset(b+(len>>1),0,sizeof(lnt)*(len>>1)); 64 return ; 65 } 66 void derivation(lnt *a,lnt *b) 67 { 68 b[len-1]=0; 69 for(int i=1;i<len;i++)b[i-1]=a[i]*i%mod; 70 return ; 71 } 72 void integral(lnt *a,lnt *b) 73 { 74 b[0]=0; 75 for(int i=1;i<len;i++)b[i]=a[i-1]*ksm(i,mod-2)%mod; 76 return ; 77 } 78 void Get_ln(lnt *a,lnt *b,int lim) 79 { 80 len=1<<lim; 81 derivation(a,B); 82 Get_inv(a,C,lim); 83 len=1<<lim; 84 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 85 NTT(B,1); 86 NTT(C,1); 87 for(int i=0;i<len;i++)B[i]=B[i]*C[i]%mod; 88 NTT(B,-1); 89 lnt INV=ksm(len,mod-2); 90 for(int i=0;i<len;i++)B[i]=B[i]*INV%mod; 91 integral(B,b); 92 return ; 93 } 94 void Get_exp(lnt *a,lnt *b,int lim) 95 { 96 if(lim==1) 97 { 98 b[0]=1; 99 return ; 100 } 101 Get_exp(a,b,lim-1); 102 len=1<<lim; 103 memset(A,0,sizeof(lnt)*len); 104 memset(B,0,sizeof(lnt)*len); 105 memset(C,0,sizeof(lnt)*len); 106 memset(D,0,sizeof(lnt)*len); 107 memset(E,0,sizeof(lnt)*len); 108 Get_ln(b,D,lim); 109 len=1<<lim; 110 for(int i=0;i<len;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(lim-1)); 111 memcpy(E,a,sizeof(lnt)*(len>>1)); 112 NTT(b,1);NTT(D,1);NTT(E,1); 113 for(int i=0;i<len;i++)D[i]=b[i]*(1ll-D[i]+E[i])%mod; 114 NTT(D,-1); 115 lnt INV=ksm(len,mod-2); 116 for(int i=0;i<len>>1;i++)b[i]=D[i]*INV%mod; 117 memset(b+(len>>1),0,sizeof(lnt)*(len>>1)); 118 return ; 119 } 120 int main() 121 { 122 scanf("%d",&n); 123 for(int i=0;i<n;i++)scanf("%lld",&f[i]); 124 int lim=0; 125 while((1<<lim)<(n<<1))lim++; 126 Get_exp(f,g,lim); 127 for(int i=0;i<n;i++)printf("%lld ",(g[i]%mod+mod)%mod); 128 return 0; 129 }