Decimal Precision Lost after 15th Digit - Bad PI

Question

Attempting to print more than 15 decimal digits of PI result in incorrect decimals printing after the 15th decimal. This despite the 30 correct decimal values being assigned and despite using long double to hold the value. The following test case clearly shows the error.

This was unexpected. If there would be any error in the digits, I would not expect to see any error until the 25th digit after exhausting the IEEE-754 significand. What rule is in play here that might explain by I can't print back the same 30 digits I just assigned to sPI below. This also affects the ability to print the representations of M_PI contained in math.h.

#include <stdio.h>

int
main (void) {

    // static PI approximation (glibc man 1.17)
    long double sPI = 3.14159265358979323846264338327;
    char strPI[]   = "3.14159265358979323846264338327";

    printf ("\n %s (strPI - string - correct)\n", strPI);
    printf (" %.29Lf (sPI - long double - INCORRECT)\n\n", sPI);

    return (0);
}

output:

3.14159265358979323846264338327 (strPI - string - correct)
3.14159265358979311599796346854 (sPI - long double - INCORRECT)
                 ^^^^^^^^^^^^^^

Presumably, this decimal error would apply to any decimal number with greater than 16 decimal precision. When printed as a string, PI prints fine (obviously), but when printed as a double -- the decimal precision breaks down after the 15th decimal. What is causing that?

Very interesting, as suggested adding an L at the end of the floating point literal did help:

 3.14159265358979323846264338327 (strPI - string - correct)
 3.14159265358979323851280895941 (sPI - long double - INCORRECT)

That provided 3 additional decimals of precision. For clarity, this was run on Linux 3.14.1 kernel, gcc 4.8.2 on an old AMD Phenom X4 9850. (AMD Turion Based Laptop & Intel P4 gives same result)

Attempting with quadmath.h and __float128 type assigned to sPI, the results were the same as for the long double. A few more digits were available, but precision still broke down on the 19th digit:

3.14159265358979323846264338327 (strPI - string - correct)
3.1415926535897932385128089594061862 (sPI - long double - INCORRECT)

Looks like `long double` is still 64-bit double precision. Not all compilers support the 80-bit extended precision and the ones that do may require a compiler switch. — Mysticial, Jun 30 '14 at 08:29
You haven't said what compiler you are using. Many compilation platforms offer 80-bit `long double` (GCC on Linux/x86, Clang on Linux/x86, …) but even more define `long double` as the same type as `double` (Visual Studio). Note: even if `long double` was the 80-bit extended precision floating-point format, you wouldn't get 30 π decimals digits out of it, more like 19-21. — Pascal Cuoq, Jun 30 '14 at 08:37
You realise that assigning a `double` or `long double` constant to anything with more precision will not make you happy, right? Use the `l` suffic for `long double` and parse the constant from a string if you're going past `long double`. — tmyklebu, Jun 30 '14 at 13:57

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

You are not storing the 'long double' value into the variable, but the default double instead. The compiler reads the floating point value, stores it as a default type, and only casts it to long double "afterwards". You can see this when you compare its value to a "normally" assigned double.

To hint the compiler to store the constant as a long double, add the modifier suffix L or l at the end of the floating point constant.

Example:

#include <stdio.h>

int main (void)
{
    // static PI approximation (glibc man 1.17)
    double sPI_d     = 3.14159265358979323846264338327;
    long double sPI  = 3.14159265358979323846264338327;
    long double sPI_L= 3.14159265358979323846264338327L;
    char strPI[]     ="3.14159265358979323846264338327";

    printf ("\n %s (strPI - string - correct)\n", strPI);
    printf (" %.29f (sPI - double)\n", sPI_d);
    printf (" %.29Lf (sPI - long double - INCORRECT)\n", sPI);
    printf (" %.29Lf (sPI - long double - BETTER)\n", sPI_L);

    return 0;
}

Output:

 3.14159265358979323846264338327 (strPI - string - correct)
 3.14159265358979311599796346854 (sPI - double)
 3.14159265358979311599796346854 (sPI - long double - INCORRECT)
 3.14159265358979323851280895941 (sPI - long double - BETTER)

See also What is the precision of long double in C++? -- you cannot expect more than around 18 significant digits for a long double.

Thanks Johnware - I will also investigate `quadmath.h` and `__float128` to see if I can get another place or two with 80 bit floating point. — David C. Rankin, Jun 30 '14 at 09:16
`quadmath.h` seems interesting! Bog standard `math.h` defines lots of digits, but it seems casting `(long double)M_PI` comes "too late", as I still got its double value. — Jongware, Jun 30 '14 at 09:22

David Heffernan · Accepted Answer · 2014-06-30T08:36:22.073

3

Double precision floating point values are represented as a 64 bit value, with finite precision of 15-16 decimal digits. It would therefore appear that long double is implemented as 64 bit double precision on your system.

In case you are not familiar with these issues, I highly recommend David Goldberg's timeless article What Every Computer Scientist Should Know About Floating-Point Arithmetic.

edited Jun 30 '14 at 08:36

answered Jun 30 '14 at 08:28

David Heffernan

601,492
42
1,072
1,490

Thank you David. That was the document I could not get to come back in my search. I was stranded reading addendums to FXT, MPFR, GnuMP, and trying to find the page that explained it, but alas. So we get 15 digits in the real world unless using a multi-precision lib. – David C. Rankin Jun 30 '14 at 08:35
2

Some systems have a `long double` that could be 80 bits, and others 128 bits. And others, like yours it would seem, have 64 bit `long double`. But types like `long double` will always have finite precision. For more precision you need specialised libraries. – David Heffernan Jun 30 '14 at 08:38
That also explains why the only way to output the digits of `PI` in the algorithms that compute it to millions of places is one digit at a time. That also was confusing at first. I had a little snippet of code that could correctly print `PI` to any decimal number desired -- but it did it one at a time. – David C. Rankin Jun 30 '14 at 08:40

Decimal Precision Lost after 15th Digit - Bad PI

2 Answers2

Linked