为什么定义 PI = 4*ATAN(1.d0)

Question

将 PI 定义为的动机是什么

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

在 Fortran 77 代码中？我了解它的工作原理，但是，原因是什么？

Answer 1

这种风格确保在为 PI 赋值时使用任何架构上可用的最大精度。

Answer 2

因为 Fortran 没有 PI 的内置常量。但是，与其手动输入数字并可能会出错或无法在给定的实现上获得最大可能的精度，不如让库为您计算结果，这样可以保证不会发生这些缺点。

这些是等效的，您有时也会看到它们：

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

Answer 3

我相信这是因为这是 pi 上最短的系列。这也意味着它是最准确的。

Gregory-Leibniz 级数 (4/1 - 4/3 + 4/5 - 4/7...) 等于 pi。

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

所以，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...

等于 Gregory-Leibniz 级数，因此等于 pi，大约 3.1415926535 8979323846 2643383279 5028841971 69399373510。

另一种使用 atan 和查找 pi 的方法是：

pi = 16*atan(1/5) - 4*atan(1/239)，但我认为这更复杂。

我希望这会有所帮助！

（说实话，我认为Gregory-Leibniz系列是基于atan的，而不是基于Gregory-Leibniz系列的4*atan(1)。换句话说，REAL证明是：

sin^2 x + cos^2 x = 1 [定理] 如果 x = pi/4 弧度，则 sin^2 x = cos^2 x，或 sin^2 x = cos^2 x = 1/2。

那么，sin x = cos x = 1/(root 2)。 tan x (sin x / cos x) = 1, atan x (1 / tan x) = 1。

所以如果 atan(x) = 1，x = pi/4，atan(1) = pi/4。 最后，4*atan(1) = pi.)

请不要让我使用 cmets - 我还是个未成年人。

Answer 4

这是因为这是将pi 计算到任意精度的精确方法。您可以简单地继续执行该函数以获得越来越高的精度，并在任何点停止以获得近似值。

相比之下，将pi 指定为常数可为您提供与最初给出的精确度一样高的精度，这可能不适合高度科学或数学的应用程序（Fortran 经常使用）。

Answer 5

这个问题比表面上看到的要多。为什么4 arctan(1)？为什么不使用任何其他表示形式，例如 3 arccos(1/2)？

这将尝试通过排除找到答案。

数学简介：在使用反三角函数如arccos、arcsin和arctan时，可以很容易地以各种方式计算 π：

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

还有很多其他的exact algebraic expressions for trigonometric values 可以在这里使用。

浮点参数 1：众所周知，有限二进制浮点表示不能表示所有实数。此类数字的一些示例是1/3, 0.97, π, sqrt(2), ...。为此，我们应该排除任何无法用数值表示反三角函数的参数的 π 数学计算。这给我们留下了参数-1,-1/2,0,1/2 和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)

浮点参数2：在二进制表示中，一个数表示为0.b_nb_n-1...b₀ x 2^m。如果反三角函数为其参数提出了最佳数值二进制近似值，我们不希望因乘法而损失精度。为此，我们应该只乘以 2 的幂。

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

注意：这在 IEEE-754 binary64 表示中可见（DOUBLE PRECISION 或 kind=REAL64 的最常见形式）。我们有

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

IEEE-754 binary32（REAL 或kind=REAL32 的最常见形式）和IEEE-754 binary128（kind=REAL128 的最常见形式）没有这种区别

实现参数：在英特尔 CPU 上，atan2 是 x86 Instruction set as FPATAN 的一部分，而其他反三角函数则派生自 atan2。一个潜在的推导可能是：

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

这可以在这些指令的汇编代码中看到（参见here）。为此，我会争论以下用法：

π = 4 arctan(1)

注意：这是一个模糊的论点。我敢肯定有些人对此有更好的看法。
有趣的阅读FPATAN： How is arctan implemented?，x87 trigonometric instructions

Fortran 论证：我们为什么要将π 近似为：

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)

而不是：

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

答案在Fortran standard。标准从不规定任何类型的REAL 都应代表IEEE-754 floating point number。 REAL 的表示取决于处理器。这意味着我可以查询selected_real_kind(33, 4931) 并期望获得binary128 floating-point number，但我可能会返回一个kind，它表示精度更高的浮点数。也许100位数，谁知道呢。在这种情况下，我上面的一串数字太短了！不能为了确定而使用this？即使那个文件也可能太短了！

有趣的事实：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"

理解为：

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

Answer 6

这听起来很像解决编译器错误的方法。或者可能是这个特定的程序依赖于这个身份是准确的，所以程序员保证了它。