chirp z-transform (評価点が等差数列を成す場合の Multipoint Evaluation)
(library/transform/chirp_z_transform.hpp)

View this file on GitHub
Last update: 2022-11-16 20:35:58+09:00
Include: #include "library/transform/chirp_z_transform.hpp"

chirp z-transform (評価点が等差数列を成す場合の Multipoint Evaluation)

問題

多項式 $\displaystyle f(x) = \sum _ {i = 0} ^ {n - 1} f _ i x ^ i$ と評価点の初項 $a$ および公比 $r$ 与えられるので、$k = 0,\ldots, m - 1$ に対する $F _ k := f(a r ^ k)$ の値を求めよ。

概要

三角数 $t _ i = \dfrac{i(i - 1)}{2}$ に関して、以下が成り立つ。

\[ki = t _ {i + k} - t _ i - t _ k.\]

即ち、次が成り立つ。

\[\begin{aligned} F _ k &= r ^ {-t _ k} \sum _ {i = 0} ^ {n - 1} (f _ i a ^ i r ^ {-t _ i}) r ^ {t _ {i + k}}. \end{aligned}\]

$x _ i := f _ {n - 1 - i} a ^ {n - i - 1} r ^ {-t _ {n - i - 1}}, y _ i := r ^ {t _ i}$ とすると

\[F _ k = r ^ {-t _ k} \sum _ {i = 0} ^ {n - 1} x _ {n - 1 - i} y _ {k + i},\]

であるから、$F _ k$ は $x$ と $y$ の畳み込みで計算できる。

参考

https://github.com/noshi91/algorithm-encyclopedia/blob/gh-pages/_algorithms/chirp-z-transform.md

Required by

Multi Variate Convolution Circular (多変数巡回畳み込み) (library/convolution/multi_variate_convolution_circular.hpp)

Verified with

Code

#ifndef SUISEN_CHIRP_Z_TRANSFORM
#define SUISEN_CHIRP_Z_TRANSFORM

#include <algorithm>
#include <vector>

#include <atcoder/convolution>

/**
 * @brief chirp z-transform ($g _ k = f(a r^k)$)
*/

namespace suisen {
    namespace internal {
        const auto default_convolution = [](const auto& a, const auto& b) { return atcoder::convolution(a, b); };

        template <typename T>
        std::vector<T> chirp_z_transform_naive(const std::vector<T> &f, T a, T r, int m) {
            const int n = f.size();
            std::vector<T> g(m);
            T pow_r = 1;
            for (int k = 0; k < m; ++k) {
                T ark = a * pow_r, pow_ark = 1;
                for (int i = 0; i < n; ++i) {
                    g[k] += f[i] * pow_ark;
                    pow_ark *= ark;
                }
                pow_r *= r;
            }
            return g;
        }
    } // namespace internal
    /**
     * @brief Calculates f(ar^k) for k=0,...,m-1 in O(M(n+m-1)+n+m) time
     */
    template <typename T, typename Convolution>
    std::vector<T> chirp_z_transform(std::vector<T> f, T a, T r, int m, Convolution&& convolution = internal::default_convolution) {
        const int n = f.size();
        std::vector<T> g(m);
        if (n == 0 or m == 0) return g;
        T pow_a = 1;
        for (int i = 0; i < n; ++i, pow_a *= a) f[i] *= pow_a;
        if (r == 0) {
            for (int i = 0; i < n; ++i) g[0] += f[i];
            for (int k = 1; k < m; ++k) g[k] += f[0];
            return g;
        }
        if (n < 60) return internal::chirp_z_transform_naive(f, a, r, m);
        const T r_inv = r.inv();

        const int l = n + m - 1;

        std::vector<T> pow_r_tri(l), pow_r_tri_inv(l);
        pow_r_tri[0] = pow_r_tri_inv[0] = 1;

        T pow_r = 1, pow_r_inv = 1;
        for (int i = 1; i < l; ++i, pow_r *= r, pow_r_inv *= r_inv) {
            pow_r_tri[i] = pow_r_tri[i - 1] * pow_r;
            pow_r_tri_inv[i] = pow_r_tri_inv[i - 1] * pow_r_inv;
        }

        std::vector<T> p(n), q(l);
        for (int i = 0; i < n; ++i) p[i] = f[i] * pow_r_tri_inv[i];
        for (int i = 0; i < l; ++i) q[i] = pow_r_tri[i];
        std::reverse(p.begin(), p.end());
        std::vector<T> pq = convolution(p, q);
        for (int k = 0; k < m; ++k) {
            g[k] = pow_r_tri_inv[k] * pq[n - 1 + k];
        }

        return g;
    }
} // namespace suisen


#endif // SUISEN_CHIRP_Z_TRANSFORM

#line 1 "library/transform/chirp_z_transform.hpp"



#include <algorithm>
#include <vector>

#include <atcoder/convolution>

/**
 * @brief chirp z-transform ($g _ k = f(a r^k)$)
*/

namespace suisen {
    namespace internal {
        const auto default_convolution = [](const auto& a, const auto& b) { return atcoder::convolution(a, b); };

        template <typename T>
        std::vector<T> chirp_z_transform_naive(const std::vector<T> &f, T a, T r, int m) {
            const int n = f.size();
            std::vector<T> g(m);
            T pow_r = 1;
            for (int k = 0; k < m; ++k) {
                T ark = a * pow_r, pow_ark = 1;
                for (int i = 0; i < n; ++i) {
                    g[k] += f[i] * pow_ark;
                    pow_ark *= ark;
                }
                pow_r *= r;
            }
            return g;
        }
    } // namespace internal
    /**
     * @brief Calculates f(ar^k) for k=0,...,m-1 in O(M(n+m-1)+n+m) time
     */
    template <typename T, typename Convolution>
    std::vector<T> chirp_z_transform(std::vector<T> f, T a, T r, int m, Convolution&& convolution = internal::default_convolution) {
        const int n = f.size();
        std::vector<T> g(m);
        if (n == 0 or m == 0) return g;
        T pow_a = 1;
        for (int i = 0; i < n; ++i, pow_a *= a) f[i] *= pow_a;
        if (r == 0) {
            for (int i = 0; i < n; ++i) g[0] += f[i];
            for (int k = 1; k < m; ++k) g[k] += f[0];
            return g;
        }
        if (n < 60) return internal::chirp_z_transform_naive(f, a, r, m);
        const T r_inv = r.inv();

        const int l = n + m - 1;

        std::vector<T> pow_r_tri(l), pow_r_tri_inv(l);
        pow_r_tri[0] = pow_r_tri_inv[0] = 1;

        T pow_r = 1, pow_r_inv = 1;
        for (int i = 1; i < l; ++i, pow_r *= r, pow_r_inv *= r_inv) {
            pow_r_tri[i] = pow_r_tri[i - 1] * pow_r;
            pow_r_tri_inv[i] = pow_r_tri_inv[i - 1] * pow_r_inv;
        }

        std::vector<T> p(n), q(l);
        for (int i = 0; i < n; ++i) p[i] = f[i] * pow_r_tri_inv[i];
        for (int i = 0; i < l; ++i) q[i] = pow_r_tri[i];
        std::reverse(p.begin(), p.end());
        std::vector<T> pq = convolution(p, q);
        for (int k = 0; k < m; ++k) {
            g[k] = pow_r_tri_inv[k] * pq[n - 1 + k];
        }

        return g;
    }
} // namespace suisen

chirp z-transform (評価点が等差数列を成す場合の Multipoint Evaluation) (library/transform/chirp_z_transform.hpp)