GerritGroot
Technical User
Hi,
I need to compute a series of factorials for high numbers, however, for very high numbers, normal factorials easily break down.
Let me firts post my code
and the outcome
The original function I got is function 1 (Func1):
In order to prevent overflow I wrote that as function 2 (Func2):
Which allows me to extend it a bit further for bigger n using the natural LOG, to avoid overflow (see code above).
However, I need to extend it even further, for higher n, upto 100000 for example without losing much precision.
Is there another mathematical trick to prevent overflow without losing much precision then the one with the LOG that I already used?
Thanks,
Gerrit
P.S.
- I don't want to get to double precision and all that, as I don't see that as a real solution
- Approximations for the factorial, like the Stirling one, don't really prevent overflow.
I need to compute a series of factorials for high numbers, however, for very high numbers, normal factorials easily break down.
Let me firts post my code
Code:
PROGRAM FactorialTest
IMPLICIT NONE
REAL :: Func1,Func2
WRITE(*,*)Func1(10,5),Func2(10,5)
WRITE(*,*)Func1(20,10),Func2(20,10)
END PROGRAM FactorialTest
REAL FUNCTION Func1(n,alfa)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n,alfa
INTEGER :: Factorial
Func1=REAL(Factorial(n))/(Factorial(alfa)*Factorial(n-alfa))
END FUNCTION Func1
REAL FUNCTION Func2(n,alfa)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n,alfa
INTEGER :: i
REAL :: sum
sum=0.
DO i=0,alfa-1,1
sum=sum+LOG((n-i)/(i+1.))
END DO
Func2=EXP(sum)
END FUNCTION Func2
RECURSIVE FUNCTION Factorial(n) RESULT(Fact)
IMPLICIT NONE
INTEGER :: Fact
INTEGER, INTENT(IN) :: n
IF(n<=0)THEN
Fact=1
ELSE
Fact=n*Factorial(n-1)
END IF
END FUNCTION Factorialand the outcome
Code:
252.000000 252.000000
11.659760 184755.900000The original function I got is function 1 (Func1):
Code:
n!/(alfa!*(n-alfa)!)In order to prevent overflow I wrote that as function 2 (Func2):
Code:
Product(from i=0 upto alfa-1) of { (n-i)/(i+1) }Which allows me to extend it a bit further for bigger n using the natural LOG, to avoid overflow (see code above).
Code:
n=10 and alfa=5 gives => 252 both functions work well
n=20 and alfa=10 gives => 184756 only the second function approximates reasonably within +/-0.1 due to the mathematical trick I used.However, I need to extend it even further, for higher n, upto 100000 for example without losing much precision.
Is there another mathematical trick to prevent overflow without losing much precision then the one with the LOG that I already used?
Thanks,
Gerrit
P.S.
- I don't want to get to double precision and all that, as I don't see that as a real solution
- Approximations for the factorial, like the Stirling one, don't really prevent overflow.