Shift of Sampling Points of Polynomial
(library/polynomial/shift_of_sampling_points.hpp)

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

Shift of Sampling Points of Polynomial

Shift of Sampling Points of Polynomial を解く実装。hos-lyric さんのコメントを解読しました。

解法

Step 1. $\displaystyle f(x) = \sum _ {i = 0} ^ {N - 1} a _ i x ^ {\underline{i}}$ を満たす列 $a=(a _ 0,\ldots,a _ {N-1})$ を求める

ラグランジュの補間公式より、$f$ は次の表示を持つ:

\[\begin{aligned} f(x) &= \sum _ {i = 0} ^ {N -1} f(i)\prod _ {j\in [N]\setminus \lbrace i \rbrace}\dfrac{x - j}{i - j} \newline &= \sum _ {i = 0} ^ {N -1} \dfrac{(-1) ^ {N-i-1} \cdot f(i)}{i!\cdot (N-1-i)!} \prod _ {j \in [N] \setminus \lbrace i \rbrace} (x - j). \end{aligned}\]

見やすさのために $p _ i = \dfrac{(-1) ^ {N-i-1} \cdot f(i)}{i!\cdot (N-1-i)!}$ とすると、次のように書くことができる。

\[f(x) = \sum _ {i = 0} ^ {N - 1} p _ i \prod _ {j \in [N] \setminus \lbrace i \rbrace} (x - j).\]

最高次の係数に注目することで $\displaystyle a _ {N - 1} = \sum _ {i = 0} ^ {N - 1} p _ i$ を得る。$f(x) - a _ {N - 1} x ^ {\underline{N - 1}}$ を計算してみると、

\[\begin{aligned} f(x) - a _ {N - 1} x ^ {\underline{N - 1}} &=\sum _ {i = 0} ^ {N - 2} p _ i \cdot (-1) \cdot (N - 1 - i) \prod _ {j \in [N - 1]\setminus\lbrace i \rbrace} (x - j) \newline &=\sum _ {i = 0} ^ {N - 2}p' _ i \prod _ {j \in [N - 1]\setminus\lbrace i \rbrace} (x - j) && (p' _ i := p _ i \cdot (-1) \cdot (N - 1 - i)) \end{aligned}\]

となり、サイズが $1$ だけ小さい問題が得られる。これを繰り返すことで、次を得る。

\[a _ {N - k} = (-1) ^ {k + 1}\sum _ {i = 0} ^ {N - k} p _ i \cdot \dfrac{(N - i - 1)!}{(N - i - k)!}.\]

これを整理して、次を得る。

\[a _ i = \sum _ {j = 0} ^ i \dfrac{f(j)}{j!} \cdot \dfrac{(-1) ^ {i - j}}{(i - j) !}.\]

これは畳み込みで高速化できる形なので、$a$ を $O(N \log N)$ 時間で計算することができた。

Step 2. $k = 0, \ldots, M-1$ に対する $f(c + k)$ の計算

下降階乗冪に関する以下の恒等式が重要。

\[(a + b) ^ {\underline{n}} = \sum _ {i = 0} ^ n \binom{n}{i} a ^ {\underline{i}} b ^ {\underline{n - i}}.\]

Step 1 で $\displaystyle f(x) = \sum _ {i = 0} ^ {N - 1} a _ i x ^ {\underline{i}}$ という表示を得ているので、これに $x = c + k$ を代入して整理すると、次のようになる。(Polynomial Taylor Shift のときと全く同じように変形する)

