如何在代码中实现阶乘函数？

Question

所以我使用泰勒级数在 fortran 90 中计算 sin(0.75) 直到某个点，所以我需要在 do while 循环中运行它（直到满足我的条件）。这意味着我需要使用阶乘，这是我的代码：

program taylor
implicit none
real :: x = 0.75  
real :: y
integer :: i = 3
do while (abs(y - sin(0.75)) > 10.00**(-7)) 
 i = i + 2
 y = x - ((x**i)/fact(i))
 print *, y
end do
end program taylor

我写 fact(i) 的地方是我需要阶乘的地方。不幸的是，Fortran 没有内在的 !功能。我将如何在这个程序中实现该功能？

谢谢。

Answer 1

以下简单函数回答了您的问题。请注意它如何返回real，而不是整数。如果性能不是问题，那么这对泰勒级数来说很好。

real function fact(n)
  integer, intent(in) :: n

  integer :: i

  if (n < 0) error stop 'factorial is singular for negative integers'
  fact = 1.0
  do i = 2, n
    fact = fact * i
  enddo
end function fact

但真正的答案是 Fortran 2008 确实具有阶乘的内在函数：the Gamma function。对于正整数n，它被定义为Gamma(n+1) == fact(n)。

（我可以想象Gamma function 是不熟悉的。它是阶乘函数的概括：Gamma(x) 是为所有复数 x 定义的，除了非正整数。定义中的偏移量是出于历史原因和你问我的问题是不必要的混淆。）

在某些情况下，您可能希望将 Gamma 函数的输出转换为整数。如果是这样，请确保通过 INT(Gamma(n+1), kind=INT64) 使用“长整数”和 USE, INTRINSIC :: ISO_Fortran_env 声明。这是防止阶乘变得非常大的预防措施。并且一如既往地注意混合模式算术！

Answer 2

这是计算 n 的另一种方法！在一行中仅使用内联函数：

product((/(i,i=1,n)/))

当然 i 必须事先声明为整数。它创建一个从 1 到 n 的数组，并获取所有组件的乘积。奖励：它甚至可以在 n = 0 时给出正确的结果。

Answer 3

您不想对泰勒级数使用阶乘函数。这意味着一遍又一遍地计算相同的术语。您应该在每次循环迭代中乘以阶乘变量。不要忘记使用real，因为整数会很快溢出。

看你同学Program For Calculating Sin Using Taylor Expansion Not Working?问题下的答案

Answer 4

你能写出给出阶乘的方程吗？ 它可能看起来像这样

PURE FUNCTION Bang(N)
IMPLICIT NONE
INTEGER, INTENT(IN) :: N
INTEGER             :: I
INTEGER             :: Bang

Bang = N
IF(N == 2) THEN
  Bang = 2
ELSEIF(N == 1) THEN
  Bang = 1
ELSEIF(N < 1) THEN
  WRITE(*,*)'Error in Bang function N=',N
  STOP
ELSE
  DO I = (N-1), 2, -1
    Bang = Bang * I
  ENDDO
ENDIF

RETURN
END FUNCTION Bang