#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); } }; 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, uint64_t maxRecurse = numeric_limits::max()) { vector> best; int bestNum = numeric_limits::max(); uint64_t numRecurse = 0; for (int KD = 1; ; ++ KD) { cerr << "KD = " << KD << "..." << endl; 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 { ++ numRecurse; if (numRecurse % 10000 == 0) cerr << "checking " << numRecurse << " times... (bestNum = " << bestNum << ")" << endl; 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 (sol[0] != mint()) { auto inv = sol[0].inverse(); for (auto &x : sol) x *= inv; } 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; cerr << "bestNum updated: " << bestNum << endl; } return true; } else { return false; } }; function rec = [&](int i, int j) { if (bestNum <= leastNum) return; if (i == K) { check(); ++ numRecurse; 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 (numRecurse > maxRecurse) return; } if (e == 0) ++ currentNum; else -- leastNum; } }; if (check()) rec(0, 0); if (numRecurse > maxRecurse) { cerr << "(suboptimal result)" << endl; break; } } if (!best.empty()) break; if (numRecurse > maxRecurse) { break; } } if (!best.empty()) { 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 x = mintToSigned(a), d = 1; for (int dd = 1; dd <= 60; ++ dd) { int xx = mintToSigned(a * dd); if (abs(x) > abs(xx)) x = xx, d = dd; } int g = gcd(abs(x), abs(d)); if (d / g < 0) g *= -1; x /= g, d /= g; stringstream ss; if (d == 1) ss << x; else ss << x << "/" << d; return ss.str(); } #define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i)) vector fact, factinv; void nCr_computeFactinv(int N) { N = min(N, mint::Mod - 1); fact.resize(N + 1); factinv.resize(N + 1); fact[0] = 1; rer(i, 1, N) fact[i] = fact[i - 1] * i; factinv[N] = fact[N].inverse(); for (int i = N; i >= 1; i --) factinv[i - 1] = factinv[i] * i; } mint nCr(int n, int r) { if (n >= mint::Mod) return nCr(n % mint::Mod, r % mint::Mod) * nCr(n / mint::Mod, r / mint::Mod); return r > n ? 0 : fact[n] * factinv[n - r] * factinv[r]; } mint catalan_number(int n) { return n == 0 ? 1 : nCr(2 * n, n) - nCr(2 * n, n - 1); } int digitsMod(const char s[], int n, int m) { int x = 0; long long pow10 = 1; for (int i = n - 1; i >= 0; i --) { if ((x += pow10 * (s[i] - '0') % m) >= m) x -= m; pow10 = pow10 * 10 % m; } return x; } template ModInt operator^(ModInt a, unsigned long long k) { ModInt r = 1; while (k) { if (k & 1) r *= a; a *= a; k >>= 1; } return r; } int main() { if (0) { vector seq; if (1) { seq = { 1,1,3,4,19,26,144,197,1171,1597,9878,13432,85216,115597,746371,1010504,6609043,8933858,59008563,79662593,530279894,715116833,790262320,454109424,458522675,507912745,683985782,349356902,641159448,914253017,666123528,256877025,283012377,259341341,985344830,99096838,997431143,359782436,337868019,275073093,823041379,299193173,618332669,623673376,788988015,348012107,960452357,132459813,484572463,876182618,173753254,441164007,288391002,290937936,415692760,211813919,329314918,283643646,565011970,845017018,390800948,209282283,329903021,700520426,397612109,649057223,946472136,142317513,639359324,720717675,983275264,507698034,651899763,234922533,86387948,840663423,815571175,144003918,212527976,763824514,542610091,452847843,921495990,979234007,906307312,656853418,297047315,193463125,69677636,636423896,799010307,909263327,147379495,774893131,945207456,461924132,658630324,333324822,851825875,36247927,769534543,614653765,155284112,370877182,608767657,471252465,903152244,327024448,257349836,431774297,714479094,114462919,350021247,139934784,801799813,539551848,41576353,802846678,659340919,672408792,575487541,368702253,575350782,619200170,562982523,429350007,150624942,945196520,6319297,894502025,374999342,42453831,221527997,3616961,580046947,341012885,276733823,579878644,670787005,162408190,615289155,917206778,768258447,411299351,832373932,18107510,835563777,233896612,738886121,393809140,644927895,625120199,336985781,516279761,519173913,217509375,754521244,424721008,850139747,438485991,920679348,219842245,442535106,128849044,269625477,753265916,440813281,928041043,42439991,348131531,381592048,654276173,47700473,43991743,548001157,730540323,887685216,932597280,801983201,142570249,471547353,382906298,745576798,286187856,510059281,818604713,507091460,180496747,177275620,665329553,200027308,550501620,759618434,249162203,361024171,165602025,282541684,651186107,424876182,377693857,567673772,859949936,529948870,981890830,720454551,978293949,769899874,526887448,623249023,930481714,794783876,957457862,564435488,398080331,34219860,458598094,396377234,642505402,706573841,542803667,314649460,302333235,577578191,692607196,984717758,51551222,571746567,181633622,552463582,538957298,753601027,442293302,778981625,146456018,791162799,991018889,504501545,952054904,231445357,40511975,512943594,225912874,965164521,550144404,476038888,439900254,975615018,195831342,689809368,815767322,935245854,688904403,993516754,439936142,668910186,937167584,657694071,25151922,724216878,897387781,13684183,263175987,781446329,542562472,747813284,240316688,653466082,475627298,753619876,473430472,816466673,551717822,299623712,258591955,913624315,779342388,313872566,912716451,968634939,645165489,119540208,327353173,70136686,329097227,867149612,458779459,300118018,5895168,924317,136386920,430435362,268428795,484792722,744629423,976335117,799648412,269161423,809591377,443320513,337653395,591911355,627672258,275133634,701549043,952539877,22665917 }; } else if (1) { nCr_computeFactinv(1000000); for (int n = 0; n < 1000; ++ n) { mint sum; for (int i = 0; i <= n; ++ i) sum += nCr(n * 3, i);//nCr(n * 7, i * 3); seq.push_back(sum); } } else { int n; string s; while (cin >> n >> s) { if (n != seq.size()) abort(); seq.push_back(digitsMod(s.c_str(), (int)s.size(), mint::Mod)); } } int N = (int)seq.size(); cerr << "N = " << N << endl; PolynomiallyRecursiveSequence prs; prs.polynomialBasisKind = PolynomiallyRecursiveSequence::MonominalBasis; auto solution = prs.findMinimalSolution(seq, 100000); 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 << ")"; if (t > 0) cout << " "; } if (hi >= 0) { cout << "a_"; if (i == 0) cout << "n"; else cout << "{n-" << i << "}"; } if (i == 0) { cout << " = "; if (K == 1) cout << "0"; } } cout << endl; mint multiplier = 1;//12; for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j) solution[i][j] *= multiplier; for (int j = D - 1; j >= 0; -- j) { int c = mintToSigned(solution[0][j]); if (c != 0) { if (c < 0) { for (int i = 0; i < K; ++ i) for (int j = 0; j < D; ++ j) solution[i][j] *= -1; } break; } } /* 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]); cout << " },\n"; } cout << "} };\n"; 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 << "vector seq(init.begin(), init.end());\n"; cout << "seq.resize(N + 1);\n"; cout << "for (int n = K - 1; n <= N; ++ n) {\n"; cout << " mint "; for (int j = 1; j < D; ++ j) { if (j != 1) cout << ", "; cout << "n" << j << " = "; if (j == 1) cout << "n"; else cout << "n" << (j - 1) << " * n1"; } cout << ";\n"; auto outputSum = [&](int i, bool negate) { int t = 0; for (int j = 0; j < D; ++ j) t += solution[i][j] != mint(); if (t == 0) { cout << "0"; return; } if (t > 1) cout << "("; int k = 0; for (int j = D - 1; j >= 0; -- j) { int c = mintToSigned(solution[i][j] * (negate ? -1 : 1)); if (c != 0) { if (k ++ > 0) { if (c < 0) { cout << " - "; c = -c; } else { cout << " + "; } } if (c < 0) { cout << "-"; c = -c; } if (j == 0) { cout << c; } else { cout << "n" << j; if (c != 1) cout << " * " << c; } } } if (t > 1) cout << ")"; }; cout << " mint sum;\n"; for (int i = 1; i < K; ++ i) { cout << " sum "; bool negative = false; for (int j = D - 1; j >= 0; -- j) { int c = mintToSigned(-solution[i][j]); if (c != 0) { negative = c < 0; break; } } cout << (negative ? "-=" : "+="); cout << " seq[n - " << i << "] * "; outputSum(i, !negative); cout << ";\n"; } cout << " seq[n] = sum / "; outputSum(0, false); cout << ";\n"; cout << "}\n"; } else { int N; while (~scanf("%d", &N)) { const int K = 8; const array init = { 1, 1, 3, 4, 19, 26, 144 }; vector seq(init.begin(), init.end()); seq.resize(N + 1); for (int n = K - 1; n <= N; ++ n) { mint n1 = n, n2 = n1 * n1, n3 = n2 * n1, n4 = n3 * n1; mint sum; sum += seq[n - 1] * (n4 * 52269863 + n3 * 305931736 + n2 * 258365537 + n1 * 276782898 - 253220885); sum -= seq[n - 2] * (n4 * 242461620 + n3 * 27603152 + n2 * 139205790 + n1 * 267856601 + 277902371); sum -= seq[n - 3] * (n4 * 317850736 + n3 * 33389931 + n2 * 29110388 + n1 * 483840402 + 262098585); sum -= seq[n - 4] * (n4 * 446541231 + n3 * 262847746 + n2 * 474399727 - n1 * 329634661 + 301985484); sum += seq[n - 5] * (n4 * 179688158 - n3 * 122495876 + n2 * 117668564 + n1 * 427842771 + 117888425); sum -= seq[n - 6] * (n4 * 319613082 - n3 * 141115369 + n2 * 359231985 + n1 * 179539766 - 342110626); sum -= seq[n - 7] * (n4 * 76661316 - n3 * 336327801 + n2 * 475890910 + n1 * 251448872 - 347026197); seq[n] = sum / (n4 * 420387432 + n3 * 79612576 - n2 * 73247243 + n1); } mint ans = seq[N]; printf("%d\n", ans.get()); } } }