// 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;
}