fork(29) download
  1. #include <cstdio>
  2. #include <algorithm>
  3. #include <cmath>
  4. //очистить выделенную память
  5. void clear(double **arr, int n)
  6. {
  7. for(int i = 0; i < n; i++)
  8. delete[] arr[i];
  9. delete[] arr;
  10. }
  11. //создать копию массива
  12. double** clone(double **arr, int n)
  13. {
  14. double **newArr = new double*[n];
  15. for(int row = 0; row < n; row++)
  16. {
  17. newArr[row] = new double[n];
  18. for(int col = 0; col < n; col++)
  19. newArr[row][col] = arr[row][col];
  20. }
  21. return newArr;
  22. }
  23. //print the matrix
  24. void show(double **matrix, int n)
  25. {
  26. for(int row = 0; row < n; row++){
  27. for(int col = 0; col < n; col++)
  28. printf("%lf\t", matrix[row][col]);
  29. printf("\n");
  30. }
  31. printf("\n");
  32. }
  33. //матричное умножение матриц
  34. double** matrix_multi(double **A, double **B, int n)
  35. {
  36. double **result = new double*[n];
  37. //заполнение нулями
  38. for(int row = 0; row < n; row++){
  39. result[row] = new double[n];
  40. for(int col = 0; col < n; col++){
  41. result[row][col] = 0;
  42. }
  43. }
  44.  
  45. for(int row = 0; row < n; row++){
  46. for(int col = 0; col < n; col++){
  47. for(int j = 0; j < n; j++){
  48. result[row][col] += A[row][j] * B[j][col];
  49. }
  50. }
  51. }
  52. return result;
  53. }
  54. //умножение матрицы на число
  55. void scalar_multi(double **m, int n, double a){
  56. for(int row = 0; row < n; row++)
  57. for(int col = 0; col < n; col++){
  58. m[row][col] *= a;
  59. }
  60. }
  61. //вычисление суммы двух квадратных матриц
  62. void sum(double **A, double **B, int n)
  63. {
  64. for(int row = 0; row < n; row++)
  65. for(int col = 0; col < n; col++)
  66. A[row][col] += B[row][col];
  67. }
  68.  
  69. //вычисление определителя
  70. double det(double **matrix, int n) //квадратная матрица размера n*n
  71. {
  72. double **B = clone(matrix, n);
  73. //приведение матрицы к верхнетреугольному виду
  74. for(int step = 0; step < n - 1; step++)
  75. for(int row = step + 1; row < n; row++)
  76. {
  77. double coeff = -B[row][step] / B[step][step]; //метод Гаусса
  78. for(int col = step; col < n; col++)
  79. B[row][col] += B[step][col] * coeff;
  80. }
  81. //Рассчитать определитель как произведение элементов главной диагонали
  82. double Det = 1;
  83. for(int i = 0; i < n; i++)
  84. Det *= B[i][i];
  85. //Очистить память
  86. clear(B, n);
  87. return Det;
  88. }
  89.  
  90. int main()
  91. {
  92. //Исходная матрица, динамический двухмерный массив
  93. int n; scanf("%d", &n);
  94. double **A = new double*[n];
  95. for(int row = 0; row < n; row++)
  96. {
  97. A[row] = new double[n];
  98. for(int col = 0; col < n; col++)
  99. scanf("%lf", &A[row][col]);
  100. }
  101.  
  102. /* Численное вычисление обратной матрицы по методу Ньютона-Шульца
  103.   1. Записать начальное приближение [Pan, Schreiber]:
  104.   1) Транспонировать данную матрицу
  105.   2) Нормировать по столбцам и строкам
  106.   2. Повторять процесс до достижения заданной точности.
  107.   */
  108.  
  109. double N1 = 0, Ninf = 0; //норма матрицы по столбцам и по строкам
  110. double **A0 = clone(A, n); //инициализация начального приближения
  111. for(size_t row = 0; row < n; row++){
  112. double colsum = 0, rowsum = 0;
  113. for(size_t col = 0; col < n; col++){
  114. rowsum += fabs(A0[row][col]);
  115. colsum += fabs(A0[col][row]);
  116. }
  117. N1 = std::max(colsum, N1);
  118. Ninf = std::max(rowsum, Ninf);
  119. }
  120. //транспонирование
  121. for(size_t row = 0; row < n - 1; row++){
  122. for(size_t col = row + 1; col < n; col++)
  123. std::swap(A0[col][row], A0[row][col]);
  124. }
  125. scalar_multi(A0, n, (1 / (N1 * Ninf))); //нормирование матрицы
  126. //инициализация удвоенной единичной матрицы нужного размера
  127. double **E2 = new double*[n];
  128. for(int row = 0; row < n; row++)
  129. {
  130. E2[row] = new double[n];
  131. for(int col = 0; col < n; col++){
  132. if(row == col)
  133. E2[row][col] = 2;
  134. else
  135. E2[row][col] = 0;
  136. }
  137. }
  138. double **inv = clone(A0, n); //A_{0}
  139. double EPS = 0.001; //погрешность
  140. if(det(A, n) != 0){ //если матрица не вырождена
  141. while(fabs(det(matrix_multi(A, inv, n), n) - 1) >= EPS) //пока |det(A * A[k](^-1)) - 1| >= EPS
  142. {
  143. double **prev = clone(inv, n); //A[k-1]
  144. inv = matrix_multi(A, prev, n); //A.(A[k-1]^(-1))
  145. scalar_multi(inv, n, -1); //-A.(A[k-1]^(-1))
  146. sum(inv, E2, n); //2E - A.(A[k-1]^(-1))
  147. inv = matrix_multi(prev, inv, n); //(A[k-1]^(-1)).(2E - A.(A[k-1]^(-1)))
  148. clear(prev, n);
  149. }
  150. //вывод матрицы на экран
  151. show(inv, n);
  152. }
  153. else
  154. printf("Impossible\n");
  155. clear(A, n);
  156. clear(E2, n);
  157. return 0;
  158. }
Success #stdin #stdout 0s 3436KB
stdin
2
1 2
2 1
stdout
-0.333067	0.666400	
0.666400	-0.333067