API tools faq
paste
Guest User

Untitled

a guest
Oct 13th, 2013
130
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 9.06 KB | None | 0 0
raw download clone embed print report
  1. // gauss.cpp
  2. // Doprowadza macierze do postaci rref + Horner gratis!
  3. // autor: Trivial, 2013-10-13
  4.  
  5.  
  6. #include <cassert>
  7. #include <cmath>
  8. #include <iostream>
  9. #include <cstdlib>
  10. #include <iomanip>
  11. #include <vector>
  12. #include <set>
  13.  
  14. using namespace std;
  15.  
  16.  
  17. int gcd(int a, int b) {
  18.     int r;
  19.     while (b != 0) {
  20.         r = a % b;
  21.         a = b;
  22.         b = r;
  23.     }
  24.     return a;
  25. }
  26.  
  27. class Rational {
  28.     int p, q;   /* p/q,   q > 0 */
  29. public:
  30.     Rational(int p, int q = 1) {
  31.         assert(q != 0);
  32.         if (q < 0) {
  33.             p = -p;
  34.             q = -q;
  35.         }
  36.         int g = gcd(abs(p), q);
  37.  
  38.         this->p = p/g;
  39.         this->q = q/g;
  40.     }
  41.     friend Rational operator+(const Rational &a, const Rational &b) {
  42.         return Rational(a.p*b.q + b.p*a.q,  a.q*b.q);
  43.     }
  44.     friend Rational operator-(const Rational &a, const Rational &b) {
  45.         return Rational(a.p*b.q - b.p*a.q,  a.q*b.q);
  46.     }
  47.     friend Rational operator*(const Rational &a, const Rational &b) {
  48.         return Rational(a.p*b.p, a.q*b.q);
  49.     }
  50.     friend Rational operator/(const Rational &a, const Rational &b) {
  51.         return Rational(a.p*b.q, a.q*b.p);
  52.     }
  53.     friend ostream & operator<<(ostream &out, const Rational &a) {
  54.         out << a.p;
  55.         if (a.q != 1) {
  56.             out << "/" << a.q;
  57.         }
  58.         return out;
  59.     }
  60.     friend bool operator==(const Rational &a, const Rational &b) {
  61.         return a.p == b.p && a.q == b.q;
  62.     }
  63.     friend bool operator!=(const Rational &a, const Rational &b) {
  64.         return a.p != b.p || a.q != b.q;
  65.     }
  66.  
  67.     friend bool operator<(const Rational &a, const Rational &b) {
  68.         return a.p*b.q < b.p*a.q;
  69.     }
  70.  
  71.     Rational operator-() const { return Rational(-p, q); }
  72.     int denom() const { return q; }
  73.     int numer() const { return p; }
  74. };
  75.  
  76. void printMatrix(Rational **A, int m, int n) {
  77.     for (int i = 0; i < m; i++) {
  78.         for (int j = 0; j < n; j++) {
  79.             cout << A[i][j] << "\t";
  80.         }
  81.         cout << "\n";
  82.     }
  83.     cout << endl;
  84. }
  85.  
  86. inline void swap(Rational &a, Rational &b) {
  87.     Rational c = a;
  88.     a = b;
  89.     b = c;
  90. }
  91.  
  92. void gauss(Rational **A, int m, int n) {
  93.     cout << "Eliminacja Gaussa dla macierzy:\n";
  94.     printMatrix(A,m,n);
  95.  
  96.     int min = (m < n) ? m : n;
  97.  
  98.     for (int k = 0; k < min; k++) {
  99.         if (A[k][k] == 0) {
  100.             for (int s = k+1; s < m; s++) {
  101.                 if (A[s][k] != 0) {
  102.                     cout << "w" << k+1 << "  <=>  w" << s+1 << "\n";
  103.                     for (int j = k; j < n; j++) {
  104.                         swap(A[k][j], A[s][j]);
  105.                     }
  106.                     break;
  107.                 }
  108.             }
  109.             if (A[k][k] == 0) {
  110.                 cout << "Macierz jest osobliwa." << endl;
  111.                 return;
  112.             }
  113.             printMatrix(A, m, n);
  114.         }
  115.  
  116.         bool hasChanged = false;
  117.         for (int i = k+1; i < m; i++) {
  118.             if (A[i][k] != 0) {
  119.                 Rational mult = A[i][k] / A[k][k];
  120.                 cout << "w" << i+1 << "  :=  w" << i+1
  121.                      << " - (" << mult << ")*w" << k+1 << "\n";
  122.                 for (int j = k; j < n; j++) {
  123.                     A[i][j] = A[i][j] - mult * A[k][j];
  124.                 }
  125.                 hasChanged = true;
  126.             }
  127.         }
  128.         if (hasChanged) {
  129.             cout << endl;
  130.             printMatrix(A, m, n);
  131.         }
  132.     }
  133.     cout << "Eliminacja zakonczona sukcesem!" << endl;
  134.     //~
  135.     //~ Rational det = 1;
  136.     //~ for (int k = 0; k < n; k++) {
  137.         //~ det = det * A[k][k];
  138.     //~ }
  139.     //~
  140.     //~ return sign*det;
  141. }
  142.  
  143.  
  144.  
  145. void rref_jordan(Rational **A, int m, int n) {
  146.     int *fs = (int*) malloc(n*sizeof(int));
  147.     int fi = 0;
  148.  
  149.     printMatrix(A, m, n);
  150.     for (int k = 0; k < n; k++) {
  151.         int r = k - fi;
  152.         if (r >= m) {
  153.             for (int c = k; c < n; c++) {
  154.                 cout << "Kolumna " << c+1 << " oznaczona jako wolna" << endl;
  155.             }
  156.             break;
  157.         }
  158.  
  159.         if (A[r][k] == 0) {
  160.             for (int s = r+1; s < m; s++) {
  161.                 if (A[s][k] != 0) {
  162.                     cout << "w" << r+1 << "  <=>  w" << s+1 << "\n";
  163.                     for (int j = k; j < n; j++) {
  164.                         swap(A[r][j], A[s][j]);
  165.                     }
  166.                     break;
  167.                 }
  168.             }
  169.             if (A[r][k] == 0) {
  170.                 fs[fi++] = k;
  171.                 cout << "Kolumna " << k+1 << " oznaczona jako wolna" << endl;
  172.                 continue;
  173.             }
  174.             printMatrix(A, m, n);
  175.         }
  176.  
  177.         if (A[r][k] != 1) {
  178.             cout << "w" << r+1 << " := " << "w" << r+1
  179.                  << "/(" << A[r][k] << ")\n" << endl;
  180.             for (int j = k+1; j < n; j++) {
  181.                 A[r][j] = A[r][j] / A[r][k];
  182.             }
  183.             for (int c = 0; c < fi; c++) {
  184.                 if (fs[c] < k) {
  185.                     A[r][fs[c]] = A[r][fs[c]] / A[r][k];
  186.                 }
  187.             }
  188.             A[r][k] = 1;
  189.             printMatrix(A, m, n);
  190.         }
  191.  
  192.         bool hasChanged = false;
  193.         for (int i = 0; i < m; i++) {
  194.             if (i == r) { continue; }
  195.             if (A[i][k] != 0) {
  196.                 Rational p = A[i][k];
  197.                 for (int j = k; j < n; j++) {
  198.                     A[i][j] = A[i][j] - p*A[r][j];
  199.                 }
  200.                 for (int c = 0; c < fi; c++) {
  201.                     int j = fs[c];
  202.                     if (j < k) {
  203.                         A[i][j] = A[i][j] - p*A[r][j];
  204.                     }
  205.                 }
  206.                 cout << "w" << i+1 << "  :=  w" << i+1
  207.                      << " - (" << p << ")*w" << r+1 << "\n";
  208.                 hasChanged = true;
  209.             }
  210.         }
  211.         if (hasChanged) {
  212.             cout << endl;
  213.             printMatrix(A, m, n);
  214.         }
  215.     }
  216.  
  217.     free(fs);
  218. }
  219.  
  220. void rref(Rational **A, int m, int n) {
  221.     int *fs = 1 + (int*) malloc((n+1)*sizeof(int));
  222.     int fi = 0;
  223.     fs[-1] = -1;
  224.  
  225.     printMatrix(A, m, n);
  226.     for (int k = 0; k < n; k++) {
  227.         int r = k - fi;
  228.         if (r >= m) {
  229.             for (int c = k; c < n; c++) {
  230.                 fs[fi++] = c;
  231.                 cout << "Kolumna " << c+1 << " oznaczona jako wolna" << endl;
  232.             }
  233.             break;
  234.         }
  235.  
  236.         if (A[r][k] == 0) {
  237.             for (int s = r+1; s < m; s++) {
  238.                 if (A[s][k] != 0) {
  239.                     cout << "w" << r+1 << "  <=>  w" << s+1 << "\n";
  240.                     for (int j = k; j < n; j++) {
  241.                         swap(A[r][j], A[s][j]);
  242.                     }
  243.                     break;
  244.                 }
  245.             }
  246.             if (A[r][k] == 0) {
  247.                 fs[fi++] = k;
  248.                 cout << "Kolumna " << k+1 << " oznaczona jako wolna" << endl;
  249.                 continue;
  250.             }
  251.             printMatrix(A, m, n);
  252.         }
  253.  
  254.         bool hasChanged = false;
  255.         for (int i = r+1; i < m; i++) {
  256.             if (A[i][k] != 0) {
  257.                 Rational p = A[i][k] / A[r][k];
  258.                 for (int j = k; j < n; j++) {
  259.                     A[i][j] = A[i][j] - p*A[r][j];
  260.                 }
  261.                 for (int c = 0; c < fi; c++) {
  262.                     int j = fs[c];
  263.                     if (j < k) {
  264.                         A[i][j] = A[i][j] - p*A[r][j];
  265.                     }
  266.                 }
  267.                 cout << "w" << i+1 << "  :=  w" << i+1
  268.                      << " - (" << p << ")*w" << r+1 << "\n";
  269.                 hasChanged = true;
  270.             }
  271.         }
  272.         if (hasChanged) {
  273.             cout << endl;
  274.             printMatrix(A, m, n);
  275.         }
  276.     }
  277.  
  278.     int fsize = fi;
  279.     for (int k = n-1; k >= 0; k--) {
  280.         if (k == fs[fi-1]) {
  281.             fi--;
  282.             continue;
  283.         }
  284.         int r = k - fi;
  285.  
  286.         if (A[r][k] != 1) {
  287.             cout << "w" << r+1 << " := " << "w" << r+1
  288.                  << "/(" << A[r][k] << ")\n" << endl;
  289.             for (int j = k+1; j < n; j++) {
  290.                 A[r][j] = A[r][j] / A[r][k];
  291.             }
  292.             A[r][k] = 1;
  293.             printMatrix(A, m, n);
  294.         }
  295.  
  296.         bool hasChanged = false;
  297.         for (int i = r-1; i >= 0; i--) {
  298.             if (A[i][k] != 0) {
  299.                 Rational p = A[i][k];
  300.                 A[i][k] = A[i][k] - p*A[r][k];
  301.                 for (int a = fi; a < fsize; a++) {
  302.                     int j = fs[a];
  303.                     A[i][j] = A[i][j] - p*A[r][j];
  304.                 }
  305.                 cout << "w" << i+1 << "  :=  w" << i+1
  306.                      << " - (" << p << ")*w" << r+1 << "\n";
  307.                 hasChanged = true;
  308.             }
  309.         }
  310.         if (hasChanged) {
  311.             cout << endl;
  312.             printMatrix(A, m, n);
  313.         }
  314.     }
  315.  
  316.  
  317.     free(fs-1);
  318. }
  319.  
  320. void horner(const Rational *as, int n, Rational x0) {
  321.     assert(n >= 0);
  322.  
  323.     Rational *bs = (Rational*) malloc((n+1)*sizeof(Rational));
  324.     Rational r = 0;
  325.  
  326.     for (int i = n; i >= 0; i--) {
  327.         cout << "\t" << as[i];
  328.     }
  329.     cout << endl;
  330.  
  331.     cout << x0;
  332.     for (int i = n; i >= 0; i--) {
  333.         cout << "\t" << r*x0;
  334.         r = r*x0 + as[i];
  335.         bs[i] = r;
  336.     }
  337.     cout << endl;
  338.  
  339.     for (int i = n; i >= 0; i--) {
  340.         cout << "\t" << bs[i];
  341.     }
  342.     cout << endl;
  343.  
  344.     free(bs);
  345. }
  346.  
  347. void divisors(vector<int> *xs, int n) {
  348.     int i;
  349.     for (i = 1; i < sqrt(n); i++) {
  350.         if (n%i == 0) {
  351.             xs->push_back(i);
  352.             xs->push_back(n/i);
  353.         }
  354.     }
  355.     if (i*i == n) {
  356.         xs->push_back(i);
  357.     }
  358. }
  359.  
  360. void printPolynomial(Rational *as, int n) {
  361.     cout << "\n\nWielomian: ";
  362.     for (int i = n; i >= 2; i--) {
  363.         if (as[i] != 0) {
  364.             cout << " ";
  365.             if (as[i] != 1) {
  366.                 cout << showpos << as[i];
  367.             }
  368.             cout << "x^" << noshowpos << i;
  369.         }
  370.     }
  371.     if (n >= 1 && as[1] != 0) {
  372.         cout << " " << showpos << as[1] << "x" << noshowpos;
  373.     }
  374.     if (as[0] != 0) {
  375.         cout << " " << showpos << as[0] << noshowpos;
  376.     }
  377.     cout << "\n" << endl;
  378. }
  379.  
  380. void tryRational(Rational *as, int n) {
  381.     assert(n >= 0);
  382.  
  383.     int m = 1;
  384.     for (int i = 0; i <= n; i++) {
  385.         m *= as[i].denom();
  386.     }
  387.     for (int i = 0; i <= n; i++) {
  388.         as[i] = as[i] * m;
  389.     }
  390.     int k = abs(as[0].numer());
  391.     for (int i = 1; i <= n && k != 1; i++) {
  392.         k = gcd(k, abs(as[i].numer()));
  393.     }
  394.     if (k != 1) {
  395.         for (int i = 0; i <= n; i++) {
  396.             as[i] = as[i].numer() / k;
  397.         }
  398.     }
  399.  
  400.     printPolynomial(as, n);
  401.  
  402.     vector<int> ns;
  403.     vector<int> ds;
  404.  
  405.     divisors(&ns, abs(as[0].numer()));
  406.     divisors(&ds, abs(as[n].numer()));
  407.  
  408.     set<Rational> cs;
  409.  
  410.     for (unsigned int j = 0; j < ds.size(); j++) {
  411.         for (unsigned int i = 0; i < ns.size(); i++) {
  412.             for (int s = 1; s > -2; s -= 2) {
  413.                 cs.insert(Rational(s*ns[i],ds[j]));
  414.             }
  415.         }
  416.     }
  417.     for (set<Rational>::iterator it = cs.begin(); it != cs.end(); it++) {
  418.         Rational x0 = *it;
  419.         cout << "dla x = " << x0 << ":" << "\n";
  420.         horner(as, n, x0);
  421.         cout << endl;
  422.     }
  423.  
  424. }
  425.  
  426.  
  427. int main() {
  428.     Rational A[4][4] = {
  429.         { 2, Rational(2,3), 3, 4,  },
  430.         { 1, 3, 1, Rational(3,2),  },
  431.         { 4, 2, 3, 5,   },
  432.         { -1, 5, 2, -3, },
  433.     };
  434.     Rational *As[4] = { A[0], A[1], A[2], A[3] };
  435.     rref(As, 4, 4); // tu można użyć różnych metod
  436.  
  437. //  Rational as[6] = { -2, 2, 1, 5, 3, 1 };
  438. //  tryRational(as, 5);
  439.  
  440.     return 0;
  441. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
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!