0

Using long double arithmetic in C++, the number 50,000,056,019,485.52605438232421875 squared yields 2,500,005,601,951,690,788,240,883,712. Meanwhile, the number 50,000,056,019,485.526050567626953125 (which differs from the first number by less than 0.001) squared yields 2,500,005,601,951,690,787,704,012,800, which differs from the first square by almost 1 billion. With the differences highlighted:

50,000,056,019,485.526054382324218750 ^ 2 = 2,500,005,601,951,690,788,240,883,712 50,000,056,019,485.526050567626953125 ^ 2 = 2,500,005,601,951,690,787,704,012,800

Long double's range in my machine goes from ~1e-300 to ~1e+300. I fully understand the inability to represent all the numbers in the range, but I didn't expect such a big difference. Could anybody shed some light on this?

kroc
  • 1
  • 7
    You've got about 18 significant digits of accuracy in all cases. That's exactly what you ought to expect from floating point arithmetic. You have to look at *relative* error rather than *absolute*. – Nate Eldredge Nov 11 '20 at 02:23
  • Floating point is hard: http://coliru.stacked-crooked.com/a/e46d6bc4009361a6 – Mooing Duck Nov 11 '20 at 02:30
  • I find the best way to understand doubles is to assume that literally every operation will randomly add or remove ~0.001% (including storing it in a variable). With this assumption in mind, you'll end up writing much better code. – Mooing Duck Nov 11 '20 at 02:33
  • @NateEldredge Thank you! Could you please elaborate a bit more on why only 18 significant digits for long doubles? – kroc Nov 11 '20 at 02:36
  • 5
    If you have an x86 machine, `long double` has a 64-bit fractional part, so at best, every number can be expected to be correct to within about 1 part in 2^64. 2^64 is about 10^19. So 18 or 19 digits is about right. – Nate Eldredge Nov 11 '20 at 02:39
  • Floating point types basically store values in scientific notation. They only store a limited number of decimals, usually too few to fully represent the desired result. A lot of rounding occurs. – François Andrieux Nov 11 '20 at 04:09
  • 2
    The computed squares differ by .54 billion, not “almost 1 billion,” and the mathematical squares would differ by .38 billion. Let x be the first number, which is about 5•10^13. The second number is x − 3.814697265625•10^−6. Squaring that produces x^2 − 2•3.814697265625•10^−6•x + (3.814697265625•10^−6)^2, which differs from x^2 by about 7.6•10^−6•x, which is about 7.6•10^−6 • 5•10^13 = 38•10^7, or about .38 billion. So the rounding error in this calculation is about .16 billion out of 25•10^26 = .16•10^9 / (25•10^26) = 6.4•10^-20, less than 1 part in 2^63. – Eric Postpischil Nov 11 '20 at 11:43

1 Answers1

1

Math

Given d as a small value compared to a

The expected difference of (a+d)*(a+d) - (a)*(a) is about 2*a*d.

Here a is about 50,000,056,000,000.0, d is about 0.00000381... and 2*a*d is about 381,000,000.0. That is within a factor of two of the difference seen using long double: 536,870,912.0.

"I didn't expect such a big difference." is only "off" less than a factor of 2. Given the squares represent rounded products, the difference of the squares seen here is reasonable.

Some details:

long double

Given the following and long double as 10-byte, a and b are consecutive long double values. They differ by 1 ULP or about 1 part in 264. Both decimal code constants are exactly saved as a, b.

long double b = 50000056019485.526054382324218750L; // 0x1.6bcc5c9f50ec355cp+45
long double a = 50000056019485.526050567626953125L; // 0x1.6bcc5c9f50ec355ap+45
// difference
long double d =              0.000003814697265625L; // 0x0.0000000000000002p+45

The square of a, b are not execrably representable as long double. Instead the reported squares are the rounded values. They differ by 2 ULP.

long double bb = 2500005601951690788240883712.0L; // 0x1.027e98e7c774ede8cp+91
long double aa = 2500005601951690787704012800.0L; // 0x1.027e98e7c774ede88p+91
// difference                       536870912.0 .... 0x0.00000000000000004p+91

When multiplying 2 floating point numbers, the error in the product is expected to be within +/-0.5 ULP. Here the ULP is 5,368,70,912 and indeed the squares are correct to within 0.5*

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256