Fortran中具有大量实数的运算

  • Z时代
  • 2024-01-10
  • 分类：问答

我编写了一个Fortran代码，该代码可以计算给定列表的ith-

排列{1,2,3,...,n}，而无需计算所有其他列表，这是n!我需要的以便找到TSP的ith-路径（旅行商问题）。

当n!大，代码给了我一些错误，我测试的第i个置换发现的是不是确切值。对于n =

10，根本没有问题，但是对于n=20，代码崩溃或找到错误的值。我认为这是由于Fortran进行大数运算（大数加法）而导致的错误。

我使用Visual Fortran

Ultimate2013。在附件中，您找到了我用于目标的子例程。WeightAdjMatRete是网络的每对结之间的距离矩阵。

    ! Fattoriale
RECURSIVE FUNCTION factorial(n) RESULT(n_factorial)
IMPLICIT NONE
REAL, INTENT(IN) :: n
REAL :: n_factorial
IF(n>0) THEN
    n_factorial=n*factorial(n-1)
ELSE
    n_factorial=1.
ENDIF
ENDFUNCTION factorial
! ith-permutazione di una lista
SUBROUTINE ith_permutazione(lista_iniziale,n,i,ith_permutation)
IMPLICIT NONE
INTEGER :: k,n
REAL :: j,f
REAL, INTENT(IN) :: i
INTEGER, DIMENSION(1:n), INTENT(IN) :: lista_iniziale
INTEGER, DIMENSION(1:n) :: lista_lavoro
INTEGER, DIMENSION(1:n), INTENT(OUT) :: ith_permutation
lista_lavoro=lista_iniziale
j=i
DO k=1,n
    f=factorial(REAL(n-k))
    ith_permutation(k)=lista_lavoro(FLOOR(j/f)+1)
    lista_lavoro=PACK(lista_lavoro,MASK=lista_lavoro/=ith_permutation(k))
    j=MOD(j,f)
ENDDO
ENDSUBROUTINE ith_permutazione
! Funzione modulo, adattata
PURE FUNCTION mood(k,modulo) RESULT(ris)
IMPLICIT NONE
INTEGER, INTENT(IN) :: k,modulo
INTEGER :: ris
IF(MOD(k,modulo)/=0) THEN
    ris=MOD(k,modulo)
ELSE
    ris=modulo
ENDIF
ENDFUNCTION mood
! Funzione quoziente, adattata
PURE FUNCTION quoziente(a,p) RESULT(ris)
IMPLICIT NONE
INTEGER, INTENT(IN) :: a,p
INTEGER :: ris
IF(MOD(a,p)/=0) THEN
    ris=(a/p)+1
ELSE
    ris=a/p
ENDIF
ENDFUNCTION quoziente
! Vettori contenenti tutti i payoff percepiti dagli agenti allo state vector attuale e quelli ad ogni sua singola permutazione
SUBROUTINE tuttipayoff(n,m,nodi,nodi_rete,sigma,bvector,MatVecSomma,VecPos,lista_iniziale,ith_permutation,lunghezze_percorso,WeightAdjMatRete,array_perceived_payoff_old,array_perceived_payoff_neg)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n,m,nodi,nodi_rete
INTEGER, DIMENSION(1:nodi), INTENT(IN) :: sigma
INTEGER, DIMENSION(1:nodi), INTENT(OUT) :: bvector
REAL, DIMENSION(1:m,1:n), INTENT(OUT) :: MatVecSomma
REAL, DIMENSION(1:m), INTENT(OUT) :: VecPos
INTEGER, DIMENSION(1:nodi_rete), INTENT(IN) :: lista_iniziale
INTEGER, DIMENSION(1:nodi_rete), INTENT(OUT) :: ith_permutation
REAL, DIMENSION(1:nodi_rete), INTENT(OUT) :: lunghezze_percorso
REAL, DIMENSION(1:nodi_rete,1:nodi_rete), INTENT(IN) :: WeightAdjMatRete
REAL, DIMENSION(1:nodi), INTENT(OUT) :: array_perceived_payoff_old,array_perceived_payoff_neg
INTEGER :: i,j,k
bvector=sigma
FORALL(i=1:nodi,bvector(i)==-1)
    bvector(i)=0
ENDFORALL
FORALL(i=1:m,j=1:n)
    MatVecSomma(i,j)=bvector(m*(j-1)+i)*(2.**REAL(n-j))
ENDFORALL
FORALL(i=1:m)
    VecPos(i)=1.+SUM(MatVecSomma(i,:))
