• Najnowsze pytania
  • Bez odpowiedzi
  • Zadaj pytanie
  • Kategorie
  • Tagi
  • Zdobyte punkty
  • Ekipa ninja
  • IRC
  • FAQ
  • Regulamin
  • Książki warte uwagi

Macierze - mnożenie przez odwrotność - dziwna liczba psująca wynik

Object Storage Arubacloud
0 głosów
275 wizyt
pytanie zadane 26 grudnia 2022 w C i C++ przez Krzysztofs1234 Użytkownik (890 p.)
edycja 27 grudnia 2022 przez Krzysztofs1234

Dzień dobry,

obecnie zajmują się przeciążaniem operatora "/" w celu łatwego napisania mnożenia macierzy A przez odwrotność macierzy B (nieformalnie: dzielenie macierzy A przez macierz B). Przetestowałem działanie funkcji na odwracalnej macierzy, pisząc Ab/Ab, więc powinienem uzyskać macierz jednostkową. Część obliczeń się zgadza, ale w pewnym etapie, np. C[1][0], element przyjmuje pewną dziwną wartość i psuje wynik. Myślałem, że to przez przybliżenia, ale te same działania wpisane w strumień wyjścia dają wynik prawidłowy, czyli w tym przypadku 0.

Uprzejmie proszę o pomoc.

 

#include <iostream>
#include <vector>
#include <complex>
using namespace std;


void createIdentifyMatrix(vector<vector<complex<long double>>> &A, const unsigned n)
{
    if(n==0)
    {
        cerr << "Wymiar takiej macierzy nie moze byc rowny 0." << endl << endl;
        return;
    }

    const vector<complex<long double>> b(n);
    for(unsigned i=0; i<n; i++)
    {
        A.push_back(b);
        A[i][i]=1;
    }
}

