Math constant PI value in C

Question

Calculating PI value is one of the complex problem and wikipedia talks about the approximations done for it and says it's difficult to calculate PI accurately.

How does C calculate PI? Does it calculate it every time or is it using a less accurate fixed value?

Airy

Using pi calculated out to only 39 decimal places would allow one to compute the circumference of the entire universe to the accuracy of less than the diameter of a hydrogen atom. 16 decimal places (approximately what you get with a double) should be enough to calculate the diameter of the Solar System with the error less than a hair-width. – Bondie

I have found this online and find it to be pretty interesting: crypto.stanford.edu/pbc/notes/pi/code.html – Offprint

Answer 1

In C Pi is defined in math.h: #define M_PI 3.14159265358979323846

Senegambia

And that value is the most accurate representation available to a double precision value. – Hervey

Not quite -- in fact a conforming C implementation may not define PI in <math.h>. POSIX specifies M_PI, but again, a conforming C implementation may not define it. (POSIX imposes some requirements that conflict with the C standard.) But you can define it that way in your own program. – Bittern

If you are unlucky with your math header, use GP's #define. – Hervey

so it's fixed value and there's no higher accuracy possible? – Airy

@Airy Not if you are dealing with the built in types no. There are libraries though dealing with larger precision numbers. – Hervey

Additional info: if you use M_PI and get error that it is undefined - this can be solved by #define _USE_MATH_DEFINES – Comment

GCC 4.8.1 (MinGW) apparently does not define M_PI or PI. – Idealistic

What does the M in M_PI stand for? – Enciso

@Enciso M stands for "math". Back in the days, all "math constants" were prefixed with M_. There was also stuff like M_E, M_LN10 etc. They never made it to the standard. – Szeged

In MSVC it requires defining _USE_MATH_DEFINES. – Bellina

Does it make sense to update the answer with this? ``` #include <math.h> #ifndef M_PI #define M_PI 3.14159265358979323846 #endif ``` – Mcelrath

@Szeged Any ideas / guesses about why they never made it to the standard? – Grubby

@KeithThompson POSIX imposes some requirements that conflict with the C standard: where is to read more about such conflicting requirements? – Grubby

@Grubby POSIX is here: pubs.opengroup.org/onlinepubs/9699919799/toc.htm I don't have a more specific reference at the moment. – Bittern

@KeithThompson About POSIX: and what about reserved identifiers? As I understand, if POSIX wants to conform to the C standard, then it shall use only reserved identifiers (i.e. accept, etc. need to be available to the user). How can you comment on that? (Here is the relevant question.) – Grubby

Answer 2

The closest thing C does to "computing π" in a way that's directly visible to applications is acos(-1) or similar. This is almost always done with polynomial/rational approximations for the function being computed (either in C, or by the FPU microcode).

