hly1204's library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub hly1204/library

:heavy_check_mark: Shift Sample Points
(shift_sample_points.hpp)

Shift via Lagrange Interpolation

Given sample points $f(0), f(1), \dots ,f(n - 1)$ of polynomial $f \in \mathbb{C}\left\lbrack x\right\rbrack$ with $\deg f \lt n$, we want to compute $f(c), f(c + 1), \dots, f(c + m - 1)$ for $c \in \mathbb{Z}, m \in \mathbb{N}$.

Recall the Lagrange interpolation formula,

\[\begin{aligned} f(x) &= \sum _ {0 \leq i \lt n}\left(f(i) \prod _ {0 \leq j \lt n \atop j \neq i} \frac{x - j}{i - j}\right) \\ &= \sum _ {0 \leq i \lt n}\left(f(i)\frac{\prod _ {0 \leq j \lt n \atop j \neq i}\left(x - j\right)}{\prod _ {0 \leq j \lt n \atop j \neq i}\left(i - j\right)}\right) \\ &= \sum _ {0 \leq i \lt n}\left(f(i)\frac{x!}{(x - n)!(x - i)}\frac{(-1)^{n - 1 - i}}{i!(n - 1 - i)!}\right) \\ &= \frac{x!}{(x - n)!}\sum _ {0 \leq i \lt n}\left(f(i)\frac{1}{x - i}\frac{(-1)^{n - 1 - i}}{i!(n - 1 - i)!}\right) \end{aligned}\]

Let

\[\begin{aligned} F(x) &:= \sum _ {0 \leq i \lt n}\frac{f(i)(-1)^{n - 1 - i}}{i!(n - 1 - i)!}x^i, \\ G(x) &:= \sum _ {i \geq 0}\frac{1}{c - (n - 1) + i}x^i \end{aligned}\]

now we have

\[\begin{aligned} \left\lbrack x^{n - 1 + t}\right\rbrack\left(F(x)G(x)\right) &= \sum _ {i = 0}^{n - 1 + t}\left(\left(\left\lbrack x^i\right\rbrack F(x)\right)\left(\left\lbrack x^{n - 1 + t - i}\right\rbrack G(x)\right)\right) \\ &= \sum _ {i = 0}^{n - 1}\left(\frac{f(i)(-1)^{n - 1 - i}}{i!(n - 1 - i)!}\frac{1}{c + t - i}\right) \\ &= \frac{(c + t - n)!}{(c + t)!} f(c + t) \end{aligned}\]

for $t = 0, 1, \dots, m - 1$. We should handle the case that $c - (n - 1) + i = 0$ for a certain $i$.

Shift via O.g.f.

Let

\[F(x) := \sum _ {i \geq 0}f(i)x^i \in \mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack\]

We have $F(x) = \frac{P(x)}{(1 - x)^n}$ where $P(x) \in \mathbb{C}\left\lbrack x\right\rbrack _ {\lt n}$.

And the negative binomial coefficients are

\[\begin{aligned} \frac{1}{(1 - x)^n} &= (1 - x)^{-n} \\ &= \sum _ {k \geq 0}\binom{-n}{k}x^k \\ &= \sum _ {k \geq 0}\frac{(-n)(-n - 1)\cdots (-n - (k - 1))}{k!}x^k \\ &= \sum _ {k \geq 0}(-1)^k\frac{(n + k - 1)!}{(n - 1)!k!}x^k \\ &= \sum _ {k \geq 0}(-1)^k\binom{n + k - 1}{k}x^k \end{aligned}\]

But we are not able to compute $(-n)(-n - 1)\cdots (-n - (k - 1))$ and $k!$ fast if $k$ is large.

Shift via Falling Factorial Polynomial

Might in another document.

Depends on

Verified with

Code

#pragma once

#include "batch_inv.hpp"
#include "binomial.hpp"
#include "middle_product.hpp"
#include "swag.hpp"
#include <cassert>
#include <functional>
#include <vector>

