Convertire codice Fortran in Python – 2

Continuo da qui.

Questione preliminare: perché non usare le liste?
Risposta: le liste sono per elementi eterogenei e sono lente e avide di memoria. Se i dati sono omogenei (p.es. tutti int o tutti float o casi simili) ci sono gli array mono- e pluri-dimensionali (matrici).

Python ha un modulo per la gestione degli arrays, ma c’è di meglio: NumPy.

Ogni trattazione di questo argomento inizia con un avvertimento: attenzione alla gestione degli indici, e continua con l’esempio Fortran “per colonna” invece del più comune “per riga”. OK, da verificare. C’è poi un altro aspetto cui prestare attenzione, anzi due. In Fortran il primo indice è 1; in Pthon 0; lo sanno tutti, è banale ma è facile fare confusione, esempio:

$ py3
>>> import numpy as np
>>> ar = np.arange(10)
>>> ar
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> ar.shape
(10,)
>>> ar.size
10
>>> ar[0]
0
>>> ar[5]
5
>>> ar[ar.size]
Traceback (most recent call last):
  File "", line 1, in 
IndexError: index 10 is out of bounds for axis 0 with size 10
>>> ar[ar.size - 1]
9

OK? Quando si ha una matrice attenzione ala gestione degli indici, dicono tutti. Da verificare in Fortran (ebbene sì! ma è semplice e c’è una piacevole sorpresa).

Gestione indici per per riga, file row_col.f

      program row_col

      integer row, col
      real m(1000,1000)

      do row = 1, 1000
          do col = 1, 1000
              m(row, col) = col + 10 * row
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      do row = 1, 1000
          do col = 1, 1000
              m(row, col) = m(row, col) / 2
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      end

Compilo ed eseguo:

$ gfortran -o row-col row-col.f
$ time ./row-col
   11.0000000       220.000000       10001.0000
   5.50000000       110.000000       5000.50000

real	0m0,016s
user	0m0,012s
sys	0m0,004s
$ time ./row-col
   11.0000000       220.000000       10001.0000
   5.50000000       110.000000       5000.50000

real	0m0,018s
user	0m0,014s
sys	0m0,004s
$ time ./row-col
   11.0000000       220.000000       10001.0000
   5.50000000       110.000000       5000.50000

real	0m0,019s
user	0m0,012s
sys	0m0,008s

Uh! per un milione di dati float, real*4 in Fortran 77, il tempo di esecuzione è troppo piccolo per essere significativo. E questa è la versione lenta, ecco quella corretta, col-row.f

      program col_row

      integer row, col
      real m(1000,1000)

      do row = 1, 1000
          do col = 1, 1000
              m(col, row) = col + 10 * row
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      do row = 1, 1000
          do col = 1, 1000
              m(col, row) = m(col, row) / 2
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      end

Compilo ed eseguo:

$ gfortran -o col-row col-row.f
$ time ./col-row
   11.0000000       2002.00000       1010.00000
   5.50000000       1001.00000       505.000000

real	0m0,018s
user	0m0,018s
sys	0m0,000s
$ time ./col-row
   11.0000000       2002.00000       1010.00000
   5.50000000       1001.00000       505.000000

real	0m0,017s
user	0m0,005s
sys	0m0,013s
$ time ./col-row
   11.0000000       2002.00000       1010.00000
   5.50000000       1001.00000       505.000000

real	0m0,017s
user	0m0,013s
sys	0m0,004s

sì, migliore! di quanto? pochissimo!

Non siamo più ai tempi del PDP-11, nemmeno a quelli del VAX suo successore e la RAM oggi è chip e abbondante. Ecco perché il Web è pieno di foto di gatti, mica come ao tempi di Lenna!

OK, questo per il Fortran, entri Python, con Numpy.

Prima versione col-row.py

#!/usr/bin/python3

import numpy as np

M = np.arange(1000*1000, dtype=float).reshape(1000,1000)
for row in range(1000):
    for col in range(1000):
        M[col, row] = col + 1 + 10 * (row + 1)

print(M[0, 0], M[1, 199], M[999, 0])

for row in range(1000):
    for col in range(1000):
        M[col, row] = M[col, row] / 2

print(M[0, 0], M[1, 199], M[999, 0])

Eseguo:

$ time py3 col-row.py
11.0 2002.0 1010.0
5.5 1001.0 505.0

real	0m0,811s
user	0m0,798s
sys	0m0,012s
$ time py3 col-row.py
11.0 2002.0 1010.0
5.5 1001.0 505.0

real	0m0,807s
user	0m0,781s
sys	0m0,020s
$ time py3 col-row.py
11.0 2002.0 1010.0
5.5 1001.0 505.0

real	0m0,798s
user	0m0,776s
sys	0m0,016s

Seconda versione row-col.py

#!/usr/bin/python3

import numpy as np

M = np.arange(1000*1000, dtype=float).reshape(1000,1000)
for row in range(1000):
    for col in range(1000):
        M[row, col] = col + 1 + 10 * (row + 1)

print(M[0, 0], M[1, 199], M[999, 0])

for row in range(1000):
    for col in range(1000):
        M[row, col] = M[row, col] / 2

print(M[0, 0], M[1, 199], M[999, 0])

$ time py3 row-col.py
11.0 220.0 10001.0
5.5 110.0 5000.5

real	0m0,771s
user	0m0,755s
sys	0m0,016s
$ time py3 row-col.py
11.0 220.0 10001.0
5.5 110.0 5000.5

real	0m0,823s
user	0m0,800s
sys	0m0,024s
$ time py3 row-col.py
11.0 220.0 10001.0
5.5 110.0 5000.5

real	0m0,794s
user	0m0,774s
sys	0m0,020s

sì, Python è lento ma ampiamente accettabile.

OK, ho tutto quello che serve per la conversione; con Python sarà semplice aggiornare e aggiungere nuove funzionalità, i grafici, ma non solo. Per esempio capita (ahemm, sempre) di dover fare la matrice trasposta, quella con righe e colonne scambiate. A me che sono vecchio vedere come fa NumPy vengono le vertigini:

