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.
Your functions are written correctly. The problem is in the main program, where you do not explicitly declare the type of integrate and factorial functions, so you have implicit typing, in which case factorial is assumed REAL and integrate is assumed INTEGER. For some reason, your compiler did not warn you about type mismatch. Mine did:
$ gfortran recurs.f90
recurs.f90:26.22:
PRINT *, integrate(2)
1
Error: Return type mismatch of function 'integrate' at (1) (INTEGER(4)/REAL(8))
recurs.f90:27.22:
PRINT *, factorial(5)
1
Error: Return type mismatch of function 'factorial' at (1) (REAL(4)/INTEGER(4))
You should change your main program to:
PROGRAM main
IMPLICIT NONE
DOUBLE PRECISION, EXTERNAL :: integrate
INTEGER, EXTERNAL :: factorial
PRINT *, factorial(5)
PRINT *, integrate(2)
END PROGRAM main
Notice the IMPLICIT NONE line. This declaration statement will disable any implicit typing, and the compiler would throw an error if not all variables and functions are explicitly declared. This is a very important line in every Fortran program, and if you had it, you would've figured out your problem yourself, because it would force you to explicitly declare everything in your program.
The output now is:
120
-19.0000000000000
as expected.
As a side note, the DOUBLE PRECISION type declaration is not as flexible as using REAL with KIND parameter specified instead, e.g. anREAL(KIND=myRealKind). See answers to this question about how to use KIND properly: Fortran 90 kind parameter.