結果
| 問題 | No.444 旨味の相乗効果 |
| コンテスト | |
| ユーザー |
 askr58
|
| 提出日時 | 2026-10-27 02:54:17 |
| 言語 | C++23 (gcc 15.3.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 64 ms / 2,500 ms |
| + 301µs | |
| コード長 | 31,474 bytes |
| 記録 | |
| コンパイル時間 | 4,226 ms |
| コンパイル使用メモリ | 343,744 KB |
| 実行使用メモリ | 9,920 KB |
| 最終ジャッジ日時 | 2026-10-27 02:54:24 |
| 合計ジャッジ時間 | 6,975 ms |
|
ジャッジサーバーID （参考情報） |
judge3_0 / judge1_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 5 |
| other | AC * 23 |
コンパイルメッセージ
In file included from /home/linuxbrew/.linuxbrew/Cellar/gcc@15/15.3.0/include/c++/15/vector:68,
from /home/linuxbrew/.linuxbrew/Cellar/gcc@15/15.3.0/include/c++/15/bits/random.h:34,
from /home/linuxbrew/.linuxbrew/Cellar/gcc@15/15.3.0/include/c++/15/random:50,
from main.cpp:19:
In member function 'constexpr void std::_Vector_base<_Tp, _Alloc>::_Vector_impl_data::_M_copy_data(const std::_Vector_base<_Tp, _Alloc>::_Vector_impl_data&) [with _Tp = atcoder::static_modint<1000000007>; _Alloc = std::allocator<atcoder::static_modint<1000000007> >]',
inlined from 'constexpr void std::_Vector_base<_Tp, _Alloc>::_Vector_impl_data::_M_swap_data(std::_Vector_base<_Tp, _Alloc>::_Vector_impl_data&) [with _Tp = atcoder::static_modint<1000000007>; _Alloc = std::allocator<atcoder::static_modint<1000000007> >]' at /home/linuxbrew/.linuxbrew/Cellar/gcc@15/15.3.0/include/c++/15/bits/stl_vector.h:133:22,
inlined from 'constexpr void std::vector< <template-parameter-1-1>, <template-parameter-1-2> >::_M_move_assign(std::vector< <template-parameter-1-1>, <template-parameter-1-2> >&&, std::true_type) [with _Tp = atcoder::static_modint<1000000007>; _Alloc = std::allocator<atcoder::static_modint<1000000007> >]' at /home/linuxbrew/.linuxbrew/Cellar/gcc@15/15.3.0/include/c++/15/bits/stl_vector.h:2268:28,
inlined from 'constexpr std::vector< <template-parameter-1-1>, <template-parameter-1-2> >& std::vector< <template-parameter-1-1>, <template-parameter-1-2> >::operator=(std::vector< <template-parameter-1-1>, <template-parameter-1-2> >&&) [with _Tp = atcoder::static_modint<1000000007>; _Alloc = std::allocator<atcoder::static_modint<1000000007> >]' at /home/linuxbrew/.linuxbrew/Cellar/gcc@15/15.3.0/include/c++/15/bits/stl_vector.h:838:16,
inlined from 'constexpr Matrix<atcoder::static_modint<1000000007> >::GaussData& Matrix<atcoder::static_modint<1000000007> >::GaussData::operator=(Matrix<atcoder::static_modint<1000000007> >::GaussD
ソースコード
// SPDX-License-Identifier: CC0-1.0
//
// [English]
// To the extent possible under law, the author(s) have dedicated all copyright
// and related and neighboring rights to this software to the public domain
// worldwide. This software is distributed without any warranty.
// See <https://creativecommons.org/publicdomain/zero/1.0/> for details.
//
// [日本語]
// 法律上可能な範囲で、著作者はこのソフトウェアに関する著作権および
// 著作隣接権を含むすべての関連する権利を、全世界でパブリックドメインに
// 提供しています。本ソフトウェアは一切の保証なく配布されます。
// 自由にコピー・改変・再配布・商用利用できます（クレジット表記は不要です）。
// 詳細は <https://creativecommons.org/publicdomain/zero/1.0/deed.ja> を参照してください。
// ※ 英語の文面および CC0 1.0 の正式な条文（英語）が優先されます。
#include <iostream>
#include <bitset>
#include <random>
#include <chrono>
#include <iomanip>
#include <set>
#include <map>
#include <queue>
#include <deque>
#include <string>
#include <stack>
#include <ranges>
#include <algorithm>
#include <vector>
using namespace std;
using ll=long long;
#include <atcoder/all>
using mint=atcoder::modint1000000007;
ostream& operator<<(ostream& os,const mint& x){
os<<x.val();
return os;
}
istream& operator>>(istream& is,mint& x){
int t;
is>>t;
x=t;
return is;
}
template <typename S,typename T>
ostream& operator<<(ostream& os,const pair<S,T>& p);
template <typename S,typename T>
istream& operator>>(istream& is,pair<S,T>& p);
template <typename T,size_t n>
ostream& operator<<(ostream& os,const array<T,n>& arr);
template <typename T,size_t n>
istream& operator>>(istream& is,array<T,n>& arr);
template <typename T>
ostream& operator<<(ostream& os,const vector<T>& vec);
template <typename T>
istream& operator>>(istream& is,vector<T>& vec);
template <typename S,typename T>
ostream& operator<<(ostream& os,const pair<S,T>& p){
os<<p.first<<" "<<p.second;
return os;
}
template <typename S,typename T>
istream& operator>>(istream& is,pair<S,T>& p){
is>>p.first>>p.second;
return is;
}
template <typename T,size_t n>
ostream& operator<<(ostream& os,const array<T,n>& arr){
for(int i=0;i<n;i++)os<<arr[i]<<(i+1==n?"":" ");
return os;
}
template <typename T,size_t n>
istream& operator>>(istream& is,array<T,n>& arr){
for(int i=0;i<n;i++)is>>arr[i];
return is;
}
template <typename T>
ostream& operator<<(ostream& os,const vector<T>& vec){
for(int i=0;i<(int)vec.size();i++)os<<vec[i]<<(i+1==(int)vec.size()?"":" ");
return os;
}
template <typename T>
istream& operator>>(istream& is,vector<T>& vec){
for(int i=0;i<(int)vec.size();i++)is>>vec[i];
return is;
}
template<class... Vecs>
void input_vec(Vecs&... vs) {
const auto n = get<0>(tie(vs...)).size();
for (size_t i = 0; i < n; ++i)
((cin >> vs[i]), ...);
}
template <class... Vecs>
void output_vec(const Vecs&... vs) {
const auto n = get<0>(tie(vs...)).size();
for (size_t i = 0; i < n; ++i) {
bool first = true;
(((cout << (exchange(first, false) ? "" : " ") << vs[i])), ...);
cout << endl;
}
}
template <typename T>
vector<T> make_unique(vector<T> vec){
ranges::sort(vec);
vec.erase(unique(vec.begin(),vec.end()),vec.end());
return vec;
}
vector<int> Iota(int n,int s=0){
vector<int> res(n);
iota(res.begin(),res.end(),s);
return res;
}
vector<ll> Iotall(int n,ll s=0){
vector<ll> res(n);
iota(res.begin(),res.end(),s);
return res;
}
void YESNO(bool f){
if(f)cout<<"Yes"<<endl;
else cout<<"No"<<endl;
}
using vl=vector<ll>;
using vvl=vector<vector<ll>>;
using vvvl=vector<vector<vector<ll>>>;
using vi=vector<int>;
using vvi=vector<vector<int>>;
using vvvi=vector<vector<vector<int>>>;
// generated by ChatGPT outside of contest time.
template<class T = mint> class Vector;
template<class T> struct LinearSolution;
// ============================================================
// Matrix
// ============================================================
template<class T = mint>
class Matrix {
public:
using value_type = T;
private:
struct GaussData {
// row-major
std::vector<T> rref;
// R = P A を満たす P
// h x h, row-major
std::vector<T> transform;
bool has_transform=false;
std::vector<int> pivot_cols;
int rank = 0;
// 正方行列の場合のみ意味を持つ
T det = T(0);
};
protected:
int h_ = 0;
int w_ = 0;
// row-major:
// a_[i * w_ + j] = A[i][j]
std::vector<T> a_;
private:
mutable std::optional<GaussData> gauss_;
void require_gauss() const {
if (!gauss_) {
throw std::logic_error(
"Matrix: call gauss_jordan() first"
);
}
}
public:
// --------------------------------------------------------
// constructors
// --------------------------------------------------------
Matrix() = default;
Matrix(int h, int w)
: h_(h),
w_(w),
a_((size_t)h * w, T(0)) {
assert(h >= 0 && w >= 0);
}
Matrix(int h, int w, const T& value)
: h_(h),
w_(w),
a_((size_t)h * w, value) {
assert(h >= 0 && w >= 0);
}
explicit Matrix(
const std::vector<std::vector<T>>& a
)
: h_((int)a.size()),
w_(a.empty() ? 0 : (int)a[0].size()),
a_((size_t)h_ * w_) {
for (int i = 0; i < h_; ++i) {
assert((int)a[i].size() == w_);
std::copy(
a[i].begin(),
a[i].end(),
(*this)[i]
);
}
}
Matrix(
std::initializer_list<
std::initializer_list<T>
> init
) {
h_ = (int)init.size();
w_ = h_
? (int)init.begin()->size()
: 0;
a_.resize((size_t)h_ * w_);
int i = 0;
for (auto row : init) {
assert((int)row.size() == w_);
std::copy(
row.begin(),
row.end(),
(*this)[i++]
);
}
}
// --------------------------------------------------------
// access
// --------------------------------------------------------
int rows() const {
return h_;
}
int cols() const {
return w_;
}
bool empty() const {
return h_ == 0 || w_ == 0;
}
// A[i] = pointer to the beginning of row i
T* operator[](int i) {
return a_.data() + (size_t)i * w_;
}
const T* operator[](int i) const {
return a_.data() + (size_t)i * w_;
}
T& operator()(int i, int j) {
return a_[(size_t)i * w_ + j];
}
const T& operator()(int i, int j) const {
return a_[(size_t)i * w_ + j];
}
T* data() {
return a_.data();
}
const T* data() const {
return a_.data();
}
const std::vector<std::vector<T>> to_vector() const{
std::vector<std::vector<T>> res(h_,std::vector<T>(w_));
for(int i=0;i<h_;i++)for(int j=0;j<w_;j++)res[i][j]=a_[i*w_+j];
return res;
}
// --------------------------------------------------------
// basic matrices
// --------------------------------------------------------
static Matrix identity(int n) {
Matrix I(n, n);
for (int i = 0; i < n; ++i) {
I[i][i] = T(1);
}
return I;
}
Matrix transpose() const {
Matrix res(w_, h_);
for (int i = 0; i < h_; ++i) {
const T* ai = (*this)[i];
for (int j = 0; j < w_; ++j) {
res[j][i] = ai[j];
}
}
return res;
}
void swap_rows(int i, int j) {
if (i == j) {
return;
}
std::swap_ranges(
(*this)[i],
(*this)[i] + w_,
(*this)[j]
);
}
// --------------------------------------------------------
// arithmetic
// --------------------------------------------------------
Matrix& operator+=(const Matrix& rhs) {
assert(
h_ == rhs.h_ &&
w_ == rhs.w_
);
for (size_t i = 0; i < a_.size(); ++i) {
a_[i] += rhs.a_[i];
}
return *this;
}
Matrix& operator-=(const Matrix& rhs) {
assert(
h_ == rhs.h_ &&
w_ == rhs.w_
);
for (size_t i = 0; i < a_.size(); ++i) {
a_[i] -= rhs.a_[i];
}
return *this;
}
Matrix& operator*=(const T& c) {
for (T& x : a_) {
x *= c;
}
return *this;
}
Matrix& operator/=(const T& c) {
T inv = T(1) / c;
for (T& x : a_) {
x *= inv;
}
return *this;
}
Matrix operator-() const {
Matrix res(h_, w_);
for (size_t i = 0; i < a_.size(); ++i) {
res.a_[i] = -a_[i];
}
return res;
}
// --------------------------------------------------------
// power
// --------------------------------------------------------
Matrix pow(long long e) const {
assert(h_ == w_);
assert(e >= 0);
Matrix base = *this;
Matrix res = identity(h_);
while (e) {
if (e & 1) {
res = res * base;
}
e >>= 1;
if (e) {
base = base * base;
}
}
return res;
}
// ========================================================
// Gauss-Jordan
// ========================================================
/*
R = P A
R : reduced row echelon form
P : invertible matrix representing row operations
*/
const Matrix& gauss_jordan(bool make_transform=false) const {
GaussData g;
g.rref = a_;
g.has_transform = make_transform;
if(make_transform){
g.transform.assign(
(size_t)h_ * h_,
T(0)
);
for (int i = 0; i < h_; ++i) {
g.transform[
(size_t)i * h_ + i
] = T(1);
}
}
int row = 0;
int swap_parity = 0;
T det_product = T(1);
for (
int col = 0;
col < w_ && row < h_;
++col
) {
// --------------------------------------------
// find pivot
// --------------------------------------------
int pivot = -1;
for (int i = row; i < h_; ++i) {
if (
g.rref[
(size_t)i * w_ + col
] != T(0)
) {
pivot = i;
break;
}
}
if (pivot == -1) {
continue;
}
// --------------------------------------------
// row swap
// --------------------------------------------
if (pivot != row) {
T* pr =
g.rref.data()
+ (size_t)pivot * w_;
T* rr =
g.rref.data()
+ (size_t)row * w_;
std::swap_ranges(
pr,
pr + w_,
rr
);
if(make_transform){
T* pt =
g.transform.data()
+ (size_t)pivot * h_;
T* rt =
g.transform.data()
+ (size_t)row * h_;
std::swap_ranges(
pt,
pt + h_,
rt
);
}
swap_parity ^= 1;
}
T* rr =
g.rref.data()
+ (size_t)row * w_;
T* rt =
g.transform.data()
+ (size_t)row * h_;
// --------------------------------------------
// normalize pivot
// --------------------------------------------
T p = rr[col];
det_product *= p;
T inv_p = T(1) / p;
/*
col より左は既に 0 なので更新不要。
*/
for (int j = col; j < w_; ++j) {
rr[j] *= inv_p;
}
rt=nullptr;
if(make_transform){
rt=g.transform.data()+row*h_;
for (int j = 0; j < h_; ++j) {
rt[j] *= inv_p;
}
}
// --------------------------------------------
// eliminate other rows
// --------------------------------------------
for (int i = 0; i < h_; ++i) {
if (i == row) {
continue;
}
T* ir =
g.rref.data()
+ (size_t)i * w_;
T factor = ir[col];
if (factor == T(0)) {
continue;
}
/*
pivot column は結果が 0 と分かっているので
直接代入する。
*/
ir[col] = T(0);
/*
col より左は current pivot row が全部 0。
よって更新する必要がない。
*/
for (
int j = col + 1;
j < w_;
++j
) {
ir[j]
-= factor * rr[j];
}
if(make_transform){
T* it =
g.transform.data()
+ (size_t)i * h_;
for (int j = 0; j < h_; ++j) {
it[j]
-= factor * rt[j];
}
}
}
g.pivot_cols.push_back(col);
++row;
}
g.rank = row;
// determinant
if (h_ == w_) {
if (g.rank < h_) {
g.det = T(0);
} else {
g.det =
swap_parity
? -det_product
: det_product;
}
}
gauss_ = std::move(g);
return *this;
}
// --------------------------------------------------------
// Gauss information
// --------------------------------------------------------
int rank() const {
require_gauss();
return gauss_->rank;
}
T det() const {
require_gauss();
assert(h_ == w_);
return gauss_->det;
}
std::vector<int> pivot_columns() const {
require_gauss();
return gauss_->pivot_cols;
}
std::vector<int> free_columns() const {
require_gauss();
std::vector<char> is_pivot(
w_,
false
);
for (int p : gauss_->pivot_cols) {
is_pivot[p] = true;
}
std::vector<int> res;
for (int j = 0; j < w_; ++j) {
if (!is_pivot[j]) {
res.push_back(j);
}
}
return res;
}
// --------------------------------------------------------
// RREF / transform
// --------------------------------------------------------
Matrix rref() const {
require_gauss();
Matrix res(h_, w_);
res.a_ = gauss_->rref;
return res;
}
/*
P such that
R = P A
*/
Matrix transform() const {
require_gauss();
if (!gauss_->has_transform) {
throw std::logic_error(
"Matrix: call gauss_jordan(true) "
"before transform()"
);
}
Matrix res(h_, h_);
res.a_ = gauss_->transform;
return res;
}
// --------------------------------------------------------
// inverse
// --------------------------------------------------------
std::optional<Matrix> inverse() const {
require_gauss();
if (!gauss_->has_transform) {
throw std::logic_error(
"Matrix: call gauss_jordan(true) "
"before transform()"
);
}
assert(h_ == w_);
if (gauss_->rank != h_) {
return std::nullopt;
}
/*
R = I = P A
therefore P = A^{-1}
*/
Matrix res(h_, h_);
res.a_ = gauss_->transform;
return res;
}
// --------------------------------------------------------
// linear algebra
// --------------------------------------------------------
std::vector<Vector<T>>
kernel_basis() const;
std::optional<LinearSolution<T>>
solve(const Vector<T>& b) const;
std::optional<LinearSolution<T>>
solve_direct(const Vector<T>& b) const;
// --------------------------------------------------------
// comparison
// --------------------------------------------------------
friend bool operator==(
const Matrix& lhs,
const Matrix& rhs
) {
return
lhs.h_ == rhs.h_ &&
lhs.w_ == rhs.w_ &&
lhs.a_ == rhs.a_;
}
friend bool operator!=(
const Matrix& lhs,
const Matrix& rhs
) {
return !(lhs == rhs);
}
};
// ============================================================
// Matrix arithmetic
// ============================================================
template<class T>
Matrix<T> operator+(
Matrix<T> lhs,
const Matrix<T>& rhs
) {
return lhs += rhs;
}
template<class T>
Matrix<T> operator-(
Matrix<T> lhs,
const Matrix<T>& rhs
) {
return lhs -= rhs;
}
/*
Matrix multiplication
i-k-j order.
C[i][j] += A[i][k] * B[k][j]
最内周では B[k][j], C[i][j] の両方が
contiguous access。
*/
template<class T>
Matrix<T> operator*(
const Matrix<T>& A,
const Matrix<T>& B
) {
assert(A.cols() == B.rows());
Matrix<T> C(
A.rows(),
B.cols()
);
for (int i = 0; i < A.rows(); ++i) {
const T* ai = A[i];
T* ci = C[i];
for (int k = 0; k < A.cols(); ++k) {
T x = ai[k];
if (x == T(0)) {
continue;
}
const T* bk = B[k];
for (int j = 0; j < B.cols(); ++j) {
ci[j] += x * bk[j];
}
}
}
return C;
}
// ------------------------------------------------------------
// scalar multiplication
// ------------------------------------------------------------
template<
class T,
class S,
std::enable_if_t<
std::is_constructible_v<T, S>,
int
> = 0
>
Matrix<T> operator*(
Matrix<T> A,
const S& c
) {
return A *= T(c);
}
template<
class S,
class T,
std::enable_if_t<
std::is_constructible_v<T, S>,
int
> = 0
>
Matrix<T> operator*(
const S& c,
Matrix<T> A
) {
return A *= T(c);
}
template<
class T,
class S,
std::enable_if_t<
std::is_constructible_v<T, S>,
int
> = 0
>
Matrix<T> operator/(
Matrix<T> A,
const S& c
) {
return A /= T(c);
}
// ============================================================
// Vector
// ============================================================
template<class T>
class Vector : public Matrix<T> {
using Base = Matrix<T>;
public:
using value_type = T;
Vector()
: Base(0, 1) {}
explicit Vector(int n)
: Base(n, 1) {}
Vector(
int n,
const T& value
)
: Base(n, 1, value) {}
explicit Vector(
const std::vector<T>& v
)
: Base((int)v.size(), 1) {
std::copy(
v.begin(),
v.end(),
this->a_.begin()
);
}
explicit Vector(
std::vector<T>&& v
)
: Base((int)v.size(), 1) {
this->a_ = std::move(v);
}
Vector(
std::initializer_list<T> init
)
: Base((int)init.size(), 1) {
std::copy(
init.begin(),
init.end(),
this->a_.begin()
);
}
int size() const {
return this->h_;
}
/*
Matrix::operator[] を隠す。
Matrix:
A[i] -> T*
Vector:
v[i] -> T&
*/
T& operator[](int i) {
return this->a_[i];
}
const T& operator[](int i) const {
return this->a_[i];
}
const std::vector<T> to_vector() const{
return this->a_;
}
Vector& operator+=(
const Vector& rhs
) {
assert(size() == rhs.size());
Base::operator+=(rhs);
return *this;
}
Vector& operator-=(
const Vector& rhs
) {
assert(size() == rhs.size());
Base::operator-=(rhs);
return *this;
}
Vector& operator*=(const T& c) {
Base::operator*=(c);
return *this;
}
Vector& operator/=(const T& c) {
Base::operator/=(c);
return *this;
}
Vector operator-() const {
Vector res(size());
for (int i = 0; i < size(); ++i) {
res[i] = -(*this)[i];
}
return res;
}
};
// ============================================================
// Vector arithmetic
// ============================================================
template<class T>
Vector<T> operator+(
Vector<T> lhs,
const Vector<T>& rhs
) {
return lhs += rhs;
}
template<class T>
Vector<T> operator-(
Vector<T> lhs,
const Vector<T>& rhs
) {
return lhs -= rhs;
}
template<
class T,
class S,
std::enable_if_t<
std::is_constructible_v<T, S>,
int
> = 0
>
Vector<T> operator*(
Vector<T> v,
const S& c
) {
return v *= T(c);
}
template<
class S,
class T,
std::enable_if_t<
std::is_constructible_v<T, S>,
int
> = 0
>
Vector<T> operator*(
const S& c,
Vector<T> v
) {
return v *= T(c);
}
template<
class T,
class S,
std::enable_if_t<
std::is_constructible_v<T, S>,
int
> = 0
>
Vector<T> operator/(
Vector<T> v,
const S& c
) {
return v /= T(c);
}
// ============================================================
// Matrix * Vector
// ============================================================
template<class T>
Vector<T> operator*(
const Matrix<T>& A,
const Vector<T>& x
) {
assert(A.cols() == x.size());
Vector<T> y(A.rows());
for (int i = 0; i < A.rows(); ++i) {
const T* ai = A[i];
T sum = T(0);
for (int j = 0; j < A.cols(); ++j) {
sum += ai[j] * x[j];
}
y[i] = sum;
}
return y;
}
// ============================================================
// dot product
// ============================================================
template<class T>
T dot(
const Vector<T>& a,
const Vector<T>& b
) {
assert(a.size() == b.size());
T res = T(0);
for (int i = 0; i < a.size(); ++i) {
res += a[i] * b[i];
}
return res;
}
// ============================================================
// General solution
// ============================================================
/*
Solution set:
x = particular
+ c_0 basis[0]
+ ...
+ c_{k-1} basis[k-1]
*/
template<class T>
struct LinearSolution {
Vector<T> particular;
std::vector<Vector<T>> basis;
bool unique() const {
return basis.empty();
}
int dimension() const {
return (int)basis.size();
}
Vector<T> get(
const std::vector<T>& coeffs
) const {
assert(
coeffs.size()
== basis.size()
);
Vector<T> x = particular;
for (
int i = 0;
i < (int)basis.size();
++i
) {
x += basis[i] * coeffs[i];
}
return x;
}
};
// ============================================================
// Kernel
// ============================================================
template<class T>
std::vector<Vector<T>>
Matrix<T>::kernel_basis() const {
require_gauss();
std::vector<char> is_pivot(
w_,
false
);
for (int p : gauss_->pivot_cols) {
is_pivot[p] = true;
}
std::vector<Vector<T>> basis;
for (
int free_col = 0;
free_col < w_;
++free_col
) {
if (is_pivot[free_col]) {
continue;
}
Vector<T> v(w_);
v[free_col] = T(1);
for (
int row = 0;
row < gauss_->rank;
++row
) {
int p =
gauss_->pivot_cols[row];
v[p] =
-gauss_->rref[
(size_t)row * w_
+ free_col
];
}
basis.push_back(
std::move(v)
);
}
return basis;
}
// ============================================================
// Solve Ax = b
// ============================================================
template<class T>
std::optional<LinearSolution<T>>
Matrix<T>::solve(
const Vector<T>& b
) const {
require_gauss();
if (!gauss_->has_transform) {
throw std::logic_error(
"Matrix: call gauss_jordan(true) "
"before transform()"
);
}
assert(b.size() == h_);
/*
R = P A
Ax = b
<=>
R x = P b
*/
std::vector<T> rhs(
h_,
T(0)
);
for (int i = 0; i < h_; ++i) {
const T* tr =
gauss_->transform.data()
+ (size_t)i * h_;
T sum = T(0);
for (int j = 0; j < h_; ++j) {
sum += tr[j] * b[j];
}
rhs[i] = sum;
}
// zero rows of R
for (
int i = gauss_->rank;
i < h_;
++i
) {
if (rhs[i] != T(0)) {
return std::nullopt;
}
}
LinearSolution<T> sol;
/*
free variables = 0
*/
sol.particular =
Vector<T>(w_);
for (
int row = 0;
row < gauss_->rank;
++row
) {
int p =
gauss_->pivot_cols[row];
sol.particular[p]
= rhs[row];
}
sol.basis =
kernel_basis();
return sol;
}
template<class T>
std::optional<LinearSolution<T>>
Matrix<T>::solve_direct(
const Vector<T>& b
) const {
assert(b.size() == h_);
/*
A と b を直接同じ行基本変形にかける。
[ A | b ]
↓
[ R | c ]
変換行列 P は明示的には作らない。
*/
std::vector<T> rref = a_;
std::vector<T> rhs(h_);
for (int i = 0; i < h_; ++i) {
rhs[i] = b[i];
}
std::vector<int> pivot_cols;
pivot_cols.reserve(std::min(h_, w_));
int row = 0;
for (
int col = 0;
col < w_ && row < h_;
++col
) {
// --------------------------------------------
// pivot search
// --------------------------------------------
int pivot = -1;
for (int i = row; i < h_; ++i) {
if (
rref[(size_t)i * w_ + col]
!= T(0)
) {
pivot = i;
break;
}
}
if (pivot == -1) {
continue;
}
// --------------------------------------------
// row swap
// --------------------------------------------
if (pivot != row) {
T* pr =
rref.data()
+ (size_t)pivot * w_;
T* rr =
rref.data()
+ (size_t)row * w_;
std::swap_ranges(
pr,
pr + w_,
rr
);
std::swap(
rhs[pivot],
rhs[row]
);
}
T* rr =
rref.data()
+ (size_t)row * w_;
// --------------------------------------------
// normalize pivot row
// --------------------------------------------
T inv_p =
T(1) / rr[col];
for (int j = col; j < w_; ++j) {
rr[j] *= inv_p;
}
rhs[row] *= inv_p;
// --------------------------------------------
// eliminate all other rows
// --------------------------------------------
for (int i = 0; i < h_; ++i) {
if (i == row) {
continue;
}
T* ir =
rref.data()
+ (size_t)i * w_;
T factor = ir[col];
if (factor == T(0)) {
continue;
}
// この要素は必ず 0 になる
ir[col] = T(0);
for (
int j = col + 1;
j < w_;
++j
) {
ir[j]
-= factor * rr[j];
}
rhs[i]
-= factor * rhs[row];
}
pivot_cols.push_back(col);
++row;
}
const int rank = row;
// --------------------------------------------
// consistency check
// --------------------------------------------
/*
rank 以降の行は
0 0 ... 0 | rhs[i]
なので rhs[i] != 0 なら解なし。
*/
for (int i = rank; i < h_; ++i) {
if (rhs[i] != T(0)) {
return std::nullopt;
}
}
// --------------------------------------------
// construct general solution
// --------------------------------------------
LinearSolution<T> sol;
/*
特殊解:
自由変数をすべて 0 にする。
*/
sol.particular =
Vector<T>(w_);
for (int i = 0; i < rank; ++i) {
int p = pivot_cols[i];
sol.particular[p] = rhs[i];
}
// --------------------------------------------
// kernel basis
// --------------------------------------------
std::vector<char> is_pivot(
w_,
false
);
for (int p : pivot_cols) {
is_pivot[p] = true;
}
sol.basis.reserve(w_ - rank);
/*
各自由変数 f について
x_f = 1
他の自由変数 = 0
と置く。
RREF の pivot 行
x_p + Σ R[p,f] x_f = 0
より
x_p = -R[p,f]
*/
for (int f = 0; f < w_; ++f) {
if (is_pivot[f]) {
continue;
}
Vector<T> v(w_);
v[f] = T(1);
for (int i = 0; i < rank; ++i) {
int p = pivot_cols[i];
const T* ri =
rref.data()
+ (size_t)i * w_;
v[p] = -ri[f];
}
sol.basis.push_back(
std::move(v)
);
}
return sol;
}
int main(){
cin.tie(nullptr);
ios::sync_with_stdio(false);
cout<<fixed<<setprecision(10);
using Mat=Matrix<mint>;
using Vec=Vector<mint>;
ll n,c;
cin>>n>>c;
vector<mint> a(n);
cin>>a;
vector<mint> b(n+1);
b[0]=1;
for(int i=0;i<n;i++){
for(int j=n-1;j>=0;j--)b[j+1]-=b[j]*a[i];
}
Vec v(n);
v[0]=1;
for(int i=1;i<n;i++){
for(int j=1;j<=i;j++)v[i]-=v[i-j]*b[j];
}
Mat mat(n,n);
for(int i=0;i<n;i++)mat[n-1][i]=-b[n-i];
for(int i=1;i<n;i++)mat[n-i-1][n-i]=1;
mat=mat.pow(c);
v=mat*v;
for(int i=0;i<n;i++)v[0]-=a[i].pow(c);
cout<<v[0]<<endl;
}