$ py3
>>> import numpy as np
>>> M = np.arange(12).reshape(3,4)
>>> M
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> M.T
array([[ 0,  4,  8],
       [ 1,  5,  9],
       [ 2,  6, 10],
       [ 3,  7, 11]])

Il manuale di NumPy è qui.

OK? vero che Python è (quasi) l’unico linguaggio in città. O almeno quello più sexy 🤩

⭕

Convertire codice Fortran in Python – 1

Si può fare ma c’è qualche aspetto da considerare.

Intanto le abitudini di una volta, complice una possibilità del Fortran non presente in Python. Nel Fortran (parlo delle versioni vecchie 77 e precedenti; il 77 è tuttora usato per quel che ne so) ci sono due tipi di sottoprogrammi: le funzioni e le subroutines. Queste ultime, affini alle funzioni C di tipo void erano di gran lunga le preferite. C’è una giustificazione: le variabili sono passate per indirizzo (reference, non so se i due termini sono rigorosamente equivalenti, questi erano quelli correnti). Esempio:

N = 8
CALL DOPPIO (N)
PRINT* N

SUBROUTINE DOPPIO (N)
N = 2 * N
END

Il risultato che ottengo è 16, cioè il valore di N modificato nella subroutine viene conservato nel main, è cambiato il valore memorizzato all’indirizzo della variabile e questo è quello che ho passato al sottoprogramma. Con Python (per quel che ne so) questo non è possibile, ecco qua:

$ py3
>>> def doppio (n):
...   n = 2 * n
...   print('in doppio:', n)
...
>>> n = 8
>>> doppio (n)
in doppio: 16
>>> print (n)
8

alla funzione doppio ho passato il valore 8 e al di fuori di questa il valore di n resta 8 indipendentemente da quanto avvenuto nel sottoprogramma. Ovvia la soluzione, la stesa da sempre, in Fortran:

N = 8
N = IDOPPIO (N)
PRINT* N

FUNTION IDOPPIO (N)
N = 2 * N
RETURN (N)
END

Notare che ho cambiato il nome alla funzione, variabili e funzioni sono di default REAL*4 tranne quelle il cui nome inizia con le lettere dell’intervallo [I..N] che sono INTEGER*2 o INTEGER*4 a seconda della versione e delle opzioni di compilazione. Caratteristico del Fortran è che la precisione viene espressa in byte e non in bit.

In Python  la funzione può ritornare una lista, torna comoda per le subroutines:

$ py3
>>> def pot(a, b, p):
...   a = a ** p
...   b = b ** p
...   return a, b
...
>>> r = pot(3, 4, 2)
>>> print(r)
(9, 16)
>>> # ovviamente
... a, b = pot(3, 4, 2)
>>> print(a)
9
>>> print(b)
16

Quindi una prima regola pratica: scrivere funzioni che ritornano i valori delle variabili che ci servono. Ma ci sono dei casi in cui questa pratica non basta, il caso di arrays e matrici di dimensioni rilevanti.

In Fortran c’è l’istruzione COMMON:

COMMON /DATI/ A, B(10), C, D
COMMON /MAT/ P(20, 20), R(10, 10, 10)

Ogni funzione o subroutine che usi A, B, C o D dovrà dichiarare il COMMON /DATI/, idem per il caso di /MAT/.

In Python si ottiene un risultato simile (meno granulare, meno raffinato ma forse meno soggetto a bug) con le variabili globali:

$ py3
>>> #a = 0
>>> b = []
>>>
>>> def popola_dati():
...   global a
...   a = 42
...   for c in range(10):
...     b.append(c ** 2)
...
>>>
>>> popola_dati()
>>> print('a =', a)
a = 42
>>> print('b =', b)
b = [0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

OK, notare che a devo dichiararla global nella funzione a differenza di b che dev’essere inizializzato nel main e diventa automaticamente globale. b è una lista, non esattamente l’equivalente degli arrays e matrici del Fortran. Le liste sono versatili ma onerose, da modificare, in un post successivo, pausa 😎

⭕

Traduzioni da un linguaggio a un altro

Un discorso che devo prendere un po’ alla lunga, da lontano, roba di prima di voi giovani, ma poi arrivo al dunque, cioè all’oggi, promesso.

In Fortran π, pi, il rapporto tra le misure della circonferenza e il raggio non è definito. Ma tutti sanno che l’arcotangente di 1 rad vale π/4.
E allora diventa tutto semplice, gli antichi (non gli antichissimi ma non cambia molto) usavano il Fortran 77 (1978).
Vediamo se ricordo bene perché possono esserci diversi modi di scrivere, tenendo presente che i tipi sono impliciti e allora

      pi1 = 4.0 * atan(1.0)
      pi2 = 4.0 * atan(1)
      pi3 = 4 * atan(1.0)
      pi4 = 4.0 * datan(1.0)
      
      print *, pi1, pi2, pi3, pi4
      end

Partito molto male: 2 errori –no uno solo in 2 varietà– correggo.

La riga pi2 = 4.0 * atan(1) è irrecuperabile, a correggerla diventa uguale alla precedente.
La riga pi4 = 4.0 * datan(1.0) è più intrigante. Intanto la D iniziale sta per double precision, quello che in C si chiama double.

      pi1 = 4.0 * atan(1.0)
      pi3 = 4 * atan(1.0)
      pi4 = 4.0 * datan(1.d0)
      pi5 = 4.0d0 * datan(1.d0)
      pi6 = 4.0d0 * datan(dble(1.0))
      
      print "(f12.10)" , pi1, pi3, pi4, pi5, pi6
      end

OK, ottengo gli stessi valori. Notare la promozione a double nella moltiplicazione per il calcolo di pi4. Notare anche la funzione dble, cast a double direbbe il C-ista.

Ma non va come previsto: otteniamo sempre lo stesso numero di cifre corrette, sia usando i real*4 che i real*8 corrispondenti a float e double del C.
Devo dichiarare il tipo quando non rientra in quelli di default. Infatti questo viene automaticamente assegnato a integer (2 bytes nel Fortran IV, 4 nel 77) se la prima lettera della variabile è nell’intervallo [I-N], real*4 per tutte le altre. Si ricorda inoltre che al di fuori delle stringhe quotate si ha una conversione automatica al maiuscolo.
Ecco:

      double precision pi4
      real * 8 pi6

      pi1 = 4.0 * atan(1.0)
      pi3 = 4 * atan(1.0)
      pi4 = 4.0 * datan(1.d0)
      pi5 = 4.0d0 * datan(1.d0)
      pi6 = 4.0d0 * datan(dble(1.0))
      
      print "(f12.10)", pi1, pi3, pi4, pi5, pi6
      end

OK 😀 Ah! dimenticavo: gli spazi non sono significativi, real*8 e real * 8 sono sinonimi, lo sarebbe pure r e a l * 8.

Lo so che nessuno usa più il Fortran, neanch’io, ma devo raccontare questo altrimenti poi non si capisce.

Neanche e, quella di Euler, dei logaritmi naturali, è predefinita. Ma facciamo come per π:

      real * 8 e2

      e1 = exp(1.0)
      e2 = dexp(1.d0)      

      print "(f12.10)", e1, e2
      end

OK 😀 notare che al posto di 1.0 posso scrivere 1. come posso scrivere .1 per 0.1; i vecchi lo facevano sempre.

Finora niente di sconvolgente ma ecco un passo ancora

      e = exp(1.0)
      pi = 4.0 * atan(1.0)
      ie2 = 2 * e
      ipi2 = 2 * pi
        
      print *, e, pi
      print *, ie2, ipi2
      end

OK, i cast a integer*4 tronca i real*4. Ma attenzione, questo avviene –ovviamente– dopo la moltiplicazione. Ovviamente è diverso da questo caso:

      e = exp(1.0)
      pi = 4.0 * atan(1.0)
      ie = e
      ipi = pi
      ie2 = 2 * ie
      ipi2 = 2 * ipi
        
      print *, e, pi
      print *, ie2, ipi2
      end

Se si vuole l’arrotondamento all’intero più vicino ci sono due modi:

      e = exp(1.0)
      ie1 = e + 0.5
      ie2 = int(e + 0.5)
      ie3 = nint(e)

      print *, ie1, ie2, ie3
      end

Siccome la conversione da real a integer avviene automaticamente per troncamento i vecchi usavano sommare 1/2 al valore da convertire. La funzione int() in questo caso è implicita per ie1 e esplicita per ie2. Esiste la funzione di arrotondamento nint() ma è usata raramente (mai).

Finora ho barato; no, non barato ma omesso di dire tutta la verità. Alcune funzioni intrinseche del Fortran sono da sempre (OK, quasi) polimorfe.
Già vista datan() ma risulta più chiaro con la funzione max():

      m1 = max(3, 1, 4, 1, 5)
      m2 = max0(3, 1, 4, 1, 5)
      m3 = max1(3.0, 1.0, 4.0, 1.0, 5.0)
      rm2 = amax0(3, 1, 4, 1, 5)
      rm3 = amax1(3.0, 1.0, 4.0, 1.0, 5.0)
      print *, m1, m2, m3
      print *, rm2, rm3
      end

Sì, conta la lettera iniziale (A per i real, D per i double) e l’eventuale numero finale (0 per integer, 1 per real). Come se non fosse abbastanza ci sono state variazioni (tante) tra le versioni 77 e 90 (e successive), RTFM (cit.).

Ci sarebbe da raccontare molto su common ma salto.
Come pure salto tutto il capitolo sulle matrici, c’è Octave visto recentemente.

Invece una cosa sempre in agguato: il passaggio di variabili ai sottoprogrammi. In Fortran ci sono le function, funzioni come quelle di C e Python e quasi tutti i linguaggi. Ma ci sono anche le subroutines –anzi sono nate prima– sottoprogrammi chiamate con la keyword call che non ritornano un valore come le funzioni ma modificano i dati passati. Tutti, perché vengono passati per indirizzo (questo è vero solo per le versioni non nuovissime ma il codice è spesso piuttosto attempato). Un esempio minimo:

*     main
      n = 5
      call foo(n, nq)
      print *, n, nq
      call foo(n, nq)
      print *, n, nq
      call foo(n, nq)
      print *, n, nq
      end
      
      subroutine foo(k, l)
      k = k + 1
      l = k ** 2
      end

Sembra facile, traduco in Python:

### versione iniziale, non funziona

def foo(k, l):
      k = k + 1
      l = k ** 2

#  main
n = 5
foo(n, nq)
print(n, nq)
foo(n, nq)
print(n, nq)
foo(n, nq)
print(n, nq)

Ehmmm… no; non basta aggiustare la sintassi delle singole righe. Il valore di l | nq dev’essere ritornato espressamente. Inoltre è inutile passarlo tra gli argomenti:

### ancora non funziona

def foo(k):
      k = k + 1
      l = k ** 2
      return l

#  main
n = 5
nq = foo(n)
print(n, nq)
nq = foo(n)
print(n, nq)
nq = foo(n)
print(n, nq)

No! non ci sono ancora: ovvio, gli passo sempre lo stesso parametro n = 5; la modifica all’interno di foo() è locale e in main non viene vista, correggo:

def foo(k):
      k = k + 1
      l = k ** 2
      return l

#  main
n = 5
nq = foo(n)
n += 1
print(n, nq)
n += 1
nq = foo(n)
print(n, nq)
n += 1
nq = foo(n)
print(n, nq)

OK! ma il codice è orrendo! Non voglio vedere niente di simile 😦

Come visto anche per un caso molto banale la traduzione automatica, riga per riga non è possibile. Ci sarebbe da dare la colpa al fortrainer che se avesse fatto le cose per bene… scrivendo la subroutine foo() così:

      subroutine foo(k, l)
      j = k + 1
      l = j ** 2
      end

la variabile k sarebbe stata solo di input –o anche intent(in) per le versioni nuove del Fortran– non sarebbe ritornata modificata e il main avrebbe dovuto essere scritto meglio.

Adesso provate –Gedankenexperiment– ad applicare questo a codice reale 👿 Ma secondo me va fatto. Prima o poi. Non necessariamente da me.

Il Fortran iniziale, I/O, FORMAT

Il titolo è fuorviante, riassumo una conversazione con un collega dei tempi antichi riguardo al post Lo sviluppo iniziale dei linguaggi di programmazione.
Naturalmente riguarda il Fortran l’unico linguaggio che c’era sul ‘puter che usavamo. Il programma postato da Knuth e Trabb Pardo (K&T) è per la versione 0 (zero), forse solo virtuale. La prima versione con cui ho avuto a che fare è la IV, anzi un’evoluzione compatibile con la 66 ma c’era codice scritto per la versione II (1958), trasferito da schede a nastri o dischi o semplicemente ricopiato dal listato.
In ogni caso il programma postato da K&T deve essere aggiornato (ehi! 50 anni!): non solo le funzioni non hanno (dalla preistoria) il nome che deve terminare con F ma non ho nemmeno un lettore di schede 😦
E non saprei nemmeno trovare le schede 😦 e la perforatrice 😦

Però il programma prevede la lettura di dati formattati in modo preciso, cosa che da terminale non è agevole (bisogna essere proprio pistini) ma si può fare un formattatore, come i millemila come si usava una volta.

Ah! un altro problema: la gestione dei file: allora si aprivano da OS, il programma lo dava per scontato e alla fine si chiudevano, sempre da OS. I comandi Prime erano rispettivamente OPEN seguito dai parametri che non ricordo ma roba tipo 1 1 (unità 5 in lettura [aggiornamento: il secondo parametro è da verificare, probamilente non come pensavo]) e CLOSE 1 (unità 5); siccome di solito si dovevano chiudere più unità c’era il comodo CLOSE ALL. Su Prime (come con il DEC) i comandi potevano essere abbreviati, nel nostro caso al posto di OPEN bastava O, era il comando più importante iniziante con quella lettera, e C stava per CLOSE.
Quindi, se ricordo bene, il mio command-file per l’esecuzione del programma FMT sarebbe stato così:

O DATI-FREE 1 1
O DATI 2 2
SEG FMT
C ALL

L’istruzione di lancio SEG FMT valeva nell’ipotesi che il programma fosse stato compilato per la versione segmentata (prevista per i programmi grossi) usata (a sproposito) abitualmente. Però se uno era razionale vedeva che bastava la versione normale, anche più veloce quindi la riga sarebbe stata RUN FMT, RUN abbreviabile in R. OK, mode nostalgia OFF.

La versione corrente usa le normali istruzioni OPEN e CLOSE, non disponibili allora. Per evidenziare il nuovo uso il minuscolo, adesso equivalente al maiuscolo, allora vietato.

Ecco fmt.f

C FMT legge in free-format e scrive formattato

      DIMENSION A(11)
      
      open(1, file='dati-free')
      READ(1, *) A
      close(1)
      
      open(2, file='dati')
      WRITE(2, 2) A
      close(2)
2     FORMAT(6F12.4)

      END

Aggiornamento: per il Prime le prime 4 unità di I/O erano predefinite e riservate per cui avrei dovuto usare numeri più alti, tipicamente 5 e 6 al posto di 1 in OPEN e 2 in WRITE.

E questo è il file di dati dati-free

1 10 100 1000 10000 98765 987654
987654.3 2.7818 3.1415 -42

Compilo ed eseguo

ottenendo il file dati

      1.0000     10.0000    100.0000   1000.0000  10000.0000  98765.0000
 987654.0000 987654.3125      2.7818      3.1415    -42.0000

Uh! visto il secondo dato della seconda riga? è 987654.3125 invece del previsto 987654.3000 😳 ma il formato REAL*4 (quello usato per le variabili inizianti con una lettera fuori dell’intervallo [I-N] riservato agli interi) garantisce 7-8 cifre significative.

Per essere proprio sicuri-sicuri si può fare la verifica, con il file verifica.f

C VERIFICA

      DIMENSION A(11)
      
      open(1, file='dati')
      READ(1, 2) A
      close(1)
2     FORMAT(6F12.4)

      DO 3 I = 1, 11
        PRINT *, A(I)
3     CONTINUE     

      END

OK, ci siamo. A questo punto torno al programma di K&P, opportunamente aggiornato (tpk.f).

C     THE TPK ALGORHTM, FORTRAN STYLE
*     aggiornato, modifiche in minuscolo

      FUN(T) = SQRT(ABS(T)) + 5.0**3
      DIMENSION A(11)
1     FORMAT(6F12.4)
      open(1, file = 'dati') 
      READ(1, 1) A
      close(1)
      DO 10 J =1, 11
      I = 11 - J
      Y = FUN(A(I+1))
      IF(400.0 - Y) 4, 8, 8
4     write(*, 5) I
5     FORMAT(I10, 10H TOO LARGE)
      GO TO 10
8     write(*, 9) I, Y
9     FORMAT(I10, F12.7)
10    CONTINUE
      STOP 52525
      end

OK 😀

C’è ancora una cosa emersa ieri: le funzioni, quando non ci sono si fanno. A dire tutta la verità i fortrainers preferivano –di gran lunga– le SUBROUTINEs ma questo sarebbe un altro rant, lunghissimissimo 😉

Era tipico di IBM fornire funzioni che non si avevano (no, niente nomi). Per gioco supponiamo non fossero definite ABS e SQRT.
La prima è banale (e lasciata come esercizio) mentre la seconda è leggermente più impegnativa. Siccome il mondo è vario da adesso in poi continuo con Python 3.x, il Basic dei nostri giorni.

Invece di limitarci alla radice quadrata conviene probabilmente essere più generali; che ne dite di usare i logaritmi? OK, la Wiki (come faremmo senza?) ci dice tutto. Illo tempore probabilmente avrei cercato sul Dwight –hey! c’è online–, c’è il capitolo Logarithmic Functions, p.130

Per l’esponenziale c’è Exponential Functions a p.125

Quindi ecco myfunc.py:

# modulo per Python 3.x; 
# per 2.x importare division from tuture

my_e = 2.71828182845904523536 

def my_log(x):
    r = (x - 1) / (x + 1)
    d = r
    v = d
    c = 1
    while d > 1e-8:
        c += 2
        t = 1
        for i in range(1, c + 1):
            t *= r
        d = 1/c * t
        v += d

    return round(2 * v, 8)

def my_exp(x):
    fact = 1
    t = x/fact
    e = 1 + t
    num = x
    c = 1
    while t > 1.e-8:
        c += 1
        fact *= c
        num *= x
        t = num / fact
        e += t
    st = str(e)
    ppoint = st.index(".")
    if ppoint < 8:
        e = round(e, 8 - ppoint)
    else:
        e = round(e, 1) 

    return e
    
def my_sqrt(x):
    return my_exp(1/2 * my_log(x))

Variabili: uso e abuso

Un tweet di un programmatore funzionale che riporta un tweet che –OK, è qui.

OK, pythonisticamente è corretto. [2] è una lista contenente l’intero 2.

In Python, come in taaanti linguaggi ‘normali‘ le variabili sono –come dice il nome– variabili e possono essere ridefinite.
Ma forse Python esagera, consideriamo il ciclo seguente:

Sembra tutto OK, vero?
Ma forse… approfondiamo:

Uhmmm… qui l’indice del ciclo viene modificato all’interno del ciclo stesso ma ciò non turba l’interprete che al giro successivo sa di dover usare l’elemento successivo della lista. La stessa cosa capita se usiamo range(), che infatti produce la lista usata dal ciclo:

OK 😀

Fin qui è tutto abbastanza scontato, tranne per i nuovi che usano solo linguaggi strettamente funzionali (un po’ li invidio) ma a me che sono vecchio è tornata in mente una cosa di tanto tempo fa. Adesso ve la racconto.
Illo tempore ogni computer aveva il suo sistema operativo, proprietario e di solito diverso da tutti gli altri. Anche i compilatori, di sicuro il Fortran e pochi altri. Quando è nato il PC IBM con il DOS di Microsoft è successo che (numero dei sistemi operativi)++, nient’altro; e peraltro tutt’altro che buono. Ma sto divagando.

Riprendo l’esempio che appalla i funzionali, in Fortran 77, quello corrente di allora si ha (uso i al posto di x per le convenzioni implicite del linguaggio):

* c0.f
      i = 1
      do 10, i = 2, 4
          print *, i
10    end do
      print *, i      
      
      end

OK, a parte s/for/do/ e * per dire ‘default’ in print tutto come previsto.

Occorre anche considerare una cosa che tutti i fortrainers davano per scontato: all’uscita del ciclo il contatore è incrementato. Ma è inaffidabile, tra l’altro varia rispetto alla versione precedente del compilatore (Fortran IV).
Vediamo se si può, come in Python, …

* c1.f
      i = 1
      do 10, i = 2, 4
          i = 2 * i
          print *, i
10    end do
      print *, i
      
      end

No, nope! come previsto il compilatore GNU segnala la violazione.

Ci può essere un caso più sottile, quello capitatomi (io ero il debugger):

* c2.f
      i = 1
      do 10, i = 2, 4
          call doppio(i)
10    end do
      print *, i
      
      end
      
      subroutine doppio(n)
      n = 2 * n
      print *, n
      end

OOPS! 😳 cos’è successo? I fortrainers lo sanno: le variabili vengono sempre passate per indirizzo a funzioni e subroutines. In questo caso l’indice del ciclo viene modificato all’interno del ciclo stesso ma il compilatore non può saperlo. È cura del programmatore –che queste cose le dovrebbe sapere– diminuire il carico di lavoro al debugger, sempre che il bug venga rilevato, possibilmente presto 😀

Al programmatore era stato senz’altro stato detto anche “cambia il nome alle variabili” che è una frase ce dev’essere interpretata e tradotta in linguaggio normale. Vediamola applicata:

* c3.f
      i = 1
      do 10, i = 2, 4
          call doppio(i)
10    end do
      print *, i
      
      end
      
      subroutine doppio(k)
      n = k
      n = 2 * n
      print *, n
      end

OK 😀 C’era un altro motivo per cambiare le variabili, cosa che oggi probabilmente ha meno senso: le variabili locali sono near pointers mentre le globali sono far pointers. Non credo di dire cose nuove, ma forse dimenticate al momento opportuno 😉

Breve paragone tra l’“overloading” in C++, C(11) e Fortran(90)

In C++ possiamo scrivere codice come questo:

#include <iostream>
 
int miafunz(int a, int b)
{
    std::cerr << "int, int\n";
    return a + b;
}
 
int miafunz(double a, int b)
{
    std::cerr << "double, int\n";
    return 2.0*a + 2.5*b;
}
 
int main()
{
    double d0 = 1.0;
    int di = 2;
    std::cout << miafunz(d0, di) << "\n";
    std::cout << miafunz(5, di) << "\n";
}

Lo compiliamo e lo eseguiamo:

Quello che accade dietro le quinte è che il compilatore in realtà genera due funzioni diverse e “capisce” qual è quella giusta da chiamare basandosi sul tipo degli argomenti.

Per il programmatore è tutto trasparente, mentre compilatore e linker hanno una gatta da pelare in più. Il problema è questo: le due funzioni sono in realtà diverse, pezzi di codice diverso, e ciascuna deve avere associato un simbolo diverso¹, anche se il programmatore usa solo il nome che ha deciso (miafunz).

Se non è definita e imposta un’ABI specifica per tutte le implementazioni del C++, ci possono essere “incomprensioni” perché ogni compilatore “decora” i nomi delle funzioni come meglio crede. (Cfr. name mangling)

Detto altrimenti: il compilatore C++ si deve inventare due nomi per le due funzioni, laddove noi usiamo un unico nome. È abbastanza ragionevole ipotizzare che i nuovi nomi siano costruiti a partire dal nome da noi scelto, secondo un certo schema² — ma non è affatto necessario che sia così. Nel caso di g++ (4.7.2),

ed è lo stesso con clang++ (3.0), per la cronaca. È buono che ci sia una certa convergenza, ma purtroppo in realtà non possiamo farci affidamento. Invece di spiare l’ELF potete dire al g++ di fermarsi alla generazione del codice assembly (opzione -S) e dare uno sguardo ai simboli usati.

Il C non ha di questi problemi

Il C non ha di questo problemi perché non permette l’overloading. Se scrivete

#include <stdio.h>

int miafunz(int a, int b)
{
    fprintf(stderr, "int, int\n");
    return a + b;
}

int miafunz(double a, int b)
{
    fprintf(stderr, "double, int\n");
    return 2.0*a + 2.5*b;
}


int main(void)
{
    double d0 = 1.0;
    int di = 2;
    printf("%d\n", miafunz(d0, di));
    printf("%d\n", miafunz(5, di));
    return 0;
}

il compilatore si ribella,

perché stiamo ridefinendo una funzione (ho usato gcc 5.2 perché evidenzia meglio gli errori, in stile clang). In C i simboli delle funzioni devono essere unici (all’interno del loro campo di visibilità³).

Potremmo scrivere

int miafunz_ii(int a, int b)
{
    fprintf(stderr, "int, int\n");
    return a + b;
}

int miafunz_di(double a, int b)
{
    fprintf(stderr, "double, int\n");
    return 2.0*a + 2.5*b;
}

ma poi starebbe a noi scegliere la funzione “giusta” per gli argomenti giusti. Fattibile, ma potenziale fonte di errori e comunque appesantisce la lettura (e la scrittura…). Un esempio più concreto è quello di certe funzioni matematiche per le quali in effetti nella libreria esistono diverse versioni a seconda del tipo di argomenti usati e restituiti; p.es. sin (double), sinf (float), sinl (long double).

Il C11 è venuto incontro a questa esigenza tramite una macro, _Generic, il cui uso però risulta limitato, cioè non dà la stessa elasticità dell’overloading del C++: in pratica è utile solo per tgmath.h (type generic math).

Riciclo un esempio da un mio post su un altro blog.

• generic.c

Il Fortran è diverso

Quando le persone sentono Fortran credono di essere tornati indietro nel tempo: è un linguaggio vecchio, ormai in disuso, è brutto, ecc… Non è così. Di sicuro è “di nicchia”. Ma comunque, a differenza di altri linguaggi di quell’epoca, è andato avanti e si è ammodernato. L’ultima versione dello standard è il Fortran 2008, ma è pianificato il Fortran 2015.

Una delle caratteristiche del Fortran moderno sono le interfacce. In pratica il Fortran ha lo stesso problema che il C(11) ha in parte risolto con _Generic; però la soluzione è molto più… generica!… elastica, elegante e (a mio immodesto avviso) intelligente.

Dettagli sintattici a parte, il programmatore non lascia al compilatore l’onere di generare i nomi delle funzioni: egli scrive le diverse funzioni come deve (p.es. miafunz_di, miafunz_ii) e specifica che miafunz è un’interfaccia a quelle funzioni (il compilatore di solito è nelle condizioni di poter dedurre da sé la funzione giusta in base al tipo degli argomenti, proprio come avviene per il C++; nei casi in cui non dovesse essere così, nell’interfaccia dovremmo specificare per intero gli argomenti e i loro tipi).

Riciclo di nuovo un esempio da un altro blog; nell’esempio my_sum è definita come interfaccia per le funzioni my_sum_r, my_sum_i, ecc…

• generic.f90
• my_sum_module.f90

---

1. A codice oggetto generato, non è sempre necessario. Per esempio nell’esempio specifico in teoria non c’è nessun bisogno che i simboli delle funzioni rimangano, perché di fatto non ci interessa esportarli. Tant’è che tali simboli possono essere rimossi, per esempio con strip. Il discorso sarebbe diverso se stessimo scrivendo del codice che deve essere linkato (una libreria statica o dinamica, per esempio).↩

2. Una soluzione ragionevole (perché intellegibile) è quella di codificare nel nome della funzione i tipi degli argomenti; vediamo infatti che miafunz(double,int) diventa _Z7miafunzdi e miafunz(int,int) diventa _Z7miafunzii.↩

3. O «confine entro cui sono visibili». È la mia traduzione approssimativa di scope. Non voglio affrontare questioni linguistiche sottili, ma la traduzione “visibilità” non mi sembra del tutto corretta. La visibilità (visibility) è la proprietà di essere visibili (la proprietà di un oggetto di poter essere visto); il termine scope fa riferimento a quanto è “estesa” questa proprietà, cioè fin dove vale o non vale per un determinato oggetto; dunque si riferisce all’“area di pertinenza” di una proprietà (in generale, e della visibilità nello specifico). La traduzione “visibilità” di per sé non veicola l’idea di “estensione”. Prendiamo per esempio la frase «the name must be unique in its scope». La traduzione «il nome deve essere unico nella sua visibilità» non corrisponde all’originale; in italiano siamo costretti ad usare una perifrasi: «il nome deve essere unico» «all’interno dell’ambito in cui è visibile», o «all’interno dell’estensione della sua visibilità», o «entro la portata della sua visibilità», o «nel suo campo di visibilità» (forse l’opzione migliore), o altre espressioni simili. Per quanto mi riguarda, in casi come questi non ho nessun problema ad usare termini inglesi nell’italiano settoriale.↩

Argomenti passati a funzione – 2

Non so se solo noi… (inteso come quelli con cui lavoravo illo tempore) ma in seguito al post Argomenti passati a funzione ho avuto una sessione di ricordi con il mio boss di allora. Riassumo.

Intanto la funzione corretta (cquad.f):

1 2 3 4 5 6 7

* quad restituisce il quadrato del parametro - corretta function quad(x) y = x y = y * y quad = y end

Verifico, con questo main (c1.f):

1 2 3 4 5 6 7 8

* versione corretta #1 x = 5 qx = quad(x) print *, ' il quadrato di', x, ' vale ', qx qx = quad(8.0) print *, ' il quadrato di 8.0 vale ', qx end

Nella funzione si copia il parametro in una variabile locale e si modifica solo questa. È così usabile come previsto nel main; qui a raccontarla tutta inizialmente avevo scritto 8 invece di 8.0 …

Però, mi fa osservare il Signor Capo, la prassi era quella di usare le subroutines, ecco squad.f:

1 2 3 4

subroutine quad(x) x = x * x end

che verifico con il main c2.f:

1 2 3 4 5 6

* versione corretta #2 x = 5 call quad(x) print *, ' il quadrato di 5.0 vale ', x end

OK, quasi. Ma questo era il modo standard di operare (almeno da noi). Notare che call quad(8.0) non avrebbe senso, non ottengo il risultato (sì un errore, invece).

Una googlata è stata molto fruttuosa: ho trovato il manuale di Prime F77 che riporta in due punti un warning a proposito di questo problema. Ma si sa che nessuno legge i manuali 👿

pagina 5-1:

e pagina 5-4

Abbiamo pensato anche a intent ma ci risulta posteriore, Fortran 90. Può darsi che qualche implementazione l’avesse già allora ma dalla ricerca con Google non ho trovato niente. Trovato invece cose interessanti:

• Fortran 77 Tutorial;
• Professional Programmer’s Guide to Fortran77;
• FORTRAN 77 Language Reference (Oracle).

Il commento di Piero Patteri, esaminato con più cura di quanto ho risposto, si riferisce a un altro aspetto ancora, interessa a qualcuno? Devo parlarne? Secondo il Boss il Fortran si usa ancora, anche se molto meno, ci sono concorrenti, tanti. Anche Matlab.

---

L’immagine lassù è del manuale del Fortran. Trovata su Ebay (in realtà KiJiJi), con questa nota: Nave scuola per migliaia di programmatori, questo volume del 1981 (riedizione di quella del 1971) ha costituito un punto di riferimento per la nascente attività di programmazione per docenti e ricercatori.
E pare sia ancora nel catalogo di Zanichelli. No, solo nel catalogo storico.

Argomenti passati a funzione

L’altro giorno ero alle prese con la filosofia della programmazione funzionale, cosa di moda ultimamente ma –scopro– vecchia ormai. E mi è tornato in mente un caso di tento tempo fa; quasi-quasi lo racconto. Anzi vado, così verifico se nel frattempo i compilatori sono migliorati. E pensa te che ho dovuto reinstallare il Fortran. Dai c’è di peggio. Per esempio la politica. Il caso è molto semplice. Pochissime righe di codice, antico; qui uso il 77 ma andrebbe bene anche il IV (1966). L’unica cosa da tener presente è che gli argomenti a una funzione sono sempre passati per indirizzo (reference). Quindi quanto si torna dalla chiamata se le variabili usate come parametri sono state modificate risulteranno modificate. Mi prendo una piccola libertà per semplificare il problema: i programmatori allora prediligevano di gran lunga le subroutine, sottoprogrammi che non ritornano esplicitamente un valore ma si affidano ai side effects; invece qui uso una function, quelle del C, c’erano da sempre anche in Fortran. Creo una piccolissima (e ridondante per visualizzare il caso in oggetto) funzione quad() che ritorna il quadrato dell’argomento che le viene passato (quad.f):

1 2 3 4 5 6 7

* quad restituisce il quadrato del parametro function quad(x) x = x * x quad = x end

Iessainou! pessima. Se uno programma così lo licenzio a calcinculo. CMQ, la prima riga è un commento, la quinta specifica il valore da ritornare, calcolata nella riga precedente, inutile, anzi dannosa. Inutile in quanto il calcolo poteva farlo sulla stessa riga che definisce il valore della funzione. Questa funzione può essere chiamata da un main, così (bug1.f):

1 2 3 4 5 6

x = 5 print *, 'calcolo il quadrato di X, con X =', x qx = quad(x) print *, ' il quadrato di', x, ' vale ', qx end

Compilo ed eseguo: Oops! anzi atz! 👿 Ma lo sapevo, x è stata modificata nella funzione. C’era (e forse c’è ancora) di peggio, consideriamo questo caso (bug2.f):

1 2 3 4 5 6

x = 5 print *, 'calcolo il quadrato di X, con X =', x qx = quad(5) print *, ' il quadrato di', x, ' vale ', qx end

Compilo, quad è già compilata, devo solo linkarla, risparmio tempo, i vecchi a queste cose erano attenti 😉 ed eseguo: OK, ecco. Come previsto. A questo punto via con il debug. Cioè no, anzi tutto come previsto. Ma allora… 👿

L’avvenire del Fortran e i suoi successori

Forse mi sono lasciato prendere la mano, non so se sia così necessario questo post. Ma siccome l’ho scritto lo pubblico 😕
Tutto è partito da qui: Why Scientists Are Still Using FORTRAN in 2014 che rimanda a Ars technica, qui: Scientific computing’s future: Can any coding language top a 1950s behemoth? dove troviamo

But almost universally, the language in which these simulation codes are written is Fortran, a relic from the 1950s.

e

Wherever you see giant simulations of the type that run for days on the world’s most massive supercomputers, you are likely to see Fortran code.

Poi dice che subito dopo viene inventato il Lisp. Che non ha successo per la sua stranezza:

Weirdness, because the prefix notation used by Lisp made expressions in that language look a lot less like normal mathematics than did math rendered in Fortran. (Fortran stood, after all, for FORmula TRANslator.) A normal chemist or engineer is far more comfortable with y = (a + b)/c than with (setf y ((/ (+ a b) c))).

Ok, l’ho sentita tante volte che mi verrebbe voglia di ignorarla. Invece no! Perché è (quasi) tutto lì. Davvero, ve ne faccio il riassunto.

(setf y (/ (+ a b) c))

Parto dalla parentesi più interna: (+ a b); questa è una lista il cui primo elemento è una funzione. Questa funzione opera sui successivi elementi della lista e li somma. Questa somma (chiamiamola A+B viene restituita: Lisp restituisce sempre l’ultima espressione calcolata.
Quindi la nostra espressione diventa:

(setf y (/ A+B c))

dove (/ A+B c) è una lista in cui il primo elemento è la funzione / che restituisce il rapporto del secondo elemento della lista per tutti i successivi.

Il valore restituito viene assegnato alla variabile y, tramite setf (setf sta per set (quote e sì, si può scrivere anche setq o set ' –non sono esattamente la stessa cosa ma il risultato è lo stesso).

Uh! visto? il Lisp è uniforme. Si comporta sempre allo stesso modo: il primo elemento della lista (chiamato CAR o first) indica come trattare il resto della lista (CDR o rest) e restituisce il risultato.
Per cui il Lisp è fatto da liste dentro liste dentro liste ancora.  Tutto lì (ok, sto semplificando). Non so voi ma a me sembra molto semplice.

Per contro il Fortran riproduce la sintassi matematica che si usa a scuola:

y = (a + b) / c

Sì, lo ammetto, siamo abituati così e per un calcolo semplice come questo meglio il Fortran. Ma per un caso non banale si vede la differenza: con il Lisp ho una concatenazione di liste e con il Fortran una serie di assegnazioni (attenzione il segno = si deve leggere “diventa“, dice il Wirth (che scrive :=)).

Poi una cosa ancora: io sono vecchio e non so nemmeno se si usano ancora ma quando frequentavo il Poli, roba anni ’70, i migliori usavano le calcolatrici HP, con RPN —reverse polish notation— quasi come il Lisp (cioè no, esattamente il contrario), nel nostro caso si sarebbe scritto

a
b
+
c
/

Con risparmio di digitalizzazione di tasti –non servono le parentesi–, si consumano meno i diti e –sopratutto– i risultati intermedi sono disponibili nello stack (mi dicono che dovrei usare il termine “pila”).

Ma tornando alla programmazione allora quando è uscito il Fortran era una cosa nuova di pakka, un compilatore per linguaggio facile per gli umani, alternativo alle criptiche istruzioni assembly. E sui calcolatori computer c’era solo il Fortran. Molto diverso da quello di oggi, istruzioni molto elementari, quasi da assembly. E tanti GOTO, spesso mascherati come IF.

Ok, fine dello sfogo, torno al post di Ars Technica secondo il quale ci sono tre candidati a sostituire il Fortran:

Haskell, bello, sintassi da matematici, funzionale (cioè odia le variabili), chissà;
Clojure, Lisp & JVM, sì c’è Java sotto; ho scritto qualcosa in passato;
Julia, nuovissimo, ne so molto poco.

Ma ce ne sono altri, secondo me: il C++, Java, Python e altri ancora. Tutti ampiamente utilizzati.
E i nuovissimi Go e Rust. Di Go sono a conoscenza di un programma funzionante, per il Web, Rust promette bene (ne parla spesso robitex su questo blog).
Inoltre proprio ieri vosto questo post: #golang: The Next Great Teaching Language. Il blog è nuovo ma promette bene. E il post è illuminante.

Invece il Fortran non è amato dai giovani (quelli che conosco), lo trovano vecchio, inusuale. E se possono lo sostituiscono con qualcosa di più pratico, per esempio Matlab (a volte basta la versione free Octave). Ma sto parlando di chi usa il computer come un mezzo complementare. E lo usa per AutoCAD, fogli di calcolo Excel e come macchina da scrivere.

Mi mancano i programmatori di una volta, quando il mondo era nuovo 😉

Precisione iii

Corto circuito

Questo post è un bellissimo corto circuito: conosco Orlando perché è un amico del GuIT che tra l’altro ho conosciuto l’anno scorso di persona al GuIT meeting di Napoli. E proprio Orlando ha commentato il post di Juhan Precisione su questo stesso blog, di qualche giorno fa su alcuni calcoli.
Questo dimostra che gli appassionati non vedono l’ora di scoprire nuove cose sul proprio argomento preferito.

Chiudo il circuito

Eccomi dunque a chiudere il circuito informatico presentrando un programma in Go che esegue le stesse operazioni eseguite in awk ed in Fortran:

package main

import (
    "fmt"
    "math/big"
)

func main() {
    s, u := big.NewInt(2), big.NewInt(1)
    for i := 1; i <= 10; i++ {
        // a = s - 1
        a := big.NewInt(0).Set(s).Sub(s, u)
        // s*a + 1 = s*(s - 1) + 1
        s.Mul(a, s).Add(s, u)
        fmt.Println(i, s)
    }
}

Tecnicamente, vorrei farvi notare che è possibile in Go concatenare le chiamate dei metodi come per esempio si legge alla riga 12 del codice sorgente, se il metodo stesso ritorna lo stesso oggetto.

Ho utilizzato il pacchetto dei grandi numeri big necessario perché dopo la sesta iterazione i numeri superano il limite massimo del tipo int64. Certo avrei potuto utilizzare il tipo uint64 ma potevo chiudere il circuito soltanto esagerando alla grande portando il ciclo fino a 10 iterazioni… ecco di seguito il risultato:

1 3
2 7
3 43
4 1807
5 3263443
6 10650056950807
7 113423713055421844361000443
8 12864938683278671740537145998360961546653259485195807
9 165506647324519964198468195444439180017513152706377497841851388766535868639572406808911988131737645185443
10 27392450308603031423410234291674686281194364367580914627947367941608692026226993634332118404582438634929548737283992369758487974306317730580753883429460344956410077034761330476016739454649828385541500213920807

A questo punto, ci sta propio bene un
Alla Prossima!
R.