API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 28th, 2017
57
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 13.19 KB | None | 0 0
raw download clone embed print report
  1. #pragma GCC optimize("O3")
  2.  
  3. #define _SCL_SECURE_NO_WARNINGS
  4. #define _CRT_SECURE_NO_WARNINGS
  5.  
  6. #include <cassert>
  7. #include <ciso646>
  8. #include <climits>
  9. #include <cmath>
  10. #include <cstdint>
  11. #include <cstdio>
  12. #include <cstdlib>
  13. #include <cstring>
  14. #include <ctime>
  15. #include <array>
  16. #include <bitset>
  17. #include <deque>
  18. #include <forward_list>
  19. #include <list>
  20. #include <map>
  21. #include <queue>
  22. #include <set>
  23. #include <stack>
  24. #include <unordered_map>
  25. #include <unordered_set>
  26. #include <vector>
  27. #include <algorithm>
  28. #include <complex>
  29. #include <functional>
  30. #include <initializer_list>
  31. #include <iterator>
  32. #include <limits>
  33. #include <locale>
  34. #include <numeric>
  35. #include <regex>
  36. #include <string>
  37. #include <utility>
  38. #include <fstream>
  39. #include <iostream>
  40. #include <sstream>
  41. #include <iomanip>
  42. using namespace std;
  43.  
  44. typedef long long ll;
  45. typedef long double ld;
  46.  
  47. #define pb push_back
  48. #define mp make_pair
  49. #define x first
  50. #define y second
  51. #define speedup cout.tie(nullptr);cin.tie(nullptr);ios_base::sync_with_stdio(false)
  52. #define coutdouble cout<<setprecision(20);cout<<fixed
  53.  
  54. #define all(v) (v).begin(),(v).end()
  55. /*------CommentInInteractive--------*/
  56. #define endl "\n"
  57. /*----------------------------------*/
  58. const int INF = 2 * 1000 * 1000 * 1000 + 322;
  59. const ll LLINF = 2LL * 1000LL * 1000LL * 1000LL * 1000LL * 1000LL * 1000LL;
  60. const int MOD = 1000 * 1000 * 1000 + 7;
  61. const int RMOD = MOD - 2;
  62. const int N = 100 + 4;
  63. const int P1M = 1336337; //Large prime ( ~1M )
  64. const int P1K = 1009; //Big prime ( ~1K )
  65. const ld EPS = 1e-9;
  66. const ld PI = 3.1415926535897932384626433832795;
  67. /*------------------------------------------------IO_OPERATORS---------------------------------------------*/
  68. template<typename T, typename U> inline ostream &operator << (ostream &_out, const pair<T, U> &_p) { _out << _p.first << ' ' << _p.second; return _out; }
  69. template<typename T, typename U> inline istream &operator >> (istream &_in, pair<T, U> &_p) { _in >> _p.first >> _p.second; return _in; }
  70. template<typename T> inline ostream &operator << (ostream &_out, const vector<T> &_v) { if (_v.empty()) return _out; _out << _v.front(); for (auto _it = ++_v.begin(); _it != _v.end(); ++_it) _out << ' ' << *_it; return _out; }
  71. template<typename T> inline istream &operator >> (istream &_in, vector<T> &_v) { for (auto &_i : _v) _in >> _i; return _in; }
  72. template<typename T> inline ostream &operator << (ostream &_out, const set<T> &_s) { if (_s.empty()) return _out; _out << *_s.begin(); for (auto _it = ++_s.begin(); _it != _s.end(); ++_it) _out << ' ' << *_it; return _out; }
  73. template<typename T> inline ostream &operator << (ostream &_out, const multiset<T> &_s) { if (_s.empty()) return _out; _out << *_s.begin(); for (auto _it = ++_s.begin(); _it != _s.end(); ++_it) _out << ' ' << *_it; return _out; }
  74. template<typename T> inline ostream &operator << (ostream &_out, const unordered_set<T> &_s) { if (_s.empty()) return _out; _out << *_s.begin(); for (auto _it = ++_s.begin(); _it != _s.end(); ++_it) _out << ' ' << *_it; return _out; }
  75. template<typename T> inline ostream &operator << (ostream &_out, const unordered_multiset<T> &_s) { if (_s.empty()) return _out; _out << *_s.begin(); for (auto _it = ++_s.begin(); _it != _s.end(); ++_it) _out << ' ' << *_it; return _out; }
  76. template<typename T, typename U> inline ostream &operator << (ostream &_out, const map<T, U> &_m) { if (_m.empty()) return _out; _out << '(' << _m.begin()->first << ": " << _m.begin()->second << ')'; for (auto _it = ++_m.begin(); _it != _m.end(); ++_it) _out << ", (" << _it->first << ": " << _it->second << ')'; return _out; }
  77. template<typename T, typename U> inline ostream &operator << (ostream &_out, const unordered_map<T, U> &_m) { if (_m.empty()) return _out; _out << '(' << _m.begin()->first << ": " << _m.begin()->second << ')'; for (auto _it = ++_m.begin(); _it != _m.end(); ++_it) _out << ", (" << _it->first << ": " << _it->second << ')'; return _out; }
  78. /*--------------------------------------------------IO_FILES-----------------------------------------------*/
  79. const char * infile =
  80. #ifdef LOCAL
  81. "input.txt"
  82. #else
  83. ""
  84. #endif //LOCAL
  85. ;
  86. const char * outfile =
  87. #ifdef LOCAL
  88. ""
  89. #else
  90. ""
  91. #endif //LOCAL
  92. ;
  93. /*----------------------------------------------------MATHS------------------------------------------------*/
  94. inline int gcd(int a, int b) { while (b) { int tmp = a%b; a = b; b = tmp; } return a; }
  95. inline int gcd(vector<int> a) { while (a.size() > 1) { int tmp1 = a.back(), tmp2 = a.back(); a.pop_back(), a.pop_back(); a.push_back(gcd(tmp1, tmp2)); } return a[0]; }
  96. inline ll lcm(int a, int b) { return (ll)a*(ll)b / gcd(a, b); }
  97. inline ll pwm(ll xx, ll pow) { ll mlt = 1; while (pow) { if (pow & 1) { mlt *= xx; xx *= xx; pow >>= 1; xx %= MOD; mlt %= MOD; } else { pow >>= 1; xx *= xx; xx %= MOD; } }return mlt; }
  98. inline ll pw(ll xx, int pow) { ll mlt = 1; while (pow) { if (pow & 1) { mlt *= xx; xx *= xx; pow >>= 1; } else { pow >>= 1; xx *= xx; } }return mlt; }
  99. inline ll rev(ll r) { return pwm(r, RMOD); }
  100. /*--------------------------------------------------------------------------------------------------------*/
  101.  
  102. #define MAX_PRINT 20
  103. #define LEN 5000
  104. #define eps 1e-17
  105.  
  106. union bit_double {
  107.     double val;
  108.     ll bits;
  109. };
  110.  
  111. inline void print_mat(double *mat, int n, int m) { // n - #columns, m - #rows
  112.     int prt_n = (n > MAX_PRINT) ? MAX_PRINT : n;
  113.     int prt_m = (m > MAX_PRINT) ? MAX_PRINT : m;
  114.     int i, j;
  115.     for (i = 0; i < prt_m; ++i, mat += n) {
  116.         for (j = 0; j < prt_n; ++j)
  117.             cout << setw(10) << mat[j] << " ";
  118.         cout << endl;
  119.     }
  120. }
  121.  
  122. inline int input_mat(double *mat, int n, int m, const char *file) {
  123.     double *end = mat + m * n;
  124.     for (; mat < end && cin >> *mat; ++mat);
  125.     return 0;
  126. }
  127.  
  128.  
  129. inline double gen_func(int i, int j) {
  130.     double res = 1. / (i + 1 + j);
  131.     //double res = i + j + 1;
  132.     return res;
  133. }
  134.  
  135. inline void gen(double *mat, int n, int m) { // n - #columns, m - #rows
  136.     int i, j;
  137.     memset(mat, 0, n * n * sizeof(double));
  138.     for (i = 0; i < m; ++i) {
  139.         for (j = 0; j < n; ++j) {
  140.             mat[i *n + j] = gen_func(i, j);
  141.         }
  142.     }
  143. }
  144.  
  145. inline void transp(double *A, int n) {
  146.     for (int i = 0; i < n; ++i) {
  147.         for (int j = i + 1; j < n; ++j) {
  148.             double tmp = A[i *n + j];
  149.             A[i * n + j] = A[j * n + i];
  150.             A[j * n + i] = tmp;
  151.         }
  152.     }
  153. }
  154.  
  155. inline void residial_cnt(double *lm, double *rm, int n, double * rs_1, double * rs_inf) {
  156.     static double norm[LEN];
  157.     const double *end = rm + n * n, *begin = rm;
  158.     double *cur = lm;
  159.     double cur_elem = -1, cur_s = 0, max_s = -1;
  160.     int i, j, k;
  161.     memset(norm, 0, n * sizeof(double));
  162.     switch (i = n & 3) {
  163.     case 3:
  164.         for (j = 0, k = 0; rm < end; ++rm, ++j) {
  165.             cur_elem += cur[j] * (*rm);
  166.             if (j == n - 1) {
  167.                 norm[k] += cur_elem;
  168.                 j = -1, ++k;
  169.                 cur_s += fabs(cur_elem), cur_elem = 0;
  170.             }
  171.         }
  172.         cur += n, rm = (double *)begin;
  173.         max_s = cur_s > max_s ? cur_s : max_s;
  174.         cur_s = 0;
  175.     case 2:
  176.         for (j = 0, k = 0; rm < end; ++rm, ++j) {
  177.             cur_elem += cur[j] * (*rm);
  178.             if (j == n - 1) {
  179.                 norm[k] += cur_elem;
  180.                 j = -1, ++k;
  181.                 cur_s += fabs(cur_elem), cur_elem = (i - 2) == k ? -1 : 0;
  182.             }
  183.         }
  184.         cur += n, rm = (double *)begin;
  185.         max_s = cur_s > max_s ? cur_s : max_s;
  186.         cur_s = 0;
  187.     case 1:
  188.         for (j = 0, k = 0; rm < end; ++rm, ++j) {
  189.             cur_elem += cur[j] * (*rm);
  190.             if (j == n - 1) {
  191.                 norm[k] += cur_elem;
  192.                 j = -1, ++k;
  193.                 cur_s += fabs(cur_elem), cur_elem = (i - 1) == k ? -1 : 0;
  194.             }
  195.         }
  196.         cur += n, rm = (double *)begin;
  197.         max_s = cur_s > max_s ? cur_s : max_s;
  198.         cur_s = 0;
  199.     }
  200.     for (; i < n;) {
  201.         // 0
  202.         for (j = 0, k = 0; rm < end; ++rm, ++j) {
  203.             cur_elem += cur[j] * (*rm);
  204.             if (j == n - 1) {
  205.                 norm[k] += cur_elem;
  206.                 j = -1, ++k;
  207.                 cur_s += fabs(cur_elem), cur_elem = i == k ? -1 : 0;
  208.             }
  209.         }
  210.         cur += n, rm = (double *)begin;
  211.         max_s = cur_s > max_s ? cur_s : max_s;
  212.         cur_s = 0; ++i;
  213.         //1
  214.         for (j = 0, k = 0; rm < end; ++rm, ++j) {
  215.             cur_elem += cur[j] * (*rm);
  216.             if (j == n - 1) {
  217.                 norm[k] += cur_elem;
  218.                 j = -1, ++k;
  219.                 cur_s += fabs(cur_elem), cur_elem = i == k ? -1 : 0;
  220.             }
  221.         }
  222.         cur += n, rm = (double *)begin;
  223.         max_s = cur_s > max_s ? cur_s : max_s;
  224.         cur_s = 0; ++i;
  225.         //2
  226.         for (j = 0, k = 0; rm < end; ++rm, ++j) {
  227.             cur_elem += cur[j] * (*rm);
  228.             if (j == n - 1) {
  229.                 norm[k] += cur_elem;
  230.                 j = -1, ++k;
  231.                 cur_s += fabs(cur_elem), cur_elem = i == k ? -1 : 0;
  232.             }
  233.         }
  234.         cur += n, rm = (double *)begin;
  235.         max_s = cur_s > max_s ? cur_s : max_s;
  236.         cur_s = 0; ++i;
  237.         //3
  238.         for (j = 0, k = 0; rm < end; ++rm, ++j) {
  239.             cur_elem += cur[j] * (*rm);
  240.             if (j == n - 1) {
  241.                 norm[k] += cur_elem;
  242.                 j = -1, ++k;
  243.                 cur_s += fabs(cur_elem), cur_elem = i == k ? -1 : 0;
  244.             }
  245.         }
  246.         cur += n, rm = (double *)begin;
  247.         max_s = cur_s > max_s ? cur_s : max_s;
  248.         cur_s = 0; ++i;
  249.     }
  250.     *rs_1 = max_s;
  251.     for (i = 1; i < n; *rs_inf = *rs_inf > norm[i] ? *rs_inf : norm[i], ++i);
  252. }
  253.  
  254. inline void mul_1(double *A, int i, int j, double alfa, double beta, int n) {
  255.     int k;
  256.     double *r1 = A + i*n, *r2 = r1 + j*n;
  257.     double elem_r1, elem_r2;
  258.     for (k = i; k < n; ++k) {
  259.         elem_r1 = r1[k], elem_r2 = r2[k];
  260.         r1[k] = r1[k] * alfa - elem_r2 * beta;
  261.         r2[k] = r2[k] * alfa + elem_r1 * beta;
  262.     }
  263. }
  264.  
  265.  
  266. inline void mul_2(double *A, int i, int j, double alfa, double beta, int n) {
  267.     int k;
  268.     double *r1 = A + i*n, *r2 = A + j*n;
  269.     double elem_r1, elem_r2;
  270.     for (k = 0; k < n; ++k) {
  271.         elem_r1 = r1[k], elem_r2 = r2[k];
  272.         r1[k] = r1[k] * alfa + elem_r2 * beta;
  273.         r2[k] = r2[k] * alfa - elem_r1 * beta;
  274.     }
  275. }
  276.  
  277. inline bool inv_U(double *U, int n) {
  278.     static double cur_row[LEN];
  279.     double *cur, *sub_r;
  280.     double factor;
  281.     int i, j, k;
  282.     cur = U + n*(n - 1);
  283.     for (i = n - 1; i >= 0; --i, cur -= n) {
  284.         memset(cur_row + i, 0, (n - i) * sizeof(double));
  285.         cur_row[i] = 1;
  286.         sub_r = U + i*n;
  287.         for (j = i + 1; j < n; ++j) {
  288.             double coef = cur[j];
  289.             sub_r += n;
  290.             for (k = j; k < n; ++k) {
  291.                 cur_row[k] -= coef * sub_r[k];
  292.             }
  293.         }
  294.         if (fabs(cur[i]) < eps) {
  295.             return false;
  296.         }
  297.         for (k = i, factor = cur[i]; k < n; ++k) {
  298.             cur_row[k] /= factor;
  299.         }
  300.         memcpy(cur + i, cur_row + i, (n - i)* sizeof(double));
  301.     }
  302.     return true;
  303. } // U = U^(-1)
  304.  
  305.  
  306. inline void QRD_upg(double *A, int n) {
  307.     int i, j;
  308.     double cos_fi, sin_fi;
  309.     double c_v, c_v2, factor;
  310.     double *cur = A;
  311.     ll *cur_ll = (ll *)A;
  312.     for (i = 0; i < n - 1; ++i, cur += n + 1, cur_ll += n + 1) {
  313.         for (j = 1; j < n - i; ++j) {
  314.             c_v = *cur, c_v2 = cur[j*n];
  315.             if (fabs(c_v2) < eps) {
  316.                 cur_ll[j*n] |= 3;
  317.                 continue;
  318.             }
  319.             factor = (c_v * c_v + c_v2 * c_v2);
  320.             sin_fi = -(cos_fi = sqrt(1. / factor));
  321.             sin_fi *= c_v2, cos_fi *= c_v;
  322.             mul_1(A, i, j, cos_fi, sin_fi, n);// A -> cur ?
  323.             if (fabs(cos_fi) < fabs(sin_fi)) { // ...ab, a = (is this sin_fi ? 1 : 0), b = (other func > 0 ? 1 : 0)
  324.                 cur[j*n] = cos_fi;
  325.                 if (sin_fi > 0) {//..01
  326.                     cur_ll[j*n] |= 1;
  327.                     cur_ll[j*n] &= ~2;
  328.                 }
  329.                 else {//...00
  330.                     cur_ll[j*n] &= ~3;
  331.                 }
  332.             }
  333.             else {
  334.                 cur[j*n] = sin_fi;
  335.                 if (cos_fi > 0) {//...11
  336.                     cur_ll[j*n] |= 3;
  337.                 }
  338.                 else {//...10
  339.                     cur_ll[j*n] |= 2;
  340.                     cur_ll[j*n] &= ~1;
  341.                 }
  342.             }
  343.         }
  344.     }
  345. }
  346.  
  347. inline void fill(double *A, double *B, int n) {
  348.     int i, j;
  349.     for (i = 0; i < n; ++i, ++B, A += n) {
  350.         for (j = i; j < n; ++j) {
  351.             B[j*n] = A[j];
  352.         }
  353.     }
  354. }
  355.  
  356. inline void mult_new(double *A, double *B, int n) {
  357.     int i, j;
  358.     bit_double cur;
  359.     double *ptr = B + n - 1;
  360.     double sin_fi, cos_fi;
  361.     for (i = n - 1; i >= 0; --i, --ptr) {
  362.         for (j = n - 1; j > i; --j) {
  363.             cur.val = ptr[j*n];
  364.             switch (cur.bits & 3) {
  365.             case 0:
  366.                 cos_fi = cur.val;
  367.                 sin_fi = -sqrt(1. - cos_fi * cos_fi);
  368.                 break;
  369.             case 1:
  370.                 cos_fi = cur.val;
  371.                 sin_fi = sqrt(1. - cos_fi * cos_fi);
  372.                 break;
  373.             case 2:
  374.                 sin_fi = cur.val;
  375.                 cos_fi = -sqrt(1. - sin_fi * sin_fi);
  376.                 break;
  377.             case 3:
  378.                 sin_fi = cur.val;
  379.                 cos_fi = sqrt(1. - sin_fi * sin_fi);
  380.                 break;
  381.             }
  382.             if (fabs(sin_fi) < eps)
  383.                 continue;
  384.             mul_2(A, i, j, cos_fi, sin_fi, n);
  385.         }
  386.     }
  387. }
  388.  
  389. inline void solve(ld tt) {
  390.     //current stats, 4K x 4K: QRD - 55 sec, inv_U - 8 sec, mult_mat - 87 sec, cnt_residial - 65 sec, overall - 215 sec
  391.     static double A[LEN * LEN], B[LEN *LEN];
  392.     int n;
  393.     double rs_1, rs_inf;
  394.     cout << setprecision(4);// cout << fixed;
  395.     cin >> n;
  396.     gen(A, n, n);
  397.     //input_mat(A, n, n, "");
  398.     cout << "/*-----A------*/\n";
  399.     print_mat(A, n, n);
  400.  
  401.     tt = clock();
  402.  
  403.     QRD_upg(A, n);
  404.    
  405.     cout << "Time: " << ((ld)clock() - tt) / CLOCKS_PER_SEC << endl;
  406.  
  407.     if (!inv_U(A, n)) {
  408.         cout << "Matrix is non - invertible!\n";
  409.         return;
  410.     }
  411.  
  412.     cout << "Time: " << ((ld)clock() - tt) / CLOCKS_PER_SEC << endl;
  413.  
  414.     fill(A, B, n);
  415.    
  416.     mult_new(B, A, n);
  417.     transp(B, n);
  418.  
  419.     cout << "/*-----A^(-1)------*/\n";
  420.     print_mat(B, n, n);
  421.  
  422.     cout << "Time: " << ((ld)clock() - tt) / CLOCKS_PER_SEC << endl;
  423.  
  424.     gen(A, n, n);
  425.     //input_mat(A, n, n, "");
  426.  
  427.     transp(A, n);
  428.     residial_cnt(A, B, n, &rs_1, &rs_inf);
  429.    
  430.     cout << "residial = " << rs_1 << " for ||.||_1, " << rs_inf << " for ||.||_inf" << endl;
  431. }
  432.  
  433. int main() {
  434.     ld tt = clock();
  435.  
  436.     if (*infile != '\0')
  437.         freopen(infile, "r", stdin);
  438.     if (*outfile != '\0')
  439.         freopen(outfile, "w", stdout);
  440.  
  441.     speedup;
  442.     //coutdouble;
  443.  
  444.     //while(true)
  445.     solve(tt);
  446.  
  447. #ifdef LOCAL
  448.     cout << "Time: " << ((ld)clock() - tt) / CLOCKS_PER_SEC << endl;
  449.     while (true);
  450. #endif // LOCAL
  451.     return 0;
  452. }
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!