Factorial

Hello

The factorial function crops up quite a lot would it be handy as an intrinsic function in ftn95? Just a thought.

program factorial_test
    implicit none
    ! Store integers as 64-bit to handle very large numbers
    integer(kind=4) :: n
    integer(kind=4) :: fact_result

    ! Test the function with n = 20
    n = 20
    fact_result = factorial(n)

    ! Output the result
    print *, 'Factorial of ', n, ' is: ', fact_result

contains

    recursive function factorial(n) result(fact)

        integer(kind=4) :: n
        integer(kind=4) :: fact

        ! Check for negative or zero values
        if (n < 0) then
            print *, 'Error: Factorial is not defined for negative numbers.'
            fact = -1  ! Error code for negative input
        else if (n == 0) then
            fact = 1  ! Factorial of 0 is 1
        else
            fact = n * factorial(n - 1)  ! Recursive case
        end if

    end function factorial

end program factorial_test

A quick code testing a large value.

Lester

Lester,

Using the 64 bit compiler, you could do this calculation in one line of fortran.

Use an implied do loop within an array constructor, and then apply the product intrinsic to the array.

program fact
implicit none
integer    :: n
integer(4) :: ifact
integer(4) :: i
real*8     :: dfact     

do n = -1, 20, 1
  ifact = max( 0, product( [(i, i=1, n, 1)] ) )         ! Max required with /Checkmate
  write(*,'(I3,A,I20)') n, '! = ' , ifact
end do
write(*,*)

do n = 20, 170, 10
  dfact = max ( 0.d0, product( [(dble(i), i=1, n, 1)] ) )
  write(*,'(I3,A,EN20.8)') n, '! = ' , dfact
end do

end program fact

The code I posted earlier is not portable to other compilers, because of a programmer error!

Code below is, which explicitly deals with negative and zero values of n.

For n .ge. 1 the implied do loop / array constructor / product intrinsic remains the same, which is the real point I was trying to demonstrate.

program fact
implicit none
integer(4) :: n
integer(4) :: ifact
integer :: i
real*8     :: dfact     

do n = -1, 20, 1
  ifact = ifactorial(n)
  write(*,'(I3,A,I20)') n, '! = ' , ifact
end do

write(*,*)

do n = 20, 170, 10
  dfact = dfactorial(n)
  write(*,'(I3,A,EN20.8)') n, '! = ' , dfact
end do

contains


function ifactorial(n)  result(sol)
integer(4), intent(in) :: n
integer(4) :: i
integer(4) :: sol
  sol = 0
  if (n .ge. 1) sol =  product( [(i, i=1, n, 1)] ) 
  if (n .eq. 0) sol =  1
end function ifactorial


function dfactorial(n)  result(sol)
integer(4), intent(in) :: n
integer :: i
real*8  :: sol
  sol = 0.d0
  if (n .ge. 1) sol =  product( [(dble(i), i=1, n, 1)] ) 
  if (n .eq. 0) sol =  1.d0
end function dfactorial

end program fact

Lester

Is there a function like this in the Fortran Standard? Our focus with respect to development has to be on conforming to the lastest standards.

Thanks for the ideas Ken, code works fine.

I am not aware of factorial being an intrinsic in the latest standard.

One thing we do have in ftn95 is the Gamma function. I know this only works on reals, so not sure how it can be flipped to give the same integer output as Ken's code?

In simple terms, Factorial(n) = Gamma(n + 1). In ftn95 the following works up to n=12, but fails for any value that is greater.

program test_gamma_factorial
implicit none
  real(kind=2) :: n = 12.0d0
  real(kind=2):: fact
  fact = gamma(n+1.0d0)
  print *, 'Factorial of ', int(n), ' is ', int(fact)
end program test_gamma_factorial

Lester

Lester, These answer both of your questions.

program q1
implicit none
  integer :: n
  real*8 :: fact
  integer :: i
  do i = 1, 25
    n = i
    fact = gamma( min( n, 20 ) + 1.0d0 )
    print *, 'Factorial of ', int(n), ' is ', int(fact,kind=4)
  end do
end program q1


program q2
implicit none
  real*8 :: fact
  integer :: i
  
  do i = 1, 175, 1
    fact = gamma( min( i, 170 ) + 1.0d0) 
    print *, 'Factorial of ', i, ' is ', fact
  end do

end program q2

Thanks for the update Ken.

In conclusion, we can say it is not required as an intrinsic function. A good discussion and helpful coding ideas.

Lester

Lester,

Looking at this again, it might be better to create an array for all factorial values the first time the function runs, and save these for the next call to the function, which then becomes a simple 'look up' problem, which executes only three lines of code.

program p
implicit none
integer i
do i = 1, 20
  print*, i, factorial(i)
end do

contains

function factorial(n) result (fact)
integer, parameter :: wp = kind(1.d0)
integer, intent(in) :: n
real(wp) :: fact
integer :: i
logical,  save :: first = .true.
real(wp), save :: fact_vals(0:165)
  ! Only do the factorial calculation once for 0! to 165! and save all
  ! values for next time the function is called. 
  if (first) then
      fact = 1.0_wp
      fact_vals(0) = fact
        do i = 1, 165, 1
          fact = fact * dble(i)
          fact_vals(i) = fact 
        end do
      first = .false.
  end if
  ! Second time the function runs values are already precomputed
  if (n .gt. 165) STOP 
  fact = fact_vals(n)
end function factorial

end program p