60

PI を次のように定義する動機は何ですか

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

Fortran 77 コード内? 仕組みは理解できたのですが、その理由は何ですか？

4

6 に答える 6

67

このスタイルにより、値を PI に割り当てるときに、任意のアーキテクチャで利用可能な最大精度が使用されるようになります。

于 2010-01-28T21:07:41.313 に答える
15

Fortran には の組み込み定数がないためPIです。ただし、数値を手動で入力して間違いを犯したり、指定された実装で可能な最大精度を得られなかったりするのではなく、ライブラリに結果を計算させることで、これらの欠点のいずれも発生しないことが保証されます。

これらは同等であり、時々それらも表示されます。

PI=DACOS(-1.D0)
PI=2.D0*DASIN(1.D0)
于 2010-01-28T21:09:06.400 に答える
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:有限の 2 進浮動小数点表現がすべての実数を表すことができないことはよく理解されています。そのような数の例は次のとおり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 binary64DOUBLE 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 binary32REAL (またはの最も一般的な形式kind=REAL32) とIEEE-754 binary128 ( の最も一般的な形式kind=REAL128)にはありません。

実装の引数:インテル CPU では、はx86 命令セットatan2の一部であり、他の逆三角関数は から派生しています。潜在的な派生物は次のとおりです。FPATANatan2

          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 標準にあります。標準では、どのような種類の もIEEE-754 浮動小数点数を表す必要があるとは決して述べていません。の表現はプロセッサに依存します。これは、照会してbinary128 浮動小数点数を取得できることを意味しますが、はるかに高い精度で浮動小数点を表す戻り値を取得する可能性があります。たぶん100桁、誰が知っていますか。この場合、上記の数字の文字列は短すぎます! 念のためにこれを使用することはできませんか？そのファイルでさえ短すぎる可能性があります。REALREALselected_real_kind(33, 4931)kind

興味深い事実 ：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

これは、コンパイラのバグの回避策のように思えます。または、この特定のプログラムはその ID が正確であることに依存しているため、プログラマーがそれを保証した可能性があります。

于 2010-01-28T21:06:14.087 に答える