title: Handling loss of precision in subtracting two doubles that are close to each other
body: <p>I have a project to do where we are to solve the matrix equation AX=B for x, given that A is a tridiagonal matrix. I did this project in C++, got the program to produce the right <code>Matrix</code> X, but when trying to report the error back to the user, <code>A*X-B</code>, I get an erroneous error!! It is due to the fact that I am subtracing <code>A*X</code> and <code>B</code>, whose entries are arbitrarily close to each other. I had two ideas on how to handle this, element-by-element: </p>
  
  <ol>
  <li>According to this article, <a href="http://en.wikipedia.org/wiki/Loss_of_significance" rel="nofollow">http://en.wikipedia.org/wiki/Loss_of_significance</a>, there could be as many as <code>-log2(1-y/x)</code> bits lost in straight subtraction <code>x-y</code>. Let's scale both <code>x</code> and <code>y</code> by <code>pow(2,bitsLost)</code>, subtract the two, and then scale them back down by dividing by <code>pow(2,bitsLost)</code></li>
  <li>Stressed so much in the numeric methods course this is for: take the arithmetic conjugate! Instead of <code>double difference = x-y;</code> use <code>double difference = (x*x-y*y)/(x+y);</code></li>
  </ol>
  
  <hr>
  
  <h2>OK, so why haven't you chose a method and moved on?</h2>
  
  <p>I tried all three methods (including straight subtraction) here: <a href="http://ideone.com/wfkEUp" rel="nofollow">http://ideone.com/wfkEUp</a> . I would like to know two things:</p>
  
  <ol>
  <li>Between the "scaling and descaling" method (for which I intentionally chose a power of two) and the arithmetic conjugate method, which one produces less error (in terms of subtracting the large numbers)?</li>
  <li>Which method is computationally more efficient? <code>/*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()*/</code></li>
  </ol>
  
  <p>Any and all help would be welcomed!!</p>
  
  <p>P.S.: All three methods seem to return the same <code>double</code> in the sample code...</p>
  
  <hr>
  
  <p><strong>Let's see some of your code</strong>
  No problem; here is my Matrix.cpp code</p>
  
  <pre><code>#include "ExceptionType.h"
  #include "Matrix.h"
  #include "MatrixArithmeticException.h"
  #include &lt;iomanip&gt;
  #include &lt;iostream&gt;
  #include &lt;vector&gt;
  
  Matrix::Matrix()
  {
      //default size for Matrix is 1 row and 1 column, whose entry is 0
      std::vector&lt;long double&gt; rowVector(1,0);
      this-&gt;matrixData.assign(1, rowVector);
  }
  
  Matrix::Matrix(const std::vector&lt;std::vector&lt;long double&gt; &gt;&amp; data)
  {
      this-&gt;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 &lt; getRowCount(); i++)
      {
          if (largestSize &lt; matrixData[i].size())
              largestSize = matrixData[i].size();
      }
      //make sure that all rowVectors are of that dimension
      for (unsigned i = 0; i &lt; getRowCount(); i++)
      {
          //if we find a rowVector where this isn't the case
          if (matrixData[i].size() &lt; 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&amp; B)
  {
      //if the sizes coincide
      if ((getRowCount() == B.getRowCount()) &amp;&amp; (getColumnCount() == B.getColumnCount()))
      {
          //declare the matrixData
          std::vector&lt;std::vector&lt;long double&gt; &gt; summedData = B.matrixData;    //since we are in the scope of the Matrix, we can access private data members
          for (unsigned i = 0; i &lt; getRowCount(); i++)
          {
              for (unsigned j = 0; j &lt; 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&amp; B)
  {
      //declare negativeB
      Matrix negativeB = B;
      //negate all entries
      for (unsigned i = 0; i &lt; negativeB.getRowCount(); i++)
      {
          for (unsigned j = 0; j &lt; negativeB.getColumnCount(); j++)
          {
              negativeB.matrixData[i][j] = 0-negativeB.matrixData[i][j];
          }
      }
      //simply add the negativeB
      try
      {
          return ((*this)+negativeB);
      }
      catch (MatrixArithmeticException&amp; mistake)
      {
          //should exit or do something similar
          std::cout &lt;&lt; mistake.what() &lt;&lt; std::endl;
      }
  }
  
  Matrix Matrix::operator*(const Matrix&amp; 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&lt;long double&gt; zeroVector(B.getColumnCount(), 0);
          std::vector&lt;std::vector&lt;long double&gt; &gt; data(getRowCount(), zeroVector);
          for (unsigned i = 0; i &lt; getRowCount(); i++)
          {
              for (unsigned j = 0; j &lt; B.getColumnCount(); j++)
              {
                  long double sum = 0; //set sum to zero
                  for (unsigned k = 0; k &lt; 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&amp; operator&lt;&lt;(std::ostream&amp; outputStream, const Matrix&amp; theMatrix)
  {
      //Here, you should use the &lt;&lt; again, just like you would for ANYTHING ELSE.
      //first, print a newline
      outputStream &lt;&lt; "\n";
      //setting precision (optional)
      outputStream.precision(11);
      for (unsigned i = 0; i &lt; theMatrix.getRowCount(); i++)
      {
          //print '['
          outputStream &lt;&lt; "[";
          //format stream(optional)
          for (unsigned j = 0; j &lt; theMatrix.getColumnCount(); j++)
          {
              //print numbers
              outputStream &lt;&lt; std::setw(17) &lt;&lt; theMatrix.matrixData[i][j];
              //print ", "
              if (j &lt; theMatrix.getColumnCount() - 1)
                  outputStream &lt;&lt; ", ";
          }
          //print ']'
          outputStream &lt;&lt; "]\n";
      }
      return outputStream;
  }
  </code></pre>