API tools faq
paste
Advertisement
prog3r

2.0 Transform recursive formula into a matrix

prog3r
Apr 22nd, 2025
486
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 12.68 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 = pair<int,vector<Single>>;
  6. using MultiZeroCoef = vector<Single>;
  7. using PolyM = vector<Multi>;
  8. map<MultiZeroCoef, PolyM> Eq;
  9. int LAST_IDX_SET = -1e9;
  10. constexpr int MOD = 1e9+7;
  11. vector<vector<int>> temp_m(1000, vector<int>(1000));
  12. ostream& operator<<(ostream& out, const vector<vector<int>>& mat) {
  13.     for (const auto &x : mat) {
  14.         for (const auto &y : x) {
  15.             out << y << " ";
  16.         }out << "\n";
  17.     }
  18.     return out;
  19. }
  20. void operator*=(vector<vector<int>> &a, const vector<vector<int>>& b) {
  21.     for (const auto &x : a) {
  22.         assert(x.size() == a[0].size());
  23.     }
  24.     for (const auto &y : b) {
  25.         assert(y.size() == b[0].size());
  26.     }
  27.     int skal_size = a[0].size();
  28.     assert(a[0].size() == b.size());
  29.     for (int i = 0; i < a.size(); i += 1) {
  30.         for (int j = 0; j < b[0].size(); j += 1) {
  31.             temp_m[i][j] = 0;
  32.             for (int k = 0; k < skal_size; k += 1) {
  33.                 temp_m[i][j] += a[i][k] * b[k][j] % MOD;
  34.             }
  35.             temp_m[i][j] %= MOD;
  36.         }
  37.     }
  38.     for (auto &x : a) {
  39.         x.resize(b[0].size());
  40.     }
  41.     for (int i = 0; i < a.size(); i += 1) {
  42.         for (int j = 0; j < a[i].size(); j += 1) {
  43.             a[i][j] = temp_m[i][j];
  44.         }
  45.     }
  46. }
  47. void operator^=(vector<vector<int>>& mat, int b) {
  48.     assert(b >= 0);
  49.     assert(mat.size() == mat[0].size());
  50.     vector<vector<int>> ret(mat.size(), vector<int>(mat.size()));
  51.     for (int i = 0; i < mat.size(); i += 1) {
  52.         ret[i][i] = 1;
  53.     }
  54.     while (b) {
  55.         if (b&1) {
  56.             ret *= mat;
  57.         }
  58.         mat *= mat;
  59.         b >>= 1;
  60.     }
  61.     mat = ret;
  62. }
  63. int query_fast(vector<vector<int>> Matrix, const int n, const vector<vector<int>>& base,
  64.                const vector<MultiZeroCoef>& list, vector<int>* f) {
  65.     if (n <= LAST_IDX_SET) {
  66.         return f->at(n);
  67.     }
  68.     Matrix ^= n - LAST_IDX_SET;
  69.     Matrix *= base;
  70.     for (int i = 0; i < Matrix.size(); i += 1) {
  71.         if (list[i] == MultiZeroCoef{{f, 0}}) {
  72.             return Matrix[i][0];
  73.         }
  74.     }
  75.     assert(false);
  76. }
  77. int query_slow(const int n, const map<MultiZeroCoef, PolyM>& slow_eq, vector<int>* f) {
  78.     // better run twice, first return value shouldn't be considered
  79.     if (n <= LAST_IDX_SET) {
  80.         return f->at(n);
  81.     }
  82.     int curr = LAST_IDX_SET+1;
  83.     for (const auto &[a, b] : slow_eq) {
  84.         for (const auto &S : a) {
  85.             S.first->resize(n+1);
  86.         }
  87.         for (const auto &M : b) {
  88.             for (const auto &S : M.second) {
  89.                 if (!(curr+S.second >= 0)) {
  90.                     throw invalid_argument("Недостаточно начальных значений");
  91.                 }
  92.             }
  93.         }
  94.     }
  95.     while (curr <= n) {
  96.         for (const auto &[a, b] : slow_eq) {
  97.             assert(a.size() == 1);
  98.             int su = 0;
  99.             for (const auto &M : b) {
  100.                 int mul = M.first;
  101.                 for (const auto &S : M.second) {
  102.                     (mul *= S.first->at(curr+S.second)) %= MOD;
  103.                 }
  104.                 (su += mul) %= MOD;
  105.             }
  106.             assert(a.front().second == 0);
  107.             a.front().first->at(curr) = su;
  108.         }
  109.         curr += 1;
  110.     }
  111.     return f->at(n);
  112. }
  113. signed main() {
  114.     vector<int> f(200), r(200), sz(200), fib(200), nech(200), I2(200), I(200);
  115.     vector<int> none(200, 0), one(200, 1);
  116.     Eq[{{&none, 0}}] = {
  117.             {1, {{&none, -1},}},
  118.     };
  119.     Eq[{{&one, 0}}] = {
  120.             {1, {{&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.                 {1, {{&sz, -1},}},
  132.                 {1, {{&sz, -2},}},
  133.                 {1, {{&one, -1},}},
  134.         };
  135.         Eq[{{&r, 0}}] = {
  136.                 {1, {{&sz, -1},}},
  137.                 {1, {{&r, -1},}},
  138.                 {1, {{&sz, -2},}},
  139.                 {1, {{&r, -2},}},
  140.         };
  141.         Eq[{{&f, 0}}] = {
  142.                 {1, {{&f, -1}}},
  143.                 {1, {{&f, -2}}},
  144.  
  145.                 {1, {{&sz, -1}, {&r, -2}}},
  146.                 {1, {{&sz, -1}, {&sz, -2}}},
  147.                 {1, {{&r, -2},}},
  148.                 {1, {{&sz, -2},}},
  149.  
  150.                 {1, {{&sz, -2}, {&r, -1}}},
  151.                 {1, {{&sz, -2}, {&sz, -1}}},
  152.                 {1, {{&r, -1},}},
  153.                 {1, {{&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.                 {1, {{&f, -1}}},
  162.                 {1, {{&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.                 {1, {{&fib, -1}}},
  174.                 {1, {{&fib, -2}}},
  175.         };
  176.         Eq[{{&f, 0}}] = {
  177.                 {1, {{&fib, 0}, {&fib, -1}}},
  178.                 {1, {{&f, -1}}},
  179.         };
  180.     };
  181.     int QUERY_TEMP; auto _102644G = [&] () -> void {
  182.         int n;
  183.         cin >> n >> QUERY_TEMP;
  184.         for (int i = 0; i < n; i += 1) {
  185.             cin >> f[i];
  186.         }
  187.         for (int i = 0; i < n; i += 1) {
  188.             int x;
  189.             cin >> x;
  190.             Eq[{{&f, 0}}].push_back({x, {{&f, -i-1}}});
  191.         }
  192.         LAST_IDX_SET = n-1;
  193.         for (int i = 0; i < n; i += 1) {
  194.             nech[i] = ((2*i-1)%MOD+MOD)%MOD;
  195.             I2[i] = i*i;
  196.         }
  197.         Eq[{{&nech, 0}}] = {
  198.                 {1, {{&nech, -1}}},
  199.                 {1, {{&one, -1}}},
  200.                 {1, {{&one, -1}}},
  201.         };
  202.         iota(I.begin(), I.end(), 0);
  203.         Eq[{{&I, 0}}] = {
  204.                 {1, {{&I, -1}}},
  205.                 {1, {{&one, -1}}},
  206.         };
  207.         Eq[{{&I2, 0}}] = {
  208.                 {1, {{&I2, -1}}},
  209.                 {1, {{&nech, 0}}},
  210.         };
  211.         int p, q, r;
  212.         cin >> p >> q >> r;
  213.         Eq[{{&f, 0}}].push_back({p, {{&one, -1}}});
  214.         Eq[{{&f, 0}}].push_back({q, {{&I, 0}}});
  215.         Eq[{{&f, 0}}].push_back({r, {{&I2, 0}}});
  216.     };
  217.     test_on_fibonacci();
  218. //    test_on_fibonacci_trees();
  219. //    another_test();
  220. //    _102644G();
  221.     if (Eq.empty()) {
  222.         throw invalid_argument("No recursive formulas given (Не даны рекуррентные уравнения)");
  223.     }
  224.     auto slow_Eq = Eq;
  225.     map<vector<int>*, string> mp;
  226.     mp[&none] = "NULL";
  227.     mp[&one] = "ONE";
  228.     mp[&sz] = "sz";
  229.     mp[&r] = "r";
  230.     mp[&f] = "f";
  231.     mp[&I] = "I";
  232.     mp[&I2] = "I2";
  233.     mp[&nech] = "nech";
  234.     auto check = [&] () -> void {
  235.         for (auto &[a, b] : Eq) {
  236.             for (auto &M : b) {
  237.                 sort(M.second.begin(), M.second.end());
  238.             }
  239.             sort(b.begin(), b.end());
  240.             if (a.size() == 1 && (a.front().first == &none || a.front().first == &one)) {
  241.                 continue;
  242.             }
  243.             PolyM nw_b;
  244.             for (const auto &M : b) {
  245.                 Multi nw_m;
  246.                 nw_m.first = M.first;
  247.                 bool bad = false;
  248.                 bool got_one = false;
  249.                 for (const auto &S : M.second) {
  250.                     if (S.first != &one) {
  251.                         nw_m.second.push_back(S);
  252.                     } else {
  253.                         got_one = true;
  254.                     }
  255.                     bad |= S.first == &none;
  256.                 }
  257.                 if (!bad) {
  258.                     if(!nw_m.second.size() && got_one) {
  259.                         nw_m.second.push_back({&one, -1});
  260.                     }
  261.                     if (nw_m.second.size()) {
  262.                         nw_b.push_back(nw_m);
  263.                     }
  264.                 }
  265.             }
  266.             b = nw_b;
  267.             for (const auto &S : a) {
  268.                 assert(S.second <= 0);
  269.             }
  270.             for (const auto &M : b) {
  271.                 for (const auto &S : M.second) {
  272.                     assert(S.second <= 0);
  273.                 }
  274.             }
  275.         }
  276.         while (true) {
  277.             set<Multi> add;
  278.             for (const auto &[a, b] : Eq) {
  279.                 for (const auto &M : b) {
  280.                     bool good = true;
  281.                     for (const auto &S : M.second) {
  282.                         good &= S.second < 0;
  283.                     }
  284.                     if (good) {
  285.                         auto nw = M;
  286.                         for (auto &S : nw.second) {
  287.                             S.second += 1;
  288.                         }
  289.                         if (!Eq.count(nw.second)) {
  290.                             add.insert(nw);
  291.                         }
  292.                     }
  293.                 }
  294.             }
  295.             for (const auto &y : add) {
  296.                 Eq[y.second] = {{1, {y.second}}};
  297.             }
  298.             if (!add.size()) {
  299.                 break;
  300.             }
  301.         }
  302.     };
  303. //    auto print = [&] () -> void {
  304. //        for (const auto &[a, b] : Eq) {
  305. //            for (const auto &S : a) {
  306. //                clog << mp[S.first] << "[" << S.second << "]";
  307. //            }
  308. //            clog << " = ";
  309. //            for (int i = 0; i < b.size(); i += 1) {
  310. //                clog << b[i].first;
  311. //                for (const auto &S : b[i].second) {
  312. //                    clog << mp[S.first] << "[" << S.second << "]";
  313. //                }
  314. //                if (i+1 != b.size()) {
  315. //                    clog << " + ";
  316. //                }
  317. //            } clog << "\n";
  318. //        }
  319. //        clog << "--------\n";
  320. //    };
  321.     while (true) {
  322.         check();
  323. //        print();
  324.         bool done = false;
  325.         for (auto &[a, b] : Eq) {
  326.             for (auto &M : b) {
  327.                 for (auto &S : M.second) {
  328.                     if (S.second == 0) {
  329.                         assert(Eq.count({S}));
  330.                         auto replace = Eq[{S}];
  331.                         S.first = &one;
  332.                         S.second = -1;
  333.                         for (auto &h: replace) {
  334.                             (h.first *= M.first) %= MOD;
  335.                         }
  336.                         for (const auto &fromM: M.second) {
  337.                             for (auto &h: replace) {
  338.                                 h.second.push_back(fromM);
  339.                             }
  340.                         }
  341.                         S.first = &none;
  342.                         S.second = -1;
  343.                         for (const auto &x: replace) {
  344.                             b.push_back(x);
  345.                         }
  346.                         done = true;
  347.                         break;
  348.                     }
  349.                 } if (done) break;
  350.             } if (done) break;
  351.         }
  352.         if (!done) {
  353.             break;
  354.         }
  355.     }
  356.     map<MultiZeroCoef, int> order;
  357.     vector<MultiZeroCoef> list;
  358.     for (const auto &[a, b] : Eq) {
  359.         assert(!order.count(a));
  360.         order[a] = order.size();
  361.         list.push_back(a);
  362.     }
  363.     vector<vector<int>> Matrix(order.size(), vector<int>(order.size()));
  364.     for (const auto &[a, b] : Eq) {
  365.         assert(order.count(a));
  366.         for (const auto &S : a) {
  367.             if (!(LAST_IDX_SET >= -S.second)) {
  368.                 throw invalid_argument("Недостаточно начальных значений");
  369.             }
  370.         }
  371.         for (const auto &M : b) {
  372.             auto from = M;
  373.             for (auto &x : from.second) {
  374.                 x.second += 1;
  375.             }
  376.             assert(order.count(from.second));
  377.             Matrix[order[a]][order[from.second]] += M.first;
  378.         }
  379.     }
  380.     vector<vector<int>> base;
  381.     for (const auto &S : list) {
  382.         int tot = 1;
  383.         for (const auto &[ptr, diff] : S) {
  384.             (tot *= ptr->at(LAST_IDX_SET+diff)) %= MOD;
  385.         }
  386.         base.push_back({tot});
  387.     }
  388.     int q = 1e3;
  389.     q = 1;
  390.     mt19937_64 mt(chrono::high_resolution_clock::now().time_since_epoch().count());
  391.     uniform_int_distribution<int> distrib(1, 1e4);
  392.     for (int ii = 0; ii < q; ii += 1) {
  393.         int n = distrib(mt);
  394.         cin >> n;
  395. //        query_slow(n, slow_Eq, &f);
  396. //        cout << query_fast(Matrix, n, base, list, &f) << "\n";
  397. //        clog << "f[" << n << "]: ";
  398. //        clog << query_fast(Matrix, n, base, list, &f) << ", "
  399. //             << query_slow(n, slow_Eq, &f) << "\n";
  400. //        assert(query_fast(Matrix, n, base, list, &f) == query_slow(n, slow_Eq, &f));
  401.         cout << query_fast(Matrix, n, base, list, &f) << "\n";
  402.     }
  403. }
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!