template<typename Tp>
inline std::vector<Tp> shift_sample_points(const std::vector<Tp> &f, Tp c, int m) {
    if (f.empty()) return std::vector<Tp>(m);
    assert(m > 0);
    const int n = f.size();
    auto &&bin  = Binomial<Tp>::get(n);
    std::vector<Tp> F(n), G(n + m - 1);
    for (int i = 0; i < n; ++i) {
        F[i] = f[i] * bin.inv_factorial(i) * bin.inv_factorial(n - 1 - i);
        if ((n - 1 - i) & 1) F[i] = -F[i];
    }
    for (int i = 0; i < n + m - 1; ++i) {
        const auto v = c + (i - (n - 1));
        // We don't care about G[i] when v = 0.
        // We assigned 1 for G[i] when v = 0 for calling batch_inv().
        G[i] = (v == 0) ? Tp(1) : v;
    }
    auto res = middle_product(batch_inv(G), F);
    SWAG<Tp, std::multiplies<>> P(std::multiplies<>{});
    // P = c!/(c-n)! = prod[-n < i <= 0](c + i)
    for (int i = -n + 1; i <= 0; ++i) P.push_back(c + i);
    // res[i] <- (c+i)!/(c+i-n)! * res[i]
    for (int i = 0; i < m; ++i) {
        if (i) P.pop_front(), P.push_back(c + i);
        const auto v = P.query().value();
        // 0 <= c+i < n iff (c+i)!/(c+i-n)! = 0
        res[i] = (v == 0) ? f[(c + i).val()] : v * res[i];
    }
    return res;
}
#line 2 "shift_sample_points.hpp"

#line 2 "batch_inv.hpp"

#include <cassert>
#include <vector>

template<typename Tp> inline std::vector<Tp> batch_inv(const std::vector<Tp> &a) {
    if (a.empty()) return {};
    const int n = a.size();
    std::vector<Tp> b(n);
    Tp v = 1;
    for (int i = 0; i < n; ++i) b[i] = v, v *= a[i];
    assert(v != 0);
    v = v.inv();
    for (int i = n - 1; i >= 0; --i) b[i] *= v, v *= a[i];
    return b;
}
#line 2 "binomial.hpp"

#include <algorithm>
#line 5 "binomial.hpp"

template<typename Tp> class Binomial {
    std::vector<Tp> factorial_, invfactorial_;

    Binomial() : factorial_{Tp(1)}, invfactorial_{Tp(1)} {}

    void preprocess(int n) {
        if (const int nn = factorial_.size(); nn < n) {
            int k = nn;
            while (k < n) k *= 2;
            k = std::min<long long>(k, Tp::mod());
            factorial_.resize(k);
            invfactorial_.resize(k);
            for (int i = nn; i < k; ++i) factorial_[i] = factorial_[i - 1] * i;
            invfactorial_.back() = factorial_.back().inv();
            for (int i = k - 2; i >= nn; --i) invfactorial_[i] = invfactorial_[i + 1] * (i + 1);
        }
    }

public:
    static const Binomial &get(int n) {
        static Binomial bin;
        bin.preprocess(n);
        return bin;
    }

    Tp binom(int n, int m) const {
        return n < m ? Tp() : factorial_[n] * invfactorial_[m] * invfactorial_[n - m];
    }
    Tp inv(int n) const { return factorial_[n - 1] * invfactorial_[n]; }
    Tp factorial(int n) const { return factorial_[n]; }
    Tp inv_factorial(int n) const { return invfactorial_[n]; }
};
#line 2 "middle_product.hpp"

#line 2 "fft.hpp"

#line 5 "fft.hpp"
#include <iterator>
#include <memory>
#line 8 "fft.hpp"

