多项式Ⅰ
多项式Ⅰ
update20181229
这篇屑文章咕咕了，我更改了一些代码写法（非递归->递归），从牛顿迭代的思想理解了一下，这里的一些思想可能比较繁琐，可以看看多项式Ⅱ
前言：本人只会瞎背板子，根本没怎么理解，而且只学了一两个最简单的（学一点更新一点
前置定义：多项式的度表示多项式的最高次数。
多项式求逆
-
定义
\[A(x)A^{-1}(x)\equiv 1\pmod {x^n} \]\(n\)是\(A(x)\)的度，这里是两个多项式进行卷积之后对一个多项式(\(x^n\))取模
-
为什么要取模？
如果不取模，除了只有常数项的多项式以外，其他多项式的逆元都是无穷级数。
-
一个取模的感性理解
对\(x^n\)取模相当于把次数大于\(n\)的项都扔掉
-
-
倍增
若求出了\(A(x)A'(x)\equiv 1\pmod {x^{\lceil\frac{n}{2}\rceil}}\)的\(A'(x)\)
\[A(x)A^{-1}(x)\equiv 1\pmod {x^n} \]\[A(x)A^{-1}(x)\equiv 1 \pmod {x^{\lceil\frac{n}{2}\rceil}} \]\[A(x)(A^{-1}(x)-A'(x))\equiv 0 \pmod {x^{\lceil\frac{n}{2}\rceil}} \]\[A^{-1}(x)-A'(x)\equiv 0 \pmod {x^{\lceil\frac{n}{2}\rceil}} \]\[(A^{-1}(x)-A'(x))^2\equiv 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^2\equiv F \pmod {2^n} \]\[T'^2\equiv F\pmod {2^{\lceil\frac{n}{2}\rceil}} \]\[(T+T')(T-T')\equiv 0 \pmod {2^{\lceil\frac{n}{2}\rceil}} \]\[(T-T')\equiv 0\pmod {2^{\lceil\frac{n}{2}\rceil}} \]\[T^2-2TT'+T'^2\equiv 0\pmod {2^n} \]\[F-2TT'+T'^2\equiv 0\pmod {2^n} \]\[T\equiv \frac{F+T'^2}{2T'}\pmod {2^n} \]或者再化简一步
\[T\equiv \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]; }
一定还会更新的...
