• 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 27 grudnia 2022 przez mokrowski Mędrzec (156,220 p.)
A debuggowałeś ten kod? Jeśli tak, to co Ci powiedział debugger?
komentarz 27 grudnia 2022 przez Krzysztofs1234 Użytkownik (890 p.)

W zasadzie to nic takiego:

komentarz 27 grudnia 2022 przez mokrowski Mędrzec (156,220 p.)
Hmm... poczytaj może jak użyć debuggera. To co pokazałeś to jego uruchomienie bez pułapek.
komentarz 27 grudnia 2022 przez Krzysztofs1234 Użytkownik (890 p.)

Które okno powinno mnie interesować? Bo na przykład Watches nie daje mi możliwości podglądu wartości przypisywanych do C[x][x]. Z kolei CPU Registers i Disassembly niewiele mi mówi.

komentarz 28 grudnia 2022 przez Oscar Nałogowiec (29,320 p.)

@Krzysztofs1234, Chodzi ci o te -1e-19? To normalne odchylenia (słowo błąd sugeruje możliwość poprawy) wyniku. Odejmujesz dwie liczby o rzędzie wielkości 1 i gdzieś na 19. pozycji za przecinkiem zostaje jakaś różnica między nimi. Tak to już jest z obliczeniami zmiennoprzecinkowymi. W zasadzie nigdy nie należy stosować z nimi operatora == bo dwie takie liczby nigdy nie są równe, chyba że po bezpośrednim przypisaniu jednej zmiennej do drugiej. Należy sprawdzać czy są dostatecznie bliskie. zwykły float ma dokładność jakiś 24 bitów czyli 16 mln - mniej więcej 8 cyfr, double  jest 2 razy dokładniejszy ok 16 cyfr. Procesory Intela używają wewnętrznie 80-bitowego formatu do obliczeń na jednostce zmiennoprzecinkowej, ale potem jest to przybliżane go 64 bitów.

Pośrednie wartości wyświetlają się z dokładnością do jakiś 7 cyfr po przecinku - ale to jest przybliżenie wyświetlania - same wartości w pamięci komputera są znacznie dokładniejsze. To przybliżenie wyświetlania ukrywa fakt, że one nie są tak naprawde tak "równe" jak to widać. Np. takie 1.66667 to gdzieś kilkanaście miejsc po przecinku nie ma końcówki 667 tylko jakieś dziwne cyferki. I jak odejmiesz od siebie dwie takie wartości to zostanie różnica tych końcowych ogonów. W porównaniu do argumentów (rzędu 1) to ten wynik to zero.

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

Więc twierdzisz, że to błąd przybliżenia spowodował taki wynik. Ja z początku też tak myślałem, ale, jak wcześniej wspomniałem, sprawdziłem, co się stanie jeśli przepiszę obliczenia do strumienia wyjścia:

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

complex<long double> x(0, 0);
x+=c(2)*c(11/36);
x+=c(10)*c(1/18);
x+=c(-7)*c(1/6);
cout << x;

Nie podważam tego, że przybliżenie może generować i generuje błędy, bo jest to rzecz naturalna, ale mam wątpliwości, czy akurat w tym przypadku tak jest, bo wykonanie powyższego kodu spowoduje wyświetlenie wartości (0, 0).

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

Hmm....

#include <iostream>
#include <complex>

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

int main() {

	std::complex<long double> x(0.0, 0.0);
	x+=c(2.0)*c(11.0/36.0);
	x+=c(10.0)*c(1.0/18.0);
	x+=c(-7.0)*c(1.0/6.0);
	std::cout << x;
}

Wiesz... 1/6 to nie to samo co 1.0/6.0.

komentarz 28 grudnia 2022 przez Krzysztofs1234 Użytkownik (890 p.)
Ale dlaczego tak jest? Czy liczby 1 i 1.0 potencjalnie nie powinny mieć jedynie zer po przecinku?

Obawiam się, że to może być ten błąd, przez który mam niepoprawne wyniki. W jaki sposób mogę te błędy zredukować, żeby jednak tam na początku wyszło mi 1?
komentarz 28 grudnia 2022 przez Oscar Nałogowiec (29,320 p.)
1 - to int.

1.0 - to float albo double.

 

Dzielenie intów do dzielenie całkowite - 1/6 to po prostu 0 (reszta z dzielenia: 1)

1.0/6.0 to ok .1(6)

W programie masz zmienne o zadeklarowanym typie.
komentarz 28 grudnia 2022 przez Krzysztofs1234 Użytkownik (890 p.)

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

Miałem więc mnożenie przez 0 i dodawanie 0, więc niedziwne, że mi 0 wyszło. :D

    complex<long double> x(0, 0);
    x+=c(2)*c(11.0/36);
    x+=c(10)*c(1.0/18);
    x+=c(-7)*c(1.0/6);
    cout << x << endl;

    std::complex<long double>y(0.0, 0.0);
    y+=c(2.0)*c(11.0/36.0);
    y+=c(10.0)*c(1.0/18.0);
    y+=c(-7.0)*c(1.0/6.0);
    std::cout << y;

Teraz wyniki są takie same, więc pytanie:

Ale dlaczego tak jest? Czy liczby 1 i 1.0 potencjalnie nie powinny mieć jedynie zer po przecinku?

wycofuję, bo było postawione na złych fundamentach.

 Pozostaję przy podejrzeniach, że jest to błąd przybliżenia, który zaburza wynik. Jeśli tak jest, co można z tym zrobić?

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

3 - to int

3.4F - to float

3.4 - to double

Tak więc:

float val = 234.4;

to operacja przypisania wartości double (234.4) do wartości float. Nastąpi utrata dokładności

float val = 234.4F;

to operacja przypisania wartości float do float. Nie nastąpi utrata dokładności

double val = 234.4;

to przypisanie wartości double do double. Nie nastąpi utrata dokładności

Konsekwencje...

float v = 0.1;
if (v != 0.1) {
    std::cout << "Ups.. jednak nie jest równe!\n";
}

Komunikat wyświetli się bo w if(...), porównujesz float (zmienna v) do double (wartość 0.1) a one nie są sobie równe.

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

Zdaję sobie sprawę ze wszystkiego, co jest napisane w powyższym komentarzu, już od lat. Moje pytania pozostają niezmienne:

Pozostaję przy podejrzeniach, że jest to błąd przybliżenia, który zaburza wynik. Jeśli tak jest, co można z tym zrobić?

 

komentarz 28 grudnia 2022 przez Oscar Nałogowiec (29,320 p.)
Nastąpi utrata dokładności w każdym przypadku, tylko trochę inna - większość ułamków o skończonej reprezentacji dziesiętnej w nie jest dokładnie reprezentowalne w zapisie binarnym. Dokładnie zapisują się ułamki typu 0.5, 0.25, 0.125 czyli 1/2. 1/4, 1/8 itp
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ź 548 wizyt
+1 głos
1 odpowiedź 714 wizyt
pytanie zadane 7 stycznia 2016 w C i C++ przez dreez Początkujący (320 p.)
0 głosów
3 odpowiedzi 677 wizyt
pytanie zadane 4 stycznia 2016 w C i C++ przez Macek Kolo Mądrala (5,480 p.)

92,696 zapytań

141,607 odpowiedzi

320,114 komentarzy

62,056 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!

...