template<typename Tp> class FftInfo {
    static Tp least_quadratic_nonresidue() {
        for (int i = 2;; ++i)
            if (Tp(i).pow((Tp::mod() - 1) / 2) == -1) return Tp(i);
    }

    const int ordlog2_;
    const Tp zeta_;
    const Tp invzeta_;
    const Tp imag_;
    const Tp invimag_;

    mutable std::vector<Tp> root_;
    mutable std::vector<Tp> invroot_;

    FftInfo()
        : ordlog2_(__builtin_ctzll(Tp::mod() - 1)),
          zeta_(least_quadratic_nonresidue().pow((Tp::mod() - 1) >> ordlog2_)),
          invzeta_(zeta_.inv()), imag_(zeta_.pow(1LL << (ordlog2_ - 2))), invimag_(-imag_),
          root_{Tp(1), imag_}, invroot_{Tp(1), invimag_} {}

public:
    static const FftInfo &get() {
        static FftInfo info;
        return info;
    }

    Tp imag() const { return imag_; }
    Tp inv_imag() const { return invimag_; }
    Tp zeta() const { return zeta_; }
    Tp inv_zeta() const { return invzeta_; }
    const std::vector<Tp> &root(int n) const {
        // [0, n)
        assert((n & (n - 1)) == 0);
        if (const int s = root_.size(); s < n) {
            root_.resize(n);
            for (int i = __builtin_ctz(s); (1 << i) < n; ++i) {
                const int j = 1 << i;
                root_[j]    = zeta_.pow(1LL << (ordlog2_ - i - 2));
                for (int k = j + 1; k < j * 2; ++k) root_[k] = root_[k - j] * root_[j];
            }
        }
        return root_;
    }
    const std::vector<Tp> &inv_root(int n) const {
        // [0, n)
        assert((n & (n - 1)) == 0);
        if (const int s = invroot_.size(); s < n) {
            invroot_.resize(n);
            for (int i = __builtin_ctz(s); (1 << i) < n; ++i) {
                const int j = 1 << i;
                invroot_[j] = invzeta_.pow(1LL << (ordlog2_ - i - 2));
                for (int k = j + 1; k < j * 2; ++k) invroot_[k] = invroot_[k - j] * invroot_[j];
            }
        }
        return invroot_;
    }
};

inline int fft_len(int n) {
    --n;
    n |= n >> 1, n |= n >> 2, n |= n >> 4, n |= n >> 8;
    return (n | n >> 16) + 1;
}

namespace detail {

template<typename Iterator> inline void
butterfly_n(Iterator a, int n,
            const std::vector<typename std::iterator_traits<Iterator>::value_type> &root) {
    assert(n > 0);
    assert((n & (n - 1)) == 0);
    const int bn = __builtin_ctz(n);
    if (bn & 1) {
        for (int i = 0; i < n / 2; ++i) {
            const auto a0 = a[i], a1 = a[i + n / 2];
            a[i] = a0 + a1, a[i + n / 2] = a0 - a1;
        }
    }
    for (int i = n >> (bn & 1); i >= 4; i /= 4) {
        const int i4 = i / 4;
        for (int k = 0; k < i4; ++k) {
            const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
            const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
            const auto a02p = a0 + a2, a02m = a0 - a2;
            const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1];
            a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p;
            a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m;
        }
        for (int j = i, m = 2; j < n; j += i, m += 2) {
            const auto r = root[m], r2 = r * r, r3 = r2 * r;
            for (int k = j; k < j + i4; ++k) {
                const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1] * r;
                const auto a2 = a[k + i4 * 2] * r2, a3 = a[k + i4 * 3] * r3;
                const auto a02p = a0 + a2, a02m = a0 - a2;
                const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1];
                a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p;
                a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m;
            }
        }
    }
}

template<typename Iterator> inline void
inv_butterfly_n(Iterator a, int n,
                const std::vector<typename std::iterator_traits<Iterator>::value_type> &root) {
    assert(n > 0);
    assert((n & (n - 1)) == 0);
    const int bn = __builtin_ctz(n);
    for (int i = 4; i <= (n >> (bn & 1)); i *= 4) {
        const int i4 = i / 4;
        for (int k = 0; k < i4; ++k) {
            const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
            const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
            const auto a01p = a0 + a1, a01m = a0 - a1;
            const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1];
            a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = a01m + a23m;
            a[k + i4 * 2] = a01p - a23p, a[k + i4 * 3] = a01m - a23m;
        }
        for (int j = i, m = 2; j < n; j += i, m += 2) {
            const auto r = root[m], r2 = r * r, r3 = r2 * r;
            for (int k = j; k < j + i4; ++k) {
                const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
                const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
                const auto a01p = a0 + a1, a01m = a0 - a1;
                const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1];
                a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = (a01m + a23m) * r;
                a[k + i4 * 2] = (a01p - a23p) * r2, a[k + i4 * 3] = (a01m - a23m) * r3;
            }
        }
    }
    if (bn & 1) {
        for (int i = 0; i < n / 2; ++i) {
            const auto a0 = a[i], a1 = a[i + n / 2];
            a[i] = a0 + a1, a[i + n / 2] = a0 - a1;
        }
    }
}

} // namespace detail

