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!

Answer 1

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.

Answer 2

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.