Exact double precision by correct rounding

Question

Although my question sounds trivial, it really is NOT. Hope you can help me.

I want to implement interval arithmetic in my .NET (C#) project. This means that every number is defined by an lower bound and an upper bound. This is helpfull for problems like

1 / 3 = 0.333333333333333 (15 significant digits)

since you would then have

1 / 3 = [ 0.33333333333333 , 0.333333333333334 ] (14 significant digits each)

, so I now FOR SURE that the right answer lays between those two numbers. Without the interval representation I would already have a rounding error with me (i.e. 0.0000000000000003).

To achieve this I wrote my own Interval type that overloads all standard operators like +-*/, etc. To make this type work correctly I need to be able to round the result of 1 / 3 in two directions. Rounding the result down will give me the lower bound for my interval, rounding the result up will give me the upper bound for my interval.

.NET has the Math.Round(double,int) method which rounds the double to int decimal places. Looks great but it can't be forced to round up/down. Math.Round(1.0/3.0,14) would round down, but the also needed up-rounding to 0.33...34 can't be achieved like this.

But there are Math.Ceil and Math.Floor you might say! Okay, those methods round to the next lower or upper integer. So if I want to round to 14 decimal places I first need to reform my result:

1 / 3 = 0.333333333333333 -> *E14 -> 33333333333333.3

So now I can call Math.Ceil and Math.Floor and get both rounded results after reforming back

33333333333333 & 33333333333334 -> /E14 -> 0.33333333333333 & 0.33333333333334

Looks great, but: Let's say my number goes near the double.MaxValue. I can't just *E14 a value near double.MaxValue since this will give me an OverflowException. So this is no solution either.

And, to top all of these facts: All this fails even harder when trying to round 0.9999999999999999999999999 (more than 15 digits) since the internal representation is already rounded to 1 before I can even start trying to round down.

I could try to somehow parse a string containing the double but this won't help since (1/3 * 3).ToString() will already print 1 instead of 0.99...9.

Decimal does not work either since I don't want that deep precision, 14 digits are enough; but I still want that double range!

In C++, where several interval arithmetic implementations exist, this problem could be solved by telling the processor dynamically to swith its roundmode to for example "always down" or "always up". I couldn't find any way to do this in .NET.

So, do you have any ideas? Thanks in advance!

I probably misunderstood your question, but can't you simply add `1/(Math.Pow(10, nDecimalPlaces))` to the `Math.Round` result? — Matteo Umili, Jan 25 '17 at 09:01
This fails when a result is 0.999999999999999 since C# rounds it to 1. I would then get [ 0.99999999999999 , 1.00000000000001 ] which is okay, but not that nice since I am losing precision on the upper bound (1 would be enough as upper bound). — selmaohneh, Jan 25 '17 at 09:04
I am not sure this is how one does [interval arithmetic](https://en.wikipedia.org/wiki/Interval_arithmetic). But when you want to round up or round down, why round to 14 _decimal_ figures? It does not make much sense for a binary type like `double`. Instead you could round to say 46 significant bits. That seems more "uniform". _Edit:_ Do not forget that a number like `0.333333333333334` (14 decimals) is not precisely representable in a binary type like `double`. If instead you use a binary "cut-off", like 46 bits or whatever, the endpoints of your intervals would be "round" in a binary sense. — Jeppe Stig Nielsen, Jan 25 '17 at 09:06
14 comes from the fact that double is pricise to up to 15 significant digits. Because the 15th might be rounded, I wanted to make sure that this digit does not count and so only take 14. — selmaohneh, Jan 25 '17 at 09:09
@JeppeStigNielsen do you have an exmaple how this would look/work with bits? I mean the rounding of 1/3 for exmaple? — selmaohneh, Jan 25 '17 at 09:11
Not sure if this works, but this question shows how to call `_control87` (which can be used to change the rounding mode) from C#: http://stackoverflow.com/questions/18811466/how-do-i-force-the-c-sharp-compiler-to-throw-an-exception-when-any-math-operatio — Simon Byrne, Jan 25 '17 at 09:11
For a `double`, the bits for `1.0/3.0` would look like this: `0 01111111101 0101 01010101 01010101 01010101 01010101 01010101 01010101` There are 64 bits. The first bit is the sign (positive number). The next eleven bits are the magnitude. Then comes fifty-two bits giving the actual digits in binary (mantissa); an initial 1 is implicit/omitted. This 1-0-1-0 pattern should go on for ever to get exactly 1/3. — Jeppe Stig Nielsen, Jan 25 '17 at 09:24
@selmaohneh Then you can check if `Math.Round` result is greater than the original number and subtract instead of adding... — Matteo Umili, Jan 25 '17 at 09:27
@PatriciaShanahan http://stackoverflow.com/questions/15330644/get-next-smallest-double-number — David Heffernan, Jan 25 '17 at 15:14

score 1 · Answer 1 · edited May 23 '17 at 12:24

Assume nextDown(x) is a function that returns the largest double that is less than x, and nextUp(x) is a function that returns the smallest double that is greater than x. See Get next smallest Double number for implementation ideas.

Where you would have rounded a lower bound result down, instead use the nextDown of the round-to-nearest result. Where you would have rounded an upper bound up, use the nextUp of the round-to-nearest result.

This method ensures the interval continues to contain the exact real number result. It introduces extra rounding error - in some cases the lower bound will be one ULP smaller than it should be, and/or the upper bound will be one ULP bigger. However, it is a minimal widening of the interval, much less widening than you would get working in decimal or by suppressing low significance bits.

score 0 · Answer 2 · answered Jan 25 '17 at 13:15

This might be more like a long comment than a real answer.

This code returns an "interval" (I just use Tuple<,>, you can use your own Interval type) based on truncating the seven least significant bits:

static Tuple<double, double> GetMinMaxIntervalBasedOnBinaryNumbersThatAreRoundOnLastSevenBits(double number)
{
  if (double.IsInfinity(number) || double.IsNaN(number))
    return Tuple.Create(number, number); // maybe treat this case differently

  var i = BitConverter.DoubleToInt64Bits(number);

  const int numberOfBitsToClear = 7; // your seven, can change this value, must be below 52
  const long precision = 1L << numberOfBitsToClear;
  const long bitMask = ~(precision - 1L);

  //truncate i
  i &= bitMask;

  return Tuple.Create(BitConverter.Int64BitsToDouble(i), BitConverter.Int64BitsToDouble(i + precision));
}

Disclaimer: I am not sure if this is useful for any purpose. In particular not sure it is useful for interval arithmetic.

With this code, GetMinMaxIntervalBasedOnBinaryNumbersThatAreRoundOnLastSevenBits(1.0 / 3.0) returns the tuple (0.333333333333329, 0.333333333333336).

This code, just like the code you ask for in your question, has the obvious "issue" that if the original value is close to (or even equal to) one of the "round" numbers we use, then the returned interval is "skewed", with the original number being close to one of the ends of the interval. For example, with input 42.0 (already round), you get out the tuple (42, 42.0000000000009).

One good thing about this code is I expect it to be extremely fast.

Thanks for the effort! Was playing around with BitConverter for a while, too. Your method fails for big inputs. 999999999999 returns 999999999999.016 for exmaples. An error 0.016 is pretty big already... — selmaohneh, Jan 25 '17 at 14:17
@selmaohneh That is what I want! No, `0.016` is ___small___ in that context! You say in your question that you do not want to multiply by `1E14` because it could overflow to infinity. That means you want to handle huge numbers. For example, you want to include the number `3.14E298`. My method will work fine with that. What do you think the error is with `3.14E298`? — Jeppe Stig Nielsen, Jan 25 '17 at 14:31

Exact double precision by correct rounding

2 Answers2