# Solution

Mar 17th, 2023 (edited)
566
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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; }
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.     }
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]; }
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.     }
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.     }
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.     }
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=;
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.     }
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.     }
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.     }
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=;
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.     }
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.     }
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.     }
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.     }
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. }