• 多项式相关操作学习笔记


    多项式相关操作学习笔记

    标签: 多项式


    说在前边

    记录一下相关的多项式操作,顺便存个模板。(多点求值之后的部分,有点写不动了。。。留坑留坑


    多项式


    定义

    给定一个环(R)((R)通常是交换环,可以是有理数、实数或者复数等等)以及一个未知数(X),则任何形同:$$a_0 + a_1X + ... +a_{n-1}X^{n-1} + a_nX^n$$的代数表达式叫做(R)上的一元多项式。其中(a_0,a_1,... ,a_n)(R)中的元素。未知数不代表任何值,但环 (R)上的所有运算都对它适用。在不至于混淆的情形下,一般将一元多项式简称为多项式。可以证明,两个多项式的和、差与积仍然是多项式,即多项式组成一个环 (R[X]),称为(R)上的(一元)多项式环。—— From Wiki


    多项式乘法

    计算两个多项式相乘时,首先使用乘法对加法的分配律将各项拆出,然后运用乘法结合律整合每一项,最后和加法一样整合同类项,就能得到乘积多项式,形式化的说:

    [P = a_0 + a_1X + ... +a_{n-1}X^{n-1} + a_nX^n \ Q = b_0 + b_1X + ... +b_{n-1}X^{n-1} + b_nX^m \ PQ = sum_{i+j=0}a_ib_j +sum_{i+j=1}a_ib_jX + ... + sum_{i+j=n+m}a_ib_jX^{n+m} ]


    离散傅里叶变换(DFT)

    定义

    对于(N)点序列(x[n]),它的(DFT)为:

    [X[k] = sum_{n = 0}^{N-1}e^{-ifrac{2pi}{N}nk}x[n] ~~k = 0,1,...,N-1]

    离散傅里叶变换的逆变换(IDFT)为:

    [x[n] = frac{1}{N}sum_{k = 0}^{N-1}e^{ifrac{2pi}{N}nk}X[k]~~n = 0,1,...,N-1 ]


    DFT与圆周卷积

    1. 圆周移位
      (f(n) = x((n+m))_NR_N(n)),成(f(n))(x(n))(m)点圆周移位序列。
      步骤:1)将x(n)以N为周期进行周期延拓;2)移位m点;3)取主值序列

    2. (x(n), y(n))进行DFT,(x(n) leftrightarrow X(k), y(n) leftrightarrow Y(k))

    [F(k) = X(k)Y(k) leftrightarrow f(n) = sum_{m=0}^{N-1}x(m)y((n-m))_NR_N(n) ]


    有限长序列的圆周卷积与线性卷积

    (x(n)):(M)点序列,(y(n)):(N)点序列

    线性卷积:$$x(n)*y(n) = sum_{m=-infty}^{infty} x(m)y(n-m)$$

    可知(f(n))的非(0)区间为([0,M+N-2])

    圆周卷积是线性卷积以(L)为周期延拓后,取([0,L-1])间的主值序列。
    (L geq N+M-1) 时,圆周卷积等价于线性卷积


    DSP角度的解释

    从DTFT到DFS,从DFS到DFT,从DFT到FFT,从一维到二维


    一些变换

    快速傅里叶变换(FFT)

    很久之前参考算导写的FFT学习笔记


    快速哈特莱变换(FHT)

    FHT是一种实序列变换


    快速数论变换(NTT)

    NTT是在整数域进行的,FFT中通过n次单位福根运算,即满足(W^n = 1)(W),而NTT,将(g^{frac{P-1}{N}}(mod~P)) 看作与 (e^{-frac{2pi i}{N}}) 等价,这里(g)是模素数(P)的原根。综上,数论变换的公式如下:

    [X_k ≡ sum_{n=0}^{N-1}x_ng^{nk}(mod~P)~~~ k = 0, 1, ...,N-1 ]

    逆变换为:

    [x_n ≡ frac{1}{N}sum_{k=0}^{N-1}X_kg^{-nk}(mod~P)~~~ n = 0, 1, ...,N-1 ]

    Code

    #include <cstdio>
    #include <algorithm>
    typedef long long ll;
    constexpr int N = 1 << 18 | 5; 
    constexpr int Mod = 998244353;
    constexpr int  G = 3;
    ll qpow(ll a, ll b) {
        ll ans = 1;
        for(; b; a = (a * a) % Mod, b >>= 1)
            if(b & 1) ans = (ans * a) % Mod;
        return ans;
    }
    int siz, rev[N];
    ll w[N];
    inline void init(int n) {
        for(siz = 1; siz < n; siz <<= 1);
        for(int i = 0; i < siz; ++i)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (siz >> 1));
        ll pr = qpow(G, (Mod-1) / siz);
        w[siz / 2] = 1;
        for(int i = 1; i < siz / 2; ++i) w[siz/2+i] = (w[siz/2+i-1]*pr) % Mod;
        for(int i = siz / 2 - 1; i >= 0; --i) w[i] = w[i << 1];
    }
    inline void DFT(ll A[], int L) {
        static ll t[N];
        int shift = __builtin_ctz(siz) - __builtin_ctz(L);
        for(int i = 0; i != L; ++i) t[rev[i] >> shift] = A[i];
        for(int d = 1; d != L; d <<= 1)
            for(int i = 0; i != L; i += d << 1)
                for(int j = 0; j != d; ++j) {
                    ll tmp = t[i + d + j] * w[d + j] % Mod;
                    t[i + d + j] = t[i + j] - tmp;
                    t[i + j] = t[i + j] + tmp;
                    if(t[i + d + j] < 0) t[i + d + j] += Mod;
                    if(t[i + j] >= Mod) t[i + j] -= Mod;
                }
        for(int i = 0; i != L; ++i) A[i] = t[i] % Mod;
    }
    inline void IDFT(ll A[], int L) {
        std::reverse(A + 1, A + L);
        DFT(A, L);
        ll Inv_L = qpow(L, Mod - 2);
        for(int i = 0; i != L; ++i) (A[i] *= Inv_L) %= Mod;
    }
    int n, m, L;
    ll a[N], b[N];
    int main() {
    	scanf("%d%d",&n,&m);
    	for(int i = 0; i <= n; ++i) scanf("%lld", &a[i]);
        for(int i = 0; i <= m; ++i) scanf("%lld", &b[i]);
        for(L = 1; L <= m+n; L <<= 1);  init(L);
        DFT(a, L); DFT(b, L);
        for(int i = 0; i <= L; ++i) a[i] = a[i] * b[i] % Mod;
        IDFT(a, L);
        for(int i = 0; i <= n+m; ++i) printf("%lld ",a[i]); puts("");
        return 0;
    }
    

    任意模数FFT

    题意:给两个多项式(f(x))(g(x))以及一个模数(p(p leq 10^9)), 求(f * g (mod ~ p))

    做法:任意模数NTT,最大的数为(p^2 max(n,m) leq 10^{23}), 所以一般选三个模数即可,求出这三个模数下的答案,然后中国剩余定理。

    设当前这一项的系数为(x),三个模数分别为(A, B, C), 那么:

    [x ≡ x_1 (mod ~ A) \ x ≡ x_2 (mod ~ B) \ x ≡ x_3 (mod ~ C) ]

    先合并前两项:

    [x_1 + k_1A = x_2 + k_2B = x \ x_1 + k_1A ≡ x_2 ~(mod ~ B) \ k_1 ≡ ~frac {x_2 - x_1} {A}~(mod ~B) ]

    及求出(x ≡ x_1 + k_1A ~(mod ~AB)), 记作(x_4 = x_1 + k_1A)

    [x_4 + k_4AB = x_3 + k_3C \ x_4 + k_4AB ≡ x_3 ~(mod ~C) \ k_4 ≡ frac{x3 - x4}{AB} ~(mod ~C) ]

    n那么(x ≡ x_4 + k_4AB ~(mod~ABC)),因为(x < ABC), 所以(x = x_4 + k_4AB)

    Code[洛谷P4245]

    #include <cstdio>
    #include <cstring>
    #include <iostream>
    #include <algorithm>
    typedef long long ll;
    constexpr int N = 1 << 18 | 5;
    constexpr ll Mod1 = 469762049;
    constexpr ll Mod2 = 998244353;
    constexpr ll Mod3 = 1004535809;
    constexpr ll M = Mod1 * Mod2;
    constexpr int  G = 3;
    using namespace std;
    inline ll Mul(ll a, ll b, ll Mod) {
        ll ans = 0; ((a %= Mod) += Mod)%= Mod;
        for(; b; a = (a + a) % Mod, b >>= 1)
            if(b & 1) ans = (ans + a) % Mod;
        return ans;
    }
    inline ll qpow(ll a, ll b, ll Mod) {
        ll ans = 1; ((a %= Mod) += Mod)%= Mod;
        for(; b; a = (a * a) % Mod, b >>= 1)
            if(b & 1) ans = (ans * a) % Mod;
        return ans;
    }
    int siz, rev[N];
    ll w[4][N];
    inline void init_w(int typ, ll Mod) {
        ll pr = qpow(G, (Mod-1) / siz, Mod);
        w[typ][siz / 2] = 1;
        for(int i = 1; i < siz / 2; ++i) w[typ][siz/2+i] = (w[typ][siz/2+i-1]*pr) % Mod;
        for(int i = siz / 2 - 1; i >= 0; --i) w[typ][i] = w[typ][i << 1];
    }
    inline void init(int n) {
        for(siz = 1; siz < n; siz <<= 1);
        for(int i = 0; i < siz; ++i)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (siz >> 1));
        init_w(1,Mod1); init_w(2,Mod2); init_w(3,Mod3);
    }
    inline void DFT(ll A[], int L, int typ, ll Mod) {
        static ll t[N];
        int shift = __builtin_ctz(siz) - __builtin_ctz(L);
        for(int i = 0; i != L; ++i) t[rev[i] >> shift] = A[i];
        for(int d = 1; d != L; d <<= 1)
            for(int i = 0; i != L; i += d << 1)
                for(int j = 0; j != d; ++j) {
                    ll tmp = t[i + d + j] * w[typ][d + j] % Mod;
                    t[i + d + j] = t[i + j] - tmp;
                    t[i + j] = t[i + j] + tmp;
                    if(t[i + d + j] < 0) t[i + d + j] += Mod;
                    if(t[i + j] >= Mod) t[i + j] -= Mod;
                }
        for(int i = 0; i != L; ++i) A[i] = t[i] % Mod;
    }
    inline void IDFT(ll A[], int L, int typ, ll Mod) {
        std::reverse(A + 1, A + L);
        DFT(A, L, typ, Mod);
        ll Inv_L = qpow(L, Mod - 2, Mod);
        for(int i = 0; i != L; ++i) (A[i] *= Inv_L) %= Mod;
    }
    
    int n, m, L, P;
    ll a[N], b[N], c[N], d[N], ans[3][N];
    
    ll CRT(ll a1, ll a2, ll a3) {
        ll k1 = Mul((a2 - a1),qpow(Mod1, Mod2-2, Mod2), Mod2);
        ll a4 = a1 + k1*Mod1;
        ll k4 = Mul((a3 - a4) , qpow(M, Mod3-2, Mod3), Mod3);
        return (a4%P + Mul(k4 , M, P)) % P;
    }
    
    int main() {
    	scanf("%d%d%d",&n,&m,&P);
    	for(int i = 0; i <= n; ++i) scanf("%lld", &a[i]);
        for(int i = 0; i <= m; ++i) scanf("%lld", &b[i]);
        for(L = 1; L <= m+n; L <<= 1);  init(L);
        
        memcpy(c, a, sizeof(a)); memcpy(d, b, sizeof(a)); DFT(c, L, 1, Mod1); DFT(d, L, 1, Mod1);
        for(int i = 0; i <= L; ++i) ans[0][i] = c[i] * d[i] % Mod1;    
        memcpy(c, a, sizeof(a)); memcpy(d, b, sizeof(a)); DFT(c, L, 2, Mod2); DFT(d, L, 2, Mod2);
        for(int i = 0; i <= L; ++i) ans[1][i] = c[i] * d[i] % Mod2;
        memcpy(c, a, sizeof(a)); memcpy(d, b, sizeof(a)); DFT(c, L, 3, Mod3); DFT(d, L, 3, Mod3);
        for(int i = 0; i <= L; ++i) ans[2][i] = c[i] * d[i] % Mod3;
        IDFT(ans[0], L, 1, Mod1); IDFT(ans[1], L, 2, Mod2); IDFT(ans[2], L, 3, Mod3);
        
        for(int i = 0; i <= n+m; ++i) {
            printf("%lld ",CRT(ans[0][i], ans[1][i], ans[2][i]));
        }puts("");
        return 0;
    }
    

    多项式的各种操作

    多项式逆元

    题意:给定多项式(A(x)),请求出一个多项式(B(x)),满足(A(x)*B(x) ≡ 1 (mod~x^n))

    做法:

    1. (n = 1)时,(A(x) ≡ c (mod~x)),这样(A^{-1}(x) ≡ c^{-1}(mod~x))
    2. (n > 1)时,设(B(x) ≡ A^{-1}(x)(mod~x)),则有$$A(x)B(x) ≡ 1 (mod~x^n)$$
      假设在 (mod~x^{lceil frac{n}{2} ceil})意义下,(A(x))的逆元是(B(x))并且我们已经求出,那么有$$A(x)B'(x) ≡ 1 (mod~x^{lceil frac{n}{2} ceil})$$
      显然$$A(x)B(x) ≡ 1 ( mod~x^{lceil frac{n}{2} ceil})$$
      相减可得$$B(x)-B'(x) ≡ 0 (mod~x^{lceil frac{n}{2} ceil})$$
      两边平方得$$B^2(x) - 2B'(x)B(x) + B'^2(x) ≡ 0 (mod~x^n)$$
      两边同乘(A(x)),移项可得$$B(x) ≡ 2B'(x) - A(x)B'^2(x) (mod ~x^n)$$

    再利用(FFT)加速运算即可。

    Code[洛谷P4238]

    int n, L;
    ll a[N], b[N], tmp[N];
    void PloyInv(ll A[], ll B[], int n) {
        if(n == 1) {
            B[0] = qpow(A[0], Mod - 2, Mod); return;
        }
        PloyInv(A, B, (n + 1) >> 1);
        int p = 1; while(p < n<<1) p <<= 1;
        copy(A, A+n, tmp); fill(tmp+n, tmp+p, 0);
        DFT(tmp, p); DFT(B, p);
        for(int i = 0; i != p; ++i) {
            B[i] = (2 - tmp[i] * B[i] % Mod) * B[i] % Mod;
            if(B[i] < 0) B[i] += Mod;
        }
        IDFT(B, p);
        fill(B + n, B + p, 0);
    }
    int main() {
        scanf("%d",&n);
        for(int i = 0; i < n; ++i) scanf("%lld", &a[i]);
        for(L = 2; L <= 2*n-2; L <<= 1); init(L);
        PloyInv(a, b, n);
        for(int i = 0; i < n; ++i) printf("%lld ",b[i]); puts("");
        return 0;
    }
    

    多项式除法

    题意:给定一(n)次多项式(A(x))和一个(m(m leq n))次多项式(B(x)),要求求出(D(x), R(x))满足$$A(x) = D(x)B(x) + R(x)$$且(degD leq degA - degB = n-m, ~degR < m)

    做法:首先,想办法消除(R(x))的影响。首先我定义一个运算

    [A^R(x)= x^n A( frac{1}{x} ) ]

    可以发现,这个操作在进行系数反转。

    现在将原式的 (x) 全部替换为 (frac{1}{x}) ,然后两边同乘(x^n),可得

    [x^nA(frac{1}{x}) = x^{n-m}D(frac{1}{x})x^mB(frac{1}{x}) + x^{n-m+1}x^{m-1}R(frac{1}{x}) \ A^R(x) = D^R(x)B^R(x) + x^{n-m+1}R^R(x) ]

    观察这个式子,(D(x))反转后,次数不会高于(n-m),而(x^{n-m+1}R^R(x))这个部分最低次高于(n-m),把上式放到(mod ~x^{n-m+1})意义下,(R(x))就会被消除,现在式子变为$$A^R(x) ≡ D{R}(x)BR(x) (mod ~x^{n-m+1})$$这样只需要求(B^R(x))的逆元即可,(D(x))可以通过回带(R(x))求得。

    void PloyDiv(ll A[], ll B[], ll D[], ll R[], int n, int m) {
        static ll A0[N], B0[N];
        int p = 1, t = n - m + 1;
        while(p < t<<1) p <<= 1;
    
        fill(A0, A0 + p, 0);
        reverse_copy(B, B+m, A0);
        PloyInv(A0, B0, t);
        fill(B0+t, B0+p, 0);
        DFT(B0, p);
        reverse_copy(A, A+n, A0);
        fill(A0 + t, A0 + p, 0);
        DFT(A0, p);
        for(int i = 0; i != p; ++i)
            A0[i] = Mul(A0[i] , B0[i]);
        IDFT(A0, p);
        reverse(A0, A0+t); copy(A0, A0+t,D);
    
        p = 1; while(p < n) p <<= 1;
        fill(A0 + t, A0 + p, 0);
        DFT(A0, p);
        copy(B, B+m, B0);
        fill(B0 + m, B0 + p, 0);
        DFT(B0, p);
        for(int i = 0; i != p; ++i) A0[i] = Mul(A0[i] , B0[i]);
        IDFT(A0, p);
        for(int i = 0; i != m; ++i) R[i] = Dec(A[i] , A0[i]);
        fill(R+m, R+p, 0);
    }
    
    int main() {
        scanf("%d%d",&n,&m); ++n; ++m;
        for(int i = 0; i < n; ++i) scanf("%lld", &a[i]);
        for(int i = 0; i < m; ++i) scanf("%lld", &b[i]);
        int LL = max(max(n, m), (n - m + 1) << 1);
        L = 1; while(L < LL) L <<= 1; init(L);
        
        PloyDiv(a, b, D, R, n, m);
        
        for(int i = 0; i < n-m+1; ++i) printf("%lld ",D[i]); puts("");
        for(int i = 0; i < m-1; ++i) printf("%lld ",R[i]); puts("");
        return 0;
    }
    

    多项式牛顿迭代

    题意:已知函数(G(z)),求一个多项式 (F(z) ~mod ~z^n) ,满足$$G(F(z)) ≡ 0 ~(mod ~z^n)$$

    做法:类似多项式求逆的做法。
    当n = 1时,(G(F(z)) ≡ 0 (mod~ z))要单独求出。现在假设已经求出$$G(F_0(z)) ≡ 0(modz^{lceil frac{n}{2} ceil})$$考虑如何将答案,扩展至 (mod~z^n)下,把(G(F(z)))(F_0(z))泰勒展开$$G(F(z))=G(F_0(z))+frac{G'(F_0(z))}{1!}(F(z)-F_0(z)) + frac{G''(F_0(z))}{2!}(F(z)-F_0(z))^2 + ...$$因为(F(z))(F_0(z))的最后 (lceil frac{n}{2} ceil) 项相同,所以 ((F(z)-F_0(z))^2) 的最低非0项次数大于 (2lceil frac{n}{2} ceil) ,所以有$$G(F(z))≡G(F_0(z))+frac{G'(F_0(z))}{1!}(F(z)-F_0(z)) (modz^n)$$因为(G(F(z)) ≡ 0 (mod z^n)), 那么有$$F(z) ≡F_0(z) - frac{G(F_0(z))}{G'(F_0(z))} ~(mod ~z^n)$$


    多项式开方

    题意:给定一个 (n-1) 次多项式 (A(x)),求一个在 (mod~x^n) 意义下的多项式 (B(x)),满足(B^2(x) ≡ A(x)~(mod~x^n))

    做法:移向可得(B^2(x) - A(x) ≡ 0 ~(mod~x^n)) ,设 (G(B(x)) = B^2(x) - A(x)), 有(G'(B(x)) = 2B(x)),带入牛顿迭代的式子里,则有:

    [B(x) ≡ B_0(x) - frac{B_0^2(x)-A(x)}{2B_0(x)} ≡ frac{B_0^2(x)+A(x)}{2B_0(x)}~(mod~x^n) ]

    注意到当 (n = 1) 时,如果题目不保证 (a_0 = 1),可能需要计算二次剩余

    Code[洛谷P5205]

    void PloySqrt(ll A[], ll B[], int n) {
        static ll T[N];
        if(n == 1) {
            B[0] = 1; return;
        }
        int p = 1, hf =  (n+1) >> 1; while(p < n<<1) p <<= 1;
        PloySqrt(A, B, hf); fill(B + hf, B + p, 0);
        PloyInv(B, T, n); fill(T + n, T + p, 0);
        DFT(T, p);
    
        fill(B + hf, B + p, 0);
        DFT(B, p >> 1);
        for(int i = 0; i != (p >> 1); ++i)
            B[i] = (B[i]*B[i]) % Mod;
        IDFT(B, p >> 1);
        for(int i = 0; i != n; ++i)
            B[i] = (A[i] + B[i]) % Mod * Inv2 % Mod;
        DFT(B, p);
        for(int i = 0; i != p; ++i)
            B[i] = (B[i] * T[i]) % Mod;
        IDFT(B, p);
    }
    int main() {
        scanf("%d",&n);
        for(int i = 0; i < n; ++i) scanf("%lld", &a[i]);
        for(L = 2; L <= n << 1; L <<= 1); init(L);
        Inv2 = qpow(2,Mod-2,Mod);
        PloySqrt(a, b, n);
        for(int i = 0; i < n; ++i) printf("%lld ",b[i]); puts("");
        return 0;
    }
    

    多项式 ln 和多项式 exp


    多项式ln

    • 定义:对于多项式 (A(x) = sum_{i geq 1}a_ix^i) 多项式ln有$$ln(1-A(x)) = -sum_{i geq 1}frac{A^i(x)}{i}$$

    • 计算
      现在有(A(x) = 1 + sum_{i geq 1} a_ix^i)
      (ln(A(x))) 求导有$$(lnA(x))' = frac{A'(x)}{A(x)}$$ 只需计算 (A(x)) 的逆元,再完成求导和积分运算即可

    Code[洛谷4725]

    void PloyD(ll A[],ll B[], int n) {
        copy(A, A+n, B);
        for(int i = 0; i < n-1; ++i) {
            B[i] = 1ll*B[i+1]*(i+1) % Mod;
        } B[n-1] = 0;
    }
    void PloyF(ll A[], ll B[], int n) {
        copy(A, A+n, B);
        for(int i = n-1; i; --i) {
            B[i] = B[i-1] * Inv[i] % Mod;
        } B[0] = 0;
    }
    void PloyLn(ll A[], ll B[], int n) {
        static ll T[N];
        int p = 1; while(p < n << 1) p <<= 1;
        PloyInv(A, T, n); fill(T+n, T+p, 0); DFT(T, p);
        PloyD(A,B,n); fill(B+n, B+p, 0); DFT(B, p);
        for(int i = 0; i != p; ++i) T[i] = (B[i] * T[i]) % Mod;
        IDFT(T, p);
        PloyF(T, B, n); fill(B+n, B+p, 0);
    }
    int main() {
        scanf("%d",&n);
        for(int i = 0; i < n; ++i) scanf("%lld", &a[i]);
        for(L = 2; L <= n << 1; L <<= 1); init(L);
        Inv[1] = 1;
        for(ll i = 2; i <= L; i++)
            Inv[i] = (Mod - Mod / i) % Mod * Inv[Mod % i] % Mod;
        PloyLn(a, b, n);
        for(int i = 0; i < n; ++i) printf("%lld ",b[i]); puts("");
        return 0;
    }
    

    多项式exp

    • 定义:对于多项式 (A(x) = sum_{i geq 1}a_ix^i) 多项式exp有$$e^{A(x)} = sum_{i geq 0}frac{A^i(x)}{i!}$$

    • 计算

    (F(x) = e^{A(x)}) ,变形后得到 $$ln ~F(x) - A(x) = 0$$ 构造函数 $$G(F(x)) = ln~F(x) - A(x) $$, (G'(F(x)) = frac{1}{F(x)}),代入牛顿迭代$$F(x) ≡ F_0(x)(1 - ln~F_0(x) + A(x)) (modx^n)$$ 当 (n = 1)(F(x) = 1)

    Code[洛谷P4726]

    void Ployexp(ll A[], ll B[], int n) {
        static ll T[N];
        if(n == 1) {
            B[0] = 1; return;
        }
        int p = 1, hf = (n + 1) >> 1; while(p < n<<1) p <<= 1;
        Ployexp(A, B, hf); fill(B + hf, B + p, 0);
        PloyLn(B, T, n);
        for(int i = 0; i != n; ++i) T[i] = (A[i] - T[i] + Mod) % Mod;
        ++ T[0]; if(T[0] >= Mod) T[0] -= Mod;
        DFT(T, p); DFT(B, p);
        for(int i = 0; i != p; ++i) B[i] = B[i] * T[i] % Mod;
        IDFT(B, p);
    }
    int main() {
        scanf("%d",&n);
        for(int i = 0; i < n; ++i) scanf("%lld", &a[i]);
        for(L = 2; L <= n << 1; L <<= 1); init(L);
        Inv[1] = 1;
        for(ll i = 2; i <= L; i++)
            Inv[i] = (Mod - Mod / i) % Mod * Inv[Mod % i] % Mod;
        Ployexp(a, b, n);
        for(int i = 0; i < n; ++i) printf("%lld ",b[i]); puts("");
        return 0;
    }
    

    多项式任意次幂

    题意:给出多项式(A(x)),求(A^k(x),k in Q)

    做法:

    1. k为整数时可以快速幂
    2. 可以利用(A^k(x) = e^{k~ln~A(x)}) 计算

    多项式三角函数

    题意:给出多项式(A(x)),求(sin A(x)~(mod~x^n)), (cos A(x)~(mod~x^n))

    做法:利用欧拉公式:$$e^{iA(x)} = cosA(x) + i sinA(x)$$与常规的多项式exp基本一致,只能用FFT,运算全部为复数,注意边界


    多项式多点求值


    多项式快速插值


    多项式复合逆


    多项式相关题目及应用

    分治FFT

    题意:给定长度为(n-1)的序列(g[1],...,g[n-1])计算(f[0],...,f[n-1])

    [f[i] = sum_{j=1}^i f[i-j]g[j] = sum_{j=0}^{i-1} f[j]g[i-j] ]

    (f[0] = 1)

    做法:CDQ+FFT。考虑已经求出([l,mid])(f[i]),计算他们堆([mid+1,r])的贡献$w[x], x in [mid+1,r] $ 有

    [ w[x] = sum_{i=l}^{mid} f[i]g[x-i] ightarrow w[x] = sum_{i=l}^{x-1} f[i]g[x-i] ]

    (设([mid+1,r])(f[i]=0)
    (k = i-l),

    [w[x] = sum_{k=0}^{x-l-1} f[k+l]g[x-k-l] = sum_{k=0}^{x-l-1} a[k]b[x-l-1-k] \ ]

    我们凑出对应的(a[i], b[i])

    [a[k] = f[k+l]\ b[x-l-1-k] = g[x-l-k] ightarrow \ b[x-l-1+k] = g[x-l+k] ightarrow \ b[k-1] = g[k] ightarrow b[k] = g[k+1] ]

    Code[洛谷P4721 分治FFT]

    int n, L;
    ll Ans[N], g[N], a[N], b[N];
    void solve(int l, int r) {
        if(l == r) return;
        int mid = (l + r) >> 1;
        solve(l,mid);
        int len = r - l + 1, p = 1;
        for(;p < len; p <<= 1); // 我们最多需要r-l+1前项!!!
        for(int i = 0; i < p; ++ i) a[i] = b[i] = 0;
        for(int i = l; i <= mid; ++ i) a[i-l] = Ans[i];
        for(int i = 0; i <= r-l; ++ i) b[i] = g[i+1];
        DFT(a,p), DFT(b,p);
        for(int i = 0; i < p; ++i) a[i] = a[i]*b[i]%Mod;
        IDFT(a,p);
        rep(i,mid+1,r) {
            Ans[i] += a[i-l-1]; 
            if(Ans[i] >= Mod) Ans[i] -= Mod;
        }
        solve(mid+1,r);
    }
    int main() {
        read(n); rep(i,1,n-1) read(g[i]);
        for(L = 2; L <= n << 1; L <<= 1); init(L);
        Ans[0] = 1;
        solve(0,n-1);
        rep(i,0,n-1) write(Ans[i]), putchar(' '); putchar('
    ');
    }
    

    参考资料

    1. DFT与圆周卷积
    2. 任意模数FFT
    3. 从DTFT到DFS,从DFS到DFT,从DFT到FFT,从一维到二维
    4. 一位神犇的Blog
    5. 另一位神犇的Blog
  • 相关阅读:
    ASP.NET性能优化篇(转载)
    Apache相关
    UVa11292 The Dragon of Loowater
    POJ2653 Pickup sticks
    POJ2155 Matrix
    POJ3009 Curling 2.0
    POJ1066 Treasure Hunt
    UVa11729 Commando War
    Ubuntu下解决压缩文件的文件名乱码问题
    HDU3415 Max Sum of MaxKsubsequence
  • 原文地址:https://www.cnblogs.com/RRRR-wys/p/10342867.html
Copyright © 2020-2023  润新知