多项式——多项式除法

最新推荐文章于 2022-08-06 20:18:03 发布

爱寂寞的时光

最新推荐文章于 2022-08-06 20:18:03 发布

阅读量801

点赞数

分类专栏：算法离散数学 # 组合数学、概率论、博弈论与数理统计等基础数学文章标签：算法

本文链接：https://blog.csdn.net/jiahonghao2002/article/details/126090079

版权

组合数学、概率论、博弈论与数理统计等基础数学同时被 3 个专栏收录

40 篇文章 2 订阅

订阅专栏

离散数学

17 篇文章 2 订阅

订阅专栏

算法

13 篇文章 1 订阅

订阅专栏

多项式——多项式除法

给定一个 $n$ 次多项式 $A (x)$ 和一个 $m$ 次多项式 $B (x)$ ，请求出多项式 $D (x)$ , $R (x)$ ，满足以下条件：

$D (x)$ 次数为 $n - m$ ， $R (x)$ 次数小于 $m$
$A (x) = D (x) B (x) + R (x)$

所有的运算在模 $998244353$ 意义下进行。

P4512 【模板】多项式除法

多项式 $A (x)$ 的从 $a_m$ 到 $a_n$ 的系数可由下面的矩阵方程表示：

$\begin{bmatrix} a_n \\ \vdots \\ a_{m+1} \\ a_{m} \end{bmatrix} = \begin{bmatrix} b_m & \dots & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots \\ \dots & \dots & b_m & 0 \\ \dots & \dots & b_{m-1} & b_m \end{bmatrix} \begin{bmatrix}d_{n-m} \\ \vdots \\ d_1 \\ d_0\end{bmatrix}$

为了使我们计算方便，我们引入一种多项式转置运算：

$A^R(x) = x^nA(x^{-1})= a_n + a_{n-1} x + \dots + a_0 x^n$

$B^R(x) = x^m B(x^{-1}) = b_m + b_{m-1} x + \dots + b_0 x^m$

$D^R(x) = x^{n-m}D(x^{-1}) = d_{n-m} + d_{n-m-1} x + \dots + d_0 x^{n-m}$

那么这个线性系统就可以写为：

$A^R(x) \equiv B^R(x)D^R(x) \mod x^{n-m+1}$

根据多项式的乘法逆元，你可以无损的恢复 $D (x)$ 的系数了：

$D^R(x) \equiv A^R(x) (B^R(x))^{-1} \mod x^{n-m+1}$

进而你可以恢复 $R (x) = A (x) - B (x) D (x)$ 。

请务必开 -O2 否则没有返回值优化。

#include <bits/stdc++.h>

using namespace std;

using ll = long long;

#ifdef LLT_DBG
#define FR freopen("in.txt", "r", stdin)
#else
#define FR
#endif

template <ll P>
ll fpow(ll a, ll b)
{
    ll res = 1;
    for (; b; b >>= 1, a = (a * a) % P)
        if (b & 1)
            res = (res * a) % P;
    return res;
}

template <ll G, ll P>
struct NTT
{
    int _n;
    int E;
    vector<int> rev;

    /**
     * @brief 构建一个 NTT 计算器
     *
     * @param n 多项式最高项数
     */
    NTT(int n)
    {
        _n = 1;
        E = 0;
        while (_n < n)
        {
            _n <<= 1;
            E++;
        }
        rev.resize(_n);
        // 逆位置对换
        for (int i = 1; i < _n; i++)
        {
            rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (E - 1));
        }
    }

    void _rNTT(vector<ll> &A, ll k)
    {
        for (int i = 0; i < _n; i++)
            if (i < rev[i])
                swap(A[i], A[rev[i]]);

        for (int e = 1; e <= E; e++)
        {
            int m = 1 << e;

            for (int i = 0; i < _n; i += m)
            {
                int hf = m / 2;
                ll g = 1;
                ll gn = fpow<P>(fpow<P>(G, (P - 1) / m), k);

                for (int j = 0; j < hf; j++)
                {
                    ll x = A[i + j];
                    ll y = (A[i + j + hf] * g) % P;
                    A[i + j] = (x + y) % P;
                    A[i + j + hf] = (x - y) % P;
                    g = (g * gn) % P;
                }
            }
        }
    }

    /**
     * @brief NTT 过程
     *
     * @param A 系数数组
     */
    void doNTT(vector<ll> &A)
    {
        _rNTT(A, 1);
    }

    /**
     * @brief NTT 逆过程
     *
     * @param A 点值数组
     */
    void doINTT(vector<ll> &A)
    {
        ll ni = fpow<P>(_n, P - 2);
        _rNTT(A, P - 2);
        for (int i = 0; i < _n; i++)
            A[i] = (A[i] * ni) % P;
    }
};