// FFT_n: A(x) |-> bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))]
template<typename Iterator> inline void fft_n(Iterator a, int n) {
    using Tp = typename std::iterator_traits<Iterator>::value_type;
    detail::butterfly_n(a, n, FftInfo<Tp>::get().root(n / 2));
}

template<typename Tp> inline void fft(std::vector<Tp> &a) { fft_n(a.begin(), a.size()); }

// IFFT_n: bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))] |-> A(x)
template<typename Iterator> inline void inv_fft_n(Iterator a, int n) {
    using Tp = typename std::iterator_traits<Iterator>::value_type;
    detail::inv_butterfly_n(a, n, FftInfo<Tp>::get().inv_root(n / 2));
    const Tp iv = Tp::mod() - (Tp::mod() - 1) / n;
    for (int i = 0; i < n; ++i) a[i] *= iv;
}

template<typename Tp> inline void inv_fft(std::vector<Tp> &a) { inv_fft_n(a.begin(), a.size()); }

// IFFT_n^T: A(x) |-> 1/n FFT_n((x^n A(x^(-1))) mod (x^n - 1))
template<typename Iterator> inline void transposed_inv_fft_n(Iterator a, int n) {
    using Tp    = typename std::iterator_traits<Iterator>::value_type;
    const Tp iv = Tp::mod() - (Tp::mod() - 1) / n;
    for (int i = 0; i < n; ++i) a[i] *= iv;
    detail::butterfly_n(a, n, FftInfo<Tp>::get().inv_root(n / 2));
}

template<typename Tp> inline void transposed_inv_fft(std::vector<Tp> &a) {
    transposed_inv_fft_n(a.begin(), a.size());
}

// FFT_n^T : FFT_n((x^n A(x^(-1))) mod (x^n - 1)) |-> n A(x)
template<typename Iterator> inline void transposed_fft_n(Iterator a, int n) {
    using Tp = typename std::iterator_traits<Iterator>::value_type;
    detail::inv_butterfly_n(a, n, FftInfo<Tp>::get().root(n / 2));
}

template<typename Tp> inline void transposed_fft(std::vector<Tp> &a) {
    transposed_fft_n(a.begin(), a.size());
}

template<typename Tp> inline std::vector<Tp> convolution_fft(std::vector<Tp> a, std::vector<Tp> b) {
    if (a.empty() || b.empty()) return {};
    const int n   = a.size();
    const int m   = b.size();
    const int len = fft_len(n + m - 1);
    a.resize(len);
    b.resize(len);
    fft(a);
    fft(b);
    for (int i = 0; i < len; ++i) a[i] *= b[i];
    inv_fft(a);
    a.resize(n + m - 1);
    return a;
}

template<typename Tp> inline std::vector<Tp> square_fft(std::vector<Tp> a) {
    if (a.empty()) return {};
    const int n   = a.size();
    const int len = fft_len(n * 2 - 1);
    a.resize(len);
    fft(a);
    for (int i = 0; i < len; ++i) a[i] *= a[i];
    inv_fft(a);
    a.resize(n * 2 - 1);
    return a;
}

template<typename Tp>
inline std::vector<Tp> convolution_naive(const std::vector<Tp> &a, const std::vector<Tp> &b) {
    if (a.empty() || b.empty()) return {};
    const int n = a.size();
    const int m = b.size();
    std::vector<Tp> res(n + m - 1);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j) res[i + j] += a[i] * b[j];
    return res;
}

template<typename Tp>
inline std::vector<Tp> convolution(const std::vector<Tp> &a, const std::vector<Tp> &b) {
    if (std::min(a.size(), b.size()) < 60) return convolution_naive(a, b);
    if (std::addressof(a) == std::addressof(b)) return square_fft(a);
    return convolution_fft(a, b);
}
#line 7 "middle_product.hpp"

// see:
// [1]: Guillaume Hanrot, Michel Quercia, Paul Zimmermann. The Middle Product Algorithm I.
// [2]: Alin Bostan, Grégoire Lecerf, Éric Schost. Tellegen's principle into practice.

template<typename Tp>
inline std::vector<Tp> middle_product_naive(const std::vector<Tp> &f, const std::vector<Tp> &g) {
    const int m = f.size();
    const int n = g.size();
    assert(m >= n);
    std::vector<Tp> res(m - n + 1);
    for (int i = n - 1; i < m; ++i)
        for (int j = i - (n - 1); j <= i; ++j) res[i - (n - 1)] += f[j] * g[i - j];
    return res;
}

