2333

時間  2019-11-11
原文   原文鏈接

\(\color{red}{\text{多項式全集之二 任長任模的FFT:}}\)

三模NTT實現任模FFT：

\(1.\)爲何要用MTT：當\(p\)不是NTT模數或者多項式長度大於模數限制時，就要使用MTT。算法

\(2.\)MTT的使用原理：咱們對初始多項式取模，那麼若是在不取模卷積狀況下，答案\(x\)不會超過\(N\times p^2\)。咱們取三個NTT模數\(p_1,p_2,p_3\)，分別作多項式乘法，獲得\(x\)分別\(mod~p_1,p_2,p_3\)的答案，經過CRT合併能夠獲得\(x~mod~p_1p_2p_3\)的答案，若是\(x<p_1p_2p_3\)那麼就能夠獲得準確的答案，再對\(p\)取模便可。優化

\(3.\)CRT合併的小優化：ui

\(step~0:\)初始式子
\[{\begin{cases}x\equiv c_1(mod~p_1)\\x\equiv c_2(mod~p_2)\\x\equiv c_3(mod~p_3)\end{cases}}\]
\(step~1:\)把一式二式合併(LL範圍內)。
\[{\begin{cases}x\equiv a(mod~p_1p_2)\\x\equiv c_3(mod~p_3)\end{cases}}\]
\(step~2:\)再次合併(不須要\(long~double\) 快速乘)。spa

\(4.\)經常使用NTT模數：code

如下模數的共同\(g=3189\)get

\(p=r\times 2^k+1\) \(k\) \(g\)
\(104857601\) \(22\) \(3\)
\(167772161\) \(25\) \(3\)
\(469762049\) \(26\) \(3\)
\(950009857\) \(21\) \(7\)
\(998244353\) \(23\) \(3\)
\(1004535809\) \(21\) \(3\)
\(2013265921\) \(27\) \(31\)
\(2281701377\) \(27\) \(3\)
\(3221225473\) \(30\) \(5\)
拆係數FFT(CFFT)實現任模FFT：

\(1.\)實現原理：運用實數FFT不取模作乘法，而後取模迴歸到整數。可是因爲偏差較大（值域是\(10^{23}\)），咱們令\(t=\sqrt{m}\)把係數\(a_i=k_it+b_i\)，對\(k_i,t_i\)交叉作四遍卷積，求出答案按係數貢獻取模加入。it

\(2.\)可按合併DFT的方法優化DFT次數。table

\(bluestein\)算法實現任長FFT：

當\(m\)不是\(2\)的冪次的時候，咱們從式子入手：
令\(X_i=a_iw_m^{\frac {i^2} 2},Y_i=w_m^{\frac{-i^2}2}\)模板

\(\begin{align*} A_k & = \sum_{j=0}^{m-1}a_jw_m^{jk}\\ & = \sum_{j=0}^{m-1}a_jw_m^{\frac{j^2+k^2-{(k-j)}^2}{2}}\\ &=w_m^{\frac {k^2} 2}\sum_{j=0}^{m-1}a_jw_m^{\frac{j^2} 2}w_m^{\frac{-{(k-j)}^2}{2}}\\ &=w_m^{\frac {k^2} 2}\sum_{j=0}^{m-1}X_jY_{k-j} \end{align*}\)class

喜聞樂見的模板：

三模NTT模板（注意：不能夠MTT回來，由於係數會取模）

