SOLVED: From Python to Fortan an easy code

1

Hello, I am trying to emulate this code written in Python but I am not able to:

  from math import factorial as fact
  def serie(n):
    return 1/(n*(n+1)*fact(n+1))

  sum([serie(i) for i in range(1,1000)])

This is my code in Modern Fortarn:

program Serie
    use iso_fortran_env
    implicit none

    integer, parameter :: MAX = 1000
    integer:: n
    real(kind=16) :: s
    real(kind=16), dimension(MAX):: fact

    ! Array of factorial
    fact(1) = 1
    do n=2, MAX
        fact(n) = fact(n-1) * n
    end do

    ! Main
    s = 0
    do n=1, MAX
        s = s + 1/(n*(n+1)*fact(n+1))
    end do

    ! Output
    print *, s
end program Serie

Can you help me?, thank you

2

Solved:

program Serie
    use iso_fortran_env
    implicit none

    integer, parameter :: MAX = 1000
    integer :: n
    real(kind=16) :: s
    real(kind=16), dimension(MAX) :: fact

    ! Array of factorial
    fact(1) = 1
    do n=2, MAX+1
        fact(n) = fact(n-1) * n
    end do

    ! Main
    s = 0
    !$omp parallel
    do n=1, MAX
        s = s + 1/(n*(n+1)*fact(n+1))
    end do
    !$omp end parallel

    ! Output
    write (*,"(f)") s
end program Serie
3

I am using the intel IFX compiler and it does not warn not shows any error with the array bounds ¿there must be a default flag?

4

use the -check option

This won’t give you a corect result

  • the loop will be entirely executed by all the threads
  • there is a race condition on s

The correct version:

    s = 0
    !$omp parallel do reduction(+:s)
    do n=1, MAX
        s = s + 1/(n*(n+1)*fact(n+1))
    end do
    !$omp end parallel do

Note that real(kind=16) arithmetic is emulated on x86 CPUs, thus quite slow.

1 Like
5

Just couple of comments. One would not normally store and compute the fact(:) array if you are only going to use each element a single time. Instead, the factorial would be computed in the loop itself. Also, the factorials grow very fast, so the individual terms in the series very quickly go to zero, and there is no need to sum all 1000 terms. Here is a short code that does the work without the fact(:) array, and tests for when the loop can be exited early.

program sumseries
   use, intrinsic :: iso_fortran_env, only: wp => real128
   implicit none
   integer, parameter :: nmax = 1000
   real(wp) :: s, fact, term
   integer :: n
   fact = 1.0_wp  ! (n+1)!
   s = 0.0_wp
   do n = 1, nmax
      fact =  fact * (n+1)
      term = 1.0_wp / ( fact * n * (n+1) )
      if ( term < spacing(s) ) exit  ! early exit
      s = s + term
   enddo
   write(*,'(a,i0,a,f0.36)') 'n=', n, ' s= ', s
end program sumseries
$ gfortran sumseries.f90  && a.out
n=29 s= .281718171540954764639712528647337493

Hopefully I have not missed a term by ±1 or something, but that is the way one might solve this task in fortran. If you want to change the precision, then just change the wp parameter appropriately.

2 Likes