#pragma GCC optimize ("O3") #pragma GCC target ("avx") #include #include #include #include #include #include #include #include #include #include #include #include #include #define _rep(_1, _2, _3, _4, name, ...) name #define rep2(i, n) rep3(i, 0, n) #define rep3(i, a, b) rep4(i, a, b, 1) #define rep4(i, a, b, c) for (int i = int(a); i < int(b); i += int(c)) #define rep(...) _rep(__VA_ARGS__, rep4, rep3, rep2, _)(__VA_ARGS__) using namespace std; using i64 = long long; using u8 = unsigned char; using u32 = unsigned; using u64 = unsigned long long; using f80 = long double; // for dense class matrix { public: using R = int; using vec = vector; using poly = vec; struct xor128 { xor128(u32 seed=521288629) : z(seed) {} u32 operator() () { return next(); } u32 operator() (u32 r) { return next() % r; } u32 next() { u32 t = x ^ (x << 11); x = y; y = z; z = w; return (w ^= (w >> 19) ^ (t ^ (t >> 8))); } u32 x = 123456789, y = 362436069, w = 88675123, z = 521288629; }; static void reset_seed(int S) { xrand = xor128(S); } private: struct matrix64 { matrix64(const matrix& rhs) : m(rhs.m_), n(rhs.n_), a(rhs.a, rhs.a + m * n) {} i64* operator [] (int i) { return a.data() + n * i; } void swap_row(int i, int j) { if (i == j) return; rep(k, n) swap((*this)[i][k], (*this)[j][k]); } void print() { printf("matrix64(%d, %d):\n", m, n); puts("["); rep(i, m) { printf(" [%d", fixed((*this)[i][0])); rep(j, 1, n) printf(", %d", fixed((*this)[i][j])); puts("],"); } puts("]"); } int m, n; vector a; }; static R gcd(R a, R b) { return b ? gcd(b, a % b) : a; } static R mul_mod(R a, R b) { return i64(a) * b % mod; } static i64 add_mod64(i64 a, i64 b) { return (a += b - lmod) < 0 ? a + lmod : a; } static i64 sub_mod64(i64 a, i64 b) { return (a -= b) < 0 ? a + lmod : a; } static R fixed(R a) { return (a %= mod) < 0 ? a + mod : a; } static R mod_inv(R a, R m=0) { R b = (m == 0) ? mod : m; R s = 1, t = 0; while (b) { R q = a / b; swap(a %= b, b); swap(s -= t * q, t); } if (a != 1) { fprintf(stderr, "Error: Modular multiplicative inverse does not exist.\n"); assert(0); } return s < 0 ? s + mod : s; } public: static void set_mod(R m) { mod = m; lmod = i64(mod) << 32; modq = min(u64(1) << 16, (u64(-1) / mod - 1) / mod); } public: matrix(int n) : matrix(n, n) {} matrix(int m, int n, bool zfill=true) : m_(m), n_(n) { alloc(); if (zfill) clear(); } matrix(const matrix& m) : matrix(m.m_, m.n_, false) { copy(m.a, m.a + m_ * n_, a); } matrix(initializer_list< initializer_list > l) : matrix(l.size(), l.begin()->size()) { int i = 0; for (auto& v : l) { int j = 0; for (auto c : v) (*this)[i][j++] = fixed(c); ++i; } } matrix& operator = (matrix&& rhs) { if (this == &rhs) return *this; this->m_ = rhs.m_; this->n_ = rhs.n_; free(); this->a = rhs.a; rhs.a = nullptr; return *this; } ~matrix() { free(); } R* operator [] (int i) { return a + i * n_; } const R* operator [] (int i) const { return a + i * n_; } void clear() { fill(a, a + m_ * n_, 0); } void swap_row(int i, int j) { if (i == j) return; rep(k, n_) swap((*this)[i][k], (*this)[j][k]); } static matrix random_matrix(int m, int n) { auto ret = matrix(m, n, false); rep(i, m) rep(j, n) ret[i][j] = xrand(mod); return ret; } static vec random_vector(int n) { auto ret = vec(n); rep(i, n) ret[i] = xrand(mod); return ret; } matrix transpose() const { auto ret = matrix(n_, m_, false); rep(i, m_) rep(j, n_) ret[j][i] = (*this)[i][j]; return ret; } matrix operator * (const matrix& rhs) const { assert(n_ == rhs.m_); auto ret = matrix(m_, rhs.n_); auto rhsT = rhs.transpose(); rep(i, m_) mul_mat_vec(rhsT, (*this)[i], ret[i]); return ret; } matrix& operator *= (const matrix& rhs) { return *this = *this * rhs; } matrix krylov(int k, vec v) const { auto ret = matrix(k + 1, n_); rep(i, k) { copy(v.begin(), v.end(), ret[i]); mul_mat_vec(*this, v.data(), v.data()); } copy(v.begin(), v.end(), ret[k]); return ret; } template static void vec_sub(i64* a, T* b, int n, R c) { rep(i, n) a[i] = sub_mod64(a[i], i64(R(b[i])) * c); } static vec solve_linear_equations_in_Z_pZ(const matrix& in) { matrix64 mat = matrix64(in); assert(mat.m + 1 == mat.n); int rank = 0; rep(i, mat.m) { int piv = -1; for (int j = mat.m - 1; j >= rank; --j) if (mat[j][i] %= mod) piv = j; if (piv < 0) continue; if (rank != piv) mat.swap_row(rank, piv); rep(k, i + 1, mat.n) mat[rank][k] %= mod; auto inv = mod_inv(mat[rank][i]); rep(j, rank + 1, mat.m) if (mat[j][i]) { vec_sub(mat[j] + i, mat[rank] + i, mat.n - i, mul_mod(mat[j][i], inv)); } ++rank; } for (int i = rank - 1; i >= 0; --i) if (mat[i][i]) { R c = mat[i][mat.m] = mul_mod(mat[i][mat.m] % mod, mod_inv(mat[i][i])); rep(j, i) mat[j][mat.m] = sub_mod64(mat[j][mat.m], i64(c) * int(mat[j][i])); } vec ret(mat.n); ret[0] = fixed(1); rep(i, rank) ret[ret.size() - 1 - i] = mat[i][mat.m]; return ret; } static vec solve_linear_equations_in_Z_nZ(const matrix& in) { matrix64 mat = matrix64(in); assert(mat.m + 1 == mat.n); int rank = 0; rep(i, mat.m) { int piv = -1; R g = mod; rep(j, rank, mat.m) if ( (mat[j][i] %= mod) && g >= 0) { R k = gcd(R(mat[j][i]), mod); g = (k == 1) ? -1 : gcd(g, k); piv = j; } if (g == mod) continue; if (g > 0) { // To Do: not optimal R g2 = gcd(mat[piv][i], mod); for (int j = rank; j < mat.n && g2 != g; ++j) if (mat[j][i]) { R g3 = gcd(mat[j][i], g2); if (g3 == g2) continue; g2 = g3; while (gcd(mat[piv][i], mod) != g2) { R q = mat[j][i] / mat[piv][i]; rep(k, i + 1, mat.n) mat[piv][k] %= mod; if (q) vec_sub(mat[j] + i, mat[piv] + i, mat.n - i, q); mat.swap_row(piv, j); } } assert(g2 == g); } else { g = 1; } if (rank != piv) mat.swap_row(rank, piv); rep(k, i + 1, mat.n) mat[rank][k] %= mod; auto inv = mod_inv(mat[rank][i] / g, mod / g); rep(j, rank + 1, mat.m) if (mat[j][i]) { vec_sub(mat[j] + i, mat[rank] + i, mat.n - i, mul_mod(mat[j][i] / g, inv)); } ++rank; } for (int i = rank - 1; i >= 0; --i) if (mat[i][i]) { R g = gcd(gcd(R(mat[i][i]), R(mat[i][mat.m] %= mod)), mod); R inv = mod_inv(mat[i][i] / g, mod / g); R c = mat[i][mat.m] = mul_mod(mat[i][mat.m] / g, inv); rep(j, i) mat[j][mat.m] = sub_mod64(mat[j][mat.m], i64(c) * int(mat[j][i])); } vec ret(mat.n); ret[0] = fixed(1); rep(i, rank) ret[ret.size() - 1 - i] = mat[i][mat.m]; return ret; } static vec linear_recurrence(matrix kry, const vec& v, bool in_Z_pZ) { int m = kry.m_, n = kry.n_; assert(m == n + 1 && n == int(v.size())); rep(i, n) kry[n][i] = (kry[n][i] ? mod - kry[n][i] : 0); kry = kry.transpose(); auto ret = in_Z_pZ ? solve_linear_equations_in_Z_pZ(kry) : solve_linear_equations_in_Z_nZ(kry); return ret; } vec pow_vec(i64 e, const vec& v, bool in_Z_pZ=false) const { assert(m_ == n_); if (mod == 1) return vec(v.size(), 0); int k = min(e, i64(m_)); auto kry = krylov(k, v); if (e <= k) return vec(kry[e], kry[e + 1]); auto poly = linear_recurrence(kry, v, in_Z_pZ); auto rem = x_pow_mod(e, poly); reverse(rem.begin(), rem.end()); tmp64.assign(v.size(), 0); rep(i, rem.size()) if (rem[i]) rep(j, v.size()) { tmp64[j] = add_mod64(tmp64[j], i64(rem[i]) * kry[i][j]); } auto ret = vec(v.size()); rep(i, v.size()) ret[i] = tmp64[i] % mod; return ret; } static poly poly_mul(const poly& f, const poly& g) { int s1 = f.size(), s2 = g.size(), s = s1 + s2 - 1; tmp64.assign(s, 0); rep(i, s1) if (f[i]) rep(j, s2) { tmp64[i + j] = add_mod64(tmp64[i + j], i64(f[i]) * int(g[j])); } poly ret(s); rep(i, s) ret[i] = tmp64[i] % mod; return ret; } static poly poly_rem(const poly& f, const poly& g) { int s1 = f.size(), s2 = g.size(); if (s1 < s2) return f; assert(g[0] == 1); tmp64.resize(s1); copy(f.begin(), f.end(), tmp64.begin()); rep(i, s1 - s2 + 1) { int c = tmp64[i] % mod; if (c) rep(j, 1, s2) tmp64[i + j] = sub_mod64(tmp64[i + j], i64(c) * int(g[j])); tmp64[i] = c; } poly ret(s2 - 1); rep(i, s2 - 1) ret[i] = tmp64[s1 - s2 + 1 + i] % mod; return ret; } static poly x_pow_mod(i64 e, const poly& f) { if (e == 0) return poly(1, 1); poly ret = poly({1, 0}); i64 mask = (i64(1) << __lg(e)) >> 1; while (mask) { ret = poly_mul(ret, ret); if (e & mask) ret.push_back(0); ret = poly_rem(ret, f); mask >>= 1; } return ret; } static void mul_mat_vec(const matrix& m, const R* v, R* res) { buff.resize(m.m_); rep(i, m.m_) buff[i] = dot(m[i], v, m.n_); copy(buff.begin(), buff.begin() + m.m_, res); } void print() const { printf("matrix(%d, %d):\n", m_, n_); puts("["); rep(i, m_) { printf(" [%d", (*this)[i][0]); rep(j, 1, n_) printf(", %d", (*this)[i][j]); puts("],"); } puts("]"); } private: void alloc() { a = new R[m_ * n_]; } void free() { delete [] a; } static int dot(const int* a, const int* b, int n) { if (n == 0) return 0; int k = min(n, modq) & ~3; int r = (k ? n % k : n), q = (k ? n / k : 0); i64 s0 = 0, s1 = 0; if (k >= 12) { rep(i, q) { u64 s = 0; rep(j, i * k, (i + 1) * k) s += i64(a[j]) * b[j]; s0 += u32(s), s1 += u32(s >> 32); } if (r) { u64 s = 0; rep(j, q * k, n) s += i64(a[j]) * b[j]; s0 += u32(s), s1 += u32(s >> 32); } } else { rep(j, n) { i64 s = i64(a[j]) * b[j]; s0 += u32(s), s1 += u32(s >> 32); } } return (((s1 % mod) << 32) + s0) % mod; } static int mod, modq; static i64 lmod; static vector buff; static vector tmp64; static xor128 xrand; int m_, n_; R* a; }; vector matrix::buff; vector matrix::tmp64; int matrix::mod, matrix::modq; i64 matrix::lmod; matrix::xor128 matrix::xrand; void solve() { const int mod = 1e9 + 7; matrix::set_mod(mod); auto pow_mod = [&] (int b, i64 e) { int ret = 1; for (; e; e >>= 1, b = i64(b) * b % mod) if (e & 1) ret = i64(ret) * b % mod; return ret; }; int A[128]; int n; i64 c; while (~scanf("%d %lld", &n, &c)) { matrix mat(n, n, true); int ans = 0; rep(i, n) { scanf("%d", &A[i]); rep(j, i + 1) mat[i][j] = A[i]; ans = (ans + mod - pow_mod(A[i], c)) % mod; } auto res = mat.pow_vec(c - 1, matrix::vec(A, A + n)); rep(i, n) ans = (ans + res[i]) % mod; printf("%d\n", ans); } } int main() { auto beg = clock(); solve(); auto end = clock(); fprintf(stderr, "%.3f sec\n", double(end - beg) / CLOCKS_PER_SEC); }