Dealing with floating point errors in .NET

问题:

I'm working on a scientific computation & visualization project in C#/.NET, and we use doubles to represent all the physical quantities. Since floating-point numbers always include a bit of rounding, we have simple methods to do equality comparisons, such as:

static double EPSILON = 1e-6;

bool ApproxEquals(double d1, double d2) {
    return Math.Abs(d1 - d2) < EPSILON;
}

Pretty standard.

However, we constantly find ourselves having to adjust the magnitude of EPSILON as we encounter situations in which the error of "equal" quantities is greater than we had anticipated. For example, if you multiply 5 large doubles together and then divide 5 times, you lose a lot of accuracy. It's gotten to the point where we can't make EPSILON too much larger or else it's going to give us false positives, but we still get false negatives as well.

In general our approach has been to look for more numerically-stable algorithms to work with, but the program is very computational and there's only so much we've been able to do.

Does anyone have any good strategies for dealing with this problem? I've looked into the Decimal type a bit, but am concerned about performance and I don't really know enough about it to know if it would solve the problem or only obscure it. I would be willing to accept a moderate performance hit (say, 2x) by going to Decimal if it would fix these problems, but performance is definitely a concern and since the code is mostly limited by floating-point arithmetic, I don't think it's an unreasonable concern. I've seen people quoting a 100x difference, which would definitely be unacceptable.

Also, switching to Decimal has other complications, such as general lack of support in the Math library, so we would have to write our own square root function, for example.

Any advice?

EDIT: by the way, the fact that I'm using a constant epsilon (instead of a relative comparison) is not the point of my question. I just put that there as an example, it's not actually a snippit of my code. Changing to a relative comparison wouldn't make a difference for the question, because the problem arises from losing precision when numbers get very big and then small again. For example, I might have a value 1000 and then I do a series of calculations on it that should result in exactly the same number, but due to loss of precision I actually have 1001. If I then go to compare those numbers, it doesn't matter much if I use a relative or absolute comparison (as long as I've defined the comparisons in a way that are meaningful to the problem and scale).

Anyway, as Mitch Wheat suggested, reordering of the algorithms did help with the problems.

回答1:

This is not a problem unique to .NET. The strategy to reduce loss of precision is to re-order calculations so that you multiply large quantities times small quantities and add/subtract similiar sized quantities (without changing the nature of the calculation, obviously).

In your example, rather than multiply 5 large quantities together and then divide by 5 large quantities, re-order to divide each large quantity by one of the divisors, and then multiply these 5 partial results.

Of interest? (if you haven't already read): What Every Computer Scientist Should Know About Floating-Point Arithmetic

回答2:

Your best answer is always better algorithms, of course. But it seems to me that if your values aren't all within a couple of orders of magnitude of 1, the using a fixed epsilon is not a good strategy. What you want to do instead is to insure that the values are equal to within some reasonable precision.

// are values equal to within 12 (or so) digits of precision?
//
bool ApproxEquals(double d1, double d2) {
    return Math.Abs(d1 - d2) < (Math.Abs(d1) * 1e-12);
}

If this were C++, then you could also pull some tricks to compare the mantissa and exponent separately, but I can't think of any way to do this safely in unmanged code.

回答3:

Due to the way real numbers are typically represented, You can do this in C (and probably in unsafe C#):

if (llabs(*(long long)&x - *(long long)&y) <= EPSILON) {
    // Close enough
}

This is obviously non-portable and probably a bad idea, but it has the significant advantage of scale-independence. That is, EPSILON can be some small constant like 1, 10 or 100 (depending on your desired tolerance), and it will handle proportional rounding errors correctly regardless of the exponent.

Disclaimer: This is my own private invention, and has not been vetted by anyone with a clue (like, say, a mathematician with a background in discrete arithmetic).