LFortran 0.48 on mac

conda でも入りました。ネタ元 https://fortran-lang.discourse.group/t/lfortran-v0-48-0-available-on-homebrew/9343

テスト

lfortran も sind の実装が、radian 用を引数をラジアン変換する形で使っているようです。 degree 専用関数で精度が最後の 1 bit まで正確になるはずなのですが、なっていません。二宮の「数値計算のつぼ」に書いてあります。

    program sine
        implicit none
        real, parameter :: pi = 4 * atan(1.0)
        real, parameter :: val(*) = [0.0, 1.0, sqrt(3.0)/2.0, 1.0/sqrt(2.0), sqrt((5.0-sqrt(5.0))/8.0), 1.0/2.0] 
        ! sin(pi  ) = 0
        ! sin(pi/2) = 1
        ! sin(pi/3) = sqrt(3)/2   
        ! sin(pi/4) = 1/sqrt(2)
        ! sin(pi/5) = sqrt((5-sqrt(5))/8)
        ! sin(pi/6) = 1/2
        real :: x, y
        integer :: i
        
        do i = 1, 6
            x = pi / i
         !   y = sin(x)
            y = sind(180.0/i)
            print '(a, i1, 2f10.7, es15.7)', 'pi/', i, y, val(i), y - val(i)   
        end do
    end program sine

(lf) [a] M1:~% ./a.out
pi/1-0.0000001 0.0000000 -8.7422777E-08
pi/2 1.0000000 1.0000000  0.0000000E+00
pi/3 0.8660254 0.8660254  5.9604645E-08
pi/4 0.7071068 0.7071068  0.0000000E+00
pi/5 0.5877852 0.5877852  0.0000000E+00
pi/6 0.5000000 0.5000000  0.0000000E+00

    program sine
        implicit none
        real, parameter :: pi = 4 * atan(1.0)
        real, parameter :: val(*) = [0.0, 1.0, sqrt(3.0)/2.0, 1.0/sqrt(2.0), sqrt((5.0-sqrt(5.0))/8.0), 1.0/2.0] 
        ! sin(pi  ) = 0
        ! sin(pi/2) = 1
        ! sin(pi/3) = sqrt(3)/2   
        ! sin(pi/4) = 1/sqrt(2)
        ! sin(pi/5) = sqrt((5-sqrt(5))/8)
        ! sin(pi/6) = 1/2
        real :: x, y
        integer :: i
        
        do i = 1, 6
            x = pi / i
            y = sin(x)
         !   y = sind(180.0/i)
            print '(a, i1, 2f10.7, es15.7)', 'pi/', i, y, val(i), y - val(i)   
        end do
    end program sine

(lf) [a] M1:~% lfortran sine.f90
pi/1-0.0000001 0.0000000 -8.7422777E-08
pi/2 1.0000000 1.0000000  0.0000000E+00
pi/3 0.8660254 0.8660254  5.9604645E-08
pi/4 0.7071068 0.7071068  0.0000000E+00
pi/5 0.5877852 0.5877852  0.0000000E+00
pi/6 0.5000000 0.5000000  0.0000000E+00

参考 intel fortran

これは degree 用三角関数の精神を理解して、正しく実装されています。

[a] M1:~% ./a.out           
pi/1 0.0000000 0.0000000  0.0000000E+00
pi/2 1.0000000 1.0000000  0.0000000E+00
pi/3 0.8660254 0.8660254  0.0000000E+00
pi/4 0.7071068 0.7071068  0.0000000E+00
pi/5 0.5877852 0.5877852  0.0000000E+00
pi/6 0.5000000 0.5000000  0.0000000E+00

gfortran

ややおかしいです。

[a] M1:~% ./a.out
pi/1-0.0000000 0.0000000 -0.0000000E+00
pi/2 1.0000000 1.0000000  0.0000000E+00
pi/3 0.8660254 0.8660254  0.0000000E+00
pi/4 0.7071068 0.7071068  0.0000000E+00
pi/5 0.5877852 0.5877852  0.0000000E+00
pi/6 0.5000000 0.5000000  0.0000000E+00

flang-new (ver. 19.1.7)

これも弧度法用のルーチンを流用しているのみ。

[a] M1:~% ./a.out
pi/1-0.0000001 0.0000000 -8.7422777E-08
pi/2 1.0000000 1.0000000  0.0000000E+00
pi/3 0.8660254 0.8660254  5.9604645E-08
pi/4 0.7071068 0.7071068  0.0000000E+00
pi/5 0.5877852 0.5877852  0.0000000E+00
pi/6 0.5000000 0.5000000  0.0000000E+00