Why define PI = 4*ATAN(1.d0)

Question

What is the motivation for defining PI as

PI=4.D0*DATAN(1.D0)

within Fortran 77 code? I understand how it works, but, what is the reasoning?

Birthplace

As an alternative, I would almost expect to see PI = 3.1415926535... etc instead – Birthplace

I found a most excellent answer to the mathematics for this equation over on the Math Stackoverflow site math.stackexchange.com/questions/1211722/… – Scopoline

If you use PI = 3.1415926535... you must add the data type suffix to get anything other than default real precision. Since you were using f66 double precision, that would be the D0 suffix. – Inulin

Note: modern Fortran doesn't require DATAN(): ATAN() is aliased to its respective single- and double-precision versions, based on the argument. – Industrial

"modern Fortran" you lost me there ;) – Birthplace

Modern Fortran is Fortran >= Fortran 2003 – Bordelon

The feature mentioned here was introduced in Fortran 77. For the last 30 years, the remaining situation where you would use datan is rather contrived and unlikely to be useful. – Inulin

Answer 1

This style ensures that the maximum precision available on ANY architecture is used when assigning a value to PI.

Pitchfork

But be warned if you follow this approach -- not all compilers and underlying maths libraries are equal in the results they return for trig functions of floating-point numbers, especially at critical points of those functions. There has recently been a long discussion of this on comp.lang.fortran where most of the Fortran experts hang out. Their conclusion -- specify the constant by pi = 3.14159...(enough digits for required precision then some for safety). – Mignonne

High Performance Mark: it would be nice if you could give a link to the comp.lang.fortran thread you mentioned! – Industrial

Unfortunately, some languages go to great lengths to ensure that the value of pi returned by this approach is not the same value that is used in computing the period of trig functions. – Kingkingbird

@supercat: Any idea why? Or, can you cite it so I can ask it as a new question? :-) – Industrial

