Operations with big real numbers in Fortran

Question

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).

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).

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

Did you try using double precision for the floating points variables in your code? — Alexander Vogt, Mar 21 '16 at 16:33
First try larger integer and real kinds http://stackoverflow.com/questions/838310/fortran-90-kind-parameter http://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4 — Vladimir F Героям слава, Mar 21 '16 at 16:33
Yes I use double precision on real numbers becouse reals permits a wide range of numbers rather Then The integers numbers. — Ilario De Vincenzo, Mar 21 '16 at 18:41
@IlarioDeVincenzo In your code, you are using single precision floating point numbers (`real :: ...`). To switch to double precision, you need to specify the corresponding `kind`, e.g., `real(kind=kind(1.d0)) :: ...` Double precision not only extends the range (exponent), but also the precision (mantissa). — Alexander Vogt, Mar 21 '16 at 18:53
@IlarioDeVincenzo: Yes, but the large real numbers may not have the required precision to store the factorial as integer. Single precision can represent all integers up to about 16 millions; double precision can represent all integers up to about 9 quadrillions (9·10^16). Factorials aren't real number; use integers. — M Oehm, Mar 21 '16 at 19:01
Even with 64-bit unsigned ints (doesFortran have unsigned ints?) you'll get only as far as 20!. You could as well use a [look-up table](https://blogs.msdn.microsoft.com/oldnewthing/20151214-00/?p=92621) instead of a recursive factorial function. — M Oehm, Mar 21 '16 at 19:05
Do you really need to compute factorials to generate a permutation ?? — , Mar 21 '16 at 20:12

score 1 · Accepted Answer · answered Mar 21 '16 at 20:48

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

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.

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-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 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))

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.

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

Thanks a lot for the answer. I find only one problem running the code you write. At certain point during iterations of my algorithm, the code interrupts because j/f+1 becomes bigger then the size of "work". How it is possible? — Ilario De Vincenzo, Mar 22 '16 at 18:36
This can only happen when `i` is equal to n! or greater. Your index `i` is zero-based. You could catch that case by making `i = mod(i, n!) ` first, but it's more likely a badly calculated input. Negative numbers will also make `j /f + 1` out of bounds. Can you debug the code and find out where bad value was passed? — M Oehm, Mar 22 '16 at 18:57
M Oehm here I show the issue. https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/621702#comment-1865326 Thank you very much — Ilario De Vincenzo, Mar 24 '16 at 12:38
I hope you don't mind if I don't debug your progarm for you. When I answered, I addressed only the factorial part of your code. Frankly, I have a hard time to understand your program, with all the long half English, half Italian variable names and with all words and lines run together without spaces. — M Oehm, Mar 24 '16 at 16:10
As for finding the problem that the permutation id exceeds n!, I can only suggest to debug and to add assertions (additional checks that ensure that the program can't be in an impossible state). You could also try to pass the program though a memory checker to rule out out-of-bound memory accesses. — M Oehm, Mar 24 '16 at 16:11

Operations with big real numbers in Fortran

1 Answers1

Linked

Related