However, an interesting issue is that computing the trigonometric functions (sin, cos, and tan) requires reduction of their argument modulo 2π. Since 2π is not a diadic rational (and not even rational), it cannot be represented in any floating point type, and thus using any approximation of the value will result in catastrophic error accumulation for large arguments (e.g. if x is 1e12, and 2*M_PI differs from 2π by ε, then fmod(x,2*M_PI) differs from the correct value of 2π by up to 1e12*ε/π times the correct value of x mod 2π. That is to say, it's completely meaningless.

A correct implementation of C's standard math library simply has a gigantic very-high-precision representation of π hard coded in its source to deal with the issue of correct argument reduction (and uses some fancy tricks to make it not-quite-so-gigantic). This is how most/all C versions of the sin/cos/tan functions work. However, certain implementations (like glibc) are known to use assembly implementations on some cpus (like x86) and don't perform correct argument reduction, leading to completely nonsensical outputs. (Incidentally, the incorrect asm usually runs about the same speed as the correct C code for small arguments.)

Anima

Do you have any references for your last paragraph? From my intuition, when you compute sin(1e16*M_PI), the argument's mantissa is completely to the left of the decimal point, and there are no fractional digits. Thus, it shouldn't matter how you do your argument reduction, because you will never recover any information that wasn't there in the first place (namely the fractional part of the mantissa). – Dd

@Fritz: The factional part of the input argument is zero -- the argument is an integer (a very large one). Argument reduction does not "recover information". Given a large x, it gives you a value y for which (mathematically) |sin(x)-sin(y)|<epsilon where epsilon is something like 1ulp, but with y in the domain of the (usually polynomial) function you use to approximate sin. – Anima

After argument reduction of a large argument, the fractional part most certainly is not zero, because pi is not rational (though technically it could happen to be zero after rounding, for very rare cases). – Anima

I was looking for the ISO C, Annex F, reference for trig functions, but apparently there's not a hard requirement for them to be the IEEE trig functions (which would require them to be correct, and correctly rounded, for all inputs). So it may be more of a QoI (quality of implementation) issue than a correctness one. But implementations based on fdlibm have argument reduction that works across the whole domain of representable values. – Anima

Okay, I understand what you are saying and thanks for the pointer to fdlibm. My line of thought was more that calculating the sine of such a large x is almost meaningless and most likely a bug in the code that is using the sine function. If x is so large that it becomes integer, then you can represent only around 6 values of x per sine period. Of course one can find a y which satisfies the condition, and that is nice, but I doubt that there is any useful application for calculating sines of integers. – Dd

More interesting than the extreme case of integer x would be an analysis of how quickly this "epsilon error" increases with x. If there is a significant error at 1000*M_PI, then I understand the problem and would wholeheartedly agree with you. – Dd

@Fritz: According to the rough analysis in my answer, the error of fmod(x,2*M_PI) vs the correct argument-reduced value grows linearly in the magnitude of x. So assuming 1ulp for x just outside the unit circle, you'd be looking at something like 1000ulp for 1000*M_PI. Near the zeros sin is close to linear with slope 1, so the # of ulp in the argument translate directly to ulp in the result. – Anima

And BTW there are interesting applications of sin(2^n) for large values of n; I don't recall what they are right off, though. – Anima

Answer 3

Just define:

#define M_PI acos(-1.0)

It should give you exact PI number that math functions are working with. So if they change PI value they are working with in tangent or cosine or sine, then your program should be always up-to-dated ;)

Tasimeter

This would evaluate quite costly fuction (acos) every time the constant is used. Hardly an effective and rational approach. – Rarely

I made this solution for the best acurate solution using #define. If someone has a problem with it then they can define a global variable: double M_PI = acos(-1.0); – Horsewhip

@JaromírAdamec Actually, any good compiler should optimize that expression into a constant, being acos a pure function from the standard library called with constant arguments; unless you compile with optimizations off, that is. – Kendall

@UnrealEagle It might be so, but it's quite controversial, as - strictly speaking - the acos() function (as any other) is not member of language proper, but the standard library.Therefore, the compiler shuld not "know" the semantics of the function. Of course, modern compilers use "instrinsic functions", but it's a bit of hacky. – Rarely

@JaromírAdamec on the contrary, the *C standard explicitly allows - or even endorses - the compiler knowing intricate details of all standard library functions! – Talie

@AnttiHaapala This might or might not be true. Concerning "modern compiler", there is usually a switch to enable using intrinsic functions (e.g. in Vistual C++ project properties: C/C++|Optimization|Enable intrinsic functions). If you happen to switch off this option, you will end up with a regular external reference to the function resolved by linker. It is also matter of question, which standard library function the compiler replaces as "intrinsic". Your opinion that "all of them" is a bit optimistic. For sure, such code is not very portable, especially not to some minor platforms. – Rarely

@JaromírAdamec I am not saying anything about the quality of a single implementation. I am just stating that a C compiler that conforms to the C specification as described in ISO 9899 is explicitly allowed to replace/inline/substitute all of the references to standard library as described in the same standard assuming that a function by the same name will have the behaviour equal to that defined in the standard. The switches that exist in various compilers exists only to make non-strictly conforming programs compile correctly in these modern implementations. – Talie

The acos function is a part of the C language proper. – Sind

Answer 4

anyway you have not a unlimited accuracy so C define a constant in this way:

#define PI 3.14159265358979323846

import math.h to use this

Winged

That means I can manage to get a higher accuracy for PI than defined by C? If i use double x = 1345*PI; But it is limited by the precision of the double variable used in the program isn't??? That means newly defined accurate value is useless?? – Airy

