2

I tried to print out the value of pi as defined in the library math.h which on my system is

# define M_PI 3.14159265358979323846 /* pi */

with this simple code:

#include <stdio.h>
#include <math.h>

int main() {

    long double a = M_PI;

    printf("%.20Lf\n", a);

    return 0;
}

what I get as output is

3.14159265358979311600

why is it different from that defined in the math.h library?

Is there a way to reproduce the one defined in the library?

My manual says the type long double extends the precision up to 24 digits. Shouldn't be enough to represent that number?

  • Got the exact same results from three online compilers - [onlinegdb (gcc 5.4.1 c99)](https://onlinegdb.com/rJ8rHyNs8), [repl.it (clang 7)](https://repl.it/repls/ThirstyUtterBooleanalgebra), and [ideone (clang 8.0)](https://ideone.com/TjrViT). Surprisingly, the latter required a standalone definition of `M_PI` - apparently their `math.h` didn't define it. Regardless - same results (3.14159265358979311600) from each of them. – Bob Jarvis - Слава Україні May 21 '20 at 11:42
  • the fact that float is accurate up to 12 decimal positions means that even if you can represent 1.7*10^-38 - 1.7*10^+38 and do calculations with these numbers, when you do try to represent such a number it will be truncated at the 12th position? in this case (for PI) there should be no difference between double and long double@EricPostpischil – Michaelangelo Meucci May 21 '20 at 11:52
  • MichaelangeloMeucci, "accurate up to 12 decimal positions" is an oversimplification (and questionable). The precision is a binary one. Example: between [1.0...2.0), the least significant binary digit is 2^-23 or 0.000000119... . Between [2.0...4.0) the least significant binary digit is 2^-22 or 0.000000238... . – chux - Reinstate Monica May 21 '20 at 12:07
  • 1
    @EricPostpischil: The `00` digits are correct; not a fabrication: the decimal expansion for `M_PI` to just beyond that point is `3.14159265358979311599796...` – Mark Dickinson May 21 '20 at 14:22
  • @MarkDickinson: Indeed, I played the odds on a 1% chance with real numbers versus a 100% chance with bad formatters and lost. – Eric Postpischil May 21 '20 at 15:00
  • `long double` does not **extend** precision. It may **have** more precision than `double`, but it cannot make new information that is not already in a `double` you assign to it. `M_PI` is a `double` constant, and the value will not change when you assign it to a `long double`. – Eric Postpischil May 21 '20 at 15:01
  • @MichaelangeloMeucci: Re “when you do try to represent such a number it will be truncated at the 12th position”: No. First, terminology. A `double` or `long double` object *represents* a certain real number. When that number is, for example, 2^1024−2^971, it is **exactly** that number. When you *convert* that number to decimal, the result is often rounded, depending on the specification, quality, and parameters of the formatter. There is no mathematical reason 2^1024−2^971 cannot be converted exactly to a decimal numeral. We just rarely do so. – Eric Postpischil May 21 '20 at 15:04
  • @MichaelangeloMeucci: Notions about how many decimal digits a binary floating-point value is good to are merely sloppy approximations of how much information they contain. As I answer [here](https://stackoverflow.com/questions/61609276/how-to-calculate-float-type-precision-and-does-it-make-sense/61614323#61614323), these descriptions of binary floating-point in terms of decimal digits are largely nonsensical. – Eric Postpischil May 21 '20 at 15:06

3 Answers3

2

The best precision for the mathematical constant π (pi) as provided by the implementation can be queried in a standard manner by calling acosl(-1);:

const long double pi = acosl(-1.0L);
printf("%.20Lf\n", pi);

Since this approach has the additional overhead of performing a computation (whereas your approach uses a compile-time constant), it is recommended to place the result in a constant variable and use that for the entire program, performing the calculation only once. The advantage is that this is portable. This will compute π with the best available precision on the target platform and requires zero change in the code when switching platforms.

Conforming implementations are forbidden by the Standard to define the constant M_PI, because implementation-defined reserved names must start with _ followed by a capital letter or another _. If you want to use it, you have to define it yourself.

DarkAtom
  • 2,589
  • 1
  • 11
  • 27
  • Why do computation when it is a constant? – Ed Heal May 21 '20 at 11:36
  • 1
    Because it is standard and portable. Different implementations support different precisions. You only have to do the calculation once and it is available for the entire program. – DarkAtom May 21 '20 at 11:39
  • Implementations are not completely forbidden from defining `M_PI`. A typical way this is handled is that a program will define a special symbol, such as, `#define _POSIX_C_SOURCE 200809L`, to request that it wants special features. It does this before including the headers. This causes the standard headers and others to modify their behavior, providing the requested features. This conforms because the standard is silent about what happens if such symbols are defined (it says the behavior is “undefined,” meaning it does not say what happens), and it does not affect strictly conforming programs. – Eric Postpischil May 21 '20 at 15:11
  • @EricPostpischil That is true, but the OP is not actually doing that, meaning that his implementation is doing something forbidden by the Standard. Of course, this is a conventional macro so nobody defines it in their programs, so it is not really a big deal. – DarkAtom May 21 '20 at 19:36
2

why is it different from that defined in the math.h library?

double constant

3.14159265358979323846 is a double. To retain long double precision, use an L. That will get us closer ...
See also @Eric Postpischil

// # define M_PI        3.14159265358979323846  /* pi */
# define M_PI_L      3.14159265358979323846L  /* pi */
//                                         ^

Is there a way to reproduce the one defined in the library?

Not quite.

Binary nature of floating point

OP's long double cannot represent every decimal encoded value exactly. Instead a binary one is used of the form: some_int*2exponent. On my machine the closest to 3.14159265358979323846 is 3.14159265358979323851...

int main() {
    //                     3.1415926535897932384626433832795‬...
    # define M_PI_L        3.14159265358979323846L  /* pi */
    long double a = M_PI_L;

    printf("%.21Lf epsilon\n", LDBL_EPSILON);
    printf("%.21Lf local epsilon\n", a - nextafterl(a, 0.0));
    printf("%.21Lf before\n", nextafterl(a, 0.0));
    printf("3.14159265358979323846L source code\n");
    printf("%.21Lf M_PI_L\n", a);
    printf("%.21Lf after\n", nextafterl(a, a*2));
}

Sample output (spaces added for clarity)

0.000000000000000000 108 epsilon
0.000000000000000000 217 local epsilon
3.141592653589793238 296 before
3.141592653589793238 46L source code
3.141592653589793238 513 M_PI_L
3.141592653589793238 730 after

My manual says the type long double extends the precision up to 24 digits. Shouldn't be enough to represent that number?

"... up to 24 digits" may be misleading. Recall floating point numbers are overwhelmingly encoded as an integer times a power of 2. Thus the precision is a binary one. Depending on the region of FP values, the corresponding decimal precision wobbles between LDBL_DIG and LDBL_DECIMAL_DIG which is 18, 21 on my machine. I would expect the same with OP.


@DarkAtom idea is good. To typically get the machine's best machine pi :

// Perform this once
const long double pi = acosl(-1.0L);

An alternative is to specify with lots of digits, more than expected on any target.

// 50 digit pi         1 2345678901234567890123456789012345678901234567890 
const long double pi = 3.1415926535897932384626433832795028841971693993751L;
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • on my system is called `# define M_PIl 3.141592653589793238462643383279502884L /* pi */` but if I try to use it with `long double a = M_PIl;` I get error: ‘M_PIl’ undeclared (first use in this function); did you mean ‘M_PI’? – Michaelangelo Meucci May 21 '20 at 11:27
  • @MichaelangeloMeucci Yes, the value of _machine_pi_ varies. Coding it in source code with lots of digits is usually a _good thing_, yet missing the 'L' has significant effects - which yours wisely does not miss. – chux - Reinstate Monica May 21 '20 at 11:31
  • @MichaelangeloMeucci Re: "I get error". `M_PI, M_PIl` should _not_ be defined in `` per the C spec, although often is there. Where did you find it? Access to it depends on compiler options. With strict spec options, I'd would expect "‘M_PIl’ undeclared...". – chux - Reinstate Monica May 21 '20 at 11:34
  • I did cat /usr/include/math.h | grep M_PI – Michaelangelo Meucci May 21 '20 at 11:37
  • @MichaelangeloMeucci See [Using M_PI with C89 standard](https://stackoverflow.com/questions/5007925/using-m-pi-with-c89-standard). Applies to C99, C11, ... too – chux - Reinstate Monica May 21 '20 at 11:47
1

Is there a way to reproduce the one defined in the library?

Since the value of π is not likely to change in the near future, it is okay to hardcode it. To do this, we may prefer to use C’s hexadecimal floating-point format because some C implementations are not good at converting decimal values to binary floating-point, and it is not clear we can rely on atanl to be as accurate as possible. Here is π to 132 bits after the radix point:

0x3.243f6a8885a308d313198a2e03707344ap0L
Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
  • Good point about some libraries and compilers having implementation/conversion weaknesses. I would expect fewer difficulties with hex. Re: "π is not likely to change" [should the value of pi change](https://medium.com/@patkerr/should-the-value-of-pi-change-f18f41599919) – chux - Reinstate Monica May 21 '20 at 15:28