fork download
copy 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3.  
  4. using ll = long long;
  5. using Matrix = vector<vector<ll>>;
  6.  
  7. // ===================================================================
  8. // This file contains a collection of Matrix Exponentiation algorithms.
  9. // Each function is ready to be used as a "black box".
  10. // Read the comments above each one to understand:
  11. //   - What it solves
  12. //   - What input it expects
  13. //   - What it returns
  14. //   - Time complexity
  15. //   - Important constraints / assumptions
  16. // ===================================================================
  17.  
  18. // ===================================================================
  19. // SECTION 1: BASIC MATRIX OPERATIONS
  20. // ===================================================================
  21.  
  22. // -------------------------------------------------------------------
  23. // matMul(A, B, mod)
  24. // -------------------------------------------------------------------
  25. // PURPOSE:
  26. //   Multiplies two matrices A and B.
  27. //
  28. // INPUT:
  29. //   A : a matrix of size (n x p), B : a matrix of size (p x m)
  30. //   mod : modulo value (e.g., 1e9+7)
  31. //
  32. // OUTPUT:
  33. //   Returns a new matrix C = A * B of size (n x m), with all entries
  34. //   reduced modulo 'mod'.
  35. //
  36. // TIME COMPLEXITY:
  37. //   O(n * p * m)  (naive triple loop)
  38. //
  39. // CONSTRAINTS / PRECONDITIONS:
  40. //   - The number of columns of A must equal the number of rows of B.
  41. //   - mod should be positive.
  42. //   - Entries should be non‑negative or already reduced; multiplication
  43. //     may overflow 64‑bit if mod is large, so 'mod' is typically < 1e9.
  44. //   - For larger matrices, consider using __int128 if needed.
  45. Matrix matMul(const Matrix& A, const Matrix& B, ll mod) {
  46.     int n = (int)A.size();
  47.     int p = (int)A[0].size();
  48.     int m = (int)B[0].size();
  49.     Matrix C(n, vector<ll>(m, 0));
  50.     for (int i = 0; i < n; ++i) {
  51.         for (int k = 0; k < p; ++k) {
  52.             if (A[i][k] == 0) continue;        // skip zeros for speed
  53.             ll aik = A[i][k];
  54.             for (int j = 0; j < m; ++j) {
  55.                 C[i][j] = (C[i][j] + aik * B[k][j]) % mod;
  56.             }
  57.         }
  58.     }
  59.     return C;
  60. }
  61.  
  62. // -------------------------------------------------------------------
  63. // matAdd(A, B, mod)
  64. // -------------------------------------------------------------------
  65. // PURPOSE:
  66. //   Adds two matrices of the same size.
  67. //
  68. // INPUT:
  69. //   A, B : matrices of size (n x m)
  70. //   mod  : modulo value
  71. //
  72. // OUTPUT:
  73. //   Returns a new matrix C = (A + B) modulo 'mod'.
  74. //
  75. // TIME COMPLEXITY:
  76. //   O(n * m)
  77. Matrix matAdd(const Matrix& A, const Matrix& B, ll mod) {
  78.     int n = (int)A.size();
  79.     int m = (int)A[0].size();
  80.     Matrix C(n, vector<ll>(m, 0));
  81.     for (int i = 0; i < n; ++i)
  82.         for (int j = 0; j < m; ++j)
  83.             C[i][j] = (A[i][j] + B[i][j]) % mod;
  84.     return C;
  85. }
  86.  
  87. // -------------------------------------------------------------------
  88. // matPow(base, exp, mod)
  89. // -------------------------------------------------------------------
  90. // PURPOSE:
  91. //   Raises a square matrix 'base' to the power 'exp' using binary
  92. //   exponentiation. This is the core function for matrix exponentiation.
  93. //
  94. // INPUT:
  95. //   base : a square matrix of size (n x n)
  96. //   exp  : exponent (non‑negative integer, can be up to 1e18)
  97. //   mod  : modulo value
  98. //
  99. // OUTPUT:
  100. //   Returns the matrix base^exp, with all entries reduced modulo 'mod'.
  101. //
  102. // TIME COMPLEXITY:
  103. //   O(n^3 * log(exp))  (each multiplication is O(n^3), repeated log exp times)
  104. //
  105. // CONSTRAINTS / PRECONDITIONS:
  106. //   - base must be square.
  107. //   - exp >= 0.
  108. //   - mod > 0.
  109. Matrix matPow(Matrix base, ll exp, ll mod) {
  110.     int n = (int)base.size();
  111.     // Initialize result as identity matrix
  112.     Matrix res(n, vector<ll>(n, 0));
  113.     for (int i = 0; i < n; ++i) res[i][i] = 1 % mod;
  114.  
  115.     while (exp > 0) {
  116.         if (exp & 1) res = matMul(res, base, mod);
  117.         base = matMul(base, base, mod);
  118.         exp >>= 1;
  119.     }
  120.     return res;
  121. }
  122.  
  123. // -------------------------------------------------------------------
  124. // matVecMul(M, v, mod)
  125. // -------------------------------------------------------------------
  126. // PURPOSE:
  127. //   Multiplies a matrix M (size n x m) by a column vector v (size m x 1).
  128. //
  129. // INPUT:
  130. //   M : matrix of size (n x m)
  131. //   v : vector of length m
  132. //   mod : modulo value
  133. //
  134. // OUTPUT:
  135. //   Returns a vector of length n = M * v (mod 'mod').
  136. //
  137. // TIME COMPLEXITY:
  138. //   O(n * m)
  139. vector<ll> matVecMul(const Matrix& M, const vector<ll>& v, ll mod) {
  140.     int n = (int)M.size();
  141.     int m = (int)M[0].size();
  142.     vector<ll> res(n, 0);
  143.     for (int i = 0; i < n; ++i) {
  144.         for (int j = 0; j < m; ++j) {
  145.             res[i] = (res[i] + M[i][j] * v[j]) % mod;
  146.         }
  147.     }
  148.     return res;
  149. }
  150.  
  151. // ===================================================================
  152. // SECTION 2: FIBONACCI AND LINEAR RECURRENCES (using Matrix)
  153. // ===================================================================
  154.  
  155. // -------------------------------------------------------------------
  156. // fib(n, mod)
  157. // -------------------------------------------------------------------
  158. // PURPOSE:
  159. //   Computes the n‑th Fibonacci number (F(0)=0, F(1)=1) modulo 'mod'.
  160. //
  161. // INPUT:
  162. //   n   : index (n >= 0, can be up to 1e18)
  163. //   mod : modulo value
  164. //
  165. // OUTPUT:
  166. //   Returns F(n) modulo 'mod'.
  167. //
  168. // TIME COMPLEXITY:
  169. //   O(log n)  (2x2 matrix exponentiation)
  170. //
  171. // CONSTRAINTS / PRECONDITIONS:
  172. //   - mod > 0.
  173. //   - Works for any non‑negative integer n.
  174. ll fib(ll n, ll mod) {
  175.     if (n == 0) return 0;
  176.     if (n == 1) return 1 % mod;
  177.     // Fibonacci transition matrix: [[1,1],[1,0]]
  178.     Matrix base = {{1 % mod, 1 % mod}, {1 % mod, 0}};
  179.     Matrix res = matPow(base, n - 1, mod);
  180.     // F(n) = res[0][0] * F(1) + res[0][1] * F(0) = res[0][0]
  181.     return res[0][0];
  182. }
  183.  
  184. // -------------------------------------------------------------------
  185. // linearRecurrence(init, coeff, n, mod)
  186. // -------------------------------------------------------------------
  187. // PURPOSE:
  188. //   Computes the n‑th term of a linear recurrence.
  189. //
  190. //   The recurrence is defined as:
  191. //     f[n] = coeff[0] * f[n-1] + coeff[1] * f[n-2] + ... + coeff[k-1] * f[n-k]
  192. //   for n >= k, with initial terms f[0], f[1], ..., f[k-1] given in 'init'.
  193. //
  194. // INPUT:
  195. //   init  : vector of length k, containing f[0] ... f[k-1]
  196. //   coeff : vector of length k, coefficients in the order shown above
  197. //   n     : index of the term to compute (0‑based, n >= 0)
  198. //   mod   : modulo value
  199. //
  200. // OUTPUT:
  201. //   Returns f[n] modulo 'mod'.
  202. //
  203. // TIME COMPLEXITY:
  204. //   O(k^3 * log n) using matrix exponentiation.
  205. //   For k up to ~50 it is acceptable; for larger k use Kitamasa (see below).
  206. //
  207. // CONSTRAINTS / PRECONDITIONS:
  208. //   - k >= 1.
  209. //   - n >= 0.
  210. //   - If n < k, the function directly returns init[n].
  211. //   - The recurrence must hold for n >= k.
  212. //   - All values should be reduced modulo 'mod'.
  213. ll linearRecurrence(const vector<ll>& init, const vector<ll>& coeff, ll n, ll mod) {
  214.     int k = (int)init.size();
  215.     if (n < k) return init[n] % mod;
  216.  
  217.     // Build companion matrix of size k x k.
  218.     // State vector: [f[t], f[t-1], ..., f[t-k+1]]^T
  219.     Matrix T(k, vector<ll>(k, 0));
  220.     for (int i = 0; i < k; ++i) T[0][i] = coeff[i] % mod;       // first row
  221.     for (int i = 1; i < k; ++i) T[i][i-1] = 1;                  // sub-diagonal
  222.  
  223.     // We need T^(n - k + 1) because we start from state at t = k-1.
  224.     Matrix Tpow = matPow(T, n - k + 1, mod);
  225.  
  226.     // Initial state at t = k-1: [f[k-1], f[k-2], ..., f[0]]^T
  227.     vector<ll> state(k);
  228.     for (int i = 0; i < k; ++i) state[i] = init[k-1-i] % mod;
  229.  
  230.     vector<ll> res = matVecMul(Tpow, state, mod);
  231.     return res[0] % mod;
  232. }
  233.  
  234. // ===================================================================
  235. // SECTION 3: ADVANCED LINEAR RECURRENCE (Kitamasa / Polynomial)
  236. // ===================================================================
  237. // These functions compute the n‑th term of a linear recurrence in
  238. // O(k^2 log n) instead of O(k^3 log n), which is useful when k is large.
  239. //
  240. // TERMINOLOGY:
  241. //   - Characteristic polynomial: derived from the recurrence.
  242. //     For recurrence f[n] = c0*f[n-1] + c1*f[n-2] + ... + c[k-1]*f[n-k],
  243. //     the characteristic polynomial is:
  244. //       P(x) = x^k - c0*x^(k-1) - c1*x^(k-2) - ... - c[k-1].
  245. //   - We compute x^n mod P(x) using binary exponentiation of polynomials.
  246. //   - Then f[n] = sum_{i=0}^{k-1} r[i] * f[i], where r is the remainder
  247. //     polynomial.
  248. // -------------------------------------------------------------------
  249.  
  250. // Helper: Multiply two polynomials modulo the characteristic polynomial.
  251. // Both polynomials have degree < k. The product is reduced using
  252. // the relation: x^k = c0*x^(k-1) + c1*x^(k-2) + ... + c[k-1].
  253. vector<ll> polyMulMod(const vector<ll>& a, const vector<ll>& b,
  254.                       const vector<ll>& coeff, ll mod) {
  255.     int k = (int)coeff.size();
  256.     vector<ll> res(2 * k - 1, 0);
  257.     // Multiply
  258.     for (int i = 0; i < k; ++i) {
  259.         if (a[i] == 0) continue;
  260.         for (int j = 0; j < k; ++j) {
  261.             res[i+j] = (res[i+j] + a[i] * b[j]) % mod;
  262.         }
  263.     }
  264.     // Reduce terms with degree >= k using the recurrence.
  265.     // We iterate from high degree down to k.
  266.     for (int deg = 2*k - 2; deg >= k; --deg) {
  267.         if (res[deg] == 0) continue;
  268.         ll coef = res[deg];
  269.         // For each i from 0 to k-1, x^deg = coeff[i] * x^(deg-1-i)
  270.         // Because x^k = coeff[0]*x^(k-1) + coeff[1]*x^(k-2) + ... + coeff[k-1]
  271.         // So x^deg = x^(deg-k) * x^k = sum_{i=0}^{k-1} coeff[i] * x^(deg-1-i)
  272.         for (int i = 0; i < k; ++i) {
  273.             res[deg - 1 - i] = (res[deg - 1 - i] + coef * coeff[i]) % mod;
  274.         }
  275.         // The term res[deg] is now eliminated.
  276.     }
  277.     // Return only the first k coefficients.
  278.     vector<ll> reduced(k);
  279.     for (int i = 0; i < k; ++i) reduced[i] = res[i] % mod;
  280.     return reduced;
  281. }
  282.  
  283. // -------------------------------------------------------------------
  284. // linearRecurrenceKitamasa(init, coeff, n, mod)
  285. // -------------------------------------------------------------------
  286. // PURPOSE:
  287. //   Computes the n‑th term of a linear recurrence using Kitamasa's
  288. //   algorithm (polynomial exponentiation) – faster than matrix method
  289. //   for large k.
  290. //
  291. // INPUT:
  292. //   init  : vector of length k, f[0] ... f[k-1]
  293. //   coeff : vector of length k, recurrence coefficients (as defined above)
  294. //   n     : index to compute (0‑based)
  295. //   mod   : modulo value
  296. //
  297. // OUTPUT:
  298. //   Returns f[n] modulo 'mod'.
  299. //
  300. // TIME COMPLEXITY:
  301. //   O(k^2 * log n)  (due to polynomial multiplication each step)
  302. //
  303. // CONSTRAINTS / PRECONDITIONS:
  304. //   - k >= 1.
  305. //   - n >= 0.
  306. //   - Works for large k (e.g., k up to 500) within time limits.
  307. //   - All arithmetic is modulo 'mod'.
  308. ll linearRecurrenceKitamasa(const vector<ll>& init, const vector<ll>& coeff,
  309.                             ll n, ll mod) {
  310.     int k = (int)init.size();
  311.     if (n < k) return init[n] % mod;
  312.  
  313.     // We want to compute x^n modulo characteristic polynomial.
  314.     // Start with polynomial representing x^1.
  315.     vector<ll> pol(k, 0);
  316.     if (k == 1) {
  317.         // For k=1, characteristic polynomial is x - c0, so x ≡ c0 (mod P).
  318.         pol[0] = coeff[0] % mod;
  319.     } else {
  320.         pol[1] = 1; // x
  321.     }
  322.     // result polynomial starts as 1 (x^0)
  323.     vector<ll> res(k, 0);
  324.     res[0] = 1;
  325.  
  326.     ll exp = n;
  327.     while (exp > 0) {
  328.         if (exp & 1) {
  329.             res = polyMulMod(res, pol, coeff, mod);
  330.         }
  331.         pol = polyMulMod(pol, pol, coeff, mod);
  332.         exp >>= 1;
  333.     }
  334.  
  335.     // Now res represents x^n mod P(x), i.e., res[i] = coefficient of x^i.
  336.     // f[n] = sum_{i=0}^{k-1} res[i] * f[i]
  337.     ll ans = 0;
  338.     for (int i = 0; i < k; ++i) {
  339.         ans = (ans + res[i] * (init[i] % mod)) % mod;
  340.     }
  341.     return ans;
  342. }
  343.  
  344. // ===================================================================
  345. // SECTION 4: GRAPH APPLICATIONS
  346. // ===================================================================
  347.  
  348. // -------------------------------------------------------------------
  349. // countWalks(adj, k, mod)
  350. // -------------------------------------------------------------------
  351. // PURPOSE:
  352. //   Given an adjacency matrix of a directed/undirected graph, computes
  353. //   the matrix (adj^k) where entry (i,j) is the number of walks of
  354. //   exactly length k from node i to node j.
  355. //
  356. // INPUT:
  357. //   adj : square matrix (n x n), where adj[i][j] = number of edges i->j
  358. //   k   : length of walks (non‑negative integer)
  359. //   mod : modulo value
  360. //
  361. // OUTPUT:
  362. //   Returns the matrix adj^k modulo 'mod'.
  363. //
  364. // TIME COMPLEXITY:
  365. //   O(n^3 * log k)
  366. //
  367. // CONSTRAINTS / PRECONDITIONS:
  368. //   - adj must be square.
  369. //   - k >= 0. For k=0, the result is the identity matrix (walk of length 0).
  370. Matrix countWalks(const Matrix& adj, ll k, ll mod) {
  371.     return matPow(adj, k, mod);
  372. }
  373.  
  374. // -------------------------------------------------------------------
  375. // walksBetween(adj, u, v, k, mod)
  376. // -------------------------------------------------------------------
  377. // PURPOSE:
  378. //   Returns the number of walks of exactly length k from node u to node v
  379. //   in the graph described by adjacency matrix 'adj'.
  380. //
  381. // INPUT:
  382. //   adj : square matrix (n x n)
  383. //   u, v: 0‑based indices of nodes
  384. //   k   : length of walk
  385. //   mod : modulo
  386. //
  387. // OUTPUT:
  388. //   Returns adj^k[u][v] modulo 'mod'.
  389. //
  390. // TIME COMPLEXITY:
  391. //   O(n^3 * log k)  (dominated by matPow)
  392. //
  393. // CONSTRAINTS / PRECONDITIONS:
  394. //   - 0 <= u,v < n.
  395. ll walksBetween(const Matrix& adj, int u, int v, ll k, ll mod) {
  396.     Matrix p = matPow(adj, k, mod);
  397.     return p[u][v] % mod;
  398. }
  399.  
  400. // ===================================================================
  401. // SECTION 5: DP WITH MATRIX EXPONENTIATION (Generic Transition)
  402. // ===================================================================
  403.  
  404. // -------------------------------------------------------------------
  405. // applyTransition(T, init, steps, mod)
  406. // -------------------------------------------------------------------
  407. // PURPOSE:
  408. //   Applies a linear transition 'steps' times to an initial state vector.
  409. //   State evolves as: state_{t+1} = T * state_t   (mod 'mod').
  410. //
  411. // INPUT:
  412. //   T     : transition matrix of size (m x m)
  413. //   init  : initial state vector of length m
  414. //   steps : number of transitions to apply (non‑negative integer)
  415. //   mod   : modulo value
  416. //
  417. // OUTPUT:
  418. //   Returns the state vector after 'steps' applications: state_steps = T^steps * init.
  419. //
  420. // TIME COMPLEXITY:
  421. //   O(m^3 * log steps)
  422. //
  423. // CONSTRAINTS / PRECONDITIONS:
  424. //   - T must be square.
  425. //   - init length must equal m.
  426. //   - steps >= 0.
  427. vector<ll> applyTransition(const Matrix& T, const vector<ll>& init,
  428.                             ll steps, ll mod) {
  429.     int m = (int)T.size();
  430.     Matrix Tpow = matPow(T, steps, mod);
  431.     return matVecMul(Tpow, init, mod);
  432. }
  433.  
  434. // ===================================================================
  435. // SECTION 6: SUM OF FIRST n TERMS OF LINEAR RECURRENCE
  436. // ===================================================================
  437.  
  438. // -------------------------------------------------------------------
  439. // sumLinearRecurrence(init, coeff, n, mod)
  440. // -------------------------------------------------------------------
  441. // PURPOSE:
  442. //   Computes the sum of the first n terms of a linear recurrence:
  443. //     S(n) = sum_{i=0}^{n-1} f[i]
  444. //   where f follows the same recurrence as defined in linearRecurrence().
  445. //
  446. // INPUT:
  447. //   init  : initial terms f[0] ... f[k-1] (length k)
  448. //   coeff : recurrence coefficients (length k)
  449. //   n     : number of terms to sum (n >= 0)
  450. //   mod   : modulo value
  451. //
  452. // OUTPUT:
  453. //   Returns S(n) modulo 'mod'.
  454. //
  455. // TIME COMPLEXITY:
  456. //   O(k^3 * log n)  using an augmented matrix of size (k+1)
  457. //
  458. // CONSTRAINTS / PRECONDITIONS:
  459. //   - n >= 0.
  460. //   - If n <= k, it computes the sum directly.
  461. //   - All arithmetic is modulo 'mod'.
  462. ll sumLinearRecurrence(const vector<ll>& init, const vector<ll>& coeff,
  463.                        ll n, ll mod) {
  464.     int k = (int)init.size();
  465.     if (n == 0) return 0;
  466.     if (n <= k) {
  467.         ll s = 0;
  468.         for (int i = 0; i < n; ++i) s = (s + init[i]) % mod;
  469.         return s;
  470.     }
  471.  
  472.     // Build augmented transition matrix of size (k+1) x (k+1)
  473.     // State: [f[t], f[t-1], ..., f[t-k+1], S(t)]^T
  474.     // where S(t) = sum_{i=0}^{t-1} f[i]
  475.     Matrix T(k+1, vector<ll>(k+1, 0));
  476.     // first row for f[t+1]
  477.     for (int i = 0; i < k; ++i) T[0][i] = coeff[i] % mod;
  478.     // shift rows
  479.     for (int i = 1; i < k; ++i) T[i][i-1] = 1;
  480.     // last row: S(t+1) = S(t) + f[t]  =>  T[k][0] = 1, T[k][k] = 1
  481.     T[k][0] = 1;
  482.     T[k][k] = 1;
  483.  
  484.     // Initial state at t = k-1:
  485.     // state[0..k-1] = f[k-1], f[k-2], ..., f[0]
  486.     vector<ll> state(k+1);
  487.     for (int i = 0; i < k; ++i) state[i] = init[k-1-i] % mod;
  488.     // S(k-1) = sum_{i=0}^{k-2} f[i]
  489.     ll sum_init = 0;
  490.     for (int i = 0; i < k-1; ++i) sum_init = (sum_init + init[i]) % mod;
  491.     state[k] = sum_init;
  492.  
  493.     // exponent = n - (k-1)
  494.     Matrix Tpow = matPow(T, n - k + 1, mod);
  495.     vector<ll> res = matVecMul(Tpow, state, mod);
  496.     return res[k] % mod;
  497. }
  498.  
  499. // ===================================================================
  500. // SECTION 7: ADDITIONAL TRICKS / PATTERNS
  501. // ===================================================================
  502.  
  503. // -------------------------------------------------------------------
  504. // (Note) Binary exponentiation for scalars:
  505. //   Already available via std::pow? But not needed, we can use
  506. //   our matPow with 1x1 matrix, or implement a fastPow function.
  507. //   We'll provide a simple fastPow for completeness.
  508. // -------------------------------------------------------------------
  509.  
  510. // fastPow(base, exp, mod): returns (base^exp) % mod.
  511. // Use this when you need scalar exponentiation.
  512. ll fastPow(ll base, ll exp, ll mod) {
  513.     base %= mod;
  514.     ll res = 1 % mod;
  515.     while (exp > 0) {
  516.         if (exp & 1) res = (res * base) % mod;
  517.         base = (base * base) % mod;
  518.         exp >>= 1;
  519.     }
  520.     return res;
  521. }
  522.  
  523. // -------------------------------------------------------------------
  524. // linearRecurrenceWithConstant(init, a, b, n, mod)
  525. // -------------------------------------------------------------------
  526. // PURPOSE:
  527. //   Computes the n‑th term of a recurrence of the form:
  528. //     f(0) = init
  529. //     f(n) = a * f(n-1) + b   for n >= 1
  530. //
  531. // INPUT:
  532. //   init : the initial value f(0)
  533. //   a    : coefficient of f(n-1)
  534. //   b    : constant term added each step
  535. //   n    : index to compute (n >= 0)
  536. //   mod  : modulo value
  537. //
  538. // OUTPUT:
  539. //   Returns f(n) modulo 'mod'.
  540. //
  541. // TIME COMPLEXITY:
  542. //   O(log n)  (using 2x2 matrix exponentiation)
  543. //
  544. // CONSTRAINTS / PRECONDITIONS:
  545. //   - n >= 0.
  546. //   - All values should be non‑negative or reduced modulo 'mod'.
  547. ll linearRecurrenceWithConstant(ll init, ll a, ll b, ll n, ll mod) {
  548.     if (n == 0) return init % mod;
  549.  
  550.     // State vector: [f(t), 1]^T
  551.     // Transition: [f(t+1)]   = [a  b] * [f(t)]
  552.     //             [   1   ]     [0  1]   [  1 ]
  553.     Matrix T = {{a % mod, b % mod}, {0, 1}};
  554.     vector<ll> state = {init % mod, 1};
  555.  
  556.     vector<ll> res = applyTransition(T, state, n, mod);
  557.     return res[0];
  558. }
  559.  
  560. // -------------------------------------------------------------------
  561. // twoSequences(a, b, n, mod)
  562. // -------------------------------------------------------------------
  563. // PURPOSE:
  564. //   Computes the n‑th terms of two sequences defined as:
  565. //     x(0) = 1, y(0) = 0
  566. //     x(n) = p * x(n-1) + q * y(n-1)
  567. //     y(n) = r * x(n-1) + s * y(n-1)
  568. //   This is a generic example; you can change the coefficients and
  569. //   initial values as needed.
  570. //
  571. // INPUT:
  572. //   p, q, r, s : coefficients of the recurrence
  573. //   n           : index to compute (n >= 0)
  574. //   mod         : modulo value
  575. //
  576. // OUTPUT:
  577. //   Returns a pair {x(n), y(n)} modulo 'mod'.
  578. //
  579. // TIME COMPLEXITY:
  580. //   O(log n)  (using 2x2 matrix exponentiation)
  581. //
  582. // CONSTRAINTS / PRECONDITIONS:
  583. //   - n >= 0.
  584. //   - Works for any integer coefficients.
  585. pair<ll, ll> twoSequences(ll p, ll q, ll r, ll s, ll n, ll mod) {
  586.     if (n == 0) return {1 % mod, 0};
  587.  
  588.     // State: [x(t), y(t)]^T
  589.     // Transition: [x(t+1)]   = [p  q] * [x(t)]
  590.     //              [y(t+1)]     [r  s]   [y(t)]
  591.     Matrix T = {{p % mod, q % mod}, {r % mod, s % mod}};
  592.     vector<ll> state = {1 % mod, 0};
  593.  
  594.     vector<ll> res = applyTransition(T, state, n, mod);
  595.     return {res[0], res[1]};
  596. }
  597.  
  598. // -------------------------------------------------------------------
  599. // sumFirstNFibonacci(n, mod)
  600. // -------------------------------------------------------------------
  601. // PURPOSE:
  602. //   Computes the sum of the first n Fibonacci numbers:
  603. //     S(n) = F(0) + F(1) + ... + F(n-1)
  604. //   where F(0)=0, F(1)=1.
  605. //
  606. // INPUT:
  607. //   n   : number of terms to sum (n >= 0)
  608. //   mod : modulo value
  609. //
  610. // OUTPUT:
  611. //   Returns S(n) modulo 'mod'.
  612. //
  613. // TIME COMPLEXITY:
  614. //   O(log n)  (using 3x3 matrix exponentiation)
  615. //
  616. // CONSTRAINTS / PRECONDITIONS:
  617. //   - n >= 0.
  618. //   - If n == 0, returns 0.
  619. ll sumFirstNFibonacci(ll n, ll mod) {
  620.     if (n == 0) return 0;
  621.  
  622.     // State: [F(t), F(t-1), S(t)]^T  where S(t) = sum_{i=0}^{t-1} F(i)
  623.     // Transition for t >= 1:
  624.     // [F(t+1)]   = [1  1  0] * [F(t)]
  625.     // [F(t)  ]     [1  0  0]   [F(t-1)]
  626.     // [S(t+1)]     [1  1  1]   [S(t)  ]
  627.     Matrix T = {{1, 1, 0}, {1, 0, 0}, {1, 1, 1}};
  628.     // Initial state at t = 1: F(1)=1, F(0)=0, S(1)=sum F(0)=0
  629.     vector<ll> state = {1 % mod, 0, 0};
  630.  
  631.     vector<ll> res = applyTransition(T, state, n - 1, mod);
  632.     return res[2];
  633. }
  634.  
  635. // -------------------------------------------------------------------
  636. // countArraysNoThreeConsecutiveEqual(M, n, mod)
  637. // -------------------------------------------------------------------
  638. // PURPOSE:
  639. //   Counts the number of arrays of length n, where each element is
  640. //   between 1 and M (inclusive), and no three consecutive elements
  641. //   are equal.
  642. //
  643. // INPUT:
  644. //   M   : maximum value of each element (M >= 1)
  645. //   n   : length of the array (n >= 1)
  646. //   mod : modulo value
  647. //
  648. // OUTPUT:
  649. //   Returns the number of valid arrays modulo 'mod'.
  650. //
  651. // TIME COMPLEXITY:
  652. //   O(M^3 * log n)  (matrix size is M x M)
  653. //
  654. // CONSTRAINTS / PRECONDITIONS:
  655. //   - Works for small M (e.g., M <= 50) because matrix is M x M.
  656. //   - If n is small (1 or 2), the answer is computed directly.
  657. //   - This is a classic example; for larger M, you need Kitamasa.
  658. ll countArraysNoThreeConsecutiveEqual(int M, ll n, ll mod) {
  659.     if (n == 1) return M % mod;
  660.     if (n == 2) return (M * M) % mod;
  661.  
  662.     // State: dp[len][last][prev] but we can compress.
  663.     // Better approach: Use a 2x2 matrix for this specific problem
  664.     // because the state can be: (number of ways ending with two equal,
  665.     // number of ways ending with two different).
  666.     // But to demonstrate large matrix, we'll keep it general.
  667.  
  668.     // For M up to 50, we can use a 2x2 matrix:
  669.     // State: [ways where last two are equal, ways where last two are different]
  670.     // Transition:
  671.     // new_equal = old_different * 1  (choose the same as last)
  672.     // new_different = old_equal * (M-1) + old_different * (M-2)
  673.     // Because from equal state, you must choose a different element (M-1 choices)
  674.     // from different state, you choose an element different from the last (M-2 choices)
  675.  
  676.     Matrix T = {{0, 1}, {M-1, M-2}};
  677.     // Initial state for length 2:
  678.     // equal ways = M (pairs like (1,1), (2,2), ...)
  679.     // different ways = M * (M-1)
  680.     vector<ll> state = {M % mod, (M * (M-1)) % mod};
  681.  
  682.     vector<ll> res = applyTransition(T, state, n - 2, mod);
  683.     return (res[0] + res[1]) % mod;
  684. }
  685.  
  686. // -------------------------------------------------------------------
  687. // matMulMinPlus(A, B)
  688. // -------------------------------------------------------------------
  689. // PURPOSE:
  690. //   Multiplies two matrices using the min-plus (or tropical) semiring:
  691. //     C[i][j] = min_k ( A[i][k] + B[k][j] )
  692. //   This is used to compute shortest paths after k steps, where
  693. //   the graph has edge weights and you want the minimum total weight
  694. //   of a walk of exactly k edges.
  695. //
  696. // INPUT:
  697. //   A, B : square matrices of the same size, containing edge weights.
  698. //          Use INF (a large number) for no edge.
  699. //
  700. // OUTPUT:
  701. //   Returns the min-plus product matrix C.
  702. //
  703. // TIME COMPLEXITY:
  704. //   O(n^3)
  705. //
  706. // CONSTRAINTS / PRECONDITIONS:
  707. //   - Matrices must be square and of the same size.
  708. //   - INF should be large enough (e.g., 4e18) to avoid overflow.
  709. const ll INF = 4e18;
  710.  
  711. Matrix matMulMinPlus(const Matrix& A, const Matrix& B) {
  712.     int n = (int)A.size();
  713.     Matrix C(n, vector<ll>(n, INF));
  714.     for (int i = 0; i < n; ++i) {
  715.         for (int k = 0; k < n; ++k) {
  716.             if (A[i][k] == INF) continue;
  717.             for (int j = 0; j < n; ++j) {
  718.                 if (B[k][j] == INF) continue;
  719.                 C[i][j] = min(C[i][j], A[i][k] + B[k][j]);
  720.             }
  721.         }
  722.     }
  723.     return C;
  724. }
  725.  
  726. // -------------------------------------------------------------------
  727. // matPowMinPlus(base, exp)
  728. // -------------------------------------------------------------------
  729. // PURPOSE:
  730. //   Raises a square matrix to the power 'exp' using min-plus
  731. //   multiplication. This computes the minimum weight of a walk of
  732. //   exactly 'exp' edges between any two nodes.
  733. //
  734. // INPUT:
  735. //   base : square matrix of edge weights (INF for no edge)
  736. //   exp  : number of edges (exp >= 0)
  737. //
  738. // OUTPUT:
  739. //   Returns base^exp under min-plus multiplication.
  740. //
  741. // TIME COMPLEXITY:
  742. //   O(n^3 * log exp)
  743. //
  744. // CONSTRAINTS / PRECONDITIONS:
  745. //   - exp >= 0.
  746. //   - For exp = 0, the result is the identity matrix for min-plus:
  747. //     C[i][i] = 0, C[i][j] = INF for i != j.
  748. Matrix matPowMinPlus(Matrix base, ll exp) {
  749.     int n = (int)base.size();
  750.     Matrix res(n, vector<ll>(n, INF));
  751.     for (int i = 0; i < n; ++i) res[i][i] = 0; // identity for min-plus
  752.  
  753.     while (exp > 0) {
  754.         if (exp & 1) res = matMulMinPlus(res, base);
  755.         base = matMulMinPlus(base, base);
  756.         exp >>= 1;
  757.     }
  758.     return res;
  759. }
  760.  
  761. // -------------------------------------------------------------------
  762. // PrecomputedPowers
  763. // -------------------------------------------------------------------
  764. // PURPOSE:
  765. //   Precomputes powers of a matrix: P[i] = base^(2^i).
  766. //   Then, for any exponent k, you can compute base^k by multiplying
  767. //   the relevant P[i] matrices.
  768. //
  769. // INPUT:
  770. //   base : square matrix
  771. //   maxExp : maximum exponent you will query (so we precompute up to log2(maxExp))
  772. //   mod : modulo value
  773. //
  774. // OUTPUT:
  775. //   The class provides a method query(exp) that returns base^exp.
  776. //
  777. // TIME COMPLEXITY:
  778. //   Precomputation: O(n^3 * log maxExp)
  779. //   Each query: O(n^3 * popcount(exp)) which is O(n^3 * log maxExp) in worst case.
  780. //
  781. // CONSTRAINTS / PRECONDITIONS:
  782. //   - base must be square.
  783. //   - maxExp >= 0.
  784. class PrecomputedPowers {
  785. private:
  786.     vector<Matrix> powers;
  787.     ll mod;
  788.  
  789. public:
  790.     PrecomputedPowers(const Matrix& base, ll maxExp, ll mod) : mod(mod) {
  791.         powers.push_back(base);
  792.         for (ll e = 2; e <= maxExp; e <<= 1) {
  793.             powers.push_back(matMul(powers.back(), powers.back(), mod));
  794.         }
  795.     }
  796.  
  797.     Matrix query(ll exp) {
  798.         int n = (int)powers[0].size();
  799.         Matrix res(n, vector<ll>(n, 0));
  800.         for (int i = 0; i < n; ++i) res[i][i] = 1 % mod;
  801.  
  802.         int bit = 0;
  803.         while (exp > 0) {
  804.             if (exp & 1) {
  805.                 res = matMul(res, powers[bit], mod);
  806.             }
  807.             exp >>= 1;
  808.             bit++;
  809.         }
  810.         return res;
  811.     }
  812. };
  813.  
  814. // -------------------------------------------------------------------
  815. // matMulSparse(A, B, mod)
  816. // -------------------------------------------------------------------
  817. // PURPOSE:
  818. //   Multiplies two sparse matrices efficiently by skipping zeros.
  819. //   A and B are represented as vector of vectors of pairs (col, value).
  820. //
  821. // INPUT:
  822. //   A, B : sparse matrices, each a vector of size n, where A[i] is a
  823. //          vector of pairs {j, value} for non-zero entries in row i.
  824. //   mod  : modulo value
  825. //
  826. // OUTPUT:
  827. //   Returns the product matrix C as a dense matrix (n x n).
  828. //
  829. // TIME COMPLEXITY:
  830. //   O( (number of non-zero in A) * (average non-zero per row in B) )
  831. //   Much faster than O(n^3) if matrices are sparse.
  832. //
  833. // CONSTRAINTS / PRECONDITIONS:
  834. //   - Matrices must be square.
  835. //   - The input format is sparse; the output is dense.
  836. Matrix matMulSparse(const vector<vector<pair<int, ll>>>& A,
  837.                     const vector<vector<pair<int, ll>>>& B, ll mod) {
  838.     int n = (int)A.size();
  839.     Matrix C(n, vector<ll>(n, 0));
  840.  
  841.     for (int i = 0; i < n; ++i) {
  842.         for (auto &p : A[i]) {
  843.             int k = p.first;
  844.             ll aik = p.second;
  845.             if (aik == 0) continue;
  846.             for (auto &q : B[k]) {
  847.                 int j = q.first;
  848.                 C[i][j] = (C[i][j] + aik * q.second) % mod;
  849.             }
  850.         }
  851.     }
  852.     return C;
  853. }
  854.  
  855. // -------------------------------------------------------------------
  856. // matrixGeometricSum(A, k, mod)
  857. // -------------------------------------------------------------------
  858. // PURPOSE:
  859. //   Computes the sum: S = A + A^2 + A^3 + ... + A^k
  860. //
  861. // INPUT:
  862. //   A   : square matrix
  863. //   k   : number of terms (k >= 1)
  864. //   mod : modulo value
  865. //
  866. // OUTPUT:
  867. //   Returns the matrix S modulo 'mod'.
  868. //
  869. // TIME COMPLEXITY:
  870. //   O(n^3 * log k)
  871. //
  872. // CONSTRAINTS / PRECONDITIONS:
  873. //   - A must be square.
  874. //   - k >= 1.
  875. Matrix matrixGeometricSum(const Matrix& A, ll k, ll mod) {
  876.     int n = (int)A.size();
  877.     // Build augmented matrix of size 2n x 2n:
  878.     // [A  I]
  879.     // [0  I]
  880.     Matrix T(2*n, vector<ll>(2*n, 0));
  881.     for (int i = 0; i < n; ++i) {
  882.         for (int j = 0; j < n; ++j) {
  883.             T[i][j] = A[i][j] % mod;          // top-left: A
  884.             T[i][j+n] = (i == j) ? 1 : 0;     // top-right: I
  885.         }
  886.         T[i+n][i+n] = 1;                      // bottom-right: I
  887.     }
  888.  
  889.     // We need sum_{i=1}^{k} A^i = (sum_{i=0}^{k} A^i) - I.
  890.     // T^(k+1) top-right block = sum_{i=0}^{k} A^i.
  891.     Matrix Tk = matPow(T, k + 1, mod);
  892.  
  893.     Matrix S(n, vector<ll>(n, 0));
  894.     for (int i = 0; i < n; ++i) {
  895.         for (int j = 0; j < n; ++j) {
  896.             S[i][j] = Tk[i][j+n] % mod;
  897.             if (i == j) {
  898.                 S[i][j] = (S[i][j] - 1 + mod) % mod; // subtract I
  899.             }
  900.         }
  901.     }
  902.     return S;
  903. }
  904.  
  905.  
  906. // -------------------------------------------------------------------
  907. // (Note) Exponentiating a matrix with non‑integer exponents? Not used.
  908. // -------------------------------------------------------------------
  909.  
  910. // ===================================================================
  911. // SECTION 8: COMMON TERMS EXPLAINED
  912. // ===================================================================
  913. //
  914. // 1. Matrix Exponentiation:
  915. //    Raising a square matrix to a power using binary exponentiation.
  916. //    Used to accelerate linear recurrences and DP transitions that can
  917. //    be expressed as repeated multiplication by a fixed matrix.
  918. //
  919. // 2. Linear Recurrence:
  920. //    A sequence where each term is a linear combination of previous terms.
  921. //    Example: Fibonacci, Tribonacci, etc.
  922. //
  923. // 3. Transition Matrix / Companion Matrix:
  924. //    A matrix that transforms the state vector from time t to t+1.
  925. //    For a recurrence of order k, the companion matrix is k x k.
  926. //
  927. // 4. State Vector:
  928. //    A column vector containing the current values needed to compute the
  929. //    next values (e.g., last k terms).
  930. //
  931. // 5. Kitamasa / Polynomial Exponentiation:
  932. //    A technique to compute the n‑th term of a linear recurrence without
  933. //    building the full matrix, by computing x^n modulo the characteristic
  934. //    polynomial. Runs in O(k^2 log n) and is better for large k.
  935. //
  936. // 6. Characteristic Polynomial:
  937. //    For recurrence f[n] = c0 f[n-1] + c1 f[n-2] + ... + c[k-1] f[n-k],
  938. //    the characteristic polynomial is P(x) = x^k - c0 x^(k-1) - ... - c[k-1].
  939. //
  940. // 7. Modulo:
  941. //    All operations are performed modulo a given number to prevent overflow
  942. //    and keep numbers small. In competitive programming, MOD is often 1e9+7
  943. //    or 998244353.
  944. //
  945. // ===================================================================
  946.  
  947. // ===================================================================
  948. // EXAMPLE USAGE (you can ignore or uncomment to test)
  949. // ===================================================================
  950.  
  951. int main() {
  952.     ios::sync_with_stdio(false);
  953.     cin.tie(nullptr);
  954.  
  955.     const ll MOD = 1000000007;
  956.  
  957.     // Fibonacci
  958.     cout << fib(10, MOD) << "\n";            // 55
  959.  
  960.     // Linear recurrence: Fibonacci (init [0,1], coeff [1,1])
  961.     vector<ll> init = {0, 1};
  962.     vector<ll> coeff = {1, 1};
  963.     cout << linearRecurrence(init, coeff, 10, MOD) << "\n";  // 55
  964.     cout << linearRecurrenceKitamasa(init, coeff, 10, MOD) << "\n"; // 55
  965.  
  966.     // Sum of first 5 Fibonacci numbers: 0+1+1+2+3 = 7
  967.     cout << sumLinearRecurrence(init, coeff, 5, MOD) << "\n"; // 7
  968.  
  969.     // Graph walks: simple directed graph with 2 nodes, edge 0->1 and 1->0
  970.     Matrix adj = {{0,1},{1,0}};
  971.     Matrix p = countWalks(adj, 3, MOD);
  972.     // p[0][1] = number of walks of length 3 from 0 to 1.
  973.     cout << p[0][1] << "\n"; // 1? Actually for two nodes, length 3 from 0 to 1: 0->1->0->1 => 1, 0->1->0->? Only one.
  974.  
  975.     return 0;
  976. }
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
55
55
55
7
1
https://ideone.com/42NGJS
language:
C++ (gcc 8.3)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!