fortran - Fortran 中的大整数阶乘函数，与 python 或 Haskell 一样高效

Question

这是我在 Fortran 中的阶乘函数。

module facmod
  implicit none

contains
  function factorial (n) result (fac)
    use FMZM 
    integer, intent(in)  :: n
    integer              :: i
    type(IM)             :: fac

    fac = 1

    if(n==0) then
       fac = 1
    elseif(n==1) then
       fac = 1
    elseif(n==2) then
       fac = 2
    elseif(n < 0) then
       write(*,*) 'Error in factorial N=', n
       stop 1
    else
       do i = 1, n
          fac = fac * i
       enddo
    endif
  end function factorial

end module facmod

program main
   use FMZM   
   use facmod, only: factorial
   implicit none
   type(IM) :: res
   integer :: n, lenr
   character (len=:), allocatable :: str
   character(len=1024) :: fmat
   print*,'enter the value of n'
   read*, n
   res = factorial(n)
   lenr = log10(TO_FM(res))+2
   allocate(character(len=lenr) :: str)
   write (fmat, "(A5,I0)") "i", lenr
   call im_form(fmat, res, str)   
   print*, trim( adjustl(str))
 end program main

我使用 FMZM 编译：

gfortran -std=f2008 fac.F90 fmlib.a -o fac

echo -e "1000" | .fac计算容易。但是，如果我给出这个echo -e "3600" | .fac，我的机器上已经出现错误：

  Error in FM.  More than       200000  type (FM), (ZM), (IM) numbers
                have been defined.  Variable  SIZE_OF_START  in file
                FMSAVE.f95  defines this value.
                Possible causes of this error and remedies:
                (1) Make sure all subroutines (also functions that do not
                    return type FM, ZM, or IM function values) have
                        CALL FM_ENTER_USER_ROUTINE
                    at the start and 
                        CALL FM_EXIT_USER_ROUTINE
                    at the end and before any other return, and all
                    functions returning an FM, ZM, or IM function value have
                        CALL FM_ENTER_USER_FUNCTION(F)
                    at the start and 
                        CALL FM_EXIT_USER_FUNCTION(F)
                    at the end and before any other return, where the actual
                    function name replaces  F  above.
                    Otherwise that routine could be leaking memory, and
                    worse, could get wrong results because of deleting some
                    FM, ZM, or IM temporary variables too soon.
                (2) Make sure all subroutines and functions declare any
                    local type FM, ZM, or IM variables as saved.  Otherwise
                    some compilers create new instances of those variables
                    with each call, leaking memory.
                    For example:
                        SUBROUTINE SUB(A,B,C,X,Y,RESULT)
                        TYPE (FM) :: A,B,C,X,Y,RESULT,ERR,TOL,H
                    Here A,B,C,X,Y,RESULT are the input variables and
                    ERR,TOL,H are local variables.  The fix is:
                        SUBROUTINE SUB(A,B,C,X,Y,RESULT)
                        TYPE (FM) :: A,B,C,X,Y,RESULT
                        TYPE (FM), SAVE :: ERR,TOL,H
                (3) Since = assignments for multiple precision variables are
                    the trigger for cleaning up temporary multiple precision
                    variables, a loop with subroutine calls that has no =
                    assignments can run out of space to store temporaries.
                    For example:
                        DO J = 1, N
                           CALL SUB(A,B,C,TO_FM(0),TO_FM(1),RESULT)
                        ENDDO
                    Most compilers will create two temporary variables with
                    each call, to hold the TO_FM values.
                    One fix is to put an assignment into the loop:
                        DO J = 1, N
                           ZERO = TO_FM(0)
                           CALL SUB(A,B,C,ZERO,TO_FM(1),RESULT)
                        ENDDO
                (4) If a routine uses allocatable type FM, ZM, or IM arrays
                    and allocates and deallocates with each call, then after
                    many calls this limit on number of variables could be 
                    exceeded, since new FM variable index numbers are
                    generated for each call to the routine.
                    A fix for this is to call FM_DEALLOCATE before actually
                    deallocating each array, so those index numbers can be
                    re-used.  For example:
                        DEALLOCATE(T)
                    becomes:
                        CALL FM_DEALLOCATE(T)
                        DEALLOCATE(T)
                (5) If none of this helps, try running this program again
                    after increasing the value of  SIZE_OF_START  and
                    re-compiling.

我错过了哪些优化或 Fortran 习语会严重影响我的表现？

例如，在 python 中，我可以对远大于 3500 的数进行阶乘：

>>> import math
>>> math.factorial(100000)

或者在 Haskell 中：

Prelude> product [1..100000]

这两种计算都不是很快，但没有错误。

如何改进我的算法或更好地使用现有库来提高 Fortran 中大整数阶乘的性能？有没有比 FMZM 更合适的大整数库？

标签： fortranbigintegerarbitrary-precisionfmzm

Answer 1

尝试这个。除了较小的外观更改外，我只是遵循了您问题中错误消息的建议：

• 添加了对 FM_ENTER_USER_FUNCTION 和 FM_EXIT_USER_FUNCTION 的调用，
• 在循环内添加了一个赋值（没有 this ii = to_im(i)，它仍然失败，但我不知道为什么，因为已经有一个赋值fac = fac * i，并且根据文档，赋值触发清理临时对象），
• factorial在主程序中重命名，因为 FMZM 中已经有一个具有此名称的函数。

使用 ifort 和 n=100000 进行测试。

module fac_mod
    implicit none
contains
    function factorial(n) result(fac)
        use FMZM
        integer, intent(in) :: n
        integer :: i
        type(IM) :: fac
        type(IM), save :: ii

        call FM_ENTER_USER_FUNCTION(fac)
        fac = to_im(1)

        if (n < 0) then
            write (*, *) "Error in factorial N=", n
            stop 1
        else if (n > 1) then
            do i = 1, n
                ii = to_im(i)
                fac = fac * ii
            end do
        end if
        call FM_EXIT_USER_FUNCTION(fac)
    end function factorial
end module fac_mod

program main
    use FMZM   
    use fac_mod, only: f=>factorial
    implicit none
    type(IM) :: res
    integer :: n, lenr
    character(:), allocatable :: str
    character(1024) :: fmat

    print *, "enter the value of n"
    read *, n
    res = f(n)
    lenr = 2 + log10(TO_FM(res))
    allocate (character(lenr) :: str)
    write (fmat, "(A5,I0)") "i", lenr
    call im_form(fmat, res, str)   
    print *, trim(adjustl(str))
end program main