This documentation is automatically generated by online-judge-tools/verification-helper
#include "library/linear_algebra/circulant_matrix.hpp"
積を高速に計算できる行列の一種。
列 $a = ( a _ 0 , \ldots , a _ { n - 1 } )$ に対して、$n\times n$ 行列 $\mathrm{Cir}(a _ 0, \ldots, a _ {n - 1})$ を以下で定義すると、これは巡回行列である。
\[\mathrm{Cir}(a _ 0, \ldots, a _ {n - 1}) \coloneqq (a _ {i - j \bmod n}) _ {0\leq i,j\lt n} =\begin{pmatrix} a _ 0 & a _ {n - 1} & a _ {n - 3} & \cdots & a _ 1\\ a _ 1 & a _ 0 & a _ {n - 2} & \cdots & a _ 2\\ a _ 2 & a _ 1 & a _ 0 & \cdots & a _ 3\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a _ {n - 1} & a _ {n - 2} & a _ {n - 3} & \cdots & a _ 0 \end{pmatrix}\]$2$ つの巡回行列 $A\coloneqq\mathrm{Cir}(a _ 0, \ldots, a _ {n - 1}),\; B \coloneqq \mathrm{Cir}(b _ 0, \ldots, b _ {n - 1})$ に対して $C\coloneqq AB$ とすれば、次が成り立つ。
\[C _ {i,j} = \sum _ {k = 0} ^ {n - 1} a _ {i - k \bmod n} b _ {k - j \bmod n} = \sum _ {k+l\equiv i-j\pmod{n}} a _ k b _ l\]右辺は $i-j\bmod n$ の値にのみ依存するから、$\displaystyle c _ i \coloneqq \sum _ {k+l\equiv i\pmod{n}} a _ k b _ l$ に対して $C=\mathrm{Cir}(c _ 0, \ldots, c _ {n - 1})$ である。これは、巡回行列の積はまた巡回行列となることを表している。更に、$c$ は $a,b$ の巡回畳み込みであるから、FFT などを用いることで $C$ は $\Theta(n\log n)$ 時間で計算可能である。
#ifndef SUISEN_CIRCULANT_MATRIX
#define SUISEN_CIRCULANT_MATRIX
#include <cassert>
#include <iostream>
#include <vector>
namespace suisen {
template <typename T>
struct CirculantMatrix {
using value_type = T;
using convolution_t = std::vector<value_type>(*)(const std::vector<value_type>&, const std::vector<value_type>&);
// empty matrix
CirculantMatrix() : CirculantMatrix(std::vector<value_type>{}) {}
/**
* +- -+
* | a[0] a[4] a[3] a[2] a[1] |
* | a[1] a[0] a[4] a[3] a[2] |
* | a[2] a[1] a[0] a[4] a[3] |
* | a[3] a[2] a[1] a[0] a[4] |
* | a[4] a[3] a[2] a[1] a[0] |
* +- -+
*/
explicit CirculantMatrix(const std::vector<value_type>& a) : _dat(a) {}
static void set_multiplication(convolution_t multiplication) {
convolve = multiplication;
}
static CirculantMatrix<value_type> e0(int n, const value_type& zero = value_type{ 0 }) {
return CirculantMatrix<value_type>{ std::vector<value_type>(n, zero) };
}
static CirculantMatrix<value_type> e1(int n, const value_type& zero = value_type{ 0 }, const value_type& one = value_type{ 1 }) {
auto dat = std::vector<value_type>(n, zero);
dat[0] = one;
return CirculantMatrix<value_type>{ dat };
}
int size() const {
return _dat.size();
}
value_type get(int i, int j) const {
const int n = size();
int k = i - j;
if (k < 0) k += n;
return _dat[k];
}
value_type operator[](const std::pair<int, int> &p) const {
return get(p.first, p.second);
}
friend CirculantMatrix<value_type> operator+(const CirculantMatrix<value_type>& mat) {
return mat;
}
friend CirculantMatrix<value_type> operator-(const CirculantMatrix<value_type>& mat) {
const int n = mat.size();
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = -mat._dat[i];
return CirculantMatrix<value_type> { std::move(res) };
}
friend CirculantMatrix<value_type> operator+(const CirculantMatrix<value_type>& lhs, const CirculantMatrix<value_type>& rhs) {
const int n = lhs.size();
assert(n == int(rhs.size()));
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = lhs._dat[i] + rhs._dat[i];
return CirculantMatrix<value_type> { std::move(res) };
}
friend CirculantMatrix<value_type> operator-(const CirculantMatrix<value_type>& lhs, const CirculantMatrix<value_type>& rhs) {
const int n = lhs.size();
assert(n == int(rhs.size()));
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = lhs._dat[i] - rhs._dat[i];
return CirculantMatrix<value_type> { std::move(res) };
}
friend CirculantMatrix<value_type> operator*(const CirculantMatrix<value_type>& lhs, const CirculantMatrix<value_type>& rhs) {
const int n = lhs.size();
assert(n == int(rhs.size()));
std::vector<value_type> res = convolve(lhs._dat, rhs._dat);
for (int i = n; i < int(res.size()); ++i) res[i - n] += res[i];
res.resize(n);
return CirculantMatrix<value_type> { std::move(res) };
}
friend std::vector<value_type> operator*(const CirculantMatrix<value_type>& mat, const std::vector<value_type>& x) {
return std::move((mat * CirculantMatrix<value_type> { x })._dat);
}
friend CirculantMatrix<value_type> operator*(const CirculantMatrix<value_type>& mat, const value_type& coef) {
const int n = mat.size();
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = coef * mat._dat[i];
return CirculantMatrix<value_type> { res };
}
friend CirculantMatrix<value_type> operator*(const value_type& coef, const CirculantMatrix<value_type>& mat) {
return mat * coef;
}
CirculantMatrix<value_type>& operator+=(const CirculantMatrix<value_type>& rhs) {
return *this = *this + rhs;
}
CirculantMatrix<value_type>& operator-=(const CirculantMatrix<value_type>& rhs) {
return *this = *this - rhs;
}
CirculantMatrix<value_type>& operator*=(const CirculantMatrix<value_type>& rhs) {
return *this = *this * rhs;
}
CirculantMatrix<value_type>& operator*=(const value_type& coef) {
return *this = *this * coef;
}
CirculantMatrix<value_type> pow(long long b) {
auto res = CirculantMatrix<value_type>::e1(size());
for (auto p = *this; b; b >>= 1) {
if (b & 1) res *= p;
p *= p;
}
return res;
}
private:
static inline convolution_t convolve{
[](const auto&, const auto&) {
std::cerr << "convolution function is not available." << std::endl;
assert(false);
return std::vector<value_type>{};
}
};
std::vector<value_type> _dat;
};
} // namespace suisen
#endif // SUISEN_CIRCULANT_MATRIX
#line 1 "library/linear_algebra/circulant_matrix.hpp"
#include <cassert>
#include <iostream>
#include <vector>
namespace suisen {
template <typename T>
struct CirculantMatrix {
using value_type = T;
using convolution_t = std::vector<value_type>(*)(const std::vector<value_type>&, const std::vector<value_type>&);
// empty matrix
CirculantMatrix() : CirculantMatrix(std::vector<value_type>{}) {}
/**
* +- -+
* | a[0] a[4] a[3] a[2] a[1] |
* | a[1] a[0] a[4] a[3] a[2] |
* | a[2] a[1] a[0] a[4] a[3] |
* | a[3] a[2] a[1] a[0] a[4] |
* | a[4] a[3] a[2] a[1] a[0] |
* +- -+
*/
explicit CirculantMatrix(const std::vector<value_type>& a) : _dat(a) {}
static void set_multiplication(convolution_t multiplication) {
convolve = multiplication;
}
static CirculantMatrix<value_type> e0(int n, const value_type& zero = value_type{ 0 }) {
return CirculantMatrix<value_type>{ std::vector<value_type>(n, zero) };
}
static CirculantMatrix<value_type> e1(int n, const value_type& zero = value_type{ 0 }, const value_type& one = value_type{ 1 }) {
auto dat = std::vector<value_type>(n, zero);
dat[0] = one;
return CirculantMatrix<value_type>{ dat };
}
int size() const {
return _dat.size();
}
value_type get(int i, int j) const {
const int n = size();
int k = i - j;
if (k < 0) k += n;
return _dat[k];
}
value_type operator[](const std::pair<int, int> &p) const {
return get(p.first, p.second);
}
friend CirculantMatrix<value_type> operator+(const CirculantMatrix<value_type>& mat) {
return mat;
}
friend CirculantMatrix<value_type> operator-(const CirculantMatrix<value_type>& mat) {
const int n = mat.size();
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = -mat._dat[i];
return CirculantMatrix<value_type> { std::move(res) };
}
friend CirculantMatrix<value_type> operator+(const CirculantMatrix<value_type>& lhs, const CirculantMatrix<value_type>& rhs) {
const int n = lhs.size();
assert(n == int(rhs.size()));
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = lhs._dat[i] + rhs._dat[i];
return CirculantMatrix<value_type> { std::move(res) };
}
friend CirculantMatrix<value_type> operator-(const CirculantMatrix<value_type>& lhs, const CirculantMatrix<value_type>& rhs) {
const int n = lhs.size();
assert(n == int(rhs.size()));
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = lhs._dat[i] - rhs._dat[i];
return CirculantMatrix<value_type> { std::move(res) };
}
friend CirculantMatrix<value_type> operator*(const CirculantMatrix<value_type>& lhs, const CirculantMatrix<value_type>& rhs) {
const int n = lhs.size();
assert(n == int(rhs.size()));
std::vector<value_type> res = convolve(lhs._dat, rhs._dat);
for (int i = n; i < int(res.size()); ++i) res[i - n] += res[i];
res.resize(n);
return CirculantMatrix<value_type> { std::move(res) };
}
friend std::vector<value_type> operator*(const CirculantMatrix<value_type>& mat, const std::vector<value_type>& x) {
return std::move((mat * CirculantMatrix<value_type> { x })._dat);
}
friend CirculantMatrix<value_type> operator*(const CirculantMatrix<value_type>& mat, const value_type& coef) {
const int n = mat.size();
std::vector<value_type> res(n);
for (int i = 0; i < n; ++i) res[i] = coef * mat._dat[i];
return CirculantMatrix<value_type> { res };
}
friend CirculantMatrix<value_type> operator*(const value_type& coef, const CirculantMatrix<value_type>& mat) {
return mat * coef;
}
CirculantMatrix<value_type>& operator+=(const CirculantMatrix<value_type>& rhs) {
return *this = *this + rhs;
}
CirculantMatrix<value_type>& operator-=(const CirculantMatrix<value_type>& rhs) {
return *this = *this - rhs;
}
CirculantMatrix<value_type>& operator*=(const CirculantMatrix<value_type>& rhs) {
return *this = *this * rhs;
}
CirculantMatrix<value_type>& operator*=(const value_type& coef) {
return *this = *this * coef;
}
CirculantMatrix<value_type> pow(long long b) {
auto res = CirculantMatrix<value_type>::e1(size());
for (auto p = *this; b; b >>= 1) {
if (b & 1) res *= p;
p *= p;
}
return res;
}
private:
static inline convolution_t convolve{
[](const auto&, const auto&) {
std::cerr << "convolution function is not available." << std::endl;
assert(false);
return std::vector<value_type>{};
}
};
std::vector<value_type> _dat;
};
} // namespace suisen