Tiny числа вместо нуля?

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

Я ввожу произвольную матрицу как таковую:

2 1 1
1 2 1
1 1 2

И получил вычисление обратного, и я получил правильный результат:

0.75 -0.25 -0.25
-0.25 0.75 -0.25
-0.25 -0.25 0.75

Но когда я попытался перемножить эти два, чтобы убедиться, что у меня есть единичная матрица, я получил:

1 5.5111512e-017 0
0 1 0
-1.11022302e-0.16 0 1

Почему я получаю эти результаты? Я бы понял, если бы я умножал странные числа, где я мог бы понять некоторые ошибки округления, но полученная сумма такова:

2 * -0.25 + 1 * 0.75 + 1 * -0.25

, что явно равно 0, а не 5.111512e-017

Если бы я вручную заставил его выполнить расчет; например:

std::cout << (2 * -0.25 + 1 * 0.75 + 1 * -0.25) << "\n";

Я получаю 0, как и ожидалось?

Все числа представлены как двойные. Вот моя перегрузка умножения:

Matrix operator*(const Matrix& A, const Matrix& B)
{
    if(A.get_cols() == B.get_rows())
    {
        Matrix temp(A.get_rows(), B.get_cols());
        for(unsigned m = 0; m < temp.get_rows(); ++m)
        {
            for(unsigned n = 0; n < temp.get_cols(); ++n)
            {
                for(unsigned i = 0; i < temp.get_cols(); ++i)
                {
                    temp(m, n) += A(m, i) * B(i, n);
                }
            }
        }

        return temp;
    }

    throw std::runtime_error("Bad Matrix Multiplication");
}

и функции доступа:

double& Matrix::operator()(unsigned r, unsigned c)
{
    return data[cols * r + c];
}

double Matrix::operator()(unsigned r, unsigned c) const
{
    return data[cols * r + c];
}

Вот функция для поиска инверсии:

Matrix Inverse(Matrix& M)
{
    if(M.rows != M.cols)
    {
        throw std::runtime_error("Matrix is not square");
    }

    int r = 0;
    int c = 0;
    Matrix augment(M.rows, M.cols*2);
    augment.copy(M);

    for(r = 0; r < M.rows; ++r)
    {
        for(c = M.cols; c < M.cols * 2; ++c)
        {
            augment(r, c) = (r == (c - M.cols) ? 1.0 : 0.0);
        }
    }

    for(int R = 0; R < augment.rows; ++R)
    {
        double n = augment(R, R);
        for(c = 0; c < augment.cols; ++c)
        {
            augment(R, c) /= n;
        }

        for(r = 0; r < augment.rows; ++r)
        {
            if(r == R) { continue; }
            double a = augment(r, R);

            for(c = 0; c < augment.cols; ++c)
            {
                augment(r, c) -= a * augment(R, c);
            }
        }
    }

    Matrix inverse(M.rows, M.cols);
    for(r = 0; r < M.rows; ++r)
    {
        for(c = M.cols; c < M.cols * 2; ++c)
        {
            inverse(r, c - M.cols) = augment(r, c);
        }
    }

    return inverse;
}
5
задан Rarge 3 June 2011 в 18:18
поделиться