Community
cancel
Showing results for
Did you mean:
Highlighted
Beginner
17 Views

## Operations with big real numbers in Fortran

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)
j=MOD(j,f)
ENDDO
ENDSUBROUTINE ith_permutazione

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

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
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), 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))
ENDFORALL
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))
ENDFORALL
array_perceived_payoff_neg(k)=(1./SUM(lunghezze_percorso))
ELSE
array_perceived_payoff_neg(k)=0.
ENDIF
ENDDO
ENDSUBROUTINE tuttipayoff```
Tags (1)
7 Replies
Highlighted
Black Belt
17 Views

## You can increase the

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

Highlighted
Beginner
17 Views

## I solved my problem, thanks

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.

Highlighted
Black Belt
17 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

Highlighted
Beginner
17 Views

## Thank you very much Jim. You

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

Highlighted
Black Belt
17 Views

## Intel Fortran does not have

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.

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

Highlighted
Black Belt
17 Views

## The program code in #3 has

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.

Highlighted
Black Belt
17 Views

## Good catch.

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