\[\begin{aligned} f(c+k) &= \sum _ {i = 0} ^ {N - 1} a _ i (c + k) ^ {\underline{i}} \newline &= \sum _ {i = 0} ^ {N - 1} a _ i \sum _ {j = 0} ^ i \binom{i}{j} c ^ {\underline{i - j}} k ^ {\underline{j}} \newline &= \sum _ {j = 0} ^ {N - 1} \dfrac{k ^ {\underline{j}}}{j!} \sum _ {i = j} ^ {N - 1} (a _ i \cdot i!)\cdot \dfrac{c ^ {\underline{i - j}}}{(i - j)!} \newline &= \sum _ {j = 0} ^ {k} \binom{k}{j} \sum _ {i = 0} ^ {N - 1 - j} (a _ {i + j} \cdot (i + j)!)\cdot \dfrac{c ^ {\underline{i}}}{i!} \newline &= \sum _ {j = 0} ^ {k} \binom{k}{j} \sum _ {i = 0} ^ {N - 1 - j} a' _ {(N - 1 - j) - i}\cdot \dfrac{c ^ {\underline{i}}}{i!} && (a' _ i := a _ {N - i - 1} \cdot (N - i - 1)!) \newline &= \sum _ {j = 0} ^ {k} \binom{k}{j} b _ {N - 1 - j} && (b _ i := \sum _ {j = 0} ^ i a' _ {i - j}\cdot \dfrac{c ^ {\underline{j}}}{j!}) \newline &= k! \sum _ {j = 0} ^ k \dfrac{b _ {N - 1 - j}}{j!} \cdot \dfrac{1}{(k - j)!}. \end{aligned}\]

途中で定義した $b$ は畳み込みの形をしているので $O(N \log N)$ 時間で計算することができる。また、$c ^ {\underline{i}}$ に関しては $i$ の昇順に計算することで必要な部分を全体 $O(N)$ 時間で列挙できる。

最後の式も畳み込みの形をしているので $O((M + N) \log (M + N))$ 時間で計算することができる。

以上より、全体 $O((M + N) \log (M + N))$ 時間で問題を解くことができた。

Depends on

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

Required by

Verified with

Code

#ifndef SUISEN_SHIFT_OF_SAMPLING_POINTS
#define SUISEN_SHIFT_OF_SAMPLING_POINTS

#include <vector>
#include <atcoder/convolution>

#include "library/math/factorial.hpp"

namespace suisen {
    template <typename mint, typename Convolve,
        std::enable_if_t<std::is_invocable_r_v<std::vector<mint>, Convolve, std::vector<mint>, std::vector<mint>>, std::nullptr_t> = nullptr>
    std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m, const Convolve &convolve) {
        const int n = ys.size();
        factorial<mint> fac(std::max(n, m));

        std::vector<mint> b = [&] {
            std::vector<mint> f(n), g(n);
            for (int i = 0; i < n; ++i) {
                f[i] = ys[i] * fac.fac_inv(i);
                g[i] = (i & 1 ? -1 : 1) * fac.fac_inv(i);
            }
            std::vector<mint> b = convolve(f, g);
            b.resize(n);
            return b;
        }();
        std::vector<mint> e = [&] {
            std::vector<mint> c(n);
            mint prd = 1;
            std::reverse(b.begin(), b.end());
            for (int i = 0; i < n; ++i) {
                b[i] *= fac.fac(n - i - 1);
                c[i] = prd * fac.fac_inv(i);
                prd *= t - i;
            }
            std::vector<mint> e = convolve(b, c);
            e.resize(n);
            return e;
        }();
        std::reverse(e.begin(), e.end());
        for (int i = 0; i < n; ++i) {
            e[i] *= fac.fac_inv(i);
        }

        std::vector<mint> f(m);
        for (int i = 0; i < m; ++i) f[i] = fac.fac_inv(i);
        std::vector<mint> res = convolve(e, f);
        res.resize(m);
        for (int i = 0; i < m; ++i) res[i] *= fac.fac(i);
        return res;
    }

    template <typename mint>
    std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m) {
        auto convolve = [&](const std::vector<mint> &f, const std::vector<mint> &g) { return atcoder::convolution(f, g); };
        return shift_of_sampling_points(ys, t, m, convolve);
    }
} // namespace suisen


#endif // SUISEN_SHIFT_OF_SAMPLING_POINTS

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



#include <vector>
#include <atcoder/convolution>

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



#include <cassert>
#line 6 "library/math/factorial.hpp"

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 8 "library/polynomial/shift_of_sampling_points.hpp"

namespace suisen {
    template <typename mint, typename Convolve,
        std::enable_if_t<std::is_invocable_r_v<std::vector<mint>, Convolve, std::vector<mint>, std::vector<mint>>, std::nullptr_t> = nullptr>
    std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m, const Convolve &convolve) {
        const int n = ys.size();
        factorial<mint> fac(std::max(n, m));

        std::vector<mint> b = [&] {
            std::vector<mint> f(n), g(n);
            for (int i = 0; i < n; ++i) {
                f[i] = ys[i] * fac.fac_inv(i);
                g[i] = (i & 1 ? -1 : 1) * fac.fac_inv(i);
            }
            std::vector<mint> b = convolve(f, g);
            b.resize(n);
            return b;
        }();
        std::vector<mint> e = [&] {
            std::vector<mint> c(n);
            mint prd = 1;
            std::reverse(b.begin(), b.end());
            for (int i = 0; i < n; ++i) {
                b[i] *= fac.fac(n - i - 1);
                c[i] = prd * fac.fac_inv(i);
                prd *= t - i;
            }
            std::vector<mint> e = convolve(b, c);
            e.resize(n);
            return e;
        }();
        std::reverse(e.begin(), e.end());
        for (int i = 0; i < n; ++i) {
            e[i] *= fac.fac_inv(i);
        }

        std::vector<mint> f(m);
        for (int i = 0; i < m; ++i) f[i] = fac.fac_inv(i);
        std::vector<mint> res = convolve(e, f);
        res.resize(m);
        for (int i = 0; i < m; ++i) res[i] *= fac.fac(i);
        return res;
    }

    template <typename mint>
    std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m) {
        auto convolve = [&](const std::vector<mint> &f, const std::vector<mint> &g) { return atcoder::convolution(f, g); };
        return shift_of_sampling_points(ys, t, m, convolve);
    }
} // namespace suisen

Shift of Sampling Points of Polynomial (library/polynomial/shift_of_sampling_points.hpp)