fortran66のブログ

fortran について書きます。Amazonのアソシエイトとして収入を得ています。

【メモ帳】cos(pi/2^n)

倍角・半角の公式から

sin(2x) = 2sin(x)cos(x)=2sin(x)sqrt(1-sin2(x))

sin2(2x) = 4sin2(x)(1-sin2(x))

t = sin2(x) と置くと

4t2-4t+sin2(2x) = 0

この二次方程式を解くと

sin(x) = sqrt(1-sqrt(1-sin2(2x))) / sqrt(2)

これから cos(x) も求まる。1±cos(2x)〜sin2(x), cos2(x) と二次方程式の解の複号とからもわかる。

x=pi から始めてこの公式を繰り返し適用すると、

sin(x/2n)=sqrt(2-sqrt(2+sqrt(2+... )))) /2 , n >=2  

cos(x/2n)=sqrt(2+sqrt(2+sqrt(2+... )))) /2, n >=1

となる。

古代人は、これを逆行する形で微小な x からはじめて、 sin(x) ~ x により初期値を決めて、有限の大きさの三角関数表を作ったとか。(正確には現代の sin / cos ではないが)

プログラム

    program cos2n
        use, intrinsic :: iso_fortran_env 
        implicit none
        integer, parameter :: kd = real64
        real(kd), parameter :: pi = 4 * atan(1.0_kd)
        real(kd) :: t, x, y
        integer :: n, k 
        !
        ! cos(pi / 2^n) = sqrt(2 + sqrt(2 + ..... sqrt(2))....)) / 2 ;  n - 1 times
        !
        do n = 1, 26
            t = 0.0_kd    
            do k = 1, n - 1
                t = sqrt(2.0_kd + t)
            end do
            t = t / 2.0_kd
            !
        !    x = pi / 2**n
        !    y = cos(x)
            y = cosd(180.0_kd/2**n)
            print *, t, y, abs(t - y) 
        end do    
    end program cos2n

実行結果

  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
  0.707106781186548       0.707106781186548       0.000000000000000E+000
  0.923879532511287       0.923879532511287       0.000000000000000E+000
  0.980785280403230       0.980785280403230       0.000000000000000E+000
  0.995184726672197       0.995184726672197       0.000000000000000E+000
  0.998795456205172       0.998795456205172       0.000000000000000E+000
  0.999698818696204       0.999698818696204       0.000000000000000E+000
  0.999924701839145       0.999924701839145       0.000000000000000E+000
  0.999981175282601       0.999981175282601       0.000000000000000E+000
  0.999995293809576       0.999995293809576       0.000000000000000E+000
  0.999998823451702       0.999998823451702       0.000000000000000E+000
  0.999999705862882       0.999999705862882       0.000000000000000E+000
  0.999999926465718       0.999999926465718       0.000000000000000E+000
  0.999999981616429       0.999999981616429       1.110223024625157E-016
  0.999999995404107       0.999999995404107       0.000000000000000E+000
  0.999999998851027       0.999999998851027       0.000000000000000E+000
  0.999999999712757       0.999999999712757       0.000000000000000E+000
  0.999999999928189       0.999999999928189       1.110223024625157E-016
  0.999999999982047       0.999999999982047       0.000000000000000E+000
  0.999999999995512       0.999999999995512       0.000000000000000E+000
  0.999999999998878       0.999999999998878       1.110223024625157E-016
  0.999999999999719       0.999999999999719       0.000000000000000E+000
  0.999999999999930       0.999999999999930       0.000000000000000E+000
  0.999999999999982       0.999999999999982       0.000000000000000E+000
  0.999999999999996       0.999999999999996       1.110223024625157E-016
  0.999999999999999       0.999999999999999       0.000000000000000E+000