ENDFORALL
DO k=1,nodi
    IF(VecPos(mood(k,m))<=factorial(REAL(nodi_rete))) THEN
        CALL ith_permutazione(lista_iniziale,nodi_rete,VecPos(mood(k,m))-1.,ith_permutation)
        FORALL(i=1:(nodi_rete-1))
            lunghezze_percorso(i)=WeightAdjMatRete(ith_permutation(i),ith_permutation(i+1))
        ENDFORALL
        lunghezze_percorso(nodi_rete)=WeightAdjMatRete(ith_permutation(nodi_rete),ith_permutation(1))
        array_perceived_payoff_old(k)=(1./SUM(lunghezze_percorso))
    ELSE
        array_perceived_payoff_old(k)=0.
    ENDIF
    IF(VecPos(mood(k,m))-SIGN(1,sigma(m*(quoziente(k,m)-1)+mood(k,m)))*2**(n-quoziente(k,m))<=factorial(REAL(nodi_rete))) THEN
        CALL ith_permutazione(lista_iniziale,nodi_rete,VecPos(mood(k,m))-SIGN(1,sigma(m*(quoziente(k,m)-1)+mood(k,m)))*2**(n-quoziente(k,m))-1.,ith_permutation)
        FORALL(i=1:(nodi_rete-1))
            lunghezze_percorso(i)=WeightAdjMatRete(ith_permutation(i),ith_permutation(i+1))
        ENDFORALL
        lunghezze_percorso(nodi_rete)=WeightAdjMatRete(ith_permutation(nodi_rete),ith_permutation(1))
        array_perceived_payoff_neg(k)=(1./SUM(lunghezze_percorso))
    ELSE
        array_perceived_payoff_neg(k)=0.
    ENDIF
ENDDO
ENDSUBROUTINE tuttipayoff

回答：

不要使用浮点数来表示阶乘；阶乘是整数的乘积，因此最好用整数表示。

阶乘增长很快，因此使用实数可能很诱人，因为实数可以表示1.0e +

30之类的巨大数字。但是浮点数仅在其大小上是精确的。它们的尾数仍然有限，它们可能很大，因为它们的指数可能很大。

32位实数可以表示精确的整数，最大约为1600万。之后，只有每个偶数整数最多可以表示3200万，而每个第四个整数最多可以表示6400万。64位整数更好，因为它们可以表示9兆位以内的精确整数。

64位整数可以进一步移动1024倍：它们可以表示2 ^ 63或大约9百万个整数（9e + 18）整数。这足以代表20 ！：

 20! = 2,432,902,008,176,640,000
2^63 = 9,223,372,036,854,775,808

Fortran允许您根据其应代表的小数位数选择一种整数：

integer, (kind=selected_int_kind(18))

使用它可以对64位整数进行计算。这将使您的阶乘最多为20!。不过，它不会比这更进一步：大多数计算机仅支持最大64位的整数，因此selected_int_kind(19)会给您一个错误。

这是程序的64位整数置换部分。注意所有类型转换如何表示地板和天花板消失。

program permute
    implicit none
    integer, parameter :: long = selected_int_kind(18)
    integer, parameter :: n = 20
    integer, dimension(1:n) :: orig
    integer, dimension(1:n) :: perm
    integer(kind=long) :: k
    do k = 1, n
        orig(k) = k
    end do
    do k = 0, 2000000000000000_long, 100000000000000_long
        call ith_perm(perm, orig, n, k)
        print *, k
        print *, perm
        print *
    end do
end program
function fact(n)
    implicit none
    integer, parameter :: long = selected_int_kind(18)
    integer(kind=long) :: fact
    integer, intent(in) :: n
    integer :: i
    fact = 1
    i = n
    do while (i > 1)
       fact = fact * i
       i = i - 1
    end do
end function fact
subroutine ith_perm(perm, orig, n, i)
    implicit none
    integer, parameter :: long = selected_int_kind(18)
    integer, intent(in) :: n
    integer(kind=long), intent(in) :: i
    integer, dimension(1:n), intent(in) :: orig
    integer, dimension(1:n), intent(out) :: perm
    integer, dimension(1:n) :: work
    integer :: k
    integer(kind=long) :: f, j
    integer(kind=long) :: fact
    work = orig
    j = i
    do k = 1, n
        f = fact(n - k)
        perm(k) = work(j / f + 1)
        work = pack(work, work /= perm(k))
        j = mod(j, f)
    end do
end subroutine ith_perm

以上是 Fortran中具有大量实数的运算 的全部内容， 来源链接： utcz.com/qa/424865.html

回到顶部