Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
28381 Discussions

Operations with big real numbers in Fortran

Ilario_De_Vincenzo
2,125 Views

Dear all, 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. I hope someone can help me. Thank you in advance.

! 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
0 Kudos
7 Replies
jimdempseyatthecove
Honored Contributor III
2,125 Views

You can increase the precision of your REAL from the default of REAL(4) to REAL(8) or REAL(16).

If (when) you exceed the precision of the mantissa (fraction) portion of the given precision, rounding or truncation occurs. For n=20 (20!) you will require a precision of 19 digits.

REAL(4) 7.22 decimal digits (fraction)
REAL(8) 15.95 decimal digits (fraction)
REAL(16) 34.02 decimal digits (fraction)

Note, not all languages support REAL(16), and those that do, may not support the same internal format (meaning the precision may differ).

Jim Dempsey

0 Kudos
Ilario_De_Vincenzo
2,125 Views

I solved my problem, thanks to your answer and to M Oehm advise, working by integers of KIND=18. 
This is the link where the problem is exposed and fixed http://stackoverflow.com/questions/36136542/operations-with-big-real-numbers-in-fortran

I have not problems when NODI_RETE=12 or less, but have problem for bigger values. I attached my code .f90 here, written in visual studio ultimate 2013, hoping someone can help me to run the program also for NODI_RETE=20, that is my goal. For modify NODI_RETE go to row 482 and change 12 to 20.

All the files needed are attached. Thanks in advance.

0 Kudos
jimdempseyatthecove
Honored Contributor III
2,125 Views

Made some minor changes. Haven't verified if code is completely correct. (for example if the REAL's calculating results below need to be changed from REAL(4) to REAL(8)).

With NODI_RETE=20

 Nodi_20_d_1.000_MinDist_1.000_MaxDist_10.000_Rete_1_M_7_T0_0.022_R_1
 Payoff=0.321 err=3209.407 Consenso=0.162 err=1618.930
 Payoff=0.314 err=0.021 Consenso=0.187 err=0.157
 Payoff=0.318 err=0.013 Consenso=0.212 err=0.131
 Payoff=0.310 err=0.025 Consenso=0.257 err=0.211
 Payoff=0.320 err=0.030 Consenso=0.290 err=0.129
 Payoff=0.311 err=0.026 Consenso=0.364 err=0.255
 Payoff=0.320 err=0.028 Consenso=0.423 err=0.163
 Payoff=0.332 err=0.038 Consenso=0.486 err=0.148
 Payoff=0.337 err=0.015 Consenso=0.588 err=0.211
 Payoff=0.322 err=0.044 Consenso=0.624 err=0.061
 Payoff=0.332 err=0.031 Consenso=0.734 err=0.175
 Payoff=0.341 err=0.027 Consenso=0.762 err=0.038
 Payoff=0.475 err=0.393 Consenso=0.937 err=0.230
 Payoff=0.485 err=0.022 Consenso=0.972 err=0.037
 Payoff=0.511 err=0.053 Consenso=0.991 err=0.020
 Payoff=0.512 err=0.002 Consenso=0.992 err=0.001
 Percorso ottimale individuato:
          17           4          15          18          13          12
          10           2          14           7          16           9
           6           5          19           8          20           3
           1          11
 con lunghezza complessiva di 66.957 contro quella minima di 34.373

 Il blocco di simulazioni richieste e' terminato!

 Vuoi proseguire? Scrivi:
 S per SI
 N per NO

Jim Dempsey

0 Kudos
Ilario_De_Vincenzo
2,125 Views

Thank you very much Jim. You solved with some minor changes my problem. Do you think are there any ways to use NODI_RETE>20 for 64 bit PC? This means to be able to calculate numbers bigger then 20! Best Regards 

 

0 Kudos
jimdempseyatthecove
Honored Contributor III
2,125 Views

Intel Fortran does not have an integer(16). Therefor to stay within the confines of the Fortran language you could change the code that uses integer(8) into real(16) and store the integers in there. Real(16) uses an exponent of 15 bits, sign is 1 bit, leaving 128-16=112 bits for precision. Or 48 bits (well 49 bits) more precision than integer(8). This would permit your factorial routine to handle up to 30!

To go above this, search the internet for multi-precision libraries, there are several available.

(here is one: https://www.researchgate.net/publication/220492964_A_Fortran-90_Based_Multiprecision_System)

As to if changing all INTEGER(kind=long_integer) into REAL(kind=long_kind) is sufficient, I will leave that up to you to test and verify.

Don't forget to change:

N=CEILING(LOG(REAL(factorial(NODI_RETE),kind=long_real))/LOG(REAL(2,kind=long_real)))

to

N=CEILING(LOG(factorial(NODI_RETE))/LOG(REAL(2,kind=long_real)))

Your interest is to avoid overflow of the precision necessary for your computation. I suggest you add some conditional compile statements to insert some asserts to detect overflows).

Jim Dempsey

0 Kudos
mecej4
Honored Contributor III
2,125 Views

The program code in #3 has the following error (at least; there may be others).

The first executable statement in BETAJSTAR uses the value of the INTENT(IN) variable index. This function is invoked in the line

DO WHILE(MAXVAL((/errore1,errore2/))>=errore_ammesso .OR. BetaJstar(A,BETA10,BETAJSTAR_TH,INDEX)<BETAJSTAR_TH)

but at this point the function argument variable index has no defined value.

0 Kudos
jimdempseyatthecove
Honored Contributor III
2,125 Views

Good catch.

Looks like a bug to me.

llario will have to correct for this. Possibly this should be M_INDEX or RETI_INDEX, I do not know.

Jim Dempsey

0 Kudos
Reply