ks.hpp

• View this file on GitHub
• Last update: 2025-01-19 15:28:01+08:00
• Include: #include "ks.hpp"
• Link: https://doi.org/10.1016/j.jsc.2009.05.004

Depends on

• FFT (fft.hpp)

Verified with

• test/convolution/2d_convolution.0.test.cpp

Code

#pragma once

#include "fft.hpp"
#include <algorithm>
#include <cassert>
#include <vector>

// returns max[0 <= j < a.size()]{a[j].size()}
template<typename Tp> inline int max_len_x_ks(const std::vector<std::vector<Tp>> &a) {
    int len = -1;
    for (int i = 0; i < (int)a.size(); ++i) len = std::max<int>(len, a[i].size());
    return len;
}

// returns a(x, x^N) where a(x,y) in R[x][y]
template<typename Tp>
inline std::vector<Tp> pack_2d_ks(const std::vector<std::vector<Tp>> &a, int N) {
    assert(N > 0);
    // y |-> x^N
    std::vector<Tp> b(N * a.size());
    for (int i = 0; i < (int)a.size(); ++i)
        for (int j = 0; j < (int)a[i].size(); ++j) b[i * N + j] += a[i][j]; // if N < a[i].size()
    return b;
}

template<typename Tp>
inline std::vector<std::vector<Tp>> unpack_2d_ks(const std::vector<Tp> &a, int N) {
    assert(N > 0);
    // x^N |-> y
    std::vector<std::vector<Tp>> b((a.size() + (N - 1)) / N, std::vector<Tp>(N));
    for (int i = 0; i < (int)((a.size() + N - 1) / N); ++i)
        for (int j = 0; j < N && i * N + j < (int)a.size(); ++j) b[i][j] = a[i * N + j];
    return b;
}

template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_naive(const std::vector<std::vector<Tp>> &a, const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = lenA + lenB - 1;
    std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(N));
    for (int i = 0; i < (int)a.size(); ++i)
        for (int j = 0; j < (int)b.size(); ++j) {
            const auto aibj = convolution(a[i], b[j]);
            for (int k = 0; k < (int)aibj.size(); ++k) ab[i + j][k] += aibj[k];
        }
    return ab;
}

// standard Kronecker substitution
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks(const std::vector<std::vector<Tp>> &a, const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = lenA + lenB - 1;
    auto ab     = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));
    ab.resize((a.size() + b.size() - 1) * N);
    return unpack_2d_ks(ab, N);
}

// see:
// [1]: David Harvey. Faster polynomial multiplication via multipoint Kronecker substitution.
//      https://doi.org/10.1016/j.jsc.2009.05.004
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks_reciprocal(const std::vector<std::vector<Tp>> &a,
                             const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = std::max(lenA, lenB);
    // original version: compute a(x, x^(-N)) b(x, x^(-N))
    // modified version (this version): compute x^(2N-2) a(x^(-1), x^N) b(x^(-1), x^N)
    // ab0 = a(x, x^N) b(x, x^N)
    auto ab0 = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));

    // returns x^(N-1) a(x^(-1), y)
    auto make_reciprocal = [](const std::vector<std::vector<Tp>> &a, int N) {
        std::vector<std::vector<Tp>> b(a.size());
        for (int i = 0; i < (int)a.size(); ++i) {
            b[i] = a[i];
            b[i].resize(N);
            std::reverse(b[i].begin(), b[i].end());
        }
        return b;
    };

    // ab1 = x^(2N-2) a(x^(-1), x^N) b(x^(-1), x^N)
    auto ab1 =
        convolution(pack_2d_ks(make_reciprocal(a, N), N), pack_2d_ks(make_reciprocal(b, N), N));
    std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(N * 2 - 1));

    /*
        example:
                 let pack: a(x,y) |-> a(x,x^3)
                      rec: a(x,y) |-> x^2 a(x^(-1),y)

                 a = 1 + 2*x + 4*x^2 + (3 + 4*x)*y
                 b = 2 + 4*x + 5*x^2 + (5 + 8*x)*y
           pack(a) = 1 + 2*x + 4*x^2 + 3*x^3 + 4*x^4
           pack(b) = 2 + 4*x + 5*x^2 + 5*x^3 + 8*x^4
      pack(rec(a)) = 4 + 2*x +   x^2 + (0 + 4*x + 3*x^2)*y
      pack(rec(b)) = 5 + 4*x + 2*x^2 + (0 + 8*x + 5*x^2)*y
                ab = 2 + 8*x + 21*x^2 + 26*x^3 + 20*x^4 +
                                       (11     + 38*x   + 67*x^2 +  52*x^3)*y +
                                                                 + (15        + 44*x   + 32*x^2)*y^2
    pack(a)pack(b) = 2 + 8*x + 21*x^2 + 37*x^3 + 58*x^4 + 67*x^5 +  67*x^6    + 44*x^7 + 32*x^8
    pack(rec(a))pack(rec(b))
                   = 20 + 26*x + 21*x^2 + 8*x^3 + 2*x^4 +
                                        + ... (overlap)
    */

    // restore ab[0]
    for (int i = 0; i < N; ++i) ab[0][i] = ab0[i];
    // ab1[0] = [x^0](x^(2N - 2) a(x^(-1), x^N) b(x^(-1), x^N))
    for (int i = 0; i < N; ++i) ab[0][(N - 1) * 2 - i] = ab1[i];
    // restore ab[1..] by subtracting the overlaped coefficients
    for (int i = 1; i < (int)(a.size() + b.size() - 1); ++i) {
        // TODO: remove redundant assignment/subtraction
        for (int j = 0; j < N * 2 - 1; ++j) {
            ab0[(i - 1) * N + j] -= ab[i - 1][j];
            ab1[(i - 1) * N + j] -= ab[i - 1][(N - 1) * 2 - j];
        }
        for (int j = 0; j < N; ++j) ab[i][j] = ab0[i * N + j];
        for (int j = 0; j < N; ++j) ab[i][(N - 1) * 2 - j] = ab1[i * N + j];
    }
    for (int i = 0; i < (int)(a.size() + b.size() - 1); ++i) ab[i].resize(lenA + lenB - 1);
    return ab;
}

