【Learning】多项式的一些东西

FFT

坑

---

NTT

将\(FFT\)中的单位复数根改成原根即可。

卡常版NTT模版

struct Mul
{   int Len;

    int wn[N], Lim;

    int rev[N];

    inline void getReverse(int * a)
    {
        static int rev[N];
        rev[0] = 0;
        for (register int i = 0; i < Len; i++)
        {
            rev[i] = (rev[i>>1] >> 1) | ((i&1) ? (Len >> 1) : 0);
            if (i < rev[i]) swap(a[i], a[rev[i]]);
        }
    }
 
    void NTT(int * a, int opt)
    {
        getReverse(a);
        
        for (register int i = 0; i < Len; i += 2)
        {
            int x = a[i], y = 1LL * wn[0] * a[i + 1] % mod;
            a[i] = Plus(x, y), a[i + 1] = Minus(x, y);
        }
        for (register int i = 0; i < Len; i += 4)
        {
            int x1 = a[i], y1 = 1LL * wn[0] * a[i + 2] % mod;
            int x2 = a[i + 1], y2 = 1LL * wn[Lim / 2] * a[i + 3] % mod;
            a[i] = Plus(x1, y1), a[i + 2] = Minus(x1, y1);
            a[i + 1] = Plus(x2, y2), a[i + 3] = Minus(x2, y2);
        }
        
        for (register int i = 8; i <= Len; i <<= 1)
        {
            int M = i >> 1;
            for (register int j = 0; j < Len; j += i)
            {
                for (register int k = 0; k < M; k += 4)
                {
                    int x1 = a[j + k], y1 = 1LL * wn[Lim / M * k] * a[j + k + M] % mod;
                    int x2 = a[j + k + 1], y2 = 1LL * wn[Lim / M * (k+1)] * a[j + k + M + 1] % mod;
                    int x3 = a[j + k + 2], y3 = 1LL * wn[Lim / M * (k+2)] * a[j + k + M + 2] % mod;
                    int x4 = a[j + k + 3], y4 = 1LL * wn[Lim / M * (k+3)] * a[j + k + M + 3] % mod;
                    a[j + k] = Plus(x1, y1), a[j + k + M] = Minus(x1, y1);
                    a[j + k + 1] = Plus(x2, y2), a[j + k + M + 1] = Minus(x2, y2);
                    a[j + k + 2] = Plus(x3, y3), a[j + k + M + 2] = Minus(x3, y3);
                    a[j + k + 3] = Plus(x4, y4), a[j + k + M + 3] = Minus(x4, y4);
                }
            }
        }
        
        if (opt == -1)
        {
            for (register int i = 0; i < Len; i++)
                a[i] = 1LL * a[i] * inv[Len] % mod;
            for(register int i = 1; i < Len; i++)
                if(i < Len-i) swap(a[i], a[Len-i]);
        }
    }

    inline void init()
    {
        Lim = Len;
        wn[0] = 1; wn[1] = power(7, (mod-1) / (Len << 1), mod);
        for (int i = 2; i < (Len << 1); i++) wn[i] = 1LL * wn[i-1] * wn[1] % mod;
    }

    inline int getLen(int l)
    {
        for (Len = 1; Len <= l; Len <<= 1);
        return Len;
    }
};

---

多项式求逆

也就是求
\[ A(x)B(x) = 1 \pmod{x^n}\]

假设我们已知
\[ A(x)G(x) = 1 \pmod{x^{\lceil{n \over 2}\rceil}} \]

两式相减得
\[ A(x)(B(x) - G(x)) = 0 \pmod{x^{\lceil{n \over 2}\rceil}}\]

\[ B(x) - G(x) = 0 \pmod{x^{\lceil{n \over 2}\rceil}}\]

\[ B^2(x) + G^2(x) - 2B(x)G(x) = 0 \pmod{x^n}\]

两边同时乘上\(A(x)\)

\[ A(x)B^2(x) + A(x)G^2(x) = 2A(x)B(x)G(x) \]

因为\(A(x)B(x) = 1\), 所以可以得到

\[ B(x) + A(x)G^2(x) = 2G(x) \]

\[ B(x) = 2G(x) - A(x)G^2(x) \]

就这样解决了

时间复杂度为\(T(n) = T(\frac{n}{2}) +O(n \log n)\)

代码