template<typename Tp>
inline std::vector<Tp> middle_product_fft(std::vector<Tp> f, std::vector<Tp> g) {
    const int m = f.size();
    const int n = g.size();
    assert(m >= n);
    std::reverse(g.begin(), g.end());
    const int len = fft_len(m);
    f.resize(len);
    g.resize(len);
    transposed_inv_fft(f);
    fft(g);
    for (int i = 0; i < len; ++i) f[i] *= g[i];
    transposed_fft(f);
    f.resize(m - n + 1);
    return f;
}

// returns (fg)_(n-1),...,(fg)_(m-1)
// f: f_0 + ... + f_(m-1)x^(m-1)
// g: g_0 + ... + g_(n-1)x^(n-1)
// requires m >= n
template<typename Tp>
inline std::vector<Tp> middle_product(const std::vector<Tp> &f, const std::vector<Tp> &g) {
    assert(f.size() >= g.size());
    if (f.size() < 60) return middle_product_naive(f, g);
    return middle_product_fft(f, g);
}
#line 2 "swag.hpp"

#line 4 "swag.hpp"
#include <cstddef>
#include <optional>
#include <stack>
#include <type_traits>
#line 9 "swag.hpp"

// see: https://www.hirzels.com/martin/papers/debs17-tutorial.pdf
// requires: Op(Op(A,B),C) = Op(A,Op(B,C))
template<typename Tp, typename Op,
         std::enable_if_t<std::is_invocable_r_v<Tp, Op, const Tp &, const Tp &>, int> = 0>
class SWAG {
public:
    Op F;
    std::stack<Tp, std::vector<Tp>> Front, Back;
    std::optional<Tp> Agg;

    explicit SWAG(Op F) : F(F) {}
    bool empty() const { return Front.empty() && Back.empty(); }
    std::size_t size() const { return Front.size() + Back.size(); }
    void push_back(const Tp &v) {
        Back.push(v);
        Agg.emplace(Agg ? F(*Agg, v) : v);
    }
    void pop_front() {
        assert(!empty());
        if (Front.empty()) {
            Front.push(Back.top());
            Back.pop();
            while (!Back.empty()) {
                Front.push(F(Back.top(), Front.top()));
                Back.pop();
            }
            Agg.reset();
        }
        Front.pop();
    }

    // returns F(...F(F(Q[0],Q[1]),Q[2]),...,Q[N-1])
    std::optional<Tp> query() const {
        if (empty()) return {};
        if (Front.empty()) return Agg;
        if (!Agg) return Front.top();
        return F(Front.top(), *Agg);
    }
};
#line 8 "shift_sample_points.hpp"
#include <functional>
#line 10 "shift_sample_points.hpp"

template<typename Tp>
inline std::vector<Tp> shift_sample_points(const std::vector<Tp> &f, Tp c, int m) {
    if (f.empty()) return std::vector<Tp>(m);
    assert(m > 0);
    const int n = f.size();
    auto &&bin  = Binomial<Tp>::get(n);
    std::vector<Tp> F(n), G(n + m - 1);
    for (int i = 0; i < n; ++i) {
        F[i] = f[i] * bin.inv_factorial(i) * bin.inv_factorial(n - 1 - i);
        if ((n - 1 - i) & 1) F[i] = -F[i];
    }
    for (int i = 0; i < n + m - 1; ++i) {
        const auto v = c + (i - (n - 1));
        // We don't care about G[i] when v = 0.
        // We assigned 1 for G[i] when v = 0 for calling batch_inv().
        G[i] = (v == 0) ? Tp(1) : v;
    }
    auto res = middle_product(batch_inv(G), F);
    SWAG<Tp, std::multiplies<>> P(std::multiplies<>{});
    // P = c!/(c-n)! = prod[-n < i <= 0](c + i)
    for (int i = -n + 1; i <= 0; ++i) P.push_back(c + i);
    // res[i] <- (c+i)!/(c+i-n)! * res[i]
    for (int i = 0; i < m; ++i) {
        if (i) P.pop_front(), P.push_back(c + i);
        const auto v = P.query().value();
        // 0 <= c+i < n iff (c+i)!/(c+i-n)! = 0
        res[i] = (v == 0) ? f[(c + i).val()] : v * res[i];
    }
    return res;
}
Back to top page