// see:
// [1]: David Harvey. Faster polynomial multiplication via multipoint Kronecker substitution.
//      https://doi.org/10.1016/j.jsc.2009.05.004
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks_negated(const std::vector<std::vector<Tp>> &a,
                          const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = std::max(lenA, lenB);
    // ab0 = a(x, x^N) b(x, x^N)
    const auto ab0 = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));

    // returns a(x, -y)
    auto make_negated = [](const std::vector<std::vector<Tp>> &a) {
        auto b = a;
        for (int i = 1; i < (int)b.size(); i += 2)
            for (int j = 0; j < (int)b[i].size(); ++j) b[i][j] = -b[i][j];
        return b;
    };

    // ab1 = a(x, -x^N) b(x, -x^N)
    const auto ab1 = convolution(pack_2d_ks(make_negated(a), N), pack_2d_ks(make_negated(b), N));

    /*
        example:
                 let pack: a(x,y) |-> a(x,x^3)

                 a = 1 + 2*x + 4*x^2 + (3 + 4*x)*y
                 b = 2 + 4*x + 5*x^2 + (5 + 8*x)*y
           pack(a) = 1 + 2*x + 4*x^2 + 3*x^3 + 4*x^4
           pack(b) = 2 + 4*x + 5*x^2 + 5*x^3 + 8*x^4
     pack(a(x,-y)) = 1 + 2*x + 4*x^2 + -(3*x^3 + 4*x^4)
     pack(b(x,-y)) = 2 + 4*x + 5*x^2 + -(5*x^3 + 8*x^4)
                ab = 2 + 8*x + 21*x^2 + 26*x^3 + 20*x^4 +
                                       (11     + 38*x   + 67*x^2 +  52*x^3)*y +
                                                                 + (15        + 44*x   + 32*x^2)*y^2
    pack(a)pack(b) = 2 + 8*x + 21*x^2 + 37*x^3 + 58*x^4 + 67*x^5 +  67*x^6    + 44*x^7 + 32*x^8
                                       (26+11)   (20+38) ...
    pack(a(x,-y))pack(b(x,-y))
                   = 2 + 8*x + 21*x^2 + 15*x^3 + -18*x^4 + ...
                                       (26-11)   (20-38) ...
    */

    std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(lenA + lenB - 1));
    for (int i = 0; i < (int)(a.size() + b.size() - 1); ++i) {
        if (i & 1) {
            for (int j = 0; j < lenA + lenB - 1; ++j)
                ab[i][j] = (ab0[i * N + j] - ab1[i * N + j]).div_by_2();
        } else {
            for (int j = 0; j < lenA + lenB - 1; ++j)
                ab[i][j] = (ab0[i * N + j] + ab1[i * N + j]).div_by_2();
        }
    }
    return ab;
}

#line 2 "ks.hpp"

#line 2 "fft.hpp"

#include <algorithm>
#include <cassert>
#include <iterator>
#include <memory>
#include <vector>

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 "ks.hpp"

