$\left(\prod_{k=0}^{m-1} f(r^k x)\right) \bmod x^N$
(library/polynomial/prod_f(r^k_x).hpp)

• View this file on GitHub
• Last update: 2023-07-09 04:04:16+09:00
• Include: #include "library/polynomial/prod_f(r^k_x).hpp"

$\left(\prod_{k=0}^{m-1} f(r^k x)\right) \bmod x^N$

$\displaystyle F(x) := \prod_{k=0}^{m-1} f(r^k x)$ とする。形式的冪級数 $f$ の $i$ 次の係数を $f_i$ と書く。

$f(x)=0$ の場合は自明なので、$f(x)\neq 0$ を仮定する。また、$f _ 0=1$ と仮定する。

$f _ 0\neq 1$ の場合は $f = c\cdot x ^ a\cdot f’\; (f’ _ 0 = 1)$ とすれば、次のようにして定数項が $1$ のケースに帰着できる。

\[F(x)=c ^ m r ^ {\frac{m(m-1)}{2}a} x ^ {\frac{m(m-1)}{2}} \prod _ {k = 0} ^ {m - 1} f'(r ^ k x).\]

$g=\log f$ とおくと、$F = \exp (\log F)$ を用いることで、$F \bmod x ^ N$ を $O(N \log N)$ 時間で計算できる。

\[\begin{aligned} \log F &= \sum _ {k = 0} ^ {m - 1} \log f(r ^ k x) \newline &= \sum _ {k = 0} ^ {m - 1} g(r ^ k x) \newline &= \sum _ {i = 0} ^ {\infty} g _ i x ^ i\sum _ {k = 0} ^ {m - 1} (r ^ i) ^ k \newline &= \sum _ {i = 0} ^ {\infty} g _ i x ^ i \cdot \begin{cases} \displaystyle m \cdot x ^ i & \text {if}\quad r ^ i = 1 \newline \displaystyle \dfrac{(r ^ m) ^ i - 1}{r ^ i - 1} & \text{otherwise} \end{cases}. \end{aligned}\]

参考

• https://maspypy.com/%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%83%BB%E5%BD%A2%E5%BC%8F%E7%9A%84%E3%81%B9%E3%81%8D%E7%B4%9A%E6%95%B0-%E9%AB%98%E9%80%9F%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%A7%E3%81%8D%E3%82%8B%E3%82%82%E3%81%AE#toc21

Depends on

• 冪乗テーブル (library/math/pow_mods.hpp)

Verified with

• test/src/polynomial/prod_f(r^k_x)/yuki2097.test.cpp

Code

#ifndef SUISEN_PROD_F_RK_X
#define SUISEN_PROD_F_RK_X

#include "library/math/pow_mods.hpp"

/**
 * @brief $\left(\prod_{k=0}^{m-1} f(r^k x)\right) \bmod x^N$
 */

namespace suisen {
    namespace internal::prod_f_rk_x {
        template <typename FPSType>
        FPSType prod_f_rk_x(FPSType f, typename FPSType::value_type r, int m, int result_size) {
            using mint = typename FPSType::value_type;
            pow_mods<mint> pow_r(r, result_size), pow_rm(r.pow(m), result_size);
            if (auto opt_sp_f = f.sparse_fps_format(15); opt_sp_f.has_value()) {
                bool all_invertible = true;
                for (int i = 1; i < result_size; ++i) {
                    if (pow_r[i] == mint{ 1 }) {
                        all_invertible = false;
                        break;
                    }
                }
                if (all_invertible) {
                    auto &sp_f = *opt_sp_f;
                    sp_f.erase(sp_f.begin());
                    FPSType g(result_size);
                    g[0] = 1;
                    for (int i = 1; i < result_size; ++i) {
                        for (auto [j, fj] : sp_f) {
                            if (j > i) break;
                            g[i] += g[i - j] * fj * (pow_r[i - j] - pow_rm[j]);
                        }
                        g[i] /= 1 - pow_r[i];
                    }
                    return g;
                }
            }
            f = f.log(result_size);
            for (int i = 1; i < result_size; ++i) f[i] *= pow_r[i] == mint{ 1 } ? mint{ m } : (pow_rm[i] - 1) / (pow_r[i] - 1);
            return f.exp(result_size);
        }
    }
    /**
     * \brief Calculates Π[k=0,m-1] f(r^k x) mod x^N in O(NlogN) time.
     * \tparam FPSType type of formal power series
     * \param f formal power series
     * \param r ratio
     * \param m the number of terms of the product
     * \param result_size N (default: size of f)
     * \return Π[k=0,m-1] f(r^k x) mod x^N 
     */
    template <typename FPSType>
    FPSType prod_f_rk_x(FPSType f, const typename FPSType::value_type r, const int m, int result_size = -1) {
        using mint = typename FPSType::value_type;
        if (result_size < 0) result_size = f.size();
        if (r == mint{ 1 }) return f.pow(m, result_size);
        if (m == 0) { FPSType res{ 1 }; res.resize(result_size); return res; }
        int z = 0;
        while (z < int(f.size()) and f[z] == mint{ 0 }) ++z;
        if (z == int(f.size()) or z >= (result_size + m - 1) / m) return FPSType(result_size, mint{ 0 });
        const mint c = f[z], d = c.pow(m) * r.pow((long long) m * (m - 1) / 2 * z);
        f >>= z, f /= c; // => f[0] = 1
        f = internal::prod_f_rk_x::prod_f_rk_x(f, r, m, result_size - z * m);
        f *= d, f <<= z * m;
        return f;
    }
} // namespace suisen

#endif // SUISEN_PROD_F_RK_X

#line 1 "library/polynomial/prod_f(r^k_x).hpp"



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



#include <vector>

namespace suisen {
    template <int base_as_int, typename mint>
    struct static_pow_mods {
        static_pow_mods() = default;
        static_pow_mods(int n) { ensure(n); }
        const mint& operator[](int i) const {
            ensure(i);
            return pows[i];
        }
        static void ensure(int n) {
            int sz = pows.size();
            if (sz > n) return;
            pows.resize(n + 1);
            for (int i = sz; i <= n; ++i) pows[i] = base * pows[i - 1];
        }
    private:
        static inline std::vector<mint> pows { 1 };
        static inline mint base = base_as_int;
        static constexpr int mod = mint::mod();
    };

    template <typename mint>
    struct pow_mods {
        pow_mods() = default;
        pow_mods(mint base, int n) : base(base) { ensure(n); }
        const mint& operator[](int i) const {
            ensure(i);
            return pows[i];
        }
        void ensure(int n) const {
            int sz = pows.size();
            if (sz > n) return;
            pows.resize(n + 1);
            for (int i = sz; i <= n; ++i) pows[i] = base * pows[i - 1];
        }
    private:
        mutable std::vector<mint> pows { 1 };
        mint base;
        static constexpr int mod = mint::mod();
    };
}


#line 5 "library/polynomial/prod_f(r^k_x).hpp"

/**
 * @brief $\left(\prod_{k=0}^{m-1} f(r^k x)\right) \bmod x^N$
 */

namespace suisen {
    namespace internal::prod_f_rk_x {
        template <typename FPSType>
        FPSType prod_f_rk_x(FPSType f, typename FPSType::value_type r, int m, int result_size) {
            using mint = typename FPSType::value_type;
            pow_mods<mint> pow_r(r, result_size), pow_rm(r.pow(m), result_size);
            if (auto opt_sp_f = f.sparse_fps_format(15); opt_sp_f.has_value()) {
                bool all_invertible = true;
                for (int i = 1; i < result_size; ++i) {
                    if (pow_r[i] == mint{ 1 }) {
                        all_invertible = false;
                        break;
                    }
                }
                if (all_invertible) {
                    auto &sp_f = *opt_sp_f;
                    sp_f.erase(sp_f.begin());
                    FPSType g(result_size);
                    g[0] = 1;
                    for (int i = 1; i < result_size; ++i) {
                        for (auto [j, fj] : sp_f) {
                            if (j > i) break;
                            g[i] += g[i - j] * fj * (pow_r[i - j] - pow_rm[j]);
                        }
                        g[i] /= 1 - pow_r[i];
                    }
                    return g;
                }
            }
            f = f.log(result_size);
            for (int i = 1; i < result_size; ++i) f[i] *= pow_r[i] == mint{ 1 } ? mint{ m } : (pow_rm[i] - 1) / (pow_r[i] - 1);
            return f.exp(result_size);
        }
    }
    /**
     * \brief Calculates Π[k=0,m-1] f(r^k x) mod x^N in O(NlogN) time.
     * \tparam FPSType type of formal power series
     * \param f formal power series
     * \param r ratio
     * \param m the number of terms of the product
     * \param result_size N (default: size of f)
     * \return Π[k=0,m-1] f(r^k x) mod x^N 
     */
    template <typename FPSType>
    FPSType prod_f_rk_x(FPSType f, const typename FPSType::value_type r, const int m, int result_size = -1) {
        using mint = typename FPSType::value_type;
        if (result_size < 0) result_size = f.size();
        if (r == mint{ 1 }) return f.pow(m, result_size);
        if (m == 0) { FPSType res{ 1 }; res.resize(result_size); return res; }
        int z = 0;
        while (z < int(f.size()) and f[z] == mint{ 0 }) ++z;
        if (z == int(f.size()) or z >= (result_size + m - 1) / m) return FPSType(result_size, mint{ 0 });
        const mint c = f[z], d = c.pow(m) * r.pow((long long) m * (m - 1) / 2 * z);
        f >>= z, f /= c; // => f[0] = 1
        f = internal::prod_f_rk_x::prod_f_rk_x(f, r, m, result_size - z * m);
        f *= d, f <<= z * m;
        return f;
    }
} // namespace suisen

Back to top page