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?

Solution 1

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

Solution 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)

Solution 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.

Solution 4

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).

Solution 5

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

View more solutions