inline void getInv(int * A, int * B, int len)
{
    static int tmp1[N], tmp2[N];
    
    B[0] = power(A[0], mod-2, mod);
    
    for (register int k = 2; k <= len; k <<= 1)
    {
        int Len = k << 1;
        
        cpy(tmp1, A, k, Len);
        cpy(tmp2, B, (k >> 1), Len);
        
        Calc.Len = Len;
        
        Calc.NTT(tmp1, 1);
        Calc.NTT(tmp2, 1);
        for (int i = 0; i < Len; i++)
            tmp1[i] = Minus(Plus(tmp2[i], tmp2[i]), 1LL * tmp1[i] * tmp2[i] % mod * tmp2[i] % mod);
        Calc.NTT(tmp1, -1);
        
        cpy(B, tmp1, k, Len);
    }
    
    return;
}

---

多项式开根

也就是求

\[ B^2(x) = A(x) \pmod{x^n} \]

假设我们已知
\[ G^2(x) = A(x) \pmod{x^{\lceil{n \over 2}\rceil}} \]

两式相减得
\[ B^2(x) - G^2(x) = 0 \pmod{x^{\lceil{n \over 2}\rceil}} \]

然后平方一下
\[ B^4(x) + G^4(x) - 2B^2(x)G^2(x) = 0 \pmod{x^n} \]

\[ B^4(x) + G^4(x) = 2B^2(x)G^2(x) \pmod{x^n} \]

配一下方
\[ B^4(x) + G^4(x) + 2B^2(x)G^2(x) = 4B^2(x)G^2(x) \pmod{x^n} \]

\[ (B^2(x) + G^2(x))^2 = (2B(x)G(x))^2 \pmod{x^n} \]

\[ B^2(x) + G^2(x) = 2B(x)G(x) \pmod{x^n} \]

因为\(A(x)B(x) = 1\), 所以可以得到

\[ B(x) = {A(x) + G^2(x) \over {2G(x)}} \]

为了方便实现，我们常把它化成

\[ B(x) = {A(x) \over 2G(x)} + {G(x) \over 2} \]

然后就解决了

时间复杂度为\(T(n) = T(\frac{n}{2}) + O(n \log n) = O(n \log n)\)

---

多项式求导

模拟即可

代码

inline void getDeri(int * a, int len)
{   for (register int i = 0; i < len; i++)
        a[i] = 1LL * a[i+1] * (i+1) % mod;
}

---

多项式积分

模拟即可

代码

inline void getInte(int * a, int len)
{   for (int i = len-1; i >= 1; i--)
        a[i] = 1LL * a[i-1] * inv[i] % mod;
    a[0] = 0;
}

---

多项式求\(\ln\)

已知多项式\(A(x)\), 求\(B(x) = \ln (A(x))\)

根据链式法则， 我们可以得到

\[B'(x) = \frac{\mathbb{d}(\ln(A(x)))}{\mathbb{d}A(x)} \frac{\mathbb{d}A(x)}{\mathbb{d}x} = \frac{A'(x)}{A(x)}\]

所以

\[B(x) = \int \frac{A'(x)}{A(x)} \mathbb{d}x\]

代码

void getLn(int * A, int len)
{
    static int tmp1[N], tmp2[N], tmp3[N];
    
    int Len = len << 1;
    
    cpy(tmp1, A, len, Len);
    cpy(tmp2, A, len, Len);
    
    getDeri(tmp1, len);
    getInv(tmp2, tmp3, len);
    
    Calc.Len = Len;
    
    Calc.NTT(tmp1, 1);
    Calc.NTT(tmp3, 1);
    for (int i = 0; i < Len; i++)
        tmp1[i] = 1LL * tmp1[i] * tmp3[i] % mod;
    Calc.NTT(tmp1, -1);
    
    memset(tmp1 + len, 0, 4 * (Len - len));
    getInte(tmp1, len);
    
    cpy(A, tmp1, len, Len);
}

时间复杂度还是\(T(n) = T(\frac{n}{2}) + O(n \log n) = O(n \log n)\)

---

牛顿迭代

我们要求\(f(x) = 0\)的根

那么可以使用泰勒公式来近似\(f(x) = 0\)的根

我们设当前的近似值是\(x_{n}\)， 我们想要得到的近似值是\(x_{n+1}\)

截取泰勒公式的线性部分： \[f(x_{n}) + f'(x_{n})(x_{n+1}-x_n) = 0\]
解方程得： \[x_{n+1} = x_n - \frac{f(x_{n})}{f'(x_{n})}\]

我们可以用这种方法来解多项式方程。

---

多项式\(\exp\)

我们要求的是 \(exp(A(x))\), 这相当于解一个方程： \(g(exp(A(x))) = \ln(exp(A(x))) - A(x) = 0\)

可以直接套用牛顿迭代法求解。

转载于:https://www.cnblogs.com/2016gdgzoi509/p/9269091.html