API tools faq
paste
Advertisement
prog3r

Transform recursive formula into a matrix

prog3r
Apr 21st, 2025
454
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 10.89 KB | None | 0 0
raw download clone embed print report
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define int long long
  4. using Single = pair<vector<int>*, int>;
  5. using Multi = vector<Single>;
  6. using PolyM = vector<Multi>;
  7. map<Multi, PolyM> Eq;
  8. int LAST_IDX_SET = -1e9;
  9. constexpr int MOD = 1e9+7;
  10. vector<vector<int>> temp_m(1000, vector<int>(1000));
  11. ostream& operator<<(ostream& out, const vector<vector<int>>& mat) {
  12.     for (const auto &x : mat) {
  13.         for (const auto &y : x) {
  14.             out << y << " ";
  15.         }out << "\n";
  16.     }
  17.     return out;
  18. }
  19. void operator*=(vector<vector<int>> &a, const vector<vector<int>>& b) {
  20.     for (const auto &x : a) {
  21.         assert(x.size() == a[0].size());
  22.     }
  23.     for (const auto &y : b) {
  24.         assert(y.size() == b[0].size());
  25.     }
  26.     int skal_size = a[0].size();
  27.     assert(a[0].size() == b.size());
  28.     for (int i = 0; i < a.size(); i += 1) {
  29.         for (int j = 0; j < b[0].size(); j += 1) {
  30.             temp_m[i][j] = 0;
  31.             for (int k = 0; k < skal_size; k += 1) {
  32.                 temp_m[i][j] += a[i][k] * b[k][j] % MOD;
  33.             }
  34.             temp_m[i][j] %= MOD;
  35.         }
  36.     }
  37.     for (auto &x : a) {
  38.         x.resize(b[0].size());
  39.     }
  40.     for (int i = 0; i < a.size(); i += 1) {
  41.         for (int j = 0; j < a[i].size(); j += 1) {
  42.             a[i][j] = temp_m[i][j];
  43.         }
  44.     }
  45. }
  46. void operator^=(vector<vector<int>>& mat, int b) {
  47.     assert(b >= 0);
  48.     assert(mat.size() == mat[0].size());
  49.     vector<vector<int>> ret(mat.size(), vector<int>(mat.size()));
  50.     for (int i = 0; i < mat.size(); i += 1) {
  51.         ret[i][i] = 1;
  52.     }
  53.     while (b) {
  54.         if (b&1) {
  55.             ret *= mat;
  56.         }
  57.         mat *= mat;
  58.         b >>= 1;
  59.     }
  60.     mat = ret;
  61. }
  62. int query_fast(vector<vector<int>> Matrix, const int n, const int mini, const vector<vector<int>>& base,
  63.                const vector<Multi>& list, vector<int>* f) {
  64.     if (n <= LAST_IDX_SET) {
  65.         return f->at(n);
  66.     }
  67.     if (!(n >= -mini)) {
  68.         throw invalid_argument("Недостаточно начальных значений");
  69.     }
  70.     Matrix ^= n - (-mini);
  71.     Matrix *= base;
  72.     for (int i = 0; i < Matrix.size(); i += 1) {
  73.         if (list[i] == Multi{{f, 0}}) {
  74.             return Matrix[i][0];
  75.         }
  76.     }
  77.     assert(false);
  78. }
  79. int query_slow(const int n, const map<Multi, PolyM>& slow_eq, vector<int>* f) {
  80.     if (n <= LAST_IDX_SET) {
  81.         return f->at(n);
  82.     }
  83.     int curr = LAST_IDX_SET+1;
  84.     for (const auto &[a, b] : slow_eq) {
  85.         for (const auto &S : a) {
  86.             S.first->resize(n+1);
  87.         }
  88.         for (const auto &M : b) {
  89.             for (const auto &S : M) {
  90.                 if (!(curr+S.second >= 0)) {
  91.                     throw invalid_argument("Недостаточно начальных значений");
  92.                 }
  93.             }
  94.         }
  95.     }
  96.     while (curr <= n) {
  97.         for (const auto &[a, b] : slow_eq) {
  98.             assert(a.size() == 1);
  99.             int su = 0;
  100.             for (const auto &M : b) {
  101.                 int mul = 1;
  102.                 for (const auto &S : M) {
  103.                     (mul *= S.first->at(curr+S.second)) %= MOD;
  104.                 }
  105.                 (su += mul) %= MOD;
  106.             }
  107.             assert(a.front().second == 0);
  108.             a.front().first->at(curr) = su;
  109.         }
  110.         curr += 1;
  111.     }
  112.     return f->at(n);
  113. }
  114. signed main() {
  115.     vector<int> f(200), r(200), sz(200), fib(200), none(200, 0), one(200, 1);
  116.     Eq[{{&none, 0}}] = {
  117.             {{&none, -1},},
  118.     };
  119.     Eq[{{&one, 0}}] = {
  120.             {{&one, -1},},
  121.     };
  122.     auto test_on_fibonacci_trees = [&] () -> void {
  123.         f[0] = 0;
  124.         f[1] = 0;
  125.         r[0] = 0;
  126.         r[1] = 0;
  127.         sz[0] = 1;
  128.         sz[1] = 1;
  129.         LAST_IDX_SET = 1;
  130.         Eq[{{&sz, 0}}] = {
  131.                 {{&sz, -1},},
  132.                 {{&sz, -2},},
  133.                 {{&one, -1},},
  134.         };
  135.         Eq[{{&r, 0}}] = {
  136.                 {{&sz, -1},},
  137.                 {{&r, -1},},
  138.                 {{&sz, -2},},
  139.                 {{&r, -2},},
  140.         };
  141.         Eq[{{&f, 0}}] = {
  142.                 {{&f, -1}},
  143.                 {{&f, -2}},
  144.  
  145.                 {{&sz, -1}, {&r, -2}},
  146.                 {{&sz, -1}, {&sz, -2}},
  147.                 {{&r, -2},},
  148.                 {{&sz, -2},},
  149.  
  150.                 {{&sz, -2}, {&r, -1}},
  151.                 {{&sz, -2}, {&sz, -1}},
  152.                 {{&r, -1},},
  153.                 {{&sz, -1},},
  154.         };
  155.     };
  156.     auto test_on_fibonacci = [&] () -> void {
  157.         f[0] = 0;
  158.         f[1] = 1;
  159.         LAST_IDX_SET = 1;
  160.         Eq[{{&f, 0}}] = {
  161.                 {{&f, -1}},
  162.                 {{&f, -2}},
  163.         };
  164.     };
  165.     auto another_test = [&] () -> void {
  166.         // f[i] = (fib[i] * fib[i-1] + f[i-1]) % MOD;
  167.         fib[0] = 0;
  168.         fib[1] = 1;
  169.         f[0] = 0;
  170.         f[1] = 0;
  171.         LAST_IDX_SET = 1;
  172.         Eq[{{&fib, 0}}] = {
  173.                 {{&fib, -1}},
  174.                 {{&fib, -2}},
  175.         };
  176.         Eq[{{&f, 0}}] = {
  177.                 {{&fib, 0}, {&fib, -1}},
  178.                 {{&f, -1}},
  179.         };
  180.     };
  181. //    test_on_fibonacci();
  182. //    test_on_fibonacci_trees();
  183.     another_test();
  184.     if (Eq.empty()) {
  185.         throw invalid_argument("No recursive formulas given (Не даны рекуррентные уравнения)");
  186.     }
  187.     auto slow_Eq = Eq;
  188.     map<vector<int>*, string> mp;
  189.     mp[&none] = "NULL";
  190.     mp[&one] = "ONE";
  191.     mp[&sz] = "sz";
  192.     mp[&r] = "r";
  193.     mp[&f] = "f";
  194.     auto check = [&] () -> void {
  195.         for (auto &[a, b] : Eq) {
  196.             for (auto &M : b) {
  197.                 sort(M.begin(), M.end());
  198.             }
  199.             sort(b.begin(), b.end());
  200.             if (a.size() == 1 && (a.front().first == &none || a.front().first == &one)) {
  201.                 continue;
  202.             }
  203.             PolyM nw_b;
  204.             for (const auto &M : b) {
  205.                 Multi nw_m;
  206.                 bool bad = false;
  207.                 bool got_one = false;
  208.                 for (const auto &S : M) {
  209.                     if (S.first != &one) {
  210.                         nw_m.push_back(S);
  211.                     } else {
  212.                         got_one = true;
  213.                     }
  214.                     bad |= S.first == &none;
  215.                 }
  216.                 if (!bad) {
  217.                     if(!nw_m.size() && got_one) {
  218.                         nw_m.push_back({&one, -1});
  219.                     }
  220.                     if (nw_m.size()) {
  221.                         nw_b.push_back(nw_m);
  222.                     }
  223.                 }
  224.             }
  225.             b = nw_b;
  226.             for (const auto &S : a) {
  227.                 assert(S.second <= 0);
  228.             }
  229.             for (const auto &M : b) {
  230.                 for (const auto &S : M) {
  231.                     assert(S.second <= 0);
  232.                 }
  233.             }
  234.         }
  235.         while (true) {
  236.             set<Multi> add;
  237.             for (const auto &[a, b] : Eq) {
  238.                 for (const auto &M : b) {
  239.                     bool good = true;
  240.                     for (const auto &S : M) {
  241.                         good &= S.second < 0;
  242.                     }
  243.                     if (good) {
  244.                         auto nw = M;
  245.                         for (auto &S : nw) {
  246.                             S.second += 1;
  247.                         }
  248.                         if (!Eq.count(nw)) {
  249.                             add.insert(nw);
  250.                         }
  251.                     }
  252.                 }
  253.             }
  254.             for (const auto &y : add) {
  255.                 Eq[y] = {y};
  256.             }
  257.             if (!add.size()) {
  258.                 break;
  259.             }
  260.         }
  261.     };
  262. //    auto print = [&] () -> void {
  263. //        for (const auto &[a, b] : Eq) {
  264. //            for (const auto &S : a) {
  265. //                clog << mp[S.first] << "[" << S.second << "]";
  266. //            }
  267. //            clog << " = ";
  268. //            for (int i = 0; i < b.size(); i += 1) {
  269. //                for (const auto &S : b[i]) {
  270. //                    clog << mp[S.first] << "[" << S.second << "]";
  271. //                }
  272. //                if (i+1 != b.size()) {
  273. //                    clog << " + ";
  274. //                }
  275. //            } clog << "\n";
  276. //        }
  277. //        clog << "--------\n";
  278. //    };
  279.     while (true) {
  280.         check();
  281. //        print();
  282.         bool done = false;
  283.         for (auto &[a, b] : Eq) {
  284.             for (auto &M : b) {
  285.                 for (auto &S : M) {
  286.                     if (S.second == 0) {
  287.                         assert(Eq.count({S}));
  288.                         auto replace = Eq[{S}];
  289.                         S.first = &one;
  290.                         S.second = -1;
  291.                         for (const auto &fromM: M) {
  292.                             for (auto &h: replace) {
  293.                                 h.push_back(fromM);
  294.                             }
  295.                         }
  296.                         S.first = &none;
  297.                         S.second = -1;
  298.                         for (const auto &x: replace) {
  299.                             b.push_back(x);
  300.                         }
  301.                         done = true;
  302.                         break;
  303.                     }
  304.                 } if (done) break;
  305.             } if (done) break;
  306.         }
  307.         if (!done) {
  308.             break;
  309.         }
  310.     }
  311.     map<Multi, int> order;
  312.     vector<Multi> list;
  313.     for (const auto &[a, b] : Eq) {
  314.         assert(!order.count(a));
  315.         order[a] = order.size();
  316.         list.push_back(a);
  317.     }
  318.     vector<vector<int>> Matrix(order.size(), vector<int>(order.size()));
  319.     int mini = 1e9;
  320.     for (const auto &[a, b] : Eq) {
  321.         assert(order.count(a));
  322.         for (const auto &S : a) {
  323.             mini = min(mini, S.second);
  324.         }
  325.         for (const auto &M : b) {
  326.             auto from = M;
  327.             for (auto &x : from) {
  328.                 x.second += 1;
  329.             }
  330.             assert(order.count(from));
  331.             Matrix[order[a]][order[from]] += 1;
  332.         }
  333.     }
  334.     if (!(LAST_IDX_SET >= -mini)) {
  335.         throw invalid_argument("Недостаточно начальных значений");
  336.     }
  337.     vector<vector<int>> base;
  338.     for (const auto &S : list) {
  339.         int tot = 1;
  340.         for (const auto &[ptr, diff] : S) {
  341.             (tot *= ptr->at(-mini+diff)) %= MOD;
  342.         }
  343.         base.push_back({tot});
  344.     }
  345.     int q = 1e3;
  346.     q = 1;
  347.     mt19937_64 mt(chrono::high_resolution_clock::now().time_since_epoch().count());
  348.     uniform_int_distribution<int> distrib(1, 1e4);
  349.     for (int ii = 0; ii < q; ii += 1) {
  350.         int n = distrib(mt);
  351.         cin >> n;
  352.         cout << query_fast(Matrix, n, mini, base, list, &f) << "\n";
  353. //        clog << "f[" << n << "]: ";
  354. //        clog << query_fast(Matrix, n, mini, base, list, &f) << ", "
  355. //             << query_slow(n, slow_Eq, &f) << "\n";
  356. //        assert(query_fast(Matrix, n, mini, base, list, &f) == query_slow(n, slow_Eq, &f));
  357.     }
  358. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!