SHARE
TWEET

afs3wq

a guest Jan 14th, 2020 98 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include<bits/stdc++.h>
  2.  
  3. using namespace std;
  4.  
  5. typedef vector < double > vector_t;
  6. typedef vector < vector < double > > matrix_t;
  7.  
  8. vector_t make_vector(size_t n) {
  9.     return vector_t(n);
  10. }
  11. matrix_t make_matrix(size_t x, size_t y) {
  12.     return matrix_t(x, vector_t(y));
  13. }
  14.  
  15. void print_vector(const vector_t &v) {
  16.     for (const auto &i : v) {
  17.         fprintf(stderr, "%.6f ", i);
  18.     }
  19.     cerr << "\n";
  20. }
  21. void print_matrix(const matrix_t &m) {
  22.     for (const auto &i : m) {
  23.         print_vector(i);
  24.     }
  25.     cerr << "\n";
  26. }
  27.  
  28. vector_t& mult(vector_t &m, double v) {
  29.     for (auto &i : m) {
  30.         i *= v;
  31.     }
  32.     return m;
  33. }
  34. vector_t& div(vector_t &m, double v) {
  35.     for (auto &i : m) {
  36.         i /= v;
  37.     }
  38.     return m;
  39. }
  40. vector_t& add(vector_t &m, double v) {
  41.     for (auto &i : m) {
  42.         i += v;
  43.     }
  44.     return m;
  45. }
  46.  
  47. vector_t& add(vector_t &a, const vector_t &b) {
  48.     for (size_t i = 0; i < a.size(); ++i) {
  49.         a[i] += b[i];
  50.     }
  51.     return a;
  52. }
  53.  
  54. matrix_t minor(const matrix_t &m, size_t a, size_t b) {
  55.     matrix_t res = make_matrix(0, 0);
  56.     for (size_t i = 0; i < m.size(); ++i) {
  57.         if (i == a) {
  58.             continue;
  59.         }
  60.         res.push_back(make_vector(0));
  61.         for (size_t j = 0; j < m.size(); ++j) {
  62.             if (j != b) {
  63.                 res[res.size() - 1].push_back(m[i][j]);
  64.             }
  65.         }
  66.     }
  67.     return res;
  68. }
  69. double det(const matrix_t &m) {
  70.     if (m.size() == 1) {
  71.         return m[0][0];
  72.     }
  73.     double res = 0;
  74.     for (size_t i = 0; i < m.size(); ++i) {
  75.         res += (i & 1 ? -1 : 1) * det(minor(m, 0, i)) * m[0][i];
  76.     }
  77.     return res;
  78. }
  79.  
  80. matrix_t inplace(const matrix_t &m, const vector_t &v, size_t n) {
  81.     matrix_t res = m;
  82.     for (size_t i = 0; i < m.size(); ++i) {
  83.         res[i][n] = v[i];
  84.     }
  85.     return res;
  86. }
  87. vector_t calculate(const matrix_t &m, const vector_t &v) {
  88.     vector_t res = make_vector(m.size());
  89.     double detrem = det(m);
  90.     for (size_t i = 0; i < m.size(); ++i) {
  91.         double t = det(inplace(m, v, i));
  92.         res[i] = detrem / t;
  93.     }
  94.     return res;
  95. }
  96.  
  97. matrix_t mult(const matrix_t &a, const matrix_t &b) {
  98.     matrix_t r = make_matrix(a.size(), b[0].size());
  99.     for (size_t i = 0; i < a.size(); ++i) {
  100.         for (size_t j = 0; j < b[0].size(); ++j) {
  101.             for (size_t l = 0; l < b.size(); ++l) {
  102.                 r[i][j] += a[i][l] * b[l][j];
  103.             }
  104.         }
  105.     }
  106.     return r;
  107. }
  108.  
  109. const matrix_t E {
  110.     { 1, 0, 0, 0 },
  111.     { 0, 1, 0, 0 },
  112.     { 0, 0, 1, 0 },
  113.     { 0, 0, 0, 1 }
  114. };
  115.  
  116. matrix_t genrevB(const matrix_t &A, size_t n) {
  117.     matrix_t res = E;
  118.     res[n] = A[n + 1];
  119.     return res;
  120. }
  121. matrix_t genB(const matrix_t &A, size_t n) {
  122.     matrix_t res = genrevB(A, n);
  123.     mult(res[n], -1);
  124.     res[n][n] = 1;
  125.     div(res[n], A[n + 1][n]);
  126.     return res;
  127. }
  128.  
  129. matrix_t& iter(matrix_t &A, size_t n) {
  130.     matrix_t B = genB(A, n), revB = genrevB(A, n);
  131.     cerr << "revB:\n";
  132.     print_matrix(revB);
  133.     cerr << "B:\n";
  134.     print_matrix(B);
  135.     return A = mult(mult(revB, A), B);
  136. }
  137.  
  138. void frubenious(matrix_t &A) {
  139.     cerr << "A:\n";
  140.     print_matrix(A);
  141.     for (int i = A.size() - 2; i >= 0; --i) {
  142.         iter(A, i);
  143.         cerr << "A:\n";
  144.         print_matrix(A);
  145.     }
  146. }
  147.  
  148. void release_lambda(matrix_t &A, double lambda) {
  149.     for (size_t i = 0; i < A.size(); ++i) {
  150.         A[i][i] -= lambda;
  151.     }
  152. }
  153.  
  154. const vector_t lambda_list { 5.811869, 4.088694, 9.482340, 2.577096 };
  155.  
  156. matrix_t get_self_vector(const matrix_t &A, size_t lambda) {
  157.     matrix_t res = make_matrix(A.size(), 1);
  158.     res[A.size() - 1][0] = lambda_list[lambda];
  159.     for (int i = A.size() - 2; i >= 0; --i) {
  160.         res[i][0] = res[i + 1][0] * lambda_list[lambda];
  161.     }
  162.     return res;
  163. }
  164.  
  165. void print_self_vectors(const matrix_t &A) {
  166.     for (size_t i = 0; i < A.size(); ++i) {
  167.         cerr << "\n";
  168.         matrix_t t = get_self_vector(A, i),
  169.                  lambda = { { lambda_list[i] } };
  170.         print_matrix(t);
  171. //      print_matrix(mult(t, lambda));
  172. //      print_matrix(mult(A, t));
  173.     }
  174. }
  175.  
  176. int main(){
  177.  
  178.     matrix_t A {
  179.         { 7.03, 0.92, 1.15, 1.135 },
  180.         { 0.92, 3.39, 1.3, 0.16 },
  181.         { 1.15, 1.3, 6.21, 2.1 },
  182.         { 1.135, 0.16, 2.1, 5.33 }
  183.     };
  184.  
  185.     frubenious(A);
  186.  
  187. //  print_vector(remans);
  188. //  print_vector(ans);
  189.  
  190.     print_self_vectors(A);
  191.  
  192. }
RAW Paste Data
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
 
Top