Обработка потери точности при вычитании двух двойных, которые находятся близко друг к другу

У меня есть проект, где мы должны решить матричное уравнение AX = B для x, учитывая, что A является трехдиагональной матрицей. Я сделал этот проект на C ++, получил программу для производства правильных Matrix X, но при попытке сообщить об ошибке пользователю, A*X-BЯ получаю ошибочную ошибку !! Это связано с тем, что я вычитаю A*X а также B, чьи записи произвольно близки друг к другу. У меня было две идеи о том, как справиться с этим, поэлементно:

1. Согласно этой статье, http://en.wikipedia.org/wiki/Loss_of_significance, может быть столько, сколько -log2(1-y/x) биты, потерянные в прямом вычитании x-y, Давайте масштабировать оба x а также y от pow(2,bitsLost), вычтите два, а затем уменьшите их, разделив на pow(2,bitsLost)
2. В курсах числовых методов это подчеркивается: возьми арифметическое сопряжение! Вместо double difference = x-y; использование double difference = (x*x-y*y)/(x+y);

Хорошо, так почему же вы не выбрали метод и пошли дальше?

Я попробовал все три метода (включая прямое вычитание) здесь: http://ideone.com/wfkEUp . Я хотел бы знать две вещи:

1. Между методом «масштабирования и удаления накипи» (для которого я намеренно выбрал степень двойки) и методом арифметического сопряжения, какой из них дает меньшую ошибку (с точки зрения вычитания больших чисел)?
2. Какой метод является вычислительно более эффективным? /*For this, I was going to say the scaling method was going to be more efficient with a linear complexity versus the seemed quadratic complexity of the conjugate method, but I don't know the complexity of log2()*/

Любая помощь будет приветствоваться!

П.С .: Все три метода, похоже, возвращают одно и то же double в примере кода …

Давайте посмотрим на ваш код
Нет проблем; вот мой код Matrix.cpp

#include "ExceptionType.h"#include "Matrix.h"#include "MatrixArithmeticException.h"#include <iomanip>
#include <iostream>
#include <vector>

Matrix::Matrix()
{
//default size for Matrix is 1 row and 1 column, whose entry is 0
std::vector<long double> rowVector(1,0);
this->matrixData.assign(1, rowVector);
}

Matrix::Matrix(const std::vector<std::vector<long double> >& data)
{
this->matrixData = data;
//validate matrixData
validateData();
}

//getter functions
//Recall that matrixData is a vector of a vector, whose elements should be accessed like matrixData[row][column].
//Each rowVector should have the same size.
unsigned Matrix::getRowCount() const { return matrixData.size(); }

unsigned Matrix::getColumnCount() const { return matrixData[0].size(); }

//matrix validator should just append zeroes into row vectors that are of smaller dimension than they should be...
void Matrix::validateData()
{
//fetch the size of the largest-dimension rowVector
unsigned largestSize = 0;
for (unsigned i = 0; i < getRowCount(); i++)
{
if (largestSize < matrixData[i].size())
largestSize = matrixData[i].size();
}
//make sure that all rowVectors are of that dimension
for (unsigned i = 0; i < getRowCount(); i++)
{
//if we find a rowVector where this isn't the case
if (matrixData[i].size() < largestSize)
{
//add zeroes to it so that it becomes the case
matrixData[i].insert(matrixData[i].end(), largestSize-matrixData[i].size(), 0);
}
}

}
//operators
//+ and - operators should check to see if the size of the first matrix is exactly the same size as that of the second matrix
Matrix Matrix::operator+(const Matrix& B)
{
//if the sizes coincide
if ((getRowCount() == B.getRowCount()) && (getColumnCount() == B.getColumnCount()))
{
//declare the matrixData
std::vector<std::vector<long double> > summedData = B.matrixData;    //since we are in the scope of the Matrix, we can access private data members
for (unsigned i = 0; i < getRowCount(); i++)
{
for (unsigned j = 0; j < getColumnCount(); j++)
{
summedData[i][j] += matrixData[i][j];   //add the elements together
}
}
//return result Matrix
return Matrix(summedData);
}
else
throw MatrixArithmeticException(DIFFERENT_DIMENSIONS);
}

Matrix Matrix::operator-(const Matrix& B)
{
//declare negativeB
Matrix negativeB = B;
//negate all entries
for (unsigned i = 0; i < negativeB.getRowCount(); i++)
{
for (unsigned j = 0; j < negativeB.getColumnCount(); j++)
{
negativeB.matrixData[i][j] = 0-negativeB.matrixData[i][j];
}
}
//simply add the negativeB
try
{
return ((*this)+negativeB);
}
catch (MatrixArithmeticException& mistake)
{
//should exit or do something similar
std::cout << mistake.what() << std::endl;
}
}

Matrix Matrix::operator*(const Matrix& B)
{
//the columnCount of the left operand must be equal to the rowCount of the right operand
if (getColumnCount() == B.getRowCount())
{
//if it is, declare data with getRowCount() rows and B.getColumnCount() columns
std::vector<long double> zeroVector(B.getColumnCount(), 0);
std::vector<std::vector<long double> > data(getRowCount(), zeroVector);
for (unsigned i = 0; i < getRowCount(); i++)
{
for (unsigned j = 0; j < B.getColumnCount(); j++)
{
long double sum = 0; //set sum to zero
for (unsigned k = 0; k < getColumnCount(); k++)
{
//add the product of matrixData[i][k] and B.matrixData[k][j] to sum
sum += (matrixData[i][k]*B.matrixData[k][j]);
}
data[i][j] = sum;   //assign the sum to data
}
}
return Matrix(data);
}
else
{
throw MatrixArithmeticException(ROW_COLUMN_MISMATCH); //dimension mismatch
}
}

std::ostream& operator<<(std::ostream& outputStream, const Matrix& theMatrix)
{
//Here, you should use the << again, just like you would for ANYTHING ELSE.
//first, print a newline
outputStream << "\n";
//setting precision (optional)
outputStream.precision(11);
for (unsigned i = 0; i < theMatrix.getRowCount(); i++)
{
//print '['
outputStream << "[";
//format stream(optional)
for (unsigned j = 0; j < theMatrix.getColumnCount(); j++)
{
//print numbers
outputStream << std::setw(17) << theMatrix.matrixData[i][j];
//print ", "if (j < theMatrix.getColumnCount() - 1)
outputStream << ", ";
}
//print ']'
outputStream << "]\n";
}
return outputStream;
}