rzikm / netqd Goto Github PK

I've been trying to implement the David H Bailey paper on Double Double and came upon this library that do just that. I stopped my implementation because I thought my results were wrong. But I discover that my implementation and yours give the same results. So we must have implemented correctly the double double as documented in the article.

However, I fail to understand what I should expect. Roughly speaking I should get about 31 digits (base 10) of precision. So when I do something like this:

3.0 + 0.000000000000000000567

I expect a result more or less

3.000000000000000000567

which is 21 digits of precision. However with your library and my unfinished implementation, I get 3.0. Almost looks like that the library is just truncating at 16 digits precision, which is basically what a double is. Since it doesn't make any sense, there is something I don't understand.

Have you come across something similar? Do you understand how the maths work and why we get such bizarre result? Is it me that is just plain stupid?

I was checking out this library, but I find that it doesn't produce results that are correct to higher precision than just double precision.

First off, here is an implementation for you of the conversion from decimal that doesn't truncate the value to double precision:

public static explicit operator DdReal(decimal value) {
	double a = (double)value;
	double b = (double)(value - (decimal)a);
	return new DdReal(a, b);
}

Using that I can successfully convert a decimal value to DdReal and back to decimal without losing any precision.

I put together this to test the correctness of the arithmetic:

private static Random _rnd = new Random();

private static decimal DecimalRnd() {
	decimal sample = 1m;
	while (sample >= 1) {
		byte[] buf = new byte[8];
		_rnd.NextBytes(buf);
		int a = BitConverter.ToInt32(buf, 0);
		int b = BitConverter.ToInt32(buf, 4);
		int c = _rnd.Next(542101087);
		sample = new Decimal(a, b, c, false, 28);
	}
	return sample;
}

public static void Test() {
	decimal a = DecimalRnd();
	decimal b = DecimalRnd();
	DdReal aa = (DdReal)a;
	DdReal bb = (DdReal)b;
	Console.WriteLine($"Decimal: {a} + {b} = {(a + b)}");
	Console.WriteLine($"DdReal:  {(decimal)aa} + {(decimal)bb} = {(decimal)(aa + bb)}");
	Console.WriteLine($"Decimal: {a} - {b} = {(a - b)}");
	Console.WriteLine($"DdReal:  {(decimal)aa} - {(decimal)bb} = {(decimal)(aa - bb)}");
}

Example output:

Decimal: 0,1814580729426223105261687684 + 0,8262370913565243172775595613 = 1,0076951642991466278037283297
DdReal:  0,1814580729426223105261687684 + 0,8262370913565243172775595613 = 1,0076951642991501004477916328
Decimal: 0,1814580729426223105261687684 - 0,8262370913565243172775595613 = -0,6447790184139020067513907929
DdReal:  0,1814580729426223105261687684 - 0,8262370913565243172775595613 = -0,6447790184139019789958151773

As you see, the result is correct to 14 significant digits, which is what you get with just double arithmetic.

Update:

After trying to port several double-double libraries from different languages, I have come to the conclusion that it can't be done reliably in C#. According to the specification the C# compiler is free to choose between 64 bit and 80 bit precision for double calculations. As the double-double implementations rely on a specific precision, it can't be done in plain C#.

Ref:
https://stackoverflow.com/questions/6683059/is-floating-point-math-consistent-in-c-can-it-be