多项式Ⅰ
update20181229
这篇屑文章咕咕了，我更改了一些代码写法（非递归->递归），从牛顿迭代的思想理解了一下，这里的一些思想可能比较繁琐，可以看看多项式Ⅱ
前言：本人只会瞎背板子，根本没怎么理解，而且只学了一两个最简单的（学一点更新一点
前置定义：多项式的度表示多项式的最高次数。
多项式求逆
-
定义
[A(x)A^{-1}(x)equiv 1pmod {x^n} ](n)是(A(x))的度，这里是两个多项式进行卷积之后对一个多项式((x^n))取模
-
为什么要取模？
如果不取模，除了只有常数项的多项式以外，其他多项式的逆元都是无穷级数。
-
一个取模的感性理解
对(x^n)取模相当于把次数大于(n)的项都扔掉
-
-
倍增
若求出了(A(x)A'(x)equiv 1pmod {x^{lceilfrac{n}{2} ceil}})的(A'(x))
[A(x)A^{-1}(x)equiv 1pmod {x^n} ][A(x)A^{-1}(x)equiv 1 pmod {x^{lceilfrac{n}{2} ceil}} ][A(x)(A^{-1}(x)-A'(x))equiv 0 pmod {x^{lceilfrac{n}{2} ceil}} ][A^{-1}(x)-A'(x)equiv 0 pmod {x^{lceilfrac{n}{2} ceil}} ][(A^{-1}(x)-A'(x))^2equiv 0 pmod {x^n} ][(A^{-1}(x))^2-2A^{-1}(x)A'(x)+A'^2(x)equiv 0 pmod {x^n} ][A^{-1}(x)-2A'(x)+A(x)A'^2(x)equiv 0 pmod {x^n} ][A^{-1}(x)equiv 2A'(x)-A(x)A'^2(x)pmod {x^n} ]上面过程以已意会为主，谢谢。正经的觉得疑惑还是从定义去理解+讨论。
-
细节
注意每次乘法时取长度取多不取少，比如那个三个多项式乘开开四倍长度。
-
Code:
#include <cstdio> #include <algorithm> #define ll long long const int N=(1<<18)+10; const int mod=998244353,G=3,Gi=332748118; #define add(x,y) ((x+y)%mod) #define Mul(x,y) (1ll*(x)*(y)%mod) int qp(int d,int k){int f=1;while(k){if(k&1)f=Mul(f,d);d=Mul(d,d);k>>=1;}return f;} int A[N],B[N],f[N],b[2][N],turn[N]; void NTT(int *a,int len,int typ) { for(int i=0;i<len;i++) if(i<turn[i]) std::swap(a[i],a[turn[i]]); for(int le=1;le<len;le<<=1) { int wn=qp(typ?G:Gi,(mod-1)/(le<<1)); for(int p=0;p<len;p+=le<<1) { int w=1; for(int i=p;i<p+le;i++,w=Mul(w,wn)) { int tx=a[i],ty=Mul(w,a[i+le]); a[i]=add(tx,ty); a[i+le]=add(tx,mod-ty); } } } if(!typ) { int inv=qp(len,mod-2); for(int i=0;i<len;i++) a[i]=Mul(a[i],inv); } } void mul(int *a,int *b,int len) { int L=-1;for(int i=1;i<len;i<<=1,++L); for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|(i&1)<<L,A[i]=B[i]=0; for(int i=0;i<len>>1;i++) A[i]=a[i],B[i]=b[i]; NTT(A,len,1),NTT(B,len,1); for(int i=0;i<len;i++) A[i]=Mul(A[i],B[i]); NTT(A,len,0); for(int i=0;i<len;i++) a[i]=A[i]; } void inv(int *a,int n) { int cur=0,len=2; b[cur][0]=qp(f[0],mod-2); while(len<=(n<<2)) { cur^=1; for(int i=0;i<len>>1;i++) b[cur][i]=add(b[cur^1][i],b[cur^1][i]); mul(b[cur^1],b[cur^1],len); mul(b[cur^1],f,len); for(int i=0;i<len;i++) b[cur][i]=add(b[cur][i],mod-b[cur^1][i]); len<<=1; } for(int i=0;i<n;i++) a[i]=b[cur][i]; } int main() { int n;scanf("%d",&n); for(int i=0;i<n;i++) scanf("%d",f+i); inv(f,n); for(int i=0;i<n;i++) printf("%d ",f[i]); return 0; }
多项式除法
[A(x)=B(x)D(x)+R(x)
]
注：(A(x))被除多项式，度为(n)，(B(x))除他的多项式，度为(m)，(D(x))商，度为(n-m)，余数为(R(x))，度不大于(m-1)
-
一种转换
(A^r(x)=x^nA(frac{1}{x}))
感性理解：把它的系数颠倒了，可以(O(n))处理
-
有什么用？
通过取模把余数扔掉。
-
-
[A(x)=B(x)D(x)+R(x) ][x^nA(x)=x^nB(x)D(x)+x^nR(x) ][A^r(x)=B^r(x)D^r(x)+x^{n-m+1}R^r(x) ][A^r(x)equiv B^r(x)D^r(x) pmod {x^{n-m+1}} ]
然后求个逆把(D^r(x))弄出来，带回去就行了
-
Code:
#include <cstdio> #include <algorithm> const int N=(1<<20)+10; const int mod=998244353,Ge=3,Gi=332748118; #define add(x,y) ((x+y)%mod) #define Mul(x,y) (1ll*(x)*(y)%mod) int qp(int d,int k){int f=1;while(k){if(k&1)f=Mul(f,d);d=Mul(d,d);k>>=1;}return f;} int F[N],G[N],b[2][N],X[N],Y[N],A[N],B[N],turn[N]; void NTT(int *a,int len,int typ) { for(int i=0;i<len;i++) if(i<turn[i]) std::swap(a[i],a[turn[i]]); for(int le=1;le<len;le<<=1) { int wn=qp(typ?Ge:Gi,(mod-1)/(le<<1)); for(int p=0;p<len;p+=le<<1) { int w=1; for(int i=p;i<p+le;i++,w=Mul(w,wn)) { int tx=a[i],ty=Mul(w,a[i+le]); a[i]=add(tx,ty); a[i+le]=add(tx,mod-ty); } } } if(!typ) { int inv=qp(len,mod-2); for(int i=0;i<len;i++) a[i]=Mul(a[i],inv); } } void polymul(int *a,int *b,int len) { int L=-1;for(int i=1;i<len;i<<=1,++L); for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|(i&1)<<L,A[i]=B[i]=0; for(int i=0;i<len>>1;i++) A[i]=a[i],B[i]=b[i]; NTT(A,len,1),NTT(B,len,1); for(int i=0;i<len;i++) A[i]=Mul(A[i],B[i]); NTT(A,len,0); for(int i=0;i<len;i++) a[i]=A[i]; } void polyinv(int *a,int n) { int cur=0,len=2; b[cur][0]=qp(a[0],mod-2); while(len<=(n<<2)) { cur^=1; for(int i=0;i<len>>1;i++) b[cur][i]=add(b[cur^1][i],b[cur^1][i]); polymul(b[cur^1],b[cur^1],len); polymul(b[cur^1],a,len); for(int i=0;i<len;i++) b[cur][i]=add(b[cur][i],mod-b[cur^1][i]); len<<=1; } for(int i=0;i<n;i++) a[i]=b[cur][i]; } void Reverse(int *a,int len) { for(int i=0;i<<1<len;i++) std::swap(a[i],a[len-i-1]); } void polydiv(int *a,int *b,int n,int m) { for(int i=0;i<n;i++) X[i]=a[i]; for(int i=0;i<m;i++) Y[i]=b[i]; Reverse(a,n),Reverse(b,m); polyinv(b,n-m+1); int len=1;while(len<=n<<2) len<<=1; polymul(a,b,len); Reverse(a,n-m+1); for(int i=n-m+1;i<len;i++) a[i]=0; polymul(Y,a,len); for(int i=0;i<m-1;i++) b[i]=add(X[i],mod-Y[i]); } int main() { int n,m;scanf("%d%d",&n,&m);++n,++m; for(int i=0;i<n;i++) scanf("%d",F+i); for(int i=0;i<m;i++) scanf("%d",G+i); polydiv(F,G,n,m); for(int i=0;i<n-m+1;i++) printf("%d ",F[i]);puts(""); for(int i=0;i<m-1;i++) printf("%d ",G[i]); return 0; }
多项式开根
-
像求逆那样进行倍增
[T^2equiv F pmod {2^n} ][T'^2equiv Fpmod {2^{lceilfrac{n}{2} ceil}} ][(T+T')(T-T')equiv 0 pmod {2^{lceilfrac{n}{2} ceil}} ][(T-T')equiv 0pmod {2^{lceilfrac{n}{2} ceil}} ][T^2-2TT'+T'^2equiv 0pmod {2^n} ][F-2TT'+T'^2equiv 0pmod {2^n} ][Tequiv frac{F+T'^2}{2T'}pmod {2^n} ]或者再化简一步
[Tequiv frac{F}{2T'}+frac{T'}{2}pmod {2^n} ] -
Code:
void polysqrt(int *a,int n) { int len=2,cur=0; sq[cur][0]=1; while(len<=(n<<2)) { cur^=1; for(int i=0;i<len>>1;i++) sq[cur][i]=mul(sq[cur^1][i],i2); for(int i=0;i<len>>1;i++) sq[cur^1][i]=add(sq[cur^1][i],sq[cur^1][i]); polyinv(sq[cur^1],len); polymul(sq[cur^1],a,len); for(int i=0;i<len;i++) sq[cur][i]=add(sq[cur][i],sq[cur^1][i]); len<<=1; } for(int i=0;i<n;i++) a[i]=sq[cur][i]; }
一定还会更新的...