#include "bits/stdc++.h" using namespace std; template struct ModInt { static const int Mod = MOD; unsigned x; ModInt() : x(0) { } ModInt(signed sig) { int sigt = sig % MOD; if (sigt < 0) sigt += MOD; x = sigt; } ModInt(signed long long sig) { int sigt = sig % MOD; if (sigt < 0) sigt += MOD; x = sigt; } int get() const { return (int)x; } ModInt &operator+=(ModInt that) { if ((x += that.x) >= MOD) x -= MOD; return *this; } ModInt &operator-=(ModInt that) { if ((x += MOD - that.x) >= MOD) x -= MOD; return *this; } ModInt &operator*=(ModInt that) { x = (unsigned long long)x * that.x % MOD; return *this; } ModInt &operator/=(ModInt that) { return *this *= that.inverse(); } ModInt operator+(ModInt that) const { return ModInt(*this) += that; } ModInt operator-(ModInt that) const { return ModInt(*this) -= that; } ModInt operator*(ModInt that) const { return ModInt(*this) *= that; } ModInt operator/(ModInt that) const { return ModInt(*this) /= that; } ModInt inverse() const { signed a = x, b = MOD, u = 1, v = 0; while (b) { signed t = a / b; a -= t * b; std::swap(a, b); u -= t * v; std::swap(u, v); } if (u < 0) u += Mod; ModInt res; res.x = (unsigned)u; return res; } bool operator==(ModInt that) const { return x == that.x; } bool operator!=(ModInt that) const { return x != that.x; } ModInt operator-() const { ModInt t; t.x = x == 0 ? 0 : Mod - x; return t; } }; typedef ModInt<1000000007> mint; struct GaussianEliminationCore { using Num = mint; using RowVec = vector; static void multiplySubtract(RowVec &x, const RowVec &y, int m, Num c) { if (c == Num()) return; for (int j = 0; j < m; ++ j) x[j] -= y[j] * c; } static Num dotProduct(const RowVec &x, const RowVec &y, int m) { Num sum = Num(); for (int j = 0; j < m; ++ j) sum += x[j] * y[j]; return sum; } vector basis; vector invDiagonal; vector order; void init(int m) { basis.assign(m, RowVec()); invDiagonal.assign(m, Num()); order.clear(); } int size() const { return (int)basis.size(); } void eliminate(RowVec &row, vector &coeffs) const { for (int i : order) { Num c = row[i] * invDiagonal[i]; coeffs[i] = c; multiplySubtract(row, basis[i], size(), c); } } void add(const RowVec &row, int k) { assert(row[k] != Num() && invDiagonal[k] == Num()); basis[k] = row; invDiagonal[k] = row[k].inverse(); order.push_back(k); } }; // idea: use falling factorial basis for polynomials struct PolynomiallyRecursiveSequence : private GaussianEliminationCore { enum PolynomalBasisKind { MonominalBasis, FallingFactorialBasis, } polynomialBasisKind = MonominalBasis; void getPolynomialBasis(vector &xs, int D, int n) { for (int j = 0; j < D; ++ j) { Num x = 1; if (polynomialBasisKind == MonominalBasis) { for (int k = 0; k < j; ++ k) x *= n; } else { for (int k = 0; k < j; ++ k) x *= n - k; } xs[j] = x; } } vector solve(const vector &seq, const vector> &indices) { int K = 0, D = 0; for (auto ix : indices) { if (K <= ix.first) K = ix.first + 1; if (D <= ix.second) D = ix.second + 1; } if (K == 0 || D == 0) return {}; int N = (int)seq.size(), m = (int)indices.size(); init(m); vector newRow(m), coeffs(m); for (int n = K - 1; n < N; ++ n) { vector ys(K), xs(D); for (int i = 0; i < K; ++ i) ys[i] = seq[n - i]; getPolynomialBasis(xs, D, n); newRow.clear(); for (auto ix : indices) newRow.push_back(ys[ix.first] * xs[ix.second]); eliminate(newRow, coeffs); for (int k = 0; k < m; ++ k) if (newRow[k] != mint()) { add(newRow, k); break; } if (order.size() == m) return {}; } int k = 0; for (; invDiagonal[k] != Num(); ++ k); vector res(m); res[k] = 1; reverse(order.begin(), order.end()); for (int i : order) { Num dp = dotProduct(res, basis[i], m); res[i] -= dp * invDiagonal[i]; } return res; } vector> findMinimalSolution(const vector &seq) { vector> best; int bestNum = numeric_limits::max(); for (int KD = 1; ; ++ KD) { for (int K = 1; K <= KD; ++ K) if (KD % K == 0) { int D = KD / K; vector enabled(KD, true); int currentNum = KD, leastNum = 0; auto check = [&]() -> bool { vector> indices; for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j) if (enabled[i * D + j]) indices.emplace_back(i, j); auto sol = solve(seq, indices); if (!sol.empty()) { if (bestNum > currentNum) { best.assign(K, vector(D)); int t = 0; for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j) if (enabled[i * D + j]) best[i][j] = sol[t ++]; bestNum = currentNum; } return true; } else { return false; } }; function rec = [&](int i, int j) { if (bestNum <= leastNum) return; if (i == K) { check(); return; } if (j == D) return rec(i + 1, 0); for (int e = 0; e < 2; ++ e) { enabled[i * D + j] = e != 0; if (e == 0) -- currentNum; else ++ leastNum; if (e == 1 || check()) { rec(i, j + 1); } if (e == 0) ++ currentNum; else -- leastNum; } }; rec(0, 0); } if (!best.empty()) break; } cerr << "bestNum = " << bestNum << endl; return best; } }; templateT gcd(T x, T y) { if (y == 0)return x; else return gcd(y, x%y); } template int mintToSigned(ModInt a) { int x = a.get(); if (x <= MOD / 2) return x; else return x - MOD; } string mintToSignedRatio(mint a, int den = 36) { int x = mintToSigned(a * den); int g = gcd(abs(x), abs(den)); if (den / g < 0) g *= -1; x /= g, den /= g; stringstream ss; if (den == 1) ss << x; else ss << x << "/" << den; return ss.str(); } int main() { if(0) { vector seq = { 1, 1, 2, 2, 8, 28, 152, 952, 7208, 62296, 605864, 6522952, 76951496, 986411272, 647501133, 653303042, 170637030, 248109503, 700583494, 619914523, 682935856, 443753916, 423068688, 507501942, 315541972, 110825117, 848156395, 798418282, 920964362, 23823302, 114894774, 279365223, 992413784, 833179437, 785518302, 524368220, 42214454, 140345871, 188150268, 808714798, 718376249, 732000901, 955005007, 139255097, 484615744, 615066955, 726914809, 856989248, 460819998, 321277105, 536397091, 555447300, 597473569, 217709372, 24981477, 143561526, 171000806, 137649694, 749333590, 700935246, 916763337, 762367836, 296796066, 236278263, 398507715, 148909632, 568524543, 926513708, 163591024, 339393165, 549241395, 548924577, 915489821, 706913104, 380913764, 993919668, 895691202, 628078606, 542382606, 735060428, 385303214, 453133962, 470556393, 439972973, 4764973, 459438929, 49172129, 93448766, 14767450, 302365655, 44994640, 637650527, 462797839, 174866371, 963824426, 761996745, 999013044, 209330964, 997280223, 561428453, 300321098, 733666220, 501835787, 614206349, 790510106, 3598758, 51698818, 133690163, 663340013, 193264830, 514261029, 953300090, 783013838, 841987389, 44621114, 487627404 }; int N = (int)seq.size(); cerr << "N = " << N << endl; PolynomiallyRecursiveSequence prs; prs.polynomialBasisKind = PolynomiallyRecursiveSequence::MonominalBasis; auto solution = prs.findMinimalSolution(seq); if (solution.empty()) { cerr << "No solution found" << endl; return 1; } int K = (int)solution.size(), D = (int)solution[0].size(); for (int n = K - 1; n < N; ++ n) { vector ys(K), xs(D); for (int i = 0; i < K; ++ i) ys[i] = seq[n - i]; prs.getPolynomialBasis(xs, D, n); mint sum; for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j) sum += solution[i][j] * ys[i] * xs[j]; if (sum != mint()) cerr << "err" << endl; } for (int i = 0; i < K; ++ i) { auto p = solution[i]; if (i == 0) { for (auto &x : p) x = -x; } int t = 0; for(auto x : p) t += x != mint(); int hi = D - 1; for (; hi >= 0 && p[hi] == mint(); -- hi); if (hi >= 0 && mintToSignedRatio(p[hi])[0] == '-') { cout << (i > 1 ? " - " : "-"); for (auto &x : p) x = -x; } else { if(i > 1) cout << " + "; } if (hi == 0 && p[hi] == 1) { } else { if (t > 1) cout << "("; auto &o = cout; bool first = true; for (int j = D - 1; j >= 0; -- j) { string c = mintToSignedRatio(p[j]); if (c != "0") { if (first && c[0] == '-') o << "-"; else if (first) o << ""; else if (c[0] == '-') o << " - "; else o << " + "; if (j != 0 && (c == "1" || c == "-1")) o << ""; else if (c[0] == '-') o << c.substr(1); else o << c; if (j == 0) o << ""; else if (j == 1) o << "n"; else { if (prs.polynomialBasisKind == PolynomiallyRecursiveSequence::MonominalBasis) { o << "n^" << j; } else { o << "n"; for (int k = 1; k < j; ++ k) o << "(n-" << k << ")"; } } first = false; } } if (first) o << "0"; if (t > 1) cout << ")"; cout << " "; } cout << "a_"; if (i == 0) cout << "n"; else cout << "{n-" << i << "}"; if (i == 0) { cout << " = "; if (K == 1) cout << "0"; } } cout << endl; cout << "const int K = " << K << ";\n"; cout << "const array init = { "; for (int j = 0; j < K - 1; ++ j) cout << (j == 0 ? "" : ", ") << seq[j].get(); cout << " };\n"; cout << "const array, K> coeffs = { {\n"; for (int i = 0; i < K; ++ i) { cout << "\t{ "; for (int j = 0; j < D; ++ j) cout << (j == 0 ? "" : ", ") << mintToSigned(solution[i][j] * 2); cout << " },\n"; } } int N; while (cin >> N) { const int K = 14; const array init = { 1, 1, 2, 2, 8, 28, 152, 952, 7208, 62296, 605864, 6522952, 76951496 }; const array, K> coeffs = { { { 4, 0, 0 }, { -18, -4, 0 }, { 6, 16, 0 }, { 18, -11, 2 }, { 116, -7, -9 }, { -138, -6, 8 }, { -64, -78, 14 }, { -90, 136, -18 }, { -632, 256, -20 }, { 892, -298, 22 }, { 1674, -354, 18 }, { -1480, 307, -16 }, { -240, 55, -3 }, { 240, -44, 2 }, } }; const mint invDen = coeffs[0][0].inverse(); vector seq(init.begin(), init.end()); seq.resize(N + 1); for (int n = K - 1; n <= N; ++ n) { mint n1 = n, n2 = n1 * n1; mint sum; for (int i = 1; i < K; ++ i) sum += (coeffs[i][0] + coeffs[i][1] * n1 + coeffs[i][2] * n2) * seq[n - i]; seq[n] = -sum * invDen; } mint ans = seq[N]; printf("%d\n", ans.get()); } }