void printMatrix(const vector<vector<complex<long double>>> &a)
{
    unsigned i, j;

    for(i=0; i<a.size(); i++)
    {
        for(j=0; j<a[i].size(); j++)
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
}

void matrixInversion(vector<vector<complex<long double>>> &A)
{
    const unsigned n=A.size();
    vector<vector<complex<long double>>> I;
    complex<long double> num;
    unsigned i, j, k;

    createIdentifyMatrix(I, n);

    for(j=0; j<n; j++)
    {
        k=j;
        for(i=j+1; i<A.size(); i++)
        {
            if(fabs(A[i][j])>fabs(A[k][j]))
            {
                k=i;
            }
        }

        if(A[k][j]==complex<long double>(0, 0))
        {
            cout << "Uklad osobliwy." << endl;
            throw;
        }

        if(A[j][j]!=A[k][j])
        {
            swap(A[j], A[k]);
            swap(I[j], I[k]);
        }

        num=A[j][j];
        A[j][j]=1;
        for(i=0; i<n; i++)
        {
            I[j][i]/=num;
            if(i==j) continue;
            A[j][i]/=num;
        }

        for(i=0; i<n; i++)
        {
            if(i==j) continue;
            num=A[i][j]/A[j][j];
            A[i][j]=0;

            for(k=0; k<=j; k++)
            {
                I[i][k]-=I[j][k]*num;
            }

            for(k=j+1; k<n; k++)
            {
                A[i][k]-=A[j][k]*num;
                I[i][k]-=I[j][k]*num;
            }
        }
    }

    A=I;
}

vector<vector<complex<long double>>> operator /(const vector<vector<complex<long double>>> &A, const vector<vector<complex<long double>>> &B)
{
    const unsigned r=A.size(), c=B.size();

    vector<vector<complex<long double>>> C(r, vector<complex<long double>>(c)), B1=B;
    matrixInversion(B1);

    unsigned i, j, k;

    for(i=0; i<r; i++)
    {
        for(j=0; j<c; j++)
        {
            cout << "C[" << i << "][" << j << "]" << endl;
            for(k=0; k<c; k++)
            {
                cout << C[i][j] << "+";
                C[i][j]+=A[k][i]*B1[j][k];
                cout << A[k][i] << "*" << B1[j][k] << "=" << C[i][j] << endl;
            }
            cout << endl;
        }
    }

    return C;
}


complex<long double> c(const long double Re, const long double Im=0) {return complex<long double>(Re, Im);}


int main()
{
    const vector<vector<vector<complex<long double>>>>
                                                      G{

                                                        {
                                                         {0,1,2,3,5},
                                                         {1,0,1,2,4},
                                                         {1,2,3,4,6},
                                                         {0,1,3,2,4}
                                                        },

                                                        {
                                                         {2,4,2,0,4},
                                                         {1,0,-1,1,2},
                                                         {0,1,3,-1,0},
                                                         {2,1,2,1,6}
                                                        },

                                                        {
                                                         {1,2,-1,2},
                                                         {2,1,0,3},
                                                         {-1,1,2,2}
                                                        },

                                                        {
                                                         {1,2,3,1},
                                                         {4,5,6,3},
                                                         {7,8,9,2}
                                                        }

                                                       },

                                                      Ch{

                                                         {
                                                          {1,1,-2,-2},
                                                          {1,2,-1,-4},
                                                          {-2,-1,6,3},
                                                          {-2,-4,3,11}
                                                         },

                                                         {
                                                          {4,2,-2},
                                                          {2,10,-7},
                                                          {-2,-7,6}
                                                         },

                                                         {
                                                          {9,-9,-6,9},
                                                          {-9,13,10,-11},
                                                          {-6,10,17,-5},
                                                          {9,-11,-5,15}
                                                         },

                                                         {
                                                          {c(4, 0), c(0, -2), c(2, -2)},
                                                          {c(0, 2), c(17, 0), c(9, 5)},
                                                          {c(2, 2), c(9, -5), c(8, 0)}
                                                         }

                                                        };

    vector<vector<complex<long double>>> Ab=Ch[1];
    printMatrix(Ab/Ab);

    return 0;
}

komentarz 28 grudnia 2022 przez mokrowski Mędrzec (156,220 p.)

@Krzysztofs1234, 

Zdaję sobie sprawę ze wszystkiego, co jest napisane w powyższym komentarzu, już od lat.

Jeśli wiesz o tym "od lat", to dlaczego dzieląc int przez int oczekujesz wyniku float? Może jednak warto cokolwiek poczytać?

https://docs.oracle.com/cd/E19957-01/800-7895/800-7895.pdf

https://archive.org/details/SterbenzFloatingPointComputation/mode/2up

https://www.amazon.com/Numerical-Recipes-3rd-Scientific-Computing/dp/0521880688

 

Rozwiązań masz wiele (a jeszcze więcej znajdziesz w literaturze):

1. Potraktować wartość "wystarczająco małą" jako 0 (sprawdzając ile wynosi epsilon danej platformy https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon )

2. Potraktować wartość "wystarczająco bliską X" jako X

3. Zastosować kontrolowane zaokrąglanie

4. Stosować specjalizowane algorytmy redukujące błąd dla obliczeń zmiennoprzecinkowych https://en.wikipedia.org/wiki/Floating-point_error_mitigation https://en.wikipedia.org/wiki/Kahan_summation_algorithm ...

5. Rezygnować z arytmetyki zmiennoprzecinkowej, na rzecz innej...

komentarz 28 grudnia 2022 przez Krzysztofs1234 Użytkownik (890 p.)

 Jeśli wiesz o tym "od lat", to dlaczego dzieląc int przez int oczekujesz wyniku float? Może jednak warto cokolwiek poczytać?

 Jak wspomniałem, była to pomyłka.

Dzięki za zwrócenie uwagi. W ostatnim kodzie, który wrzuciłem, nie zwróciłem uwagi, że przecież jest to dzielenie całkowite, przepisując po prostu ułamki w wyniku jak leci. :P

Z literatury będzie mi trudno skorzystać z uwagi na niewystarczającą znajomość angielskiego, żeby przeczytać całość ze zrozumieniem. Przez to nawet część napisana językiem matematyki może być bez kontekstu niezrozumiała. Mimo tego, dzięki za polecenie.

A może... czytałeś coś traktującego o tym samym, ale napisane (czy przetłumaczone) w języku polskim?

komentarz 28 grudnia 2022 przez Krzysztofs1234 Użytkownik (890 p.)

Im szybciej pogodzisz się z faktem że na większość zagadnień nie znajdziesz literatury w języku polskim, tym lepiej.

Może nie dasz wiary, ale jestem tego świadom od lat. :D

Dzięki za polecenia. Jestem też pod wrażeniem ilości materiału, jaki pochłonąłeś na ten temat. Skoro są publikacje na temat po kilkaset stron, to zdaje się to być nieco szersze zagadnienie niż z początku mi się zdawało.

komentarz 28 grudnia 2022 przez mokrowski Mędrzec (156,220 p.)
Poprzednio polecane pozycje przeczytałem uważnie. Czyli:

https://docs.oracle.com/cd/E19957-01/800-7895/800-7895.pdf

https://archive.org/details/SterbenzFloatingPointComputation/mode/2up

https://www.amazon.com/Numerical-Recipes-3rd-Scientific-Computing/dp/0521880688

Polskie zweryfikowałem tylko pod względem "jest ok, nie ma błędów, uproszczenia nie prowadzą na manowce". BTW... książki z computer science, traktuje się często wybiórczo poszukując określonych rozwiązań i algorytmów. Rzadko czyta się je z wypiekami na twarzy jak kryminał (od deski do deski)  z nadzieją na znalezienie informacji kto zabił...

Zaloguj lub zarejestruj się, aby odpowiedzieć na to pytanie.

Podobne pytania

0 głosów
1 odpowiedź 547 wizyt
+1 głos
1 odpowiedź 712 wizyt
pytanie zadane 7 stycznia 2016 w C i C++ przez dreez Początkujący (320 p.)
0 głosów
3 odpowiedzi 670 wizyt
pytanie zadane 4 stycznia 2016 w C i C++ przez Macek Kolo Mądrala (5,480 p.)

92,680 zapytań

141,583 odpowiedzi

320,068 komentarzy

62,041 pasjonatów

Motyw:

Akcja Pajacyk

Pajacyk od wielu lat dożywia dzieci. Pomóż klikając w zielony brzuszek na stronie. Dziękujemy! ♡

Oto polecana książka warta uwagi.
Pełną listę książek znajdziesz tutaj.

Akademia Sekuraka

Niedawno wystartował dodruk tej świetnej, rozchwytywanej książki (około 940 stron). Mamy dla Was kod: pasja (wpiszcie go w koszyku), dzięki któremu otrzymujemy 10% zniżki - dziękujemy zaprzyjaźnionej ekipie Sekuraka za taki bonus dla Pasjonatów! Książka to pierwszy tom z serii o ITsec, który łagodnie wprowadzi w świat bezpieczeństwa IT każdą osobę - warto, polecamy!

...