Определитель

[function]

Здравствуйте! Написал программу, приводящую матрицу к диагональной элементарными преобразованиями. На некоторых матрицах она работает корректно, а на других -- нет. Например, если ввести такую:

30 20 15 12
20 15 12 15
15 12 15 20
12 15 20 30

программа выдаёт следующее:

30.000000  0.000000  0.000000  0.000000 
 0.000000  0.000000  0.000000  0.000000 
 0.000000  0.000000  -38280596832649256.000000  0.000000 
 0.000000  0.000000  0.000000  -10.349020 

Что интересно, определитель она считает верно, хотя видно в данном случае, что на диагонали есть нуль (а детерминант равен 2639). Помогите, пожалуйста, найти ошибку.

   
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

void det(double **m, int *k, int n, double *D);

void det(double **m, int *k, int n, double *D)
{
    
    int i, j, p, u, s;
    double x, y;
    
    (*D)=1;
    
    for(j=0; j<n; j++)
    { 
        for(i=0; i<n; i++)
        {  
          if((m[i][j]!=0)&&(j<=i))
          { 
              for(p=0; p<n; p++)
              { 
                  if(p!=i)
                  {
                    y=m[p][j]/m[i][j];
                    for(u=0; u<n; u++)
                    {
                        m[p][u]-=y*m[i][u];
                    }
                  }
                  
              }
             
             for(s=0; s<n; s++)
             {
               if(j!=i)
               {
                 x=m[j][s];
                 m[j][s]=m[i][s];
                 m[i][s]=x;
                 (*D)=(*D)*(-1);
               }
             }
          
         
          }
          
          
          
        }
    }
    
  

[Xaositect]

А как Вы делаете вывод матрицы? Может, это не нуль на диагонали, а маленькое число?

[function]

Xaositect в сообщении #1194794 писал(а):

А как Вы делаете вывод матрицы? Может, это не нуль на диагонали, а маленькое число?

for(i=0; i<n; i++)
  {
    for(j=0; j<k[i]; j++)
    {
      fprintf(f," %lf ",m[i][j]);
      if(j==k[i]-1) fprintf(f,"\n");
    }
  }

-- 23.02.2017, 15:29 --

Откуда может взяться это маленькое число?

[SomePupil]

function в сообщении #1194797 писал(а):

Откуда может взяться это маленькое число?

Цифр после запятой может не хватить $-$ даблы не мелочатся.

[function]

Как это можно исправить?

[arseniiv]

Заменить lf на другой спецификатор формата — например, e (или E, отличается регистром).

Справка по printf и спецификаторам формата.

[function]

arseniiv в сообщении #1194828 писал(а):

Заменить lf на другой спецификатор формата — например, e (или E, отличается регистром).

Теперь, конечно, они отображаются, но проблема не ушла: почему эти числа возникают, непонятно. За массив я не вылезаю, вроде, да и деления на ноль нет...

[Dmitriy40]

function
Мне думается дело в конечной точности деления. Т.е. код

y=m[p][j]/m[i][j]; y*m[i][u];

не обязательно (не для всех возможных аргументов) восстановит точное исходное значение m[p][j] при u=j, будет какая-то погрешность. И эта погрешность вылезет на следующей итерации внешних циклов и вместо точного равенства будет лишь приближённое. А потом будет деление на эту маленькую разницу - и получите большое число, которого на самом деле быть не должно.

[mustitz]

Напрашивается такой вариант: (1) выполнить все действия на калькуляторе (2) вывести определители после каждой итерации (3) сравнить результаты.

[Dmitriy40]

Калькуляторы тоже не всегда правильно округляют, так что это не панацея.
Если входные данные в виде целых или рациональных чисел, то выходом было бы проводить все операции тоже в рациональных числах (с числителем и знаменателем дроби), при этом результат будет точным. Для действительных чисел решение в общем виде без переделки алгоритма вряд ли существует (впрочем могу и не знать).

[function]

Dmitriy40 в сообщении #1194837 писал(а):

function
Мне думается дело в конечной точности деления. Т.е. код

y=m[p][j]/m[i][j]; y*m[i][u];

не обязательно (не для всех возможных аргументов) восстановит точное исходное значение m[p][j] при u=j, будет какая-то погрешность. И эта погрешность вылезет на следующей итерации внешних циклов и вместо точного равенства будет лишь приближённое. А потом будет деление на эту маленькую разницу - и получите большое число, которого на самом деле быть не должно.

Спасибо Вам большое! Я просто обнулил соответствующие позиции ($u=j$), и всё получилось!

[Dmitriy40]

function
Я хотел сразу предложить этот метод, но потом подумал что погрешность при вычислении остальных элементов (при $u\ne j$) так не уберётся. И при других исходных данных может снова вылезти. Так что не считайте это окончательным решением.

[arseniiv]

Известна ещё модификация этого метода с перестановкой строк так, чтобы делилось всегда на как можно большее по модулю число (кажется), и в результате выдаются кроме матрицы (если что-то нужно считать дальше) ещё и новые индексы строк.

[Dmitriy40]

Да, известна, но такая модификация лишь облегчает проблему потери точности, не исключая её полностью.

[arseniiv]

Ну, это понятно. Я думаю, ваш совет насчёт рациональных чисел, если в все элементы матрицы предполагаются небольшими целыми, самый лучший. Правда, наверняка они здесь небольшие и целые просто для простоты (хотя ещё лучше не угадывать за ТС).