Advertisement
MagicWinnie

Solution

Mar 17th, 2023 (edited)
566
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 21.39 KB | None | 0 0
  1. #include <iostream>
  2. #include <iomanip>
  3. #include <cmath>
  4. #define USE_GNUPLOT 1
  5. #define GNUPLOT_NAME "gnuplot -persist"
  6.  
  7. using namespace std;
  8.  
  9. // Epsilon for comparing double values
  10. const double EPS = 1E-10;
  11.  
  12. using namespace std;
  13.  
  14. // The basic class for a general matrix containing double values.
  15. // Not using templates as there is no need for it.
  16. // In maths matrices are used with floating point values.
  17. class Matrix
  18. {
  19. protected:
  20.     // Size of matrix.
  21.     size_t rows = 0, cols = 0;
  22.     // 2D dynamic array for storing the matrix.
  23.     double **matrix;
  24.  
  25. public:
  26.     // Constructor.
  27.     Matrix() {}
  28.     Matrix(size_t rows, size_t cols)
  29.     {
  30.         this->rows = rows;
  31.         this->cols = cols;
  32.  
  33.         matrix = new double *[rows];
  34.         for (size_t i = 0; i < rows; i++)
  35.             matrix[i] = new double[cols];
  36.  
  37.         // Contains zeros as default values.
  38.         for (size_t i = 0; i < rows; i++)
  39.             for (size_t j = 0; j < cols; j++)
  40.                 matrix[i][j] = 0.0;
  41.     }
  42.     // Getter for number of rows.
  43.     size_t get_rows() const { return rows; }
  44.     // Getter for number of columns.
  45.     size_t get_cols() const { return cols; }
  46.     // Cin operator overload.
  47.     friend istream &operator>>(istream &in, Matrix &m)
  48.     {
  49.         for (size_t i = 0; i < m.rows; i++)
  50.             for (size_t j = 0; j < m.cols; j++)
  51.                 in >> m[i][j];
  52.         return in;
  53.     }
  54.     // Cout operator overload.
  55.     friend ostream &operator<<(ostream &out, Matrix &m)
  56.     {
  57.         for (size_t i = 0; i < m.rows; i++)
  58.         {
  59.             for (size_t j = 0; j < m.cols; j++)
  60.             {
  61.                 // Correct output for small negative values.
  62.                 out << ((abs(m[i][j]) < EPS and signbit(m[i][j])) ? 0.0 : m[i][j]);
  63.                 if (j != m.cols - 1)
  64.                     cout << ' ';
  65.             }
  66.             out << '\n';
  67.         }
  68.         return out;
  69.     }
  70.     // Method for accessing the values by [].
  71.     double *operator[](size_t i) { return matrix[i]; }
  72.     const double *operator[](size_t i) const { return matrix[i]; }
  73.     // = operator overloading.
  74.     Matrix &operator=(const Matrix &old)
  75.     {
  76.         if (this == &old)
  77.             return *this;
  78.  
  79.         if (rows == 0 || cols == 0)
  80.         {
  81.             rows = old.rows;
  82.             cols = old.cols;
  83.  
  84.             matrix = new double *[rows];
  85.             for (size_t i = 0; i < rows; i++)
  86.                 matrix[i] = new double[cols];
  87.         }
  88.  
  89.         if (rows != old.rows or cols != old.cols)
  90.             throw runtime_error("[ERROR] [operator=]: the dimensional problem occurred");
  91.  
  92.         for (size_t i = 0; i < rows; i++)
  93.             for (size_t j = 0; j < cols; j++)
  94.                 matrix[i][j] = old.matrix[i][j];
  95.  
  96.         return *this;
  97.     }
  98.     // + operator overloading.
  99.     Matrix operator+(const Matrix &B) const
  100.     {
  101.         if (rows != B.rows or cols != B.cols)
  102.             throw runtime_error("[ERROR] [operator+]: the dimensional problem occurred");
  103.         Matrix C(rows, cols);
  104.         for (size_t i = 0; i < rows; i++)
  105.             for (size_t j = 0; j < cols; j++)
  106.             {
  107.                 C.matrix[i][j] = matrix[i][j] + B[i][j];
  108.                 if (abs(C.matrix[i][j]) < EPS)
  109.                     C.matrix[i][j] = 0;
  110.             }
  111.         return C;
  112.     }
  113.     // - operator overloading.
  114.     Matrix operator-(const Matrix &B) const
  115.     {
  116.         if (rows != B.rows or cols != B.cols)
  117.             throw runtime_error("[ERROR] [operator-]: the dimensional problem occurred");
  118.         Matrix C(rows, cols);
  119.         for (size_t i = 0; i < rows; i++)
  120.             for (size_t j = 0; j < cols; j++)
  121.             {
  122.                 C.matrix[i][j] = matrix[i][j] - B[i][j];
  123.                 if (abs(C.matrix[i][j]) < EPS)
  124.                     C.matrix[i][j] = 0;
  125.             }
  126.         return C;
  127.     }
  128.     // * operator overloading.
  129.     Matrix operator*(const Matrix &B) const
  130.     {
  131.         if (cols != B.rows)
  132.             throw runtime_error("[ERROR] [operator*]: the dimensional problem occurred");
  133.         Matrix C(rows, B.cols);
  134.         for (size_t i = 0; i < rows; i++)
  135.             for (size_t j = 0; j < B.cols; j++)
  136.                 for (size_t k = 0; k < cols; k++)
  137.                 {
  138.                     C.matrix[i][j] += matrix[i][k] * B[k][j];
  139.                     if (abs(C.matrix[i][j]) < EPS)
  140.                         C.matrix[i][j] = 0;
  141.                 }
  142.         return C;
  143.     }
  144.     // Method to transpose the matrix.
  145.     Matrix T() const
  146.     {
  147.         Matrix B(cols, rows);
  148.         for (size_t i = 0; i < rows; i++)
  149.             for (size_t j = 0; j < cols; j++)
  150.                 B.matrix[j][i] = matrix[i][j];
  151.         return B;
  152.     }
  153. };
  154.  
  155. // Class for a square matrix that inherites the general matrix class.
  156. class SquareMatrix : public Matrix
  157. {
  158. protected:
  159.     // Method to get the inverse of the matrix.
  160.     friend SquareMatrix inverse_worker(SquareMatrix A, bool debug_info);
  161.     // Method to get the determinant of the matrix.
  162.     friend double determinant_worker(SquareMatrix A, bool debug_info);
  163.  
  164. public:
  165.     // Constructor.
  166.     SquareMatrix() {}
  167.     SquareMatrix(size_t size) : Matrix(size, size) {}
  168.     SquareMatrix(Matrix matrix) : SquareMatrix(matrix.get_cols())
  169.     {
  170.         if (matrix.get_cols() != matrix.get_rows())
  171.             throw runtime_error("[ERROR] [SquareMatrix(Matrix)]: the dimensional problem occurred");
  172.         for (int i = 0; i < matrix.get_rows(); i++)
  173.             for (int j = 0; j < matrix.get_cols(); j++)
  174.                 this->matrix[i][j] = matrix[i][j];
  175.     }
  176.     using Matrix::operator=;
  177.     // + operator overloading using downcasting and upcasting.
  178.     SquareMatrix operator+(const SquareMatrix &B) const
  179.     {
  180.         const Matrix *upcast_A = this;
  181.         const Matrix *upcast_B = &B;
  182.  
  183.         const Matrix upcast_C = (*upcast_A) + (*upcast_B);
  184.         const Matrix *upcast_C_ptr = &upcast_C;
  185.  
  186.         const SquareMatrix *C_ptr = (SquareMatrix *)(upcast_C_ptr);
  187.  
  188.         SquareMatrix C(cols);
  189.         C = *C_ptr;
  190.  
  191.         return C;
  192.     }
  193.     // - operator overloading using downcasting and upcasting.
  194.     SquareMatrix operator-(const SquareMatrix &B) const
  195.     {
  196.         const Matrix *upcast_A = this;
  197.         const Matrix *upcast_B = &B;
  198.  
  199.         const Matrix upcast_C = (*upcast_A) - (*upcast_B);
  200.         const Matrix *upcast_C_ptr = &upcast_C;
  201.  
  202.         const SquareMatrix *C_ptr = (SquareMatrix *)(upcast_C_ptr);
  203.  
  204.         SquareMatrix C(cols);
  205.         C = *C_ptr;
  206.  
  207.         return C;
  208.     }
  209.     // * operator overloading using downcasting and upcasting.
  210.     SquareMatrix operator*(const SquareMatrix &B) const
  211.     {
  212.         const Matrix *upcast_A = this;
  213.         const Matrix *upcast_B = &B;
  214.  
  215.         const Matrix upcast_C = (*upcast_A) * (*upcast_B);
  216.         const Matrix *upcast_C_ptr = &upcast_C;
  217.  
  218.         const SquareMatrix *C_ptr = (SquareMatrix *)(upcast_C_ptr);
  219.  
  220.         SquareMatrix C(cols);
  221.         C = *C_ptr;
  222.         return C;
  223.     }
  224.     // Transpose method overloading using downcasting and upcasting.
  225.     SquareMatrix T() const
  226.     {
  227.         const Matrix *upcast_A = this;
  228.  
  229.         const Matrix upcast_C = upcast_A->T();
  230.         const Matrix *upcast_C_ptr = &upcast_C;
  231.  
  232.         const SquareMatrix *C_ptr = (SquareMatrix *)(upcast_C_ptr);
  233.  
  234.         SquareMatrix C(cols);
  235.         C = *C_ptr;
  236.  
  237.         return C;
  238.     }
  239.     // Method to get the determinant of the matrix.
  240.     double determinant(bool debug_info) { return determinant_worker(*this, debug_info); }
  241.     // Method to get the inverse of the matrix.
  242.     SquareMatrix inverse(bool debug_info) { return inverse_worker(*this, debug_info); }
  243. };
  244.  
  245. // Class for an identity matrix that inherites the square matrix class.
  246. class IdentityMatrix : public SquareMatrix
  247. {
  248. public:
  249.     // Constructor.
  250.     IdentityMatrix(size_t size) : SquareMatrix(size)
  251.     {
  252.         // Main diagonal has only 1s.
  253.         for (size_t i = 0; i < size; i++)
  254.             matrix[i][i] = 1.0;
  255.     }
  256.     IdentityMatrix(Matrix matrix) : IdentityMatrix(matrix.get_cols())
  257.     {
  258.         if (matrix.get_cols() != matrix.get_rows())
  259.             throw runtime_error("[ERROR] [IdentityMatrix(Matrix)]: the dimensional problem occurred");
  260.         for (int i = 0; i < matrix.get_rows(); i++)
  261.             for (int j = 0; j < matrix.get_cols(); j++)
  262.                 this->matrix[i][j] = matrix[i][j];
  263.     }
  264. };
  265.  
  266. // Class for an elimination matrix that inherites the identity matrix class.
  267. class EliminationMatrix : public IdentityMatrix
  268. {
  269. public:
  270.     EliminationMatrix(size_t size, size_t i, size_t j, double val) : IdentityMatrix(size)
  271.     {
  272.         // Elimination matrix is an identity matrix, but with a special value at ij.
  273.         matrix[i][j] = val * -1.0;
  274.     }
  275.     EliminationMatrix(Matrix matrix) : IdentityMatrix(matrix.get_cols())
  276.     {
  277.         if (matrix.get_cols() != matrix.get_rows())
  278.             throw runtime_error("[ERROR] [EliminationMatrix(Matrix)]: the dimensional problem occurred");
  279.         for (int i = 0; i < matrix.get_rows(); i++)
  280.             for (int j = 0; j < matrix.get_cols(); j++)
  281.                 this->matrix[i][j] = matrix[i][j];
  282.     }
  283. };
  284.  
  285. // Class for a permutation matrix that inherites the identity matrix class.
  286. class PermutationMatrix : public IdentityMatrix
  287. {
  288. public:
  289.     PermutationMatrix(size_t size, size_t i, size_t j) : IdentityMatrix(size)
  290.     {
  291.         // Permutation matrix is a matrix that has rows i and j exchanged.
  292.         swap(matrix[i], matrix[j]);
  293.     }
  294.     PermutationMatrix(Matrix matrix) : IdentityMatrix(matrix.get_cols())
  295.     {
  296.         if (matrix.get_cols() != matrix.get_rows())
  297.             throw runtime_error("[ERROR] [PermutationMatrix(Matrix)]: the dimensional problem occurred");
  298.         for (int i = 0; i < matrix.get_rows(); i++)
  299.             for (int j = 0; j < matrix.get_cols(); j++)
  300.                 this->matrix[i][j] = matrix[i][j];
  301.     }
  302. };
  303.  
  304. // Class for a column vector that inherites the general matrix class.
  305. class ColumnVector : public Matrix
  306. {
  307. public:
  308.     // Constructor.
  309.     ColumnVector() {}
  310.     ColumnVector(size_t size) : Matrix(size, 1) {}
  311.     ColumnVector(Matrix matrix) : ColumnVector(matrix.get_rows())
  312.     {
  313.         if (matrix.get_cols() != 1)
  314.             throw runtime_error("[ERROR] [ColumnVector(Matrix)]: the dimensional problem occurred");
  315.         for (int i = 0; i < matrix.get_rows(); i++)
  316.             this->matrix[i][0] = matrix[i][0];
  317.     }
  318.     using Matrix::operator=;
  319.     // + operator overloading using downcasting and upcasting.
  320.     ColumnVector operator+(const ColumnVector &B) const
  321.     {
  322.         const Matrix *upcast_A = this;
  323.         const Matrix *upcast_B = &B;
  324.  
  325.         const Matrix upcast_C = (*upcast_A) + (*upcast_B);
  326.         const Matrix *upcast_C_ptr = &upcast_C;
  327.  
  328.         const ColumnVector *C_ptr = (ColumnVector *)(upcast_C_ptr);
  329.  
  330.         ColumnVector C(cols);
  331.         C = *C_ptr;
  332.  
  333.         return C;
  334.     }
  335.     // - operator overloading using downcasting and upcasting.
  336.     ColumnVector operator-(const ColumnVector &B) const
  337.     {
  338.         const Matrix *upcast_A = this;
  339.         const Matrix *upcast_B = &B;
  340.  
  341.         const Matrix upcast_C = (*upcast_A) - (*upcast_B);
  342.         const Matrix *upcast_C_ptr = &upcast_C;
  343.  
  344.         const ColumnVector *C_ptr = (ColumnVector *)(upcast_C_ptr);
  345.  
  346.         ColumnVector C(cols);
  347.         C = *C_ptr;
  348.  
  349.         return C;
  350.     }
  351.     // * operator overloading using downcasting and upcasting.
  352.     ColumnVector operator*(const ColumnVector &B) const
  353.     {
  354.         const Matrix *upcast_A = this;
  355.         const Matrix *upcast_B = &B;
  356.  
  357.         const Matrix upcast_C = (*upcast_A) * (*upcast_B);
  358.         const Matrix *upcast_C_ptr = &upcast_C;
  359.  
  360.         const ColumnVector *C_ptr = (ColumnVector *)(upcast_C_ptr);
  361.  
  362.         ColumnVector C(cols);
  363.         C = *C_ptr;
  364.  
  365.         return C;
  366.     }
  367.     // Transpose method overloading using downcasting and upcasting.
  368.     ColumnVector T() const
  369.     {
  370.         const Matrix *upcast_A = this;
  371.  
  372.         const Matrix upcast_C = upcast_A->T();
  373.         const Matrix *upcast_C_ptr = &upcast_C;
  374.  
  375.         const ColumnVector *C_ptr = (ColumnVector *)(upcast_C_ptr);
  376.  
  377.         ColumnVector C(cols);
  378.         C = *C_ptr;
  379.  
  380.         return C;
  381.     }
  382.     // Method to get norm of vector.
  383.     double norm()
  384.     {
  385.         double result = 0.0;
  386.         for (size_t i = 0; i < rows; i++)
  387.             result += matrix[i][0] * matrix[i][0];
  388.         return sqrt(result);
  389.     }
  390. };
  391.  
  392. // Class for storing an augmented matrix with different types of left and right matrices.
  393. template <typename L, typename R>
  394. class AugmentedMatrix
  395. {
  396. protected:
  397.     L matrixLeft;
  398.     R matrixRight;
  399.  
  400. public:
  401.     // Constructor.
  402.     AugmentedMatrix(const AugmentedMatrix<L, R> &old)
  403.     {
  404.         matrixLeft = old.matrixLeft;
  405.         matrixRight = old.matrixRight;
  406.     }
  407.     AugmentedMatrix(L A, R B)
  408.     {
  409.         if (A.get_rows() != B.get_rows())
  410.             throw runtime_error("[ERROR] [AugmentedMatrix(L, R)]: the dimensional problem occurred");
  411.         matrixLeft = A;
  412.         matrixRight = B;
  413.     }
  414.     // Getter for left matrix.
  415.     L &getLeft() { return matrixLeft; }
  416.     // Getter for right matrix.
  417.     R &getRight() { return matrixRight; }
  418.     AugmentedMatrix<L, R> &operator=(const AugmentedMatrix<L, R> &old)
  419.     {
  420.         if (this == &old)
  421.             return *this;
  422.  
  423.         matrixLeft = old.matrixLeft;
  424.         matrixRight = old.matrixRight;
  425.  
  426.         return *this;
  427.     }
  428.     // Cout operator overloading.
  429.     friend ostream &operator<<(ostream &out, AugmentedMatrix &m)
  430.     {
  431.         out << m.matrixLeft << m.matrixRight;
  432.         return out;
  433.     }
  434.     // Method to execute forward elimination.
  435.     friend AugmentedMatrix<L, R> ForwardElimination(AugmentedMatrix<L, R> U, bool debug_info, int &no_perms)
  436.     {
  437.         // This will be our final augmented matrix stroing the upper triangle.
  438.         size_t curr_col = 0;
  439.         size_t rows = U.matrixLeft.get_rows();
  440.         // Iterate through the rows.
  441.         for (size_t i = 0; i < rows; i++)
  442.         {
  443.             int row_with_max_pivot = -1;
  444.             // Find the maximum pivot.
  445.             double max_pivot = U.matrixLeft[i][curr_col];
  446.             for (size_t j = i + 1; j < rows; j++)
  447.             {
  448.                 if (abs(U.matrixLeft[j][curr_col]) < EPS)
  449.                     continue;
  450.                 if (abs(U.matrixLeft[j][curr_col]) > abs(max_pivot))
  451.                 {
  452.                     row_with_max_pivot = j;
  453.                     max_pivot = U.matrixLeft[j][curr_col];
  454.                 }
  455.             }
  456.             // If there is a maximum pivot, then we should swap row with it with the current row.
  457.             if (row_with_max_pivot != -1)
  458.             {
  459.                 if (debug_info)
  460.                     cout << "step: permutation\n";
  461.                 PermutationMatrix P(rows, i, row_with_max_pivot);
  462.                 // To exchange the rows we multiply by permuatation matrix.
  463.                 U.matrixLeft = P * U.matrixLeft;
  464.                 U.matrixRight = P * U.matrixRight;
  465.                 if (debug_info)
  466.                     cout << U.matrixLeft;
  467.                 no_perms++;
  468.             }
  469.             // Go through the next rows and eliminate the values in the column under the pivot.
  470.             for (size_t j = i + 1; j < rows; j++)
  471.             {
  472.                 if (abs(U.matrixLeft[j][curr_col]) < EPS or abs(U.matrixLeft[i][curr_col]) < EPS)
  473.                     continue;
  474.                 if (debug_info)
  475.                     cout << "step: elimination\n";
  476.                 double elimK = U.matrixLeft[j][curr_col] / U.matrixLeft[i][curr_col];
  477.                 if (abs(elimK) < EPS)
  478.                     elimK = 0;
  479.                 EliminationMatrix E(rows, j, i, elimK);
  480.                 // Tho eliminate the value we multiply by elimination matrix.
  481.                 U.matrixLeft = E * U.matrixLeft;
  482.                 U.matrixRight = E * U.matrixRight;
  483.                 if (debug_info)
  484.                     cout << U.matrixLeft;
  485.             }
  486.             curr_col++;
  487.         }
  488.         return U;
  489.     }
  490.     friend AugmentedMatrix<L, R> BackwardElimination(AugmentedMatrix<L, R> U, bool debug_info)
  491.     {
  492.         int curr_col = U.matrixLeft.get_cols() - 1;
  493.         int rows = U.matrixLeft.get_rows();
  494.         for (int i = rows - 1; i >= 0; i--)
  495.         {
  496.             for (int j = i - 1; j >= 0; j--)
  497.             {
  498.                 if (abs(U.matrixLeft[j][curr_col]) < EPS || abs(U.matrixLeft[i][curr_col]) < EPS)
  499.                     continue;
  500.                 EliminationMatrix E(rows, j, i, U.matrixLeft[j][curr_col] / U.matrixLeft[i][curr_col]);
  501.                 U.matrixLeft = E * U.matrixLeft;
  502.                 U.matrixRight = E * U.matrixRight;
  503.                 if (debug_info)
  504.                     cout << U;
  505.             }
  506.             curr_col--;
  507.         }
  508.         return U;
  509.     }
  510.     friend ColumnVector solve(AugmentedMatrix<SquareMatrix, ColumnVector> A, bool debug_info);
  511. };
  512.  
  513. // Method to calculate the determinant.
  514. double determinant_worker(SquareMatrix A, bool debug_info)
  515. {
  516.     AugmentedMatrix<SquareMatrix, SquareMatrix> aug(A, A);
  517.     // The determinant can be calculated as a product
  518.     // of values on main diagonal after forward elimination.
  519.     int no_perms = 0;
  520.     SquareMatrix U = ForwardElimination(aug, debug_info, no_perms).getLeft();
  521.     double res = (no_perms % 2) ? -1.0 : 1.0;
  522.     for (size_t i = 0; i < A.get_cols(); i++)
  523.         res *= U[i][i];
  524.     return res;
  525. }
  526.  
  527. SquareMatrix inverse_worker(SquareMatrix A, bool debug_info)
  528. {
  529.     IdentityMatrix I(A.get_cols());
  530.     // SquareMatrix I_matrix = I;
  531.     if (debug_info)
  532.         cout << "step #0: Augmented Matrix\n";
  533.     AugmentedMatrix<SquareMatrix, SquareMatrix> aug(A, I);
  534.     if (debug_info)
  535.         cout << aug;
  536.  
  537.     if (debug_info)
  538.         cout << "Direct way:\n";
  539.     int no_perms = 0;
  540.     AugmentedMatrix<SquareMatrix, SquareMatrix> B = ForwardElimination(aug, debug_info, no_perms);
  541.  
  542.     if (debug_info)
  543.         cout << "Way back:\n";
  544.     AugmentedMatrix<SquareMatrix, SquareMatrix> C = BackwardElimination(B, debug_info);
  545.  
  546.     if (debug_info)
  547.         cout << "Diagonal normalization:\n";
  548.  
  549.     for (int i = 0; i < C.getLeft().get_rows(); i++)
  550.     {
  551.         double pivot = C.getLeft()[i][i];
  552.         for (int j = 0; j < C.getLeft().get_cols(); j++)
  553.             C.getLeft()[i][j] /= pivot;
  554.         for (int j = 0; j < C.getRight().get_cols(); j++)
  555.             C.getRight()[i][j] /= pivot;
  556.     }
  557.     if (debug_info)
  558.         cout << C;
  559.  
  560.     return C.getRight();
  561. }
  562.  
  563. ColumnVector solve(AugmentedMatrix<SquareMatrix, ColumnVector> A, bool debug_info)
  564. {
  565.     if (debug_info)
  566.     {
  567.         cout << "step #0:\n";
  568.         cout << A;
  569.     }
  570.  
  571.     int no_perms = 0;
  572.     AugmentedMatrix<SquareMatrix, ColumnVector> B = ForwardElimination(A, debug_info, no_perms);
  573.     AugmentedMatrix<SquareMatrix, ColumnVector> C = BackwardElimination(B, debug_info);
  574.  
  575.     if (debug_info)
  576.         cout << "Diagonal normalization:\n";
  577.     for (int i = 0; i < C.getLeft().get_rows(); i++)
  578.     {
  579.         double pivot = C.getLeft()[i][i];
  580.         if (abs(pivot) < EPS)
  581.             continue;
  582.         for (int j = 0; j < C.getLeft().get_cols(); j++)
  583.             C.getLeft()[i][j] /= pivot;
  584.         for (int j = 0; j < C.getRight().get_cols(); j++)
  585.             C.getRight()[i][j] /= pivot;
  586.     }
  587.     if (debug_info)
  588.         cout << C;
  589.  
  590.     return A.getRight();
  591. }
  592.  
  593. int main(void)
  594. {
  595.     cout << fixed << setprecision(4);
  596.  
  597.     int m;
  598.     cin >> m;
  599.     double *t = new double[m];
  600.     ColumnVector b(m);
  601.     for (int i = 0; i < m; i++)
  602.     {
  603.         double t_i, b_i;
  604.         cin >> t_i >> b_i;
  605.         t[i] = t_i;
  606.         b[i][0] = b_i;
  607.     }
  608.     int n;
  609.     cin >> n;
  610.  
  611.     Matrix A(m, n + 1);
  612.     for (int i = 0; i <= n; i++)
  613.         for (int j = 0; j < m; j++)
  614.             A[j][i] = pow(t[j], (double)i);
  615.  
  616.     cout << "A:\n";
  617.     cout << A;
  618.  
  619.     Matrix A_T = A.T();
  620.     SquareMatrix A_1 = A_T * A;
  621.  
  622.     cout << "A_T*A:\n";
  623.     cout << A_1;
  624.  
  625.     SquareMatrix A_2 = A_1.inverse(false);
  626.  
  627.     cout << "(A_T*A)^-1:\n";
  628.     cout << A_2;
  629.  
  630.     ColumnVector A_3 = A_T * b;
  631.  
  632.     cout << "A_T*b:\n";
  633.     cout << A_3;
  634.  
  635.     ColumnVector A_4 = Matrix(A_2) * A_3;
  636.  
  637.     cout << "x~:\n";
  638.     cout << A_4;
  639.  
  640. #if (defined(WIN32) || defined(_WIN32)) && USE_GNUPLOT
  641.     FILE *pipe = _popen(GNUPLOT_NAME, "w");
  642. #elif USE_GNUPLOT
  643.     FILE *pipe = popen(GNUPLOT_NAME, "w");
  644. #endif
  645. #if USE_GNUPLOT
  646.     fprintf(pipe, "%s\n", "set terminal png");
  647.     fprintf(pipe, "%s\n", "set output 'output.png'");
  648.     fprintf(pipe, "%s\n", "set title \"Least Squares Approximation\"");
  649.     fprintf(pipe, "%s\n", "set key noautotitle");
  650.     fprintf(pipe, "%s\n", "set autoscale xy");
  651.     fprintf(pipe, "%s\n", "set offsets 0.05, 0.05, 0.05, 0.05");
  652.     string func;
  653.     for (int i = 0; i <= n; i++)
  654.     {
  655.         if (A_4[i][0] < 0 and i != 0)
  656.             func = func.substr(0, func.size() - 1);
  657.         func += to_string(A_4[i][0]);
  658.         func += '*';
  659.         func += "x**";
  660.         func += to_string(i);
  661.         if (i != n)
  662.             func += '+';
  663.     }
  664.     cout << func << endl;
  665.     fprintf(pipe, "plot %s lw 3, '-' w p pt 7 ps 2\n", func.c_str());
  666.     for (int i = 0; i < m; i++)
  667.         fprintf(pipe, "%lf %lf\n", t[i], b[i][0]);
  668.     fprintf(pipe, "%s\n", "e");
  669.     fflush(pipe);
  670. #endif
  671. #if (defined(WIN32) || defined(_WIN32)) && USE_GNUPLOT
  672.     _pclose(pipe);
  673. #elif USE_GNUPLOT
  674.     pclose(pipe);
  675. #endif
  676.  
  677.     delete[] t;
  678.  
  679.     return 0;
  680. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement