60

将 PI 定义为的动机是什么

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

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

4

6 回答 6

67

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

于 2010-01-28T21:07:41.313 回答
15

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

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

PI=DACOS(-1.D0)
PI=2.D0*DASIN(1.D0)
于 2010-01-28T21:09:06.400 回答
14

我相信这是因为这是 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。）

请不要给我留言——我还是个未成年人。

于 2013-11-17T02:48:52.620 回答
9

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

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

于 2010-01-28T21:06:43.710 回答
9

这个问题比表面上看到的要多。为什么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) = ...

存在许多其他可以在此处使用的三角函数值的精确代数表达式。

浮点参数 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 n b n-1 ...b 0 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 PRECISIONor形式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（or 的最常见形式REAL）kind=REAL32和IEEE-754 binary128（or 的最常见形式kind=REAL128）中不存在这种差异

实现参数：在英特尔 CPU 上，atan2是x86 指令集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)))

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

π = 4 arctan(1)

注意：这是一个模糊的论点。我敢肯定有人对此有更好的看法。
有趣的继续阅读FPATAN： arctan 是如何实现的？, x87 三角函数指令

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 标准。该标准从未声明REAL任何类型的 a 都应表示IEEE-754 浮点数。的表示REAL取决于处理器。这意味着我可以查询selected_real_kind(33, 4931)并期望获得一个binary128 浮点数，但我可能会得到一个kind表示精度更高的浮点数的返回值。也许100位数，谁知道呢。在这种情况下，我上面的一串数字太短了！不能仅仅为了确定而使用它吗？即使那个文件也可能太短了！

有趣的事实 ：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
于 2018-03-21T20:51:19.953 回答
-7

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

于 2010-01-28T21:06:14.087 回答