Polynomial Taylor Shift
(library/polynomial/polynomial_taylor_shift.hpp)

View this file on GitHub
Last update: 2024-01-30 20:57:42+09:00
Include: #include "library/polynomial/polynomial_taylor_shift.hpp"

Polynomial Taylor Shift

$\mathbb{F} _ p$ 上の $n$ 次多項式 $\displaystyle f(x) = \sum _ {i = 0} ^ {n} a _ i x ^ i$ および $c\in \mathbb{F} _ p$ に対して、多項式 $g(x) = f(x + c)$ を計算します。

アルゴリズム

\[\begin{aligned} f(x + c) &=\sum _ {i = 0} ^ {n} a _ i (x + c) ^ i \\ &=\sum _ {i = 0} ^ {n} \sum _{j = 0} ^ i a _ i \cdot \binom{i}{j}\cdot c ^ {i - j} \cdot x ^ j \\ &=\sum _ {j = 0} ^ {n} \dfrac{x ^ j}{j!} \sum _ {i = j} ^ {n} (a _ i \cdot i!) \cdot \dfrac{c ^ {i - j}}{(i - j)!} \\ &=\sum _ {j = 0} ^ {n} \dfrac{x ^ j}{j!} \sum _ {k = 0} ^ {n - j} (a _ {j + k} \cdot (j + k)!) \cdot \dfrac{c ^ k}{k!} \\ &=\sum _ {j = 0} ^ {n} \dfrac{x ^ j}{j!} \sum _ {k = 0} ^ {n - j} b _ {(n - j) - k} \cdot \dfrac{c ^ k}{k!} \quad (b _ i := a _ {n - i} \cdot (n - i)! )\\ &=\sum _ {j = 0} ^ {n} \dfrac{x ^ {n - j}}{(n - j)!} \sum _ {k = 0} ^ {j} b _ {j - k} \cdot \dfrac{c ^ k}{k!} \end{aligned}\]

で、各 $j = 0, \ldots, n$ に対する $\displaystyle \sum _ {k = 0} ^ {j} b _ {j - k} \cdot \dfrac{c ^ k}{k!}$ は高速フーリエ変換等を用いた畳み込みにより $\Theta(n \log n)$ 時間で計算できます。

従って、$0!,\ldots,n!$ の乗法逆元の前計算の下で、$f(x + c)$ も $\Theta(n \log n)$ 時間で計算できます。

なお、$0!,\ldots,n!$ の乗法逆元が存在すること、即ち $n \lt p$ を仮定していることに注意します。

Depends on

階乗テーブル (library/math/factorial.hpp)

Verified with

Code

#ifndef SUISEN_POLYNOMIAL_TAYLOR_SHIFT
#define SUISEN_POLYNOMIAL_TAYLOR_SHIFT

#include <algorithm>

#include "library/math/factorial.hpp"

namespace suisen {
    // return f(x + c) 
    template <typename FPSType, typename T>
    FPSType translate(const FPSType& f, const T c) {
        int d = f.deg();
        if (d < 0) return FPSType{ 0 };
        using mint = typename FPSType::value_type;
        factorial<mint> fac(d);
        FPSType expc(d + 1), g(d + 1);
        mint p = 1;
        for (int i = 0; i <= d; ++i, p *= c) {
            expc[i] = p * fac.fac_inv(i);
            g[d - i] = f[i] * fac(i);
        }
        g *= expc, g.resize(d + 1);
        for (int i = 0; i <= d; ++i) g[i] *= fac.fac_inv(d - i);
        std::reverse(g.begin(), g.end());
        return g;
    }
} // namespace suisen

#endif // SUISEN_POLYNOMIAL_TAYLOR_SHIFT

#line 1 "library/polynomial/polynomial_taylor_shift.hpp"



#include <algorithm>

#line 1 "library/math/factorial.hpp"



#include <cassert>
#include <vector>

namespace suisen {
    template <typename T, typename U = T>
    struct factorial {
        factorial() = default;
        factorial(int n) { ensure(n); }

        static void ensure(const int n) {
            int sz = _fac.size();
            if (n + 1 <= sz) return;
            int new_size = std::max(n + 1, sz * 2);
            _fac.resize(new_size), _fac_inv.resize(new_size);
            for (int i = sz; i < new_size; ++i) _fac[i] = _fac[i - 1] * i;
            _fac_inv[new_size - 1] = U(1) / _fac[new_size - 1];
            for (int i = new_size - 1; i > sz; --i) _fac_inv[i - 1] = _fac_inv[i] * i;
        }

        T fac(const int i) {
            ensure(i);
            return _fac[i];
        }
        T operator()(int i) {
            return fac(i);
        }
        U fac_inv(const int i) {
            ensure(i);
            return _fac_inv[i];
        }
        U binom(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[r] * _fac_inv[n - r];
        }
        template <typename ...Ds, std::enable_if_t<std::conjunction_v<std::is_integral<Ds>...>, std::nullptr_t> = nullptr>
        U polynom(const int n, const Ds& ...ds) {
            if (n < 0) return 0;
            ensure(n);
            int sumd = 0;
            U res = _fac[n];
            for (int d : { ds... }) {
                if (d < 0 or d > n) return 0;
                sumd += d;
                res *= _fac_inv[d];
            }
            if (sumd > n) return 0;
            res *= _fac_inv[n - sumd];
            return res;
        }
        U perm(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[n - r];
        }
    private:
        static std::vector<T> _fac;
        static std::vector<U> _fac_inv;
    };
    template <typename T, typename U>
    std::vector<T> factorial<T, U>::_fac{ 1 };
    template <typename T, typename U>
    std::vector<U> factorial<T, U>::_fac_inv{ 1 };
} // namespace suisen


#line 7 "library/polynomial/polynomial_taylor_shift.hpp"

namespace suisen {
    // return f(x + c) 
    template <typename FPSType, typename T>
    FPSType translate(const FPSType& f, const T c) {
        int d = f.deg();
        if (d < 0) return FPSType{ 0 };
        using mint = typename FPSType::value_type;
        factorial<mint> fac(d);
        FPSType expc(d + 1), g(d + 1);
        mint p = 1;
        for (int i = 0; i <= d; ++i, p *= c) {
            expc[i] = p * fac.fac_inv(i);
            g[d - i] = f[i] * fac(i);
        }
        g *= expc, g.resize(d + 1);
        for (int i = 0; i <= d; ++i) g[i] *= fac.fac_inv(d - i);
        std::reverse(g.begin(), g.end());
        return g;
    }
} // namespace suisen

Polynomial Taylor Shift (library/polynomial/polynomial_taylor_shift.hpp)