c语言怎么表示pi atan1,fortran-为什么定义PI = 4 * ATAN（1.d0）

这个问题不仅仅令人眼前一亮。 为什么是sin(pi) is never zero？ 为什么没有其他表示形式，例如REAL？

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

数学简介：当使用反三角函数(例如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：众所周知，有限的二进制浮点表示不能表示所有实数。 这样的数字的一些示例是sin(pi) is never zero。为此，我们应该排除π的任何数学计算，其中不能用数字表示反三角函数的参数。 这留下了争论REAL和selected_real_kind(33, 4931)。

π = 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.bnbn-1 ... b0 x 2m。 如果逆三角函数为其参数提供了最佳的数值二进制近似值，则我们不想因乘法而损失精度。 为此，我们只能乘以2的幂。

π = 4 arctan(1) = 2 arcsin(1)

= 2 arccos(0)

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

注意：这在IEEE-754 binary64表示形式(sin(pi) is never zero或REAL的最常见形式)中可见。 那里我们有

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

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

IEEE-754二进制32(最常见的sin(pi) is never zero或REAL)和IEEE-754二进制128(最常见的selected_real_kind(33, 4931))没有这种区别。

模糊实现的论点：从这一点来看，一切都取决于反三角函数的实现。 有时sin(pi) is never zero和REAL源自selected_real_kind(33, 4931)和kind作为

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

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

或更具体地从数字角度来看：

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

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

此外，sin(pi) is never zero是设置为sin(pi) is never zero的x86指令的一部分，而其他不是。 为此，我将争论以下用法：

π = 4 arctan(1)

超过所有其他。

注意：这是一个模糊的论点。 我敢肯定有人对此有更好的看法。

Fortran参数：为什么我们应该将sin(pi) is never zero近似为：

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标准。 该标准从未声明任何sin(pi) is never zero应该代表IEEE-754浮点数。 REAL的表示形式取决于处理器。 这意味着我可以查询selected_real_kind(33, 4931)并希望获得二进制128浮点数，但是我可能会得到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