// returns max[0 <= j < a.size()]{a[j].size()}
template<typename Tp> inline int max_len_x_ks(const std::vector<std::vector<Tp>> &a) {
    int len = -1;
    for (int i = 0; i < (int)a.size(); ++i) len = std::max<int>(len, a[i].size());
    return len;
}

// returns a(x, x^N) where a(x,y) in R[x][y]
template<typename Tp>
inline std::vector<Tp> pack_2d_ks(const std::vector<std::vector<Tp>> &a, int N) {
    assert(N > 0);
    // y |-> x^N
    std::vector<Tp> b(N * a.size());
    for (int i = 0; i < (int)a.size(); ++i)
        for (int j = 0; j < (int)a[i].size(); ++j) b[i * N + j] += a[i][j]; // if N < a[i].size()
    return b;
}

template<typename Tp>
inline std::vector<std::vector<Tp>> unpack_2d_ks(const std::vector<Tp> &a, int N) {
    assert(N > 0);
    // x^N |-> y
    std::vector<std::vector<Tp>> b((a.size() + (N - 1)) / N, std::vector<Tp>(N));
    for (int i = 0; i < (int)((a.size() + N - 1) / N); ++i)
        for (int j = 0; j < N && i * N + j < (int)a.size(); ++j) b[i][j] = a[i * N + j];
    return b;
}

template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_naive(const std::vector<std::vector<Tp>> &a, const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = lenA + lenB - 1;
    std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(N));
    for (int i = 0; i < (int)a.size(); ++i)
        for (int j = 0; j < (int)b.size(); ++j) {
            const auto aibj = convolution(a[i], b[j]);
            for (int k = 0; k < (int)aibj.size(); ++k) ab[i + j][k] += aibj[k];
        }
    return ab;
}

// standard Kronecker substitution
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks(const std::vector<std::vector<Tp>> &a, const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = lenA + lenB - 1;
    auto ab     = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));
    ab.resize((a.size() + b.size() - 1) * N);
    return unpack_2d_ks(ab, N);
}

// see:
// [1]: David Harvey. Faster polynomial multiplication via multipoint Kronecker substitution.
//      https://doi.org/10.1016/j.jsc.2009.05.004
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks_reciprocal(const std::vector<std::vector<Tp>> &a,
                             const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = std::max(lenA, lenB);
    // original version: compute a(x, x^(-N)) b(x, x^(-N))
    // modified version (this version): compute x^(2N-2) a(x^(-1), x^N) b(x^(-1), x^N)
    // ab0 = a(x, x^N) b(x, x^N)
    auto ab0 = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));

    // returns x^(N-1) a(x^(-1), y)
    auto make_reciprocal = [](const std::vector<std::vector<Tp>> &a, int N) {
        std::vector<std::vector<Tp>> b(a.size());
        for (int i = 0; i < (int)a.size(); ++i) {
            b[i] = a[i];
            b[i].resize(N);
            std::reverse(b[i].begin(), b[i].end());
        }
        return b;
    };

    // ab1 = x^(2N-2) a(x^(-1), x^N) b(x^(-1), x^N)
    auto ab1 =
        convolution(pack_2d_ks(make_reciprocal(a, N), N), pack_2d_ks(make_reciprocal(b, N), N));
    std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(N * 2 - 1));

    /*
        example:
                 let pack: a(x,y) |-> a(x,x^3)
                      rec: a(x,y) |-> x^2 a(x^(-1),y)

                 a = 1 + 2*x + 4*x^2 + (3 + 4*x)*y
                 b = 2 + 4*x + 5*x^2 + (5 + 8*x)*y
           pack(a) = 1 + 2*x + 4*x^2 + 3*x^3 + 4*x^4
           pack(b) = 2 + 4*x + 5*x^2 + 5*x^3 + 8*x^4
      pack(rec(a)) = 4 + 2*x +   x^2 + (0 + 4*x + 3*x^2)*y
      pack(rec(b)) = 5 + 4*x + 2*x^2 + (0 + 8*x + 5*x^2)*y
                ab = 2 + 8*x + 21*x^2 + 26*x^3 + 20*x^4 +
                                       (11     + 38*x   + 67*x^2 +  52*x^3)*y +
                                                                 + (15        + 44*x   + 32*x^2)*y^2
    pack(a)pack(b) = 2 + 8*x + 21*x^2 + 37*x^3 + 58*x^4 + 67*x^5 +  67*x^6    + 44*x^7 + 32*x^8
    pack(rec(a))pack(rec(b))
                   = 20 + 26*x + 21*x^2 + 8*x^3 + 2*x^4 +
                                        + ... (overlap)
    */

    // restore ab[0]
    for (int i = 0; i < N; ++i) ab[0][i] = ab0[i];
    // ab1[0] = [x^0](x^(2N - 2) a(x^(-1), x^N) b(x^(-1), x^N))
    for (int i = 0; i < N; ++i) ab[0][(N - 1) * 2 - i] = ab1[i];
    // restore ab[1..] by subtracting the overlaped coefficients
    for (int i = 1; i < (int)(a.size() + b.size() - 1); ++i) {
        // TODO: remove redundant assignment/subtraction
        for (int j = 0; j < N * 2 - 1; ++j) {
            ab0[(i - 1) * N + j] -= ab[i - 1][j];
            ab1[(i - 1) * N + j] -= ab[i - 1][(N - 1) * 2 - j];
        }
        for (int j = 0; j < N; ++j) ab[i][j] = ab0[i * N + j];
        for (int j = 0; j < N; ++j) ab[i][(N - 1) * 2 - j] = ab1[i * N + j];
    }
    for (int i = 0; i < (int)(a.size() + b.size() - 1); ++i) ab[i].resize(lenA + lenB - 1);
    return ab;
}

