zoukankan      html  css  js  c++  java
  • 多项式不全家桶

    多项式不全家桶(暂不更新)

    前置知识需要FFT和NTT。

    多项式求逆

    板子题链接
    求逆指的是给定一个多项式 (F(x)),你需要求出一个多项式 (G(x)),使其满足 (F(x) * G(x)equiv1pmod {x^n})

    考虑为什么要在 (pmod {x^n}) 情况下来求解,其实多项式可以近似理解为级数,显然级数的项数是无穷的,所以我们用 (pmod {x^n}) 来使得这个多项式的项数有穷。(大致理解就行)

    首先大致的思路是分治解决问题。

    显然的如果多项式只有 (1) 项,那么解就是其逆元。

    我们考虑从原式在 mod (x^{lceilfrac n 2 ceil}) 的推到 mod (x^ n) 的情 况:

    [F(x)* G_1(x)equiv1pmod {x^{lceilfrac n 2 ceil}} ]

    其中 (G_1(x)) 指在 (pmod {x^{lceilfrac n 2 ceil}}) 的情况下的解,(G(x)) 指在 (pmod {x^n}) 的情况下的解。

    那么也有:

    [F(x)* G(x)equiv1pmod {x^{lceilfrac n 2 ceil}} ]

    两个式子相减可以得到:

    [G_1(x)-G(x)equiv0pmod{x^{lceilfrac n 2 ceil}} ]

    平方得 :

    [G_1^{2}(x)-2G_1(x)G(x)+G^2(x)equiv0pmod{x^n} ]

    两边同乘多项式 (F(x)) 得:

    [F(x) * G_1^{2}(x) - 2G_1(x) + G(x)equiv0pmod{x^n} ]

    移项显然得:

    [G(x)equiv(2-F(x)G_1(x))G_1(x)pmod{x^n} ]

    然后你发现右边的东西已经可以计算了。

    复杂度非常显然的 (O(nlogn))

    但是需要注意几个问题,我会在代码里写出来。

    代码
    #include <bits/stdc++.h>
    using namespace std;
    const int N = 4e5 + 5, mod = 998244353, g1 = 3, g2 = (mod + 1) / 3;
    int n, a[N], b[N], d[N], A[N], B[N];
    inline int qpow(int a, int b, int res = 1) {
        for(; b; b >>= 1, a = 1ll * a * a % mod) 
            if(b & 1) res = 1ll * res * a % mod;
        return res;
    }
    void NTT(int *a, int flag, int len) {
        for(register int i = 0; i < len; ++i) if(i < d[i]) swap(a[i], a[d[i]]);
        for(register int l = 2, m = 1; l <= len; l <<= 1, m <<= 1) {
            int w1 = qpow(flag == 1 ? g1 : g2, (mod - 1) / l);
            for(register int *p = a; p != a + len; p += l) {
                int w = 1;
                for(register int i = 0; i < m; ++i, w = 1ll*w*w1%mod) {
                    int t = 1ll * w * p[i + m] % mod;
                    p[i + m] = (p[i] - t + mod) % mod, p[i] = (p[i] + t) % mod;
                }
            }
        }
    }
    void niyuan(int *a, int *b, int n) {
        memset(A, 0, n * 16), memset(B, 0, n * 16);//如果需要多次计算的话,别忘了清零
        b[0] = qpow(a[0], mod - 2);
        for(register int l = 2, m = 1; l <= n << 2; l <<= 1, m <<= 1) {//注意我们需要跑到 l < n * 4,即4倍
            for(register int i = 0; i < m; ++i) A[i] = a[i], B[i] = b[i];
            for(register int i = 0; i < l; ++i) d[i] = (d[i>>1] >> 1) | ((i & 1) ? m : 0);
            NTT(A, 1, l), NTT(B, 1, l);
            for(register int i = 0; i < l; ++i) b[i] = ((2ll - 1ll*A[i]*B[i]%mod) * B[i] % mod + mod) % mod;
            NTT(b, -1, l);
            int ni = qpow(l, mod - 2);
            for(register int i = 0; i < l; ++i) b[i] = 1ll * b[i] * ni % mod;
            for(register int i = m; i < l; ++i) b[i] = 0;//为什么每次需要清零,主要是由推导的式子决定的。
        }
    }
    int main() {
        scanf("%d", &n);
        for(register int i = 0; i < n; ++i) scanf("%d", &a[i]);
        niyuan(a, b, n);
        for(register int i = 0; i < n; ++i) printf("%d ", b[i]);
        return 0;
    }
    

    多项式开根

    板子题链接

    多项式开根是已知多项式 (A(x)) ,在 (pmod {x^n}) 的条件下,求得一个多项式 (B(x)),使得 (B^2(x) = A(x) equiv1pmod {x^n})

    我们依旧采用分治的思想进行解决。
    一般我们习惯设 (B(x)) 的常数项为 1,这样主要是为了方便计算。

    我们考虑从原式在 mod (x^{lceilfrac n 2 ceil}) 的推到 mod (x^ n) 的情 况:

    [B^2(x)equiv A(x) (mod x^{n}) ]

    设多项式 (B'^2(x)) 为在 ((mod x^{frac{n}{2}})) 的情况下的解,
    那么:

    [B'^2(x)equiv A(x) (mod x^{frac{n}{2}}) ]

    然后两个式子相减得到 :

    [B(x)-B'(x)equiv 0(mod x^{frac{n}{2}}) ]

    然后平方得 :

    [A(x)-2B(x)B'(x)+B'^2(x)equiv 0(mod x^n) ]

    最后移项得 :

    [B(x)equiv dfrac{A(x)+B'^2(x)}{2B'(x)} (mod x^n) ]

    非常显然的右边的式子可以用求逆和乘法解决掉。

    复杂度 (O(nlogn))

    代码
    #include <bits/stdc++.h>
    using namespace std;
    const int N = 4e5 + 5, mod = 998244353, g1 = 3, g2 = (mod + 1) / 3, ni2 = (mod + 1) / 2;
    int n, a[N], b[N], d[N], A[N], B[N], inv[N], c[N];
    inline int qpow(int a, int b, int res = 1) {
        for(; b; b >>= 1, a = 1ll * a * a % mod) 
            if(b & 1) res = 1ll * res * a % mod;
        return res;
    }
    void NTT(int *a, int flag, int len) {
        for(register int i = 0; i < len; ++i) if(i < d[i]) swap(a[i], a[d[i]]);
        for(register int l = 2, m = 1; l <= len; l <<= 1, m <<= 1) {
            int w1 = qpow(flag == 1 ? g1 : g2, (mod - 1) / l);
            for(register int *p = a; p != a + len; p += l) {
                int w = 1;
                for(register int i = 0; i < m; ++i, w = 1ll*w*w1%mod) {
                    int t = 1ll * w * p[i + m] % mod;
                    p[i + m] = (p[i] - t + mod) % mod, p[i] = (p[i] + t) % mod;
                }
            }
        }
    }
    void niyuan(int *a, int *b, int n) {
        memset(A, 0, n * 16), memset(B, 0, n * 16);
        b[0] = qpow(a[0], mod - 2);
        for(register int l = 2, m = 1; l < n << 2; l <<= 1, m <<= 1) {
            for(register int i = 0; i < m; ++i) A[i] = a[i], B[i] = b[i];
            for(register int i = 0; i < l; ++i) d[i] = (d[i>>1] >> 1) | ((i & 1) ? m : 0);
            NTT(A, 1, l), NTT(B, 1, l);
            for(register int i = 0; i < l; ++i) b[i] = ((2ll - 1ll*A[i]*B[i]%mod) * B[i] % mod + mod) % mod;
            NTT(b, -1, l);
            int ni = qpow(l, mod - 2);
            for(register int i = 0; i < l; ++i) b[i] = 1ll * b[i] * ni % mod;
            for(register int i = m; i < l; ++i) b[i] = 0;
        }
    }
    void kaigen(int *a, int *b, int n) {
        memset(c, 0, n * 16), memset(inv, 0, n * 16);
        b[0] = 1;//默认
        for(register int l = 2, m = 1; l < n << 2; l <<= 1, m <<= 1) {
            for(register int i = 0; i < m; ++i) c[i] = a[i];
            niyuan(b, inv, m);
            for(register int i = 0; i < l; ++i) d[i] = (d[i>>1] >> 1) | ((i&1) ? m : 0);
            NTT(c, 1, l), NTT(inv, 1, l);
            for(register int i = 0; i < l; ++i) c[i] = 1ll * c[i] * inv[i] % mod;
            NTT(c, -1, l);
            int ni = qpow(l, mod - 2);
            for(register int i = 0; i < l; ++i) c[i] = 1ll * c[i] * ni % mod;
            for(register int i = 0; i < m; ++i) b[i] = 1ll * ni2 * ((b[i] + c[i])%mod) % mod;
        }
    }
    int main() {
        scanf("%d", &n);
        for(register int i = 0; i < n; ++i) scanf("%d", &a[i]);
        kaigen(a, b, n);
        for(register int i = 0; i < n; ++i) printf("%d ", b[i]);
        return 0;
    }
    

    多项式求导与积分

    求导 : 多项式求导的话一定会使它的次数降低。
    多项式求导可以近似理解为对于多项式的每一项分别求导,幂函数的求导公式就是 (f'(x^n) = n* x^{n-1})
    所以多项式求导的代码也就大致出来了:

    void qiudao(int *a, int *b, int n) {
         for(register int i = 1; i < n; ++i) b[i-1] = 1ll * i * a[i] % mod; b[n-1] = 0;
    }
    

    积分
    反正我只知道积分可以大概理解为是求导的逆运算,而且我也不会算,所以直接告诉你代码吧(其实可以从求逆的角度理解):

    void jifen(int *a, int *b, int n) {
         for(register int i = 1; i < n; ++i) b[i] = 1ll * a[i-1] * qpow(i, mod - 2) % mod; b[0] = 0;
    }
    

    注意积分不能把常数项求出来,所以我们默认为0。

    多项式对数函数(多项式 ln)

    已知多项式 (A(x)),我们需要求出一个多项式使其满足 (B(x) = lnA(x))
    下面开始推式子,首先我们对两边分别求导:

    [B(x) = lnA(x) ]

    [B'(x) = frac{A'(x)}{A(x)} ]

    右边的式子我们可以理解为复合函数的求导展开。
    然后我们发现右边的式子很可做。
    然后我们就可以得到 (B'(x)),最后在积分一下就得到了 (B(x))

    代码
    #include <bits/stdc++.h>
    using namespace std;
    const int N = 4e5 + 5, mod = 998244353, g1 = 3, g2 = (mod + 1) / 3;
    int n, a[N], b[N], A[N], B[N], c[N], inv[N], d[N];
    int qpow(int a, int b, int res = 1) {
        for(; b; b >>= 1, a = 1ll * a * a % mod) 
            if(b & 1) res = 1ll * res * a % mod;
        return res;
    }
    void NTT(int *a, int flag, int len) {
        for(register int i = 0; i < len; ++i) if(i < d[i]) swap(a[i], a[d[i]]);
        for(register int l = 2, m = 1; l <= len; l <<= 1, m <<= 1) {
            int w1 = qpow(flag == 1 ? g1 : g2, (mod-1) / l);
            for(register int *p = a; p != a + len; p += l) {
                int w = 1;
                for(register int i = 0; i < m; ++i, w = 1ll*w*w1%mod) {
                    int t = 1ll * w * p[i + m] % mod;
                    p[i + m] = (p[i] - t + mod) % mod;
                    p[i] = (p[i] + t) % mod;
                }
            }
        }
    }
    void qiuni(int *a, int *b, int n) {
        memset(A, 0, n * 16), memset(B, 0, n * 16);
        b[0] = qpow(a[0], mod - 2);
        for(register int l = 2, m = 1; l < n << 2; l <<= 1, m <<= 1) {
            for(register int i = 0; i < m; ++i) A[i] = a[i], B[i] = b[i];
            for(register int i = 0; i < l; ++i) d[i] = (d[i>>1] >> 1) | ((i&1) ? m : 0);
            NTT(A, 1, l), NTT(B, 1, l);
            for(register int i = 0; i < l; ++i) b[i] = ((2ll - 1ll*A[i]*B[i]%mod) * B[i] % mod + mod) % mod;
            NTT(b, -1, l);
            int ni = qpow(l, mod - 2);
            for(register int i = 0; i < l; ++i) b[i] = 1ll * b[i] * ni % mod;
            for(register int i = m; i < l; ++i) b[i] = 0;
        }
    }
    void qiudao(int *a, int *b, int n) {
        for(register int i = 1; i < n; ++i) b[i-1] = 1ll * i * a[i] % mod; b[n-1] = 0;
    }
    void jifen(int *a, int *b, int n) {
        for(register int i = 1; i < n; ++i) b[i] = 1ll * a[i-1] * qpow(i, mod - 2) % mod; b[0] = 0;
    }
    void loyin(int *a, int *b, int n) {
        memset(c, 0, n * 16), memset(inv, 0, n * 16);
        qiudao(a, c, n), qiuni(a, inv, n);
        int len = 1;
        while(len <= n*2) len <<= 1;
        for(register int i = 0; i < len; ++i) d[i] = (d[i>>1] >> 1) | ((i&1) ? len >> 1 : 0);
        NTT(c, 1, len), NTT(inv, 1, len);
        for(register int i = 0; i < len; ++i) c[i] = 1ll * c[i] * inv[i] % mod;
        NTT(c, -1, len);
        int ni = qpow(len, mod - 2);
        for(register int i = 0; i < len; ++i) c[i] = 1ll * c[i] * ni % mod;
        jifen(c, b, n);
    }
    int main() {
        scanf("%d", &n);
        for(register int i = 0; i < n; ++i) scanf("%d", &a[i]);
        loyin(a, b, n);
        for(register int i = 0; i < n; ++i) printf("%d ", b[i]); 
        return 0;
    }
    

    分治FFT

    板子题链接

    分治FFT主要是解决这样一类自己会转移到自己的问题,比如给定多项式 (A(x)),让你求出 (B_i = sum^{i}_{j = 1} B_{i-j} * A_j),边界 (B_0 = 1)

    这就是分治FFT主要解决的问题,主要的思想是利用CDQ分治的思想,每次计算左半变对于右半边的共享,这样复杂度是 (O(nlog^2n))的。

    其实其主要思想就是cdq分治,每次做一个卷积(乘法),好像不知道怎么解释啊,自己试一下,当然主要的细节问题是边界问题,就是每次做多项式乘法的边界问题,这个在代码里标注。

    代码
    #include <bits/stdc++.h>
    using namespace std;
    const int N = 4e5 + 5, mod = 998244353, g1 = 3, g2 = (mod + 1) / 3;
    int n, a[N], b[N], d[N], cdq1[N], cdq2[N];
    int qpow(int a, int b, int res = 1) {
        for(; b; b >>= 1, a = 1ll * a * a % mod) 
            if(b & 1) res = 1ll * res * a % mod;
        return res;
    }
    void NTT(int *a, int flag, int len) {
        for(register int i = 0; i < len; ++i) if(i < d[i]) swap(a[i], a[d[i]]);
        for(register int l = 2, m = 1; l <= len; l <<= 1, m <<= 1) {
            int w1 = qpow(flag == 1 ? g1 : g2, (mod-1) / l);
            for(register int *p = a; p != a + len; p += l) {
                int w = 1;
                for(register int i = 0; i < m; ++i, w = 1ll*w*w1%mod) {
                    int t = 1ll * w * p[i + m] % mod;
                    p[i + m] = (p[i] - t + mod) % mod;
                    p[i] = (p[i] + t) % mod;
                }
            }
        }
    }
    void fenzhi(int l, int r) {
        if(l == r) return;
        int mid = (l + r) / 2, len = r - l + 1;
        fenzhi(l, mid);
        int *A = cdq1, *B = cdq2;
        memset(A, 0, len * 8), memset(B, 0, len * 8);//注意每次清空数组
        for(register int i = l; i <= mid; ++i) A[i-l] = b[i];
        for(register int i = 0; i < r-l+1; ++i) B[i] = a[i+1];
        int n = 1;
        while(n <= len) n <<= 1;//之所以不用到 n <= len << 1,是因为这是一个循环卷积。
        for(register int i = 0; i < n; ++i) d[i] = (d[i>>1] >> 1) | ((i&1) ? n >> 1 : 0);
        NTT(A, 1, n), NTT(B, 1, n);
        for(register int i = 0; i < n; ++i) A[i] = 1ll * A[i] * B[i] % mod;
        NTT(A, -1, n);
        int ni = qpow(n, mod - 2);
        for(register int i = 0; i < n; ++i) A[i] = 1ll * A[i] * ni % mod;
        for(register int i = mid + 1; i <= r; ++i) b[i] = (b[i] + A[i-l-1]) % mod;//注意取值与边界
        fenzhi(mid + 1, r);
    }
    int main() {
        scanf("%d", &n), b[0] = 1;
        for(register int i = 1; i < n; ++i) scanf("%d", &a[i]);
        fenzhi(0, n-1);
        for(register int i = 0; i < n; ++i) printf("%d ", b[i]); 
        return 0;
    }
    

    多项式指数函数(多项式 exp)

    我为什么要先讲一个分治FFT,主要是因为我的exp是写的分治FFT,复杂度是 (O(nlog^2n)),虽然复杂度比正解的 (O(nlogn)) 要慢一点,但是代码好写,比正解的代码码量少一半,而且因为正解的常数过于巨大,如果卡常不好的话,我这个更快,洛谷上的提交记录显示,我的代码是比大多数人的正解都快的,所以我选择写分治FFT。

    显然 exp 是给你一个多项式 (A(x)),让你求得一个多项式 (B(x)),使其满足 (e^{A(x)} = B(x))

    显然我们接着推式子:
    我们依旧对于左右两边求导:

    [B'(x) = A'(x) * e^{A(x)} ]

    替换得到 :

    [B'(x) = A'(x) * B(x) ]

    我们直接积分求解得:

    [B(x) = int A'(x) * B(x) ]

    然后我们可能觉得右边的式子没有什么思路,但是我们发现这个式子与分治FFT非常像,只不过是带了一个积分,但是多项式的积分公式也很好求,所以我们可以选择采用分治FFT的思路,注释在代码里。

    代码
    #include <bits/stdc++.h>
    using namespace std;
    const int N = 4e5 + 5, mod = 998244353, g1 = 3, g2 = (mod + 1) / 3;
    int n, a[N], b[N], d[N], cdq1[N], cdq2[N];
    int qpow(int a, int b, int res = 1) {
        for(; b; b >>= 1, a = 1ll * a * a % mod) 
            if(b & 1) res = 1ll * res * a % mod;
        return res;
    }
    void NTT(int *a, int flag, int len) {
        for(register int i = 0; i < len; ++i) if(i < d[i]) swap(a[i], a[d[i]]);
        for(register int l = 2, m = 1; l <= len; l <<= 1, m <<= 1) {
            int w1 = qpow(flag == 1 ? g1 : g2, (mod-1) / l);
            for(register int *p = a; p != a + len; p += l) {
                int w = 1;
                for(register int i = 0; i < m; ++i, w = 1ll*w*w1%mod) {
                    int t = 1ll * w * p[i + m] % mod;
                    p[i + m] = (p[i] - t + mod) % mod;
                    p[i] = (p[i] + t) % mod;
                }
            }
        }
    }
    void fenzhi(int l, int r) {
        if(l == r) {
            if(l) b[l] = 1ll * b[l] * qpow(l, mod - 2) % mod;//积分公式,乘一个逆元
            else b[l] = 1;
            return;
        }
        int mid = (l + r) / 2, len = r - l + 1;
        fenzhi(l, mid);
        int *A = cdq1, *B = cdq2;
        int n = 1;
        while(n < len) n <<= 1;
        memset(A, 0, n * 4), memset(B, 0, n * 4);//其他的与分治FFT无任何区别
        for(register int i = 0; i < n; ++i) d[i] = (d[i>>1] >> 1) | ((i&1) ? n >> 1 : 0);
        for(register int i = l; i <= mid; ++i) A[i-l] = b[i];
        for(register int i = 0; i < len; ++i) B[i] = a[i];
        NTT(A, 1, n), NTT(B, 1, n);
        for(register int i = 0; i < n; ++i) A[i] = 1ll * A[i] * B[i] % mod;
        NTT(A, -1, n);
        int ni = qpow(n, mod - 2);
        for(register int i = 0; i < n; ++i) A[i] = 1ll * A[i] * ni % mod;
        for(register int i = mid+1; i <= r; ++i) b[i] = (b[i] + A[i-l-1]) % mod;
        fenzhi(mid + 1, r);
    }
    int main() {
        scanf("%d", &n);
        for(register int i = 0; i < n; ++i) scanf("%d", &a[i]);
        for(register int i = 1; i < n; ++i) a[i-1] = 1ll * a[i] * i % mod; a[n-1] = 0;//求导
        fenzhi(0, n-1);
        for(register int i = 0; i < n; ++i) printf("%d ", b[i]);
        return 0;
    }
    
  • 相关阅读:
    快速整理sql表结构到wiki
    mac subline批量处理行
    iphone7忘记手机屏幕密码
    docker 常用命令
    初窥响应式布局
    用jquery写的一个图片轮播插件
    javascript中的对象和创建对象的主要模式
    用户注册界面(带js特效)
    用javascript实现简易留言板
    用javascript实现的购物车实例
  • 原文地址:https://www.cnblogs.com/longdie/p/duoxiangshi.html
Copyright © 2011-2022 走看看