namespace MTT{
    typedef long long LL;
    int n, m;
    LL p, mod;
    const LL p1 = 998244353;
    const LL p2 = 1004535809;
    const LL p3 = 104857601;
    const int g = 3189;
    LL a[300005], b[300005], c[300005], cpa[300005], cpb[300005];
    LL c3[300005], c1[300005], c2[300005];
    LL qpow(LL a, LL b, LL mod) {
        LL ans = 1;
        while(b) {
            if(b & 1)   ans = ans * a % mod;
            a = a * a % mod;
            b >>= 1;
        }
        return ans;
    }
    const LL inv12 = qpow(p1, p2 - 2, p2);
    const LL inv123 = qpow(p1 * p2 % p3, p3 - 2, p3);
    struct p_l_e{
        int wz[300005];
        void MTT(LL *a, int N, int op) {
            for(int i = 0; i < N; i++)
                if(i < wz[i]) swap(a[i], a[wz[i]]);
            for(int le = 2; le <= N; le <<= 1) {
                int mid = le >> 1;
                LL wn = qpow(g, (mod - 1) / le, mod);
                if(op == -1) wn = qpow(wn, mod - 2, mod);
                for(int i = 0; i < N ;i += le) {
                    LL w = 1, x, y;
                    for(int j = 0; j < mid; j++) {
                        x = a[i + j];
                        y = a[i + j + mid] * w % mod;
                        a[i + j] = (x + y) % mod;
                        a[i + j + mid] = (x - y + mod) % mod;
                        w = w * wn % mod;
                    }
                }
            }
        }
        void mult(LL *a, LL *b, LL *c, int M) {
            int N = 1, len = 0;
            while(N < M) N <<= 1, len++;
            for(int i = 0; i < N; i++)
                wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
            MTT(a, N, 1); MTT(b, N, 1);
            for(int i = 0; i < N; i++)  c[i] = a[i] * b[i] % mod;
            MTT(c, N, -1);
            LL t = qpow(N, mod - 2, mod);
            for(int i = 0; i < N; i++)  c[i] = c[i] * t % mod;
        }

    }PLE;
    LL CRT(LL c1, LL c2, LL c3) {
        LL x = (c1 + p1 * ((c2 - c1 + p2) % p2 * inv12 % p2));
        LL y = (x % p + p1 * p2 % p * ((c3 - x % p3 + p3) % p3 * inv123 % p3) % p) % p;
        return y;
    }
    void merge(LL *c1, LL *c2, LL *c3, LL *c, int N) {
        for(int i = 0; i < N; i++)
            c[i] = CRT(c1[i], c2[i], c3[i]);
        return;
    }
    void main() {
        scanf("%d%d%lld", &n, &m, &p); n++; m++;
        for(int i = 0; i < n; i++)  scanf("%lld", &a[i]);
        for(int i = 0; i < m; i++)  scanf("%lld", &b[i]);
        mod = p1; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c1, n + m - 1);
        mod = p2; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c2, n + m - 1);
        mod = p3; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c3, n + m - 1);
        merge(c1, c2, c3, c, n + m - 1);
        for(int i = 0; i < n + m - 1; i++)  printf("%lld ", (c[i] % p + p) % p);
        return;
    }
}

拆係數FFT模板（注意：相同係數的兩項能夠合併一塊兒IDFT。採用共軛優化法，只進行四次DFT）

namespace CFFT{
    typedef long long LL;
    int n, m, p ,sqrp; 
    int a[300005], b[300005];
    const long double pi = acos(-1);
    struct cp{
        long double x, y;
        cp() {x = y = 0;}
        cp(long double X,long double Y) {x = X; y = Y; }
        cp conj() {return (cp) {x, -y};}
    }ka[300005], kb[300005], ta[300005], tb[300005], kk[300005], kt[300005], tt[300005], c[300005], I(0, 1), d[300005];
     
    cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
    cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
    cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
    cp operator* (const cp &a, long double b) {return (cp){a.x * b, a.y * b};}
    cp operator/ (const cp &a, long double b) {return (cp){a.x / b, a.y / b};}  
    struct p_l_e{
        int wz[300005];
        void FFT(cp *a, int N, int op){
            for(int i = 0; i < N; i++)
                if(i < wz[i])   swap(a[i], a[wz[i]]);
            for(int le = 2; le <= N; le <<= 1){
                int mid = le >> 1;
                cp x, y, w, wn = (cp){cos(op * 2 * pi / le), sin(op * 2 * pi / le)};
                for(int i = 0; i < N; i += le){
                    w = (cp){1, 0};
                    for(int j = 0; j < mid; j++){
                        x = a[i + j];
                        y = a[i + j + mid] * w;
                        a[i + j] = x + y;
                        a[i + j + mid] = x - y;
                        w = w * wn;
                    }
                }
            } 
        }
        void D_FFT(cp *a, cp *b, int N, int op){
            for(int i = 0; i < N; i++)  d[i] = a[i] + I * b[i];
            FFT(d, N, op);
            d[N] = d[0];
            if(op == 1){
                for(int i = 0; i < N; i++){
                    a[i] = (d[i] + d[N - i].conj()) / 2;
                    b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;
                }
            } else {
                for(int i = 0; i < N; i++){
                    a[i] = cp(d[i].x, 0);
                    b[i] = cp(d[i].y, 0);
                }
            }
            d[N] = cp(0, 0);
        }
        void mult(int *a, int *b, int M){
            int N = 1, len = 0;
            while(N < M) N <<= 1, len++;
            for(int i = 0; i < N; i++)
                wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
            for(int i = 0; i < N; i++){
                ka[i].x = a[i] >> 15;
                kb[i].x = b[i] >> 15;
                ta[i].x = a[i] & 32767;
                tb[i].x = b[i] & 32767;
            }
            D_FFT(ta, ka, N, 1); D_FFT(tb, kb, N, 1);
            for(int i = 0; i < N; i++){
                kk[i] = ka[i] * kb[i];
                kt[i] = ka[i] * tb[i] + ta[i] * kb[i];
                tt[i] = ta[i] * tb[i];
            }
            D_FFT(tt, kk, N, -1); FFT(kt, N, -1);
            for(int i = 0; i < N; i++){
                tt[i] = tt[i] / N;
                kt[i] = kt[i] / N;
                kk[i] = kk[i] / N;
            }
        }
    }PLE;
    void main() {
        scanf("%d%d%d", &n, &m, &p); n++; m++;
        for(int i = 0; i < n; i++)  scanf("%d", &a[i]),a[i] = a[i] % p;
        for(int i = 0; i < m; i++)  scanf("%d", &b[i]),b[i] = b[i] % p;
        PLE.mult(a, b, n + m - 1);
        for(int i = 0; i < n + m - 1; i++)
            printf("%lld ",(((((LL)round(kk[i].x)) % p) << 30) + ((((LL)round(kt[i].x)) % p) << 15) + ((LL)round(tt[i].x)) % p) % p);
    }
}

\(blue\_stein\)模板：

struct polynie {
    CP getw(int m, int k, int op) {
        return CP(cos(2 * pi * k / m), op * sin(2 * pi * k / m));
    }
    int wz[MAXN];
    CP A[MAXN], B[MAXN], C[MAXN];
    void FFT(CP *a, int N, int op) {
        rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);
        for(int l = 2; l <= N; l <<= 1) {
            int mid = l >> 1;
            CP x, y, w, wn = CP(cos(pi / mid), sin(op * pi / mid));
            for(int i = 0; i < N; i += l) {
                w = CP(1, 0);
                rop(j, 0, mid) {
                    x = a[i + j];
                    y = w * a[i + j + mid];
                    a[i + j] = x + y;
                    a[i + j + mid] = x - y;
                    w = w * wn;
                }
            }
        }
    }
    void mult(CP *a, CP *b, CP *c, int M) {
        int N = 1, len = 0;
        while(N < M) N <<= 1, len++;
        rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
        FFT(a, N, 1); FFT(b, N, 1);
        rop(i, 0, N) c[i] = a[i] * b[i];
        FFT(c, N, -1);
        rop(i, 0, N) c[i].x = c[i].x / N, c[i].y = c[i].y / N;
    }
    void blue_stein(CP *a, int M, int op) {
        int M2 = M << 1;
        memset(A, 0, sizeof(A));
        memset(B, 0, sizeof(B));
        memset(C, 0, sizeof(C));
        rop(i, 0, M) A[i] = a[i] * getw(M2, 1ll * i * i % M2, op);
        rop(i, 0, M2) B[i] = getw(M2, 1ll * (i - M) * (i - M) % M2, -op);
        mult(A, B, C, M2 + M - 1);
        rop(i, 0, M) a[i] = C[i + M] * getw(M2, 1ll * i * i % M2, op);
        if(op == -1) rop(i, 0, M) a[i].x = a[i].x / M, a[i].y = a[i].y / M;
    }
}PLE;

\(\color{red}{\text{多項式全集之三 多項式求逆與除法:}}\)

多項式求逆：

\(1.\)問題描述：

已知\(F(x)\)，且\(F(x)G(x)\equiv 1 (mod~x^n)\)，求\(G(x)\)

\(2.\)推導過程：

\begin{align}
B(x)&\equiv F(x)^{-1}&(mod~x^{\lceil \frac n 2 \rceil})\
因爲
G(x)&\equiv F(x)^{-1}& (mod~x^n)\
因此
G(x)&\equiv F(x)^{-1}& (mod~x^{\lceil \frac n 2 \rceil})\
G(x)& \equiv B(x)&(mod~x^{\lceil \frac n 2 \rceil})\
G(x)-B(x)& \equiv 0&(mod~x^{\lceil \frac n 2 \rceil})\
\end{align}
兩邊平方，得：

