Very small number turns negative in Fortran

user_14515953 0
Viewed 416

I'm doing a program in Fortran90 which do a sum from i=1 to i=n where nis given. The sum is sum_{i=1}^{i=n}1/(i*(i+1)*(i+2)). This sum converges to 0.25. This is the code:

PROGRAM main

INTEGER n(4)
DOUBLE PRECISION s(4)
INTEGER i

  
OPEN(11,FILE='input')
OPEN(12,FILE='output')
DO i=1,4
  READ(11,*) n(i)
END DO

PRINT*,n
CALL suma(n,s)
PRINT*, s

END

 
SUBROUTINE suma(n,s)
INTEGER n(4),j,k
DOUBLE PRECISION s(4),add
s=0 
DO k=1,4
  DO j=1,n(k)
    add=1./(j*(j+1)*(j+2))
    s(k)=s(k)+add 
  END DO 
END DO
  
END SUBROUTINE

input

178   
1586     
18232 
142705

The output file is now empty, I need to code it. I'm just printing the results, which are: 0.249984481688 0.249999400246 0.248687836759 0.247565846142

The problem comes with the variable add. When j is bigger, add turns negative, and the sum doesn't converge well. How can I fix it?

fortran
2 Answers

The problem is an integer overflow. 142705142706142707 is a number that is too large for a 4-byte integer.

What happens then is that the number overflows and loops back to negative numbers.

As @albert said in his comment, one solution is to convert it to double precision every step of the way: ((1.d0/j) / (j+1)) / (j+2). That way, it is calculating with floating point values.

Another option would be to use 8-byte integers:

integer, parameter :: int64 = selected_int_kind(17)
integer(kind=int64) :: j

You should be very careful with your calculations, though. Finer is not always better. I recommend that you look at how floating point arithmetic is performed by a computer, and what issues this can create. See for example here on wikipedia.

user_46218230

This is likely a better way to achieve what you want. I did remove the IO. The output from the program is

% gfortran -o z a.f90 && ./z
   178 0.249984481688392
  1586 0.249999801599584
 18232 0.249999998496064
142705 0.249999999975453

program main

  implicit none  ! Never write a program without this statement

  integer, parameter :: knd = kind(1.d0) ! double precision kind

  integer n(4)
  real(knd) s(4)
  integer i

  n = [178, 1586, 18232, 142705]
  call suma(n, s)
  do i = 1, 4
    print '(I6,F18.15)', n(i), s(i)
  end do

  contains
     !
     ! Recursively, sum a(j+1) = j * a(j) / (j + 1)
     !
     subroutine suma(n, s)
        integer, intent(in) :: n(4)
        real(knd), intent(out) :: s(4)
        real(knd) aj
        integer j, k
        s = 0 
        do k = 1, 4
           aj = 1 / real(1 * 2 * 3, knd)    ! a(1)
           do j = 1, n(k)
              s(k) = s(k) + aj
              aj = j * aj / (j + 3)
           end do 
        end do
     end subroutine

end program main
user_98925340
Related