@jvriesem: If someone tries to compute, e.g. sin(31416), the correct mathematical value should be sin(31416 minus (10000 times mathematical pi), and not sin(31416 - 10000 times the nearest floating-point value to mathematical pi). – Kingkingbird

@jvriesem: Personally, I doubt that there are very many situations where someone who calls e.g. sin(x) when x is outside the range +/-pi/2 would not be equally happy with the sine of any value whose closest representation is x. Nonetheless, the Java language runtime adds extra code to its trig function implementations which, when asked to compute sin(x) for x near e.g. 2π will first subtract a precisely-representable value which is close to 2π and then subtract the difference between that value and mathematical 2π. – Kingkingbird

@jvriesem: If one imagines a floating-point system with five decimal digits of precision asked to compute sin(3.1416), it might (noting that π is about 3.1415 + 0.000092653 + 0.00000000058979), subtract 3.1415 from the value (yielding 0.0001), and then subtract 0.000092653 from that, yielding 0.000007347, and then subtract 0.00000000058979 from that, yielding 0.0000073464, of which it would then take the sine, yielding -0.0000073464. – Kingkingbird

@jvriesem: This result would be the five-digit value closest to the mathematical sine of the exact value 3.1416, but if 3.1416 is the result of a computation whose value, if computed perfectly without rounding, could have been anywhere between 3.14155 and 3.14165, then the perfectly-computed sine of a perfectly-computed value could be anywhere between -0.000057346 and 0.000042653, and any value in that range would be essentially as good as any other. – Kingkingbird

Answer 2

Because Fortran does not have a built-in constant for PI. But rather than typing in the number manually and potentially making a mistake or not getting the maximum possible precision on the given implementation, letting the library calculate the result for you guarantees that neither of those downsides happen.

These are equivalent and you'll sometimes see them too:

PI=DACOS(-1.D0)
PI=2.D0*DASIN(1.D0)

Exocarp

Note that in Fortran 77 and later generic names ACOS and ASIN are preferred to the specific names DACOS and DASIN. – Shockheaded

Answer 3

I believe it's because this is the shortest series on pi. That also means it's the MOST ACCURATE.

The Gregory-Leibniz series (4/1 - 4/3 + 4/5 - 4/7...) equals pi.

atan(x) = x^1/1 - x^3/3 + x^5/5 - x^7/7...

So, atan(1) = 1/1 - 1/3 + 1/5 - 1/7 + 1/9... 4 * atan(1) = 4/1 - 4/3 + 4/5 - 4/7 + 4/9...

That equals the Gregory-Leibniz series, and therefore equals pi, approximately 3.1415926535 8979323846 2643383279 5028841971 69399373510.

Another way to use atan and find pi is:

pi = 16*atan(1/5) - 4*atan(1/239), but I think that's more complicated.

I hope this helps!

(To be honest, I think the Gregory-Leibniz series was based on atan, not 4*atan(1) based on the Gregory-Leibniz series. In other words, the REAL proof is:

sin^2 x + cos^2 x = 1 [Theorem] If x = pi/4 radians, sin^2 x = cos^2 x, or sin^2 x = cos^2 x = 1/2.

Then, sin x = cos x = 1/(root 2). tan x (sin x / cos x) = 1, atan x (1 / tan x) = 1.

So if atan(x) = 1, x = pi/4, and atan(1) = pi/4. Finally, 4*atan(1) = pi.)

Please don't load me with comments-I'm still a pre-teen.

Typewriting

I dont understand how you get from atan x (1/tan x) = 1 to atan(x) = 1 – Scopoline

remark that with the formula pi = 16*atan(1/5) - 4*atan(1/239) is numerically not a good choice. Both 1/5 and 1/239 cannot be represented accurately using floating-point numbers. Hence, the atan will already introduce errors due to these floating-point approximations. – Light

Answer 4

There is more to this question than meets the eye. Why 4 arctan(1)? Why not any other representation such as 3 arccos(1/2)?

This will try to find an answer by exclusion.

Mathematical intro: When using inverse trigonometric functions such as arccos, arcsin and arctan, one can easily compute π in various ways:

π = 4 arctan(1) = arccos(-1) = 2 arcsin(1) = 3 arccos(1/2) = 6 arcsin(1/2)
  = 3 arcsin(sqrt(3)/2) = 4 arcsin(sqrt(2)/2) = ...

There exist many other exact algebraic expressions for trigonometric values that could be used here.

Floating-point argument 1: it is well understood that a finite binary floating-point representation cannot represent all real numbers. Some examples of such numbers are 1/3, 0.97, π, sqrt(2), .... To this end, we should exclude any mathematical computation of π where the argument to the inverse trigonometric functions cannot be represented numerically. This leaves us the arguments -1,-1/2,0,1/2 and 1.

π = 4 arctan(1) = 2 arcsin(1)
   = 3 arccos(1/2) = 6 arcsin(1/2)
   = 2 arccos(0)
   = 3/2 arccos(-1/2) = -6 arcsin(-1/2)
   = -4 arctan(-1) = arccos(-1) = -2 arcsin(-1)

Floating-point argument 2: in the binary representation, a number is represented as 0.b_nb_n-1...b₀ x 2^m. If the inverse trigonometric function came up with the best numeric binary approximation for its argument, we do not want to lose precision by multiplication. To this end we should only multiply with powers of 2.

π = 4 arctan(1) = 2 arcsin(1)
  = 2 arccos(0)
  = -4 arctan(-1) = arccos(-1) = -2 arcsin(-1)

Note: this is visible in the IEEE-754 binary64 representation (the most common form of DOUBLE PRECISION or kind=REAL64). There we have

write(*,'(F26.20)') 4.0d0*atan(1.0d0) -> "    3.14159265358979311600"
write(*,'(F26.20)') 3.0d0*acos(0.5d0) -> "    3.14159265358979356009"

This difference is not there in IEEE-754 binary32 (the most common form of REAL or kind=REAL32) and IEEE-754 binary128 (the most common form of kind=REAL128)

Implementation argument: On intel CPU, the atan2 is part of the x86 Instruction set as FPATAN while the other inverse trigonometric function are derived from atan2. A potential derivation could be:

          mathematically         numerically
ACOS(x) = ATAN2(SQRT(1-x*x),1) = ATAN2(SQRT((1+x)*(1-x)),1)
ASIN(x) = ATAN2(1,SQRT(1-x*x)) = ATAN2(1,SQRT((1+x)*(1-x)))

This can be seen in the assembly code of these instructions (See here). To this end I would argue the usage of:

π = 4 arctan(1)

Note: this is a fuzzy argument. I'm certain there are people with better opinions on this.
Interesting reads on FPATAN: How is arctan implemented?, x87 trigonometric instructions

The Fortran argument: why should we approximate π as :

integer, parameter :: sp = selected_real_kind(6, 37)
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: qp = selected_real_kind(33, 4931)

real(kind=sp), parameter :: pi_sp = 4.0_sp*atan2(1.0_sp,1.0_sp)
real(kind=dp), parameter :: pi_dp = 4.0_dp*atan2(1.0_dp,1.0_dp)
real(kind=qp), parameter :: pi_qp = 4.0_qp*atan2(1.0_qp,1.0_qp)

and not :

real(kind=sp), parameter :: pi_sp = 3.14159265358979323846264338327950288_sp
real(kind=dp), parameter :: pi_dp = 3.14159265358979323846264338327950288_dp
real(kind=qp), parameter :: pi_qp = 3.14159265358979323846264338327950288_qp

The answer lays in the Fortran standard. The standard never states that a REAL of any kind should represent an IEEE-754 floating point number. The representation of REAL is processor dependent. This implies that I could inquire selected_real_kind(33, 4931) and expect to obtain a binary128 floating-point number, but I might get a kind returned that represents a floating-point with much higher accuracy. Maybe 100 digits, who knows. In this case, my above string of numbers is to short! One cannot use this just to be sure? Even that file could be too short!

Interesting fact : sin(pi) is never zero

write(*,'(F17.11)') sin(pi_sp) => "   -0.00000008742"
write(*,'(F26.20)') sin(pi_dp) => "    0.00000000000000012246"
write(*,'(F44.38)') sin(pi_qp) => "    0.00000000000000000000000000000000008672"

which is understood as:

pi = 4 ATAN2(1,1) = π + δ
SIN(pi) = SIN(pi - π) = SIN(δ) ≈ δ

---

program print_pi
! use iso_fortran_env, sp=>real32, dp=>real64, qp=>real128

  integer, parameter :: sp = selected_real_kind(6, 37)
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: qp = selected_real_kind(33, 4931)

  real(kind=sp), parameter :: pi_sp = 3.14159265358979323846264338327950288_sp
  real(kind=dp), parameter :: pi_dp = 3.14159265358979323846264338327950288_dp
  real(kind=qp), parameter :: pi_qp = 3.14159265358979323846264338327950288_qp
  
  write(*,'("SP "A17)') "3.14159265358..."
  write(*,'(F17.11)') pi_sp
  write(*,'(F17.11)')        acos(-1.0_sp)
  write(*,'(F17.11)') 2.0_sp*asin( 1.0_sp)
  write(*,'(F17.11)') 4.0_sp*atan2(1.0_sp,1.0_sp)
  write(*,'(F17.11)') 3.0_sp*acos(0.5_sp)
  write(*,'(F17.11)') 6.0_sp*asin(0.5_sp)

  write(*,'("DP "A26)') "3.14159265358979323846..."
  write(*,'(F26.20)') pi_dp
  write(*,'(F26.20)')        acos(-1.0_dp)
  write(*,'(F26.20)') 2.0_dp*asin( 1.0_dp)
  write(*,'(F26.20)') 4.0_dp*atan2(1.0_dp,1.0_dp)
  write(*,'(F26.20)') 3.0_dp*acos(0.5_dp)
  write(*,'(F26.20)') 6.0_dp*asin(0.5_dp)

  write(*,'("QP "A44)') "3.14159265358979323846264338327950288419..."
  write(*,'(F44.38)') pi_qp
  write(*,'(F44.38)')        acos(-1.0_qp)
  write(*,'(F44.38)') 2.0_qp*asin( 1.0_qp)
  write(*,'(F44.38)') 4.0_qp*atan2(1.0_qp,1.0_qp)
  write(*,'(F44.38)') 3.0_qp*acos(0.5_qp)
  write(*,'(F44.38)') 6.0_qp*asin(0.5_qp)

  write(*,'(F17.11)') sin(pi_sp)
  write(*,'(F26.20)') sin(pi_dp)
  write(*,'(F44.38)') sin(pi_qp)

end program print_pi

Light

Along with the argument sin(pi) /= 0 I have always been reluctant to assume that acos(-1d0) would be as accurate as 4*atan(1d0), although it may be, depending on the implementation. – Inulin

I believe there is a proven algorithm for sin() to produce a correctly rounded value for large arguments, although it comes with considerable additional execution time. I haven't seen the reference. Some IBM libraries gave you this by default. Any high quality implementation will involve some simulation of extra precision which avoids significant degradation of accuracy up to arguments as large as could be useful in practice, say +-20 Pi, or the (smaller) point at which the internal m387 implementation changes abruptly from returning an accurate value to returning the argument. – Inulin

By the way, the internal Pi constant in the m387 firmware is advertised as carrying 66 bits precision. Actually, this is nothing special in the implementation; the correctly rounded 66 bit precision value has 2 low order zero bits. A high quality C compiler will give you this value for a long double constant with 21 digits. – Inulin

Answer 5

It's because this is an exact way to compute pi to arbitrary precision. You can simply continue executing the function to get greater and greater precision and stop at any point to have an approximation.

By contrast, specifying pi as a constant provides you with exactly as much precision as was originally given, which may not be appropriate for highly scientific or mathematical applications (as Fortran is frequently used with).

Allister

Answer 6

That sounds an awful lot like a work-around for a compiler bug. Or it could be that this particular program depends on that identity being exact, and so the programmer made it guaranteed.

Degradation

It's actually a very common way of setting the value of PI -- not only in Fortran, but in other languages as well. (See above comments.) – Industrial