partially yes:D... "double x = 1345*PI" by this work you loss some accuracy, because PI is accurate as it can. if you want more accuracy you should implement YOUR OWN structure and store the result in an array (Like BigInteger in java). OK? – Winged

Actually M_PI is not in C. It's part of the XSI extension option in POSIX. – Anima

@R.., are you sure it's XSI-only? The docs don't have the normal [XSI] marker they have on XSI-only stuff (see e.g. who). – Sardonic

Yes, see here: pubs.opengroup.org/onlinepubs/9699919799/basedefs/math.h.html All of the M_* constants are [XSI> tagged. – Anima

Ah, you were looking at SUSv2, which is not POSIX. This was before the UNIX/POSIX specification merger, at which time the whole SUS was XSI. Beginning with SUSv3, the Unix and POSIX standards merged, with the XSI option being the parts of Unix that were not deemed sufficiently useful or universal to mandate for all POSIX systems to support. Later in SUSv4 (POSIX 2008), a number of the more-useful XSI options were moved to the base (POSIX) standard, and the less-useful ones are gradually being marked obsolescent, so eventually with some luck the XSI option will cease to exist... – Anima

Answer 5

Depending on the library you are using the standard GNU C Predefined Mathematical Constants are here... https://www.gnu.org/software/libc/manual/html_node/Mathematical-Constants.html

You already have them so why redefine them? Your system desktop calculators probably have them and are even more accurate so you could but just be sure you're not conflicting with existing defined ones to save on compile warnings as they tend to get defaults for things like that. Enjoy!

Effulgent

Answer 6

You could make a function that calculates PI, by doing one of multiple possible infinite sum calculations(I wouldn't recommend something like atan() as those functions probably use pi in one form or another themselves). That said I would not recommend manually calculating pi. No one needs pi to be any more accurate than math.h provides with M_PI, in fact, it would probably be best for you to use just the first few digits of it instead as that would allow you to have slightly faster code. Unless you want to calculate the size of the universe down to the centimeter why bother?

If you still want to calculate it here's how.

(Method from 3blue1brown The Wallis product for pi, proved geometrically)

double calculate_pi(int accuracy){
     double result = 1;
     int a = 2;
     int b = 1;

     for(int i = 0; i < accuracy; i ++){
          result = (a/(double)b) * result;
          if(a < b){
               a = a + 2;
          }
          else if(b < a){
               b = b + 2;
          }
     }

     return result * 2;
}

December

Answer 7

I don't know exactly how C calculates PI directly as I'm more familiar with C++ than C; however, you could either have a predefined C macro or const such as:

#define PI 3.14159265359.....
const float PI = 3.14159265359.....
const double PI = 3.14159265359.....
/* If your machine,os & compiler supports the long double */
const long double PI = 3.14159265359..... 

or you could calculate it with either of these two formulas:

#define M_PI acos(-1.0);
#define M_PI (4.0 * atan(1.0)); // tan(pi/4) = 1 or acos(-1)

IMHO I'm not 100% certain but I think atan() is cheaper than acos().

Offprint

@RyanHaining That was a pure typo on my part thank you for pointing that out. I made the appropriate corrections. – Offprint

It's still not valid C – Maggoty

@RyanHaining Hmm okay; I know you can do this in C++ I wasn't sure if C permitted it. It's been too many years since I worked in C. Is there an equivalent without using; #define CONSTANT numberHere? – Offprint

"All the expressions in an initializer for an object that has static or thread storage duration shall be constant expressions or string literals" (C2011, 6.7.9/4). It's not clear to me what characteristics you do hope to have in an alternative that does not rely on a macro, but that's one of the main things you're up against. – Stander

@JohnBollinger Makes sense; my C is very rusty as it's been close to 20 years since I've worked in it and the fact that I've mostly used C++. I had originally answered this question with a C++ perspective in mind... I've updated the answer to accommodate for this... – Offprint

@FrancisCugler you don't have to #define in your first example, that one works find as const double since it is a constant expression – Maggoty

@RyanHaining Oh okay I couldn't remember correctly about constness in C. – Offprint