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

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?

回答1:

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

回答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)

回答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.

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

回答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)

fuzzy implementation argument: from this point, everything depends a bit on the implementation of the inverse trigonometric functions. Sometimes arccos and arcsin are derived from atan2 and atan2 as

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

or more specifically from a numeric perspective:

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

Furthermore, atan2 is part of the x86 Instruction set as FPATAN while the others are not. To this end I would argue the usage of:

π = 4 arctan(1)

over all others.

Note: this is a fuzzy argument. I'm certain there are people with better opinions on this.

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 a 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

回答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.