subroutine show_precision(p, m) implicit none double precision p, m double precision p0 p0 = sqrt(p**2 + m**2) print *, p, p0 - p, m**2/(p0 + p) end program lose_precision implicit none call show_precision(1D3, 1D0) call show_precision(1D6, 1D0) call show_precision(1D9, 1D0) call show_precision(1D12, 1D0) call show_precision(1D15, 1D0) end