因爲\([G(x)-B(x)]^2\)的第\(k<n\)項爲
\[\sum_{i=0}^k[g_ix^i-b_ix^i][g_{k-i}x^{k-i}-b_{k-i}x^{k-i}]\]
\(i,k-i\)必定有一項\(<\frac n 2\)，因此
\begin{align}
[G(x)-B(x)]^2&\equiv 0 &(mod~x^n)\
G^2(x)+B^2(x)-2G(x)B(x)&\equiv 0&(mod~x^n)
\end{align}
兩邊同乘\(A(x)\)，得：
\begin{align}
A(x)G^2(x)+A(x)B^2(x)-2A(x)G(x)B(x)&\equiv 0& (mod~x^n)\
G(x)+A(x)B^2(x)-2B(x)&\equiv 0& (mod~x^n)\
G(x)&\equiv 2B(x)-A(x)B^2(x)&(mod~x^n)\
G(x)&\equiv B(x)[2-A(x)B(x)]&(mod~x^n)\
\end{align}

多項式除法：

\(1.\)問題描述：

已知一個\(n\)次多項式\(A(x)\)，一個\(m\)次多項式\(B(x)\)，且\(A(x)=B(x)C(x)+D(x)\)，求\(n-m\)次多項式\(C(x)\)，\(<m\)次多項式\(D(x)\)。

\(2.\)推導過程：

由\(A(x) = B(x)C(x)+D(x)\)得：
\begin{align}
A(\frac 1 x)&=B(\frac 1 x)C(\frac 1 x)+D(\frac 1 x)\
x^nA(\frac 1 x)&=x^nB(\frac 1 x)C(\frac 1 x)+x^nD(\frac 1 x)\
x^nA(\frac 1 x)&=x^mB(\frac 1 x)x^{n-m}C(\frac 1 x)+x^{m-1}x^{n-m+1}D(\frac 1 x)\
A_r(x)&=B_r(x)C_r(x)+x^{n-m+1}D(x)\
A_r(x)&=B_r(x)C_r(x)&(mod ;x^{n-m+1})\
B_r(x)&=A_r^{-1}(x)C_r(x)&(mod ;x^{n-m+1})\
\end{align}
求逆可得\(B_r(x)\)，再反轉得\(B(x)\)，而後乘\(C(x)\)去減\(A(x)\)得\(D(x)\).

喜聞樂見的模板：
namespace INV{
    typedef long long LL;
    int n, a[300005], b[300005];
    const int mod = 998244353;
    const int g = 3189;
    int qpow(int a, int b){
        int ans = 1;
        while(b){
            if(b & 1)   ans = 1ll * ans * a % mod;
            a = 1ll * a * a % mod;
            b >>= 1;
        }
        return ans;
    }
    struct p_l_e{
        int wz[300005], i_c[300005];
        void NTT(int *a, int N, int op){
            for(int i = 0; i < N; i++)
                if(i < wz[i]) swap(a[i], a[wz[i]]);
            for(int le = 2; le <= N; le <<= 1){
                int mid = le >> 1, wn = qpow(g, (mod - 1) / le);
                if(op == -1) wn = qpow(wn, mod - 2);
                for(int i = 0; i < N; i += le){
                    LL w = 1; int x, y;
                    for(int j = 0; j < mid; j++){
                        x = a[i + j];
                        y = w * a[i + j + mid] % mod;
                        a[i + j] = (x + y) % mod;
                        a[i + j + mid] = (x - y + mod) % mod;
                        w = w * wn % mod;
                    }
                }
            }
        }
        int init(int M){
            int N = 1, len = 0;
            while(N < M) N <<= 1, len++;    
            for(int i = 0; i < N; i++)
                wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
            return N;
        }
        void INV(int *a, int *b, int deg){
            if(deg == 1){b[0] = qpow(a[0], mod - 2); return;}
            INV(a, b, (deg + 1) >> 1);
            int N = init(deg + deg - 1);
            for(int i = 0; i < deg; i++) i_c[i] = a[i];
            for(int i = deg; i < N; i++) i_c[i] = 0;
            NTT(b, N, 1);NTT(i_c, N, 1);
            for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * (2 - 1ll * b[i] * i_c[i] % mod + mod) % mod;
            NTT(b, N, -1);
            int t = qpow(N, mod - 2);
            for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;
            for(int i = deg; i < N; i++) b[i] = 0;
        }
    }PLE;
    
