! 20110512
! 川畑拓也
!
! 0.01 を 10000 回足す計算を 3 通りの方法で行う
!
! 1) 単精度実数を用いる
! 2) 倍精度実数を用いる
! 3) 単精度実数を用いるが, 補助変数を導入
!
! 単精度でも, 補助変数を導入すれば精度が向上す
! ることを確かめる

program semidble

implicit none

!--- 変数の宣言 ---!

real(4) :: R1     ! 補助変数 R
real(4) :: R2     ! 補助変数 R
real(4) :: R3     ! 補助変数 R
real(4) :: S1     ! 部分和 S
real(4) :: S2     ! 部分和 S
real(4) :: T1     ! 補助変数 T
real(4) :: T2     ! 補助変数 T
real(4) :: Esum   ! 単精度での部分和
real(8) :: Dsum   ! 倍精度での部分和
integer :: ix     ! do ループ作業用変数


!--- 変数の初期化 ---!

R1 = 0.0e0
R2 = 0.0e0
R3 = 0.0e0
S1 = 0.0e0
S2 = 0.0e0
T1 = 0.0e0
T2 = 0.0e0
Esum = 0.0e0
Dsum = 0.0d0


!--- 計算開始 ---!

do ix = 1, 10000

   R2 = R1 + 0.010e0
   T1 = S1
   S2 = S1 + R2
   T2 = S2 - T1
   R3 = R2 - T2
   R1 = R3
   S1 = S2
   T1 = T2

   Esum = Esum + 0.010e0

   Dsum = Dsum + 0.010d0

   write(*,*) 'R1 =', R1
   pause 
enddo


!--- 結果の出力 ---!

write(*,*) 'real     = ', Esum
write(*,*) 'dble     = ', Dsum 
write(*,*) 'semidble = ', S1


end program semidble
