recursion - Fortran functions returning unexpected types and values

Question

I'm working on a project that needs to implement few numerical methods in Fortran. For this, I need to write some recursive functions. Here is my code.

!     
! File:   main.F95
!

RECURSIVE FUNCTION integrate(n) RESULT(rv)
    IMPLICIT NONE
    DOUBLE PRECISION :: rv
    INTEGER, INTENT(IN) :: n
    DOUBLE PRECISION, PARAMETER :: minusone = -1.0
    IF (n == 1) THEN
        rv = 10 !exp(minusone)
        RETURN
    ELSE
        rv = 1 - (n * integrate(n - 1))
        RETURN
    END IF
END FUNCTION integrate

RECURSIVE FUNCTION factorial(n) RESULT(res)
    INTEGER res, n
    IF (n .EQ. 0) THEN
        res = 1
    ELSE
        res = n * factorial(n - 1)
    END IF
END

PROGRAM main
    DOUBLE PRECISION :: rv1
    PRINT *, factorial(5)
    PRINT *, integrate(2)

    !READ *, rv1

END PROGRAM main

For this program the output is:

         NaN
       1

If I change the order of the print statements (line 30 & 31), the output will be:

         1
-19.000000

Output should be (for the original print statement order):

  120  
  -19 

I took the factorial function from the Wikipedia Fortran 95 language features page.

• Compiler : gfortran 4.5.3 with Cygwin
• IDE: Netbeans 7.0.1
• Platform: Windows 7

Answer 1

As one of the comments mentions, a better solution is to put your subroutines and functions into a module, then use that module from your main program. This will make the interface of those procedures known to the caller -- "explicit" in Fortran terminology. Not only will the compiler correctly handle the type of the function, it will be able to check type-agreement between the arguments in the call and the arguments in the callee ("dummy arguments") for consistency.

If you use as many debugging options as possible the compiler will help you find mistakes. With gfortran, try: -O2 -fimplicit-none -Wall -Wline-truncation -Wcharacter-truncation -Wsurprising -Waliasing -Wimplicit-interface -Wunused-parameter -fwhole-file -fcheck=all -std=f2008 -pedantic -fbacktrace

module factorial_procs

   IMPLICIT NONE

contains

   RECURSIVE FUNCTION integrate(n) RESULT(rv)
       DOUBLE PRECISION :: rv
       INTEGER, INTENT(IN) :: n

       IF (n == 1) THEN
           rv = 10
           RETURN
       ELSE
           rv = 1 - (n * integrate(n - 1))
           RETURN
       END IF
   END FUNCTION integrate

   RECURSIVE FUNCTION factorial(n) RESULT(res)
       INTEGER res, n
       IF (n .EQ. 0) THEN
           res = 1
       ELSE
           res = n * factorial(n - 1)
       END IF
   END

end module factorial_procs

PROGRAM main

    use factorial_procs

    implicit none

    PRINT *, factorial(5)
    PRINT *, integrate(2)

END PROGRAM main

You'll probably find that you can only calculate factorials of very small integers by straight forward multiplication using regular integers. One fix is to use a larger integer type, e.g.,

integer, parameter :: BigInt_K = selected_int_kind (18)

Just as you could modernize and use selected_real_kind instead of Double Precision.