Double comparison (again)

Question

I have a "golden source" algorithm that is completely accurate, i.e. the output is accurate to double precision.

I want to compare another algorithm and determine whether it is also accurate to double precision.

For example, the following three numbers might be output:

A = 8.5534733167898555463e-05
X = 8.5534733167898640989e-05
Y = 8.553473316788089652e-05
sig:1 23456789012345

A is the golden source. I can see that A and X are equal to 15 significant digits, whilst Y differs in the 13th significant digit. Hence, Y is not equal to the others to double precision, whilst A and X might be equal to double precision.

I saw this question, but I'm not sure exactly how to apply it. If I use n=1, it signals X to be equal, and Y not to be, which seems reasonable but is it correct? I was hoping it would work with n=0, but that signals both X and Y to be unequal (maybe that's correct).

Answer 1

I really think the key question is: What does it mean to be accurate to double precision?

And I think there's a simple answer: If r is the real answer, and d is the double representation, then d-r<=d'-r for all other possible doubles, d'.

This means that my 2^nd algorithm can only be accurate to double precision if its output is equal to the output from the golden source or the real solution lies exactly half way between the algorithm outputs.

---

My 2^nd algorithm is not that accurate.

Now the question evolves to: What accuracy do I require? I have established in the exchange with Debasish that absolute differences are not a suitable metric. I suspect that relative differences are better but not quite what I want.

I settled on the following algorithm to report the "difference", which considers equality to 15 significant figures to be equality:

double doubles_differ(double a, double b) {
  const int significant_figures = 15;

  // This test will catch most cases
  if (std::abs(a - b) < pow(.1, significant_figures)
                        * std::max(std::abs(a), std::abs(b)))
    return 0;

  // Because we are at the edge of double precision, sometimes a case that is
  // actually equal slips through the above test. The following should be more
  // robust, but a lot slower.
  std::stringstream ss_a;
  std::stringstream ss_b;
  ss_a << std::setprecision(significant_figures - 1) << std::scientific << a;
  ss_b << std::setprecision(significant_figures - 1) << std::scientific << b;
  std::string s_a = ss_a.str();
  std::string s_b = ss_b.str();
  if (s_a == s_b)
    return 0;

  // Finally, return the difference scaled to unity.
  // Format: "7.70612131004268e-013" = significand, followed by "e±nnn".
  if (s_a.substr(s_a.length() - 5, 5) != s_b.substr(s_b.length() - 5, 5))
    throw std::runtime_error("Big diff");

  std::stringstream ss_convert_a(s_a.substr(0, s_a.length() - 6));
  std::stringstream ss_convert_b(s_b.substr(0, s_b.length() - 6));
  double a_converted;
  double b_converted;
  if (!(ss_convert_a >> a_converted))
    throw std::runtime_error("Could not convert a");
  if (!(ss_convert_b >> b_converted))
    throw std::runtime_error("Could not convert b");

  return std::abs(a_converted - b_converted);
}

Finally, some sample data:

a     = 9.9399367132570751e-016
b     = 9.9399367132633209e-016
index   1 2345678901234567
diff  = 6.3007377093526884e-012
relative diff = 6.2835472672107563e-013

Note that the exponent on the diff indicates the number of significant figures of accuracy (in this case, 12). I'm not sure what (if anything) the relative difference indicates.

There should be a better way to do this without strings...

Answer 2

Double precision numbers are stored internally in mantissa and exponent form. Equality checkin double precision floating point number is of no use, as they may be very very close even not exactly same. So, you need to define a threshold. For example, you define epsilon as say 0.000000001 or so (precision depending on your accuracy tolerance). Then, if a, b are two double numbers, check if abs(a-b) < epsilon to accept as equality or almost near equality.