    void main(){
        scanf("%d", &n);
        for(int i = 0; i < n; i++)  scanf("%d", &a[i]);
        PLE.INV(a, b, n);
        for(int i = 0; i < n; i++)  printf("%d ",b[i]);
    }
}

\(\color{red}{\text{多項式全集之四 多項式ln與exp:}}\)

多項式ln：

\(1.\)作法：

設
\[G(x) =\ln F(x)\]
兩邊求導得
\[G'(x)=\frac {F'(x)} {F(x)}\]
積分回去便可。

\(2.\)應用：
\[e^{F(x)}=\sum_{k \ge 0}\frac{F^k(x)}{k!}=G(x)\]
\[F(x) = \ln G(x)\]
這個的組合意義是：無序組合。

設\(F(x)\)，\(f_i\)表示一些東西，那麼這些東西有序組合的方案數爲
\[F^0(x) + F^1(x)+F^2(x)+\cdots=\frac 1 {1 - F(x)}\]
而無序組成的方案數爲：
\[\frac{F^0(x)}{0!} + \frac {F^1(x)}{1!}+\frac{F^2(x)}{2!}+\cdots=e^{F(x)}\]
若是無序組合方案數好求，那麼求\(\ln\)就能獲得\(F(x)\)。

例題

多項式\(\exp\)：
喜聞樂見的代碼：

多項式\(\ln\):

namespace PLE_ln{
  struct polyme {
      int li[SZ], wz[SZ];
      void NTT(int *a, int N, int op) {
          rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);
          for(int l = 2; l <= N; l <<= 1) {
              int mid = l >> 1;
              int x, y, w, wn = qpow(g, (mod - 1) / l);
              if(op) wn = qpow(wn, mod - 2);
              for(int i = 0; i < N; i += l) {
                  w = 1;
                  for(int j = 0; j < mid; ++j) {
                      x = a[i + j]; y = 1ll * w * a[i + j + mid] % mod;
                      a[i + j] = (x + y) % mod;
                      a[i + j + mid] = (x - y + mod) % mod;
                      w = 1ll * w * wn % mod;
                  }
              }
          }
      }
      void qd(int *a, int *b, int n) {
          rop(i, 0, n) b[i] = 1ll * a[i + 1] * (i + 1) % mod;
      }
      void jf(int *a, int *b, int n) {
          rop(i, 1, n) b[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;
      }
      void mult(int *a, int *b, int *c, int M) {
          int N = 1, len = 0;
          while(N < M) N <<= 1, len ++;
          rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
          NTT(a, N, 0); NTT(b, N, 0);
          rop(i, 0, N) c[i] = 1ll * a[i] * b[i] % mod;
          NTT(c, N, 1);
          int t = qpow(N, mod - 2);
          rop(i, 0, N) c[i] = 1ll * c[i] * t % mod;
      }
      void inv(int *a, int *b, int deg) {
          if(deg == 1) {b[0] = qpow(a[0], mod - 2) % mod; return;}
          inv(a, b, (deg + 1) >> 1);
          rop(i, 0, deg) li[i] = a[i];
          int N = 1, len = 0;
          while(N < deg + deg - 1) N <<= 1, len ++;
          rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
          rop(i, deg, N) li[i] = 0;
          NTT(li, N, 0); NTT(b, N, 0);
          rop(i, 0, N) b[i] = 1ll * b[i] * (2 - 1ll * li[i] * b[i] % mod + mod) % mod;
          NTT(b, N, 1);
          int t = qpow(N, mod - 2);
          for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;
          rop(i, deg, N) b[i] = 0;
      }
  }PLE;
  int a[SZ], da[SZ], ia[SZ], dla[SZ], la[SZ], n;

  void main() {

      scanf("%d", &n);
      rop(i, 0, n) scanf("%d", &a[i]);
      PLE.qd(a, da, n);
      PLE.inv(a, ia, n);
      PLE.mult(ia, da, dla, n + n - 1);
      PLE.jf(dla, la, n);
      rop(i, 0, n) printf("%d ", la[i]);
  }
}

\(\color{red}{\text{多項式全集之五 多項式快速冪與開根:}}\)

多項式快速冪方法一：

\(\color{red}{\text{多項式全集之六 多項式快速插值與多點求值:}}\)

\(\color{red}{\text{多項式全集之七 拉格朗日插值：}}\)

每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程