template <ll G, ll P>
struct Poly
{
    vector<ll> vec;
    int _n;

    Poly(int n) : vec(n), _n(n)
    {
    }
    Poly(ll an[], int n) : vec(n), _n(n)
    {
        for (int i = 0; i < n; i++)
            vec[i] = an[i];
    }

    void clip(int n)
    {
        vec.resize(n);
        _n = n;
    }

    Poly<G, P> operator+(const Poly<G, P> &o) const
    {
        Poly<G, P> p(max(_n, o._n));

        for (int i = 0; i < p._n; i++)
        {
            p.vec[i] = ((i < _n ? vec[i] : 0) + (i < o._n ? o.vec[i] : 0)) % P;
        }
        return p;
    }

    Poly<G, P> operator-(const Poly<G, P> &o) const
    {
        Poly<G, P> p(max(_n, o._n));

        for (int i = 0; i < p._n; i++)
        {
            p.vec[i] = ((i < _n ? vec[i] : 0) - (i < o._n ? o.vec[i] : 0)) % P;
        }

        return p;
    }

    Poly<G, P> operator*(ll v) const
    {
        Poly p(_n);

        for (int i = 0; i < _n; i++)
            p.vec[i] = (v * vec[i]) % P;

        return p;
    }

    Poly<G, P> operator*(const Poly<G, P> &o) const
    {
        NTT<G, P> ntt(_n + o._n - 1);
        Poly<G, P> p(ntt._n);

        vector<ll> expA(p._n);
        vector<ll> expB(p._n);

        for (int i = 0; i < _n; i++)
            expA[i] = vec[i];
        for (int i = 0; i < o._n; i++)
            expB[i] = o.vec[i];

        ntt.doNTT(expA);
        ntt.doNTT(expB);

        for (int i = 0; i < p._n; i++)
            p.vec[i] = (expA[i] * expB[i]) % P;

        ntt.doINTT(p.vec);

        p.clip(_n + o._n - 1);

        return p;
    }

    Poly<G, P> inv() const
    {
        Poly g(1);
        g.vec[0] = fpow<P>(vec[0], P - 2);

        for (int i = 1; i < _n;)
        {
            i <<= 1;
            Poly vk(i);
            for (int j = 0; j < i; j++)
                vk.vec[j] = (j < _n ? vec[j] : 0);
            g = g * 2 - g * g * vk;
            g.clip(i);
        }
        g.clip(_n);
        return g;
    }

    Poly<G, P> transpose() const
    {
        Poly p(_n);
        for (int i = 0; i < _n; i++)
            p.vec[i] = vec[i];
        reverse(p.vec.begin(), p.vec.end());
        return p;
    }

    pair<Poly<G, P>, Poly<G, P>> operator/(const Poly &o) const
    {
        Poly<G, P> ar = this->transpose();
        Poly<G, P> br = o.transpose();
        br.clip(_n - o._n + 1);
        Poly<G, P> dr = ar * br.inv();
        dr.clip(_n - o._n + 1);
        Poly<G, P> d = dr.transpose();
        Poly<G, P> r = (*this) - o * d;
        r.clip(o._n);

        return make_pair(d, r);
    }
};

ll A[100005];
ll B[100005];
const int mod = 998244353;

void solve()
{
    int n, m;
    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);
    }

    Poly<3, mod> ap(A, n + 1);
    Poly<3, mod> bp(B, m + 1);

    pair<Poly<3, mod>, Poly<3, mod>> dr = ap / bp;

    for (int i = 0; i < n - m + 1; i++)
    {
        printf("%lld ", (dr.first.vec[i] + mod) % mod);
    }

    printf("\n");

    for (int i = 0; i < m; i++)
    {
        printf("%lld ", (dr.second.vec[i] + mod) % mod);
    }
}

int main()
{
    FR;
    solve();
    return 0;
}