[Notes][多項式]雑記・多項式アルゴリズム-多項式求逆多項式取模多項式開根...

19030 ワード
ノーツ 多項式
多項式
いくつかの単項式の加算からなる代数式を多項式形と呼ぶ.例えば、f(x)=Σni=0 aixi f(x)=Σi=0 n a i x i、deg f(x)deg
せいせいかんすう
Σ∞i=0 aixiΣi=0∞a i x i生成関数は親関数とも呼ばれ、多項式アルゴリズムと関連して移行を最適化する役割を果たすことが多い.
たこうしきアルゴリズム
加減法
多項式加減算法はf(x)=Σni=0 aixi f(x)=Σi=0 na i x i,g(x)=Σni=0 bixi g(x)=Σi=0 nb i x i
(f±g)(x)=Σni=0(ai±bi)xi(f±g)(x)=Σi=0 n(a i±b i)x i
コード実装も比較的簡単です
乗算＃ジョウサン＃
多項式乗算はすべての多項式アルゴリズムの基礎であり、関数を生成するボリューム演算でもある.
f(x)=Σni=0 aixif(x)=Σi=0 na i x i,g(x)=Σni=0 bixi g(x)=Σi=0 nb i x i
則(f×g)(x)=∑2ni=0∑ij=0aj×bi−jxi ( f × g ) ( x ) = ∑ i = 0 2 n ∑ j = 0 i a j × b i − j x i
素朴な多項式乗算はO(n 2)O(n 2)であるが,FFT(高速フーリエ変換)は複数本の特性を用いてO(nlogn)O(n log小さな範囲で暴力を振るうことが多い.
inline void FFT(E *a,int r){
  for(int i=0;iif(rev[i]>i) swap(a[rev[i]],a[i]);
  for(int i=1;i1){
    E wn(cos(M_PI/i),r*sin(M_PI/i));
    for(int j=0;j1)){
      E w(1,0);
      for(int k=0;k*wn){
    E x=a[j+k],y=w*a[j+k+i];
    a[j+k]=x+y; a[j+k+i]=x-y;
      }
    }
  }
  if(r==-1) for(int i=0;i

生成関数を用いたカウント問題は往々にして1つの質量数をモデリングし,この質量数をkと表すことができる×2n+1 k × 2 n+1であり,kが奇数nで多項式の度数より大きい場合,この素数の原根を複数本の代わりに用いて多項式の係数対1の素数モデリングを実現することができ,このアルゴリズムがNTT(高速数論変換)であり,このモジュールをNTTモジュールと呼ぶ.
inline void Pre(int n){
  num=n;
  int g=Pow(3,(P-1)/num);
  w[0][0]=w[1][0]=1; for(int i=1;i0][i]=1LL*w[0][i-1]*g%P;
  for(int i=1;i1][i]=w[0][num-i];
}

inline void NTT(int *a,int n,int r){
  for(int i=1;iif(rev[i]>i) swap(a[i],a[rev[i]]);
  for(int i=1;i1)
    for(int j=0;j1)
      for(int k=0;kint x=a[j+k],y=1LL*a[j+k+i]*w[r][num/(i<<1)*k]%P;
        a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
      }
  if(!r) for(int i=0,inv=Pow(n,P-2);i1LL*a[i]*inv%P;
}

しかし、一部の狂った問題のモデル数はNTTモデルではなく、この時はCRT合併が必要です.
多項式求逆
例えばBZOJ 4555
最後にf(x)=f(x)g(x)+1 f(x)=f(x)g(x)+1が得られるが,これはn番目の項が大きくないため,f(x)f(x)を直接求めればよい非整列線形繰返しである.
f(x)=11−g(x)f(x)=1−g(x)を容易に得ることができ,h(x)=1−g(x)=1−g(x)=1−g(x)=1−g(x)とすれば,f(x)=(x)−1 f(x)=(h(x))−1となるので,h(x)h(x)の逆元求逆元を求めるとネット上での説明も多く,A(x)A(x)を求める逆元を取って、B(x)B(x)を求めて、A(x)B(x)≡1(modx 2 k)A(x)B(x)≡1(modx 2 k)A(x)≡1(modx 2 k)
G(x)G(x),A(x)G(x)≡1(modxn)A(x)G(x)≡1(modxn)を求めると
A(x)G(x)≡1(modxn) A ( x ) G ( x ) ≡ 1 ( mod x n )
A(x)G(x)−1≡0(modxn) A ( x ) G ( x ) − 1 ≡ 0 ( mod x n )
(A(x)G(x)−1)2≡0(modx2n) ( A ( x ) G ( x ) − 1 ) 2 ≡ 0 ( mod x 2 n )
A2(x)G2(x)−2A(x)G(x)+1≡0(modx2n) A 2 ( x ) G 2 ( x ) − 2 A ( x ) G ( x ) + 1 ≡ 0 ( mod x 2 n )
A(x)(2G(x)−A(x)G2(x))≡1(modx2n) A ( x ) ( 2 G ( x ) − A ( x ) G 2 ( x ) ) ≡ 1 ( mod x 2 n )
n=0 n=0の時、少し求めます
A(x)A(x)定数項の乗算逆元でよい
では、O(nlogn)O(n log⁡n)の複雑さでB(x)B(x)を2 G(x)−A(x)G 2(x)2 G(x)−A(x)G 2(x)G 2(x)を次の層演算に入れることで、総複雑度T(n)=T(n/2)+O(nlogn)=O(nlogn)T(n)/2)+O(n)=T(n/2)+O(n log⁡n)=O(n log⁡n)
void Inv(int *a,int *b,int n){
    if(n==1){
        b[0]=Pow(a[0],P-2); return ;
    }
    Inv(a,b,n>>1);
    for(int i=0;i0;
    int L=0; while(!(n>>L&1)) L++;
    for(int i=1;i1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<1,1); FFT(b,n<<1,1);
    for(int i=0;i1);i++)
        tmp[i]=(1LL*b[i]*2+P-1LL*tmp[i]*b[i]%P*b[i]%P)%P;
    FFT(tmp,n<<1,0);
    for(int i=0;i0;
}

適用すると、いくつかの現行の繰返し式を解くか、バーヌーリ数の前処理バーヌーリ数を求めることになります.http://blog.csdn.net/coldef/article/details/72876762
学習:http://blog.miskcoo.com/2015/05/polynomial-inverse
たじゅうしきかいせん
A(x)A(x)に対して、B(x)B(x)を求め、B 2(x)≡A(x)(modxn)B 2(x)≡A(x)(modxn)を同様に倍増させる方法である.G(x)G(x)が求められ、G 2(x)≡A(x)(modxn)G 2(x)≡A(x)(modxn)では
B2(x)≡G2(x)(modxn) B 2 ( x ) ≡ G 2 ( x ) ( mod x n )
[B2(x)−G2(x)]2≡0(modx2n) [ B 2 ( x ) − G 2 ( x ) ] 2 ≡ 0 ( mod x 2 n )
B4(x)−2B2(x)G2(x)+G4(x)≡0(modx2n) B 4 ( x ) − 2 B 2 ( x ) G 2 ( x ) + G 4 ( x ) ≡ 0 ( mod x 2 n )
B4(x)+2B2(x)G2(x)+G4(x)≡4B2(x)G2(x)(modx2n) B 4 ( x ) + 2 B 2 ( x ) G 2 ( x ) + G 4 ( x ) ≡ 4 B 2 ( x ) G 2 ( x ) ( mod x 2 n )
[B2(x)+G2(x)]2≡[2B(x)G(x)]2(modx2n) [ B 2 ( x ) + G 2 ( x ) ] 2 ≡ [ 2 B ( x ) G ( x ) ] 2 ( mod x 2 n )
B2(x)+G2(x)≡2B(x)G(x)(modx2n) B 2 ( x ) + G 2 ( x ) ≡ 2 B ( x ) G ( x ) ( mod x 2 n )
A(x)+G2(x)≡2B(x)G(x)(mod22n) A ( x ) + G 2 ( x ) ≡ 2 B ( x ) G ( x ) ( mod 2 2 n )
B(x)≡A(x)+G2(x)2G(x)(mod22n) B ( x ) ≡ A ( x ) + G 2 ( x ) 2 G ( x ) ( mod 2 2 n )
n=0 n=0の場合
B(x)B(x)の定数項が1(なぜか、Mancheryが言った)
時間の複雑さ
T(n)=T(n/2)+多項式求逆の複雑さ+O(nlogn)=O(nlogn)T(n)=T(n/2)+多項式求逆の複雑さ+O(nlogn)=O(nlogn)
void Sqrt(int *a,int *b,int n){
    if(n==1) return void(b[0]=1);
    Sqrt(a,b,n>>1);
    memset(invb,0,sizeof(int)*n); Inv(b,invb,n);
    int L=0; while(!(n>>L&1)) L++;
    for(int i=1;i1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<for(int i=0;i0;
    FFT(invb,n<<1,1); FFT(tmp,n<<1,1); 
    for(int i=0,inv2=P+1>>1;i1);i++) tmp[i]=1LL*tmp[i]*inv2%P*invb[i]%P;
    FFT(tmp,n<<1,0);
    for(int i=0,inv2=P+1>>1;i1LL*b[i]*inv2+tmp[i])%P;
}

学習:http://www.cnblogs.com/acha/p/6472798.html
たじゅうしきかたとり
この型取りは係数型取りではなく多項式型多項式である
UPD:暴力的多項式は型取りが簡単で、型多項式をM(x)=Σni=0 aixi M(x)=Σi=0 na i x iとすると、xn=Σn−1 i=0 aianxi x n=Σi=0 n−1 a a n x iが得られ、n n nより大きいすべての項目を持ち込めばよい
A(x)A(x)に対して、B(x)B(x)に対して、degA(x)=n、degB(x)=m(m)=m(m≦n)deg⁡A(x)=n、deg⁡B(x)=m(m≦n)D(x)D(x)D(x)R(x)、R(x)R(x)を求め、A(x)=D(x)B(x)+R(x)+R(x)+R(x)+R(x)+R(x)=x)A(x(x)=D(x)B(x)+R(x)を満足し、多項式加減算乗乗乗数で乗算すると、多項式加減算乗乗乗乗乗乗乗乗乗R(x)R(x)を得て1 x 1 xを多項式に持ち込む
A(1x)=D(1x)B(1x)+R(1x) A ( 1 x ) = D ( 1 x ) B ( 1 x ) + R ( 1 x )
両側にxn x nを乗じる
xnA(1x)=xn−mD(1x)xmB(1x)+xnR(1x) x n A ( 1 x ) = x n − m D ( 1 x ) x m B ( 1 x ) + x n R ( 1 x )
fR(x)=xdegf(x)f(1 x)f R(x)=x deg⁡f(x)f(1 x)
則
AR(x)=DR(x)BR(x)+xn−degR(x)RR(x) A R ( x ) = D R ( x ) B R ( x ) + x n − deg ⁡ R ( x ) R R ( x )
degR(x)では
AR(x)≡DR(x)BR(x)(modxn−m+1) A R ( x ) ≡ D R ( x ) B R ( x ) ( mod x n − m + 1 )
AR(x)BR(x)≡DR(x)(modxn−m+1) A R ( x ) B R ( x ) ≡ D R ( x ) ( mod x n − m + 1 )
このように多項式は逆を求めてD(x)D(x)を求めることができて、持って帰ってプラスして減らしてR(x)R(x)を得ることができます
inline void Div(int *a,int n,int *b,int m){
  static int tmp[N],A[N],B[N];
  for(int i=0;i1];
  for(int i=0;i1];
  int nn=1,d=n-m+1; for(;nn1;nn<<=1);
  for(int i=n;i0;
  for(int i=d;i0;
  for(int i=0;i0;
  Inv(t,B,nn);
  for(int i=d;i0;
  Mul(A,B,max(n,d));
  for(int i=d;i<=n<<1;i++) A[i]=0;
  for(int i=0;iif(i>d-i-1) swap(A[i],A[d-i-1]);
  for(int i=0;imax(d,m));
  for(int i=0;i

学習:http://blog.miskcoo.com/2015/05/polynomial-division
大体これらのアルゴリズムです
还有多項式はexpなどを求めます...まだできません...先に穴を掘ります
多項式exp
金策の2015年の国家合宿チーム論文を参考にすることができます
B(x)=eA(x) B ( x ) = e A ( x )
両側同導
B′(x)=eA(x)A′(x) B ′ ( x ) = e A ( x ) A ′ ( x )
B′(x)=B(x)A′(x) B ′ ( x ) = B ( x ) A ′ ( x )
各項目に対する係数
(i+1)bi+1=∑j+k=ibjak+1(k+1) ( i + 1 ) b i + 1 = ∑ j + k = i b j a k + 1 ( k + 1 )
つまり
bi=1i∑j=1ijajbi−j b i = 1 i ∑ j = 1 i j a j b i − j
分治FFT
ニュートン反復は~~
void exp(int *a,int *b,int l,int r){
  if(l==r){
    add(b[l],1LL*inv[l]*a[l]%P); return ;
  }
  int mid=l+r>>1;
  exp(a,b,l,mid); 
  static int tmpa[N],tmpb[N];
  for(int i=l;i<=mid;i++) tmpa[i-l]=b[i];
  for(int i=1;i<=r-l;i++) tmpb[i]=a[i];
  int m,L=0;
  for(m=1;m1;m<<=1,L++); m<<=1;
  for(int i=1;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<m,1); NTT(tmpb,m,1);
  for(int i=0;i<m;i++) tmpa[i]=1LL*tmpa[i]*tmpb[i]%P;
  NTT(tmpa,m,0);
  for(int i=mid+1;i<=r;i++) add(b[i],1LL*inv[i]*tmpa[i-l]%P);
  for(int i=0;i<m;i++) tmpa[i]=tmpb[i]=0;
  exp(a,b,mid+1,r);
}

適用
生成関数の畳み込み演算による最適化遷移とかよくある多項式の多点評価とかhttp://blog.miskcoo.com/2015/05/polynomial-multipoint-eval-and-interpolation私は実現しました...しかし大きな定数を持っていて、走るのは暴力に及ばない...
NOI 2017 D 13 T 3は、多項式型取りを用いて線形繰返し数列のn番目の項(nが大きい)、すなわちf(x)=f(x)g(x)+1 f(x)=f(x)g(x)+1のn番目の項を解決し、
具体的にはフォーク姉さんの論文とPicksのブログを見てhttp://www.docin.com/p-724323397.html http://picks.logdown.com/posts/197262-polynomial-divisioncodechef rngはテンプレートの問題で、上で大神のコードを見ることができます
他の用途...慣例はまず穴を掘る
私は能力の制限のため、コードはただ参考に供して、もし疑問があるならば、私と討論することを歓迎します
拡張ユークリッドアルゴリズム
マトリックス高速べき乗を用いてフィボナッチ数列を求める