! 桁落ちの実演とその回避
! 回避方法は[伊理, 藤野: 数値計算の常識]に載ってた.
! 
program ketaochi

  implicit none

  integer i,n
  real sum_i, sum_m, sum_c, a !sum_iは単なる加算, sum_mは乗算, 
                              !sum_cは桁落ち回避法を使ったもの

  a = 3.1415926535
  n = 1000000

  do i=0,n-1
     sum_i = sum_i + a
  end do
  sum_m = a*real(n)
  call res(a, n, sum_c)
  print *, (sum_m-sum_i)/sum_m, sum_m, sum_i
  print *, (sum_m-sum_c)/sum_m, sum_m, sum_c
end program ketaochi

subroutine res(a, n, sum)
  
  integer, intent(in) :: n
  real, intent(in) :: a
  real, intent(out) :: sum
  
  real t, s
  integer i

  sum = 0.0
  r = 0.0
  do i=0,n-1
     r = r+a
     t = sum
     sum = sum + r
     t = sum - t
     r = r-t
  end do
end subroutine res
