Fortran中具有大量实数的运算 [英] Operations with big real numbers in Fortran

问题描述

我编写了一个Fortran代码，该代码计算给定列表 {1,2,3，...，n} 的ith-排列，而不计算所有其他列，这是 n！我需要它来找到TSP的ith路径（旅行推销员问题）。

I wrote a Fortran code that calculates the ith-permutation of a given list {1,2,3,...,n}, without computing all the others, that are n! I needed that in order to find the ith-path of the TSP (Travelling salesman problem).

当 n！大时，代码给了我一些错误，我测试了发现的ith-置换不是确切值。对于n = 10，根本没有问题，但是对于 n = 20 ，代码崩溃或找到错误的值。我认为这是由于Fortran进行大数运算（大数加法）而导致的错误。

When n! is big, the code gives me some error and I tested that the ith-permutation found is not the exact value. For n=10, there are not problems at all, but for n=20, the code crashes or wrong values are found. I think this is due to errors that Fortran makes operating with big numbers (sums of big numbers).

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

I use Visual Fortran Ultimate 2013. In attached you find the subroutine I use for my goal. WeightAdjMatRete is the distance matrix between each pair of knots of the network.

    ! 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

推荐答案

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

Don't use floating-point numbers to represent factorials; factorials are products of integers and are therefore best represented as integers.

构造函数快速增长，因此使用实数可能很诱人，因为实数可以表示1.0之类的大数字。 e + 30。但是浮点数仅在与幅度有关时才是精确的。它们的尾数仍然有限，它们可能很大，因为它们的指数可能很大。

Factorials grow big fast, so it may be tempting to use reals, because reals can represent huge numbers like 1.0e+30. But floating-point numbers are precise only with relation to their magnitude; their mantissa still has a limited size, they can be huge because their exponents may be huge.

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

A 32-bit real can represent exact integers up to about 16 million. After that, only every even integer can be represented up to 32 million and every fourth integer up to 64 million. 64-bit integers are better, because they can represent exact integers up to 9 quadrillion.

64位整数可以扩展1024倍：它们可以表示2 ^ 63或大约9百万个（9e + 18）整数。足够代表20个！：

64-bit integers can go 1024 times further: They can represent 2^63 or about 9 quintillion (9e+18) integers. That is enough to represent 20!:

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

Fortran允许您选择一种基于整数的类型

Fortran allows you to select a kind of integer based on the decimal places it should be able to represent:

integer, (kind=selected_int_kind(18))

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

Use this to do your calculations with 64-bit integers. This will give you factorials up to 20!. It won't go further than that, though: Most machines support only integers up to 64 bit, so selected_int_kind(19) will give you an error.

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

Here's the permutation part of your program with 64-bit integers. Note how all the type conversions ald floors and ceilings disappear.

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中具有大量实数的运算的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！