// see:
// [1]: David Harvey. Faster polynomial multiplication via multipoint Kronecker substitution.
//      https://doi.org/10.1016/j.jsc.2009.05.004
template<typename Tp> inline std::vector<std::vector<Tp>>
convolution_2d_ks_negated(const std::vector<std::vector<Tp>> &a,
                          const std::vector<std::vector<Tp>> &b) {
    if (a.empty() || b.empty()) return {};
    const int lenA = max_len_x_ks(a);
    const int lenB = max_len_x_ks(b);
    if (lenA == 0 || lenB == 0) return std::vector<std::vector<Tp>>(a.size() + b.size() - 1);
    const int N = std::max(lenA, lenB);
    // ab0 = a(x, x^N) b(x, x^N)
    const auto ab0 = convolution(pack_2d_ks(a, N), pack_2d_ks(b, N));

    // returns a(x, -y)
    auto make_negated = [](const std::vector<std::vector<Tp>> &a) {
        auto b = a;
        for (int i = 1; i < (int)b.size(); i += 2)
            for (int j = 0; j < (int)b[i].size(); ++j) b[i][j] = -b[i][j];
        return b;
    };

    // ab1 = a(x, -x^N) b(x, -x^N)
    const auto ab1 = convolution(pack_2d_ks(make_negated(a), N), pack_2d_ks(make_negated(b), N));

    /*
        example:
                 let pack: a(x,y) |-> a(x,x^3)

                 a = 1 + 2*x + 4*x^2 + (3 + 4*x)*y
                 b = 2 + 4*x + 5*x^2 + (5 + 8*x)*y
           pack(a) = 1 + 2*x + 4*x^2 + 3*x^3 + 4*x^4
           pack(b) = 2 + 4*x + 5*x^2 + 5*x^3 + 8*x^4
     pack(a(x,-y)) = 1 + 2*x + 4*x^2 + -(3*x^3 + 4*x^4)
     pack(b(x,-y)) = 2 + 4*x + 5*x^2 + -(5*x^3 + 8*x^4)
                ab = 2 + 8*x + 21*x^2 + 26*x^3 + 20*x^4 +
                                       (11     + 38*x   + 67*x^2 +  52*x^3)*y +
                                                                 + (15        + 44*x   + 32*x^2)*y^2
    pack(a)pack(b) = 2 + 8*x + 21*x^2 + 37*x^3 + 58*x^4 + 67*x^5 +  67*x^6    + 44*x^7 + 32*x^8
                                       (26+11)   (20+38) ...
    pack(a(x,-y))pack(b(x,-y))
                   = 2 + 8*x + 21*x^2 + 15*x^3 + -18*x^4 + ...
                                       (26-11)   (20-38) ...
    */

    std::vector<std::vector<Tp>> ab(a.size() + b.size() - 1, std::vector<Tp>(lenA + lenB - 1));
    for (int i = 0; i < (int)(a.size() + b.size() - 1); ++i) {
        if (i & 1) {
            for (int j = 0; j < lenA + lenB - 1; ++j)
                ab[i][j] = (ab0[i * N + j] - ab1[i * N + j]).div_by_2();
        } else {
            for (int j = 0; j < lenA + lenB - 1; ++j)
                ab[i][j] = (ab0[i * N + j] + ab1[i * N + j]).div_by_2();
        }
    }
    return ab;
}

Back to top page