Attenzione alle particolarità

Qualche tempo fa io e un mio giovane collega siamo stati vittime di una disattenzione su una specifica particolare del Fortran. Colpa mia che non ho avvisato del tranello o colpa del comunque bravo P che sul Fortran era nuovo? Adesso vi racconto e giudicherete voi.

Supponiamo di avere una matrice bidimensionale, anzi faccio un caso pratico: la Tavola Pitagorica. A proposito lo sapete che si chiama così solo in italiano? Nel resto del mondo diventa qualcosa come tabella delle moltiplicazioni, ma sono fuori tema.

L’esempio è banale, eccolo, ma attenti al quirk

*234567
      program tavpit

      integer tp(10, 10)

      do 20 j = 1, 10
        do 10 i = 1, 10
            tp(i, j) = i * j
   10   end do
   20 end do
      write(*, *) ' Tavola Pitagorica'
      write(*, 90) tp

   90 format(10(1x, i3))

      end 

Compilo e eseguo

Tutto OK, come previsto 😀

Ma avete visto la particolarità? Dai ve la dico io: a differenza di tutti gli altri linguaggi di programmazione il Fortran ordina gli elementi di un array per colonna (column major order) invece che per riga (row major order).

Nell’esempio si vede come l’indice esterno i varia più velocemente dell’indice interno j. Esattamente all’opposto del C o Python o Basic o quello che volete voi. Nel caso in oggetto non ci sarebbero problemi a usare gli indici invertiti ma anticamente con le macchine lente e poca memoria il tempo di esecuzione sarebbe variato enormemente.

A noi è capitato di far scrivere la matrice dal Fortran e leggerla (conservando gli stessi nomi) con Python con il risultato che veniva ribaltata rispetto alla diagonale principale. E siccome la matrice non era simmetrica come nell’esempio c’è stato un momento di panico.

Colpa di chi? Mia o di P?
O piuttosto del Fortran, sì dev’essere così.

e Nepero e Eulero e il Fortran

Sto leggendo (quasi) un libro bellissimissimo di cui farò una recensione come si deve quando lo finisco. Ma per intanto ecco un’idea che mi ha suggerito, a p.145: e di Nepero e/o di Eulero e come è definito e come calcolarlo.

Per chi fosse curioso al riguardo c’è questo post: 2,7182818284590, numero di Nepero o di Eulero?

OK, tornando a noi, cioè a e ho pensato che è un’ottima occasione, unica, per il Fortran. Adesso vi spiego.

Si è sparsa in giro la voce che io sono un fortrainer, vero: quando io ho cominciato a lavorare sui computer normali c’era solo quel compilatore. Recentemente ho ricevuto due segnalazioni sul Fortran (non riesco a ritrovarle) sul fatto che questo produca il codice più veloce. Può darsi e forse posso anche dire qualcosa al riguardo.

Non appartengo alla prima generazione di programmatori Fortran, sono troppo giovane. I Primi Programmatori Fortran (PPF) hanno come ancestori i Veri Programmatori (VP) quelli che programmavano solo in assembler (si dovrebbe dire Assembly ma la parola “Assembly” si usa solo nella frase “si dovrebbe dire Assembly”).

La verità è che i VP si estinsero quasi completamente (con l’eccezione di piccole nicchie) alla comparsa del FORTRAN e della primissima generazione di proto PPF provenienti dai VP con l’inserimento di nuove unità.
Era decisamente più facile la scrittura e la manutenzione il codice scritto in FORTRAN e allora l’evoluzione…

Quando io sono entrato in gioco i VP si erano estinti (almeno localmente, sì lo so che l’assembler del PDP11 è stato per lungo tempo tra le domande d’esame temute per gli info, anche quando il PDP11 non c’era più). I PPF ormai erano la specie dominante, non scrivevano ma supervisionavano, c’era ormai la seconda generazione, il FORTRAN si era evoluto, c’era stato il IV (1966). IBM che allora era quasi egemone (c’era anche Digital, certo) usciva ogni tanto con un aggiornamento indicato da una lettera e diventava quello lo standard di fatto. Nel ’78 è uscita la versione classica e definitiva (quasi) del Fortran, il 77. Io avevo una formazione IV ma poi mi sono adeguato al nuovo standard, i VPF (veri programmatori Fortran) no! Hanno continuato a tirar dritto per la retta via. Come in ambito scientifico si è continuato a usare il latino per tutto il ‘700 (e oltre).

Tutta questa premessa ha un senso: sì uso un GOTO ma per rimanere fedele allo spirito del tempo (lo Zeitgeist (chissà se si scrive davvero così) di allora). Sarebbe facile farne a meno ma –dai!– quando ci vuole ci vuole.

Ecco allora calcE.

Calcolo di e

      program calcE

      real*8 e, n, delta, eprec

      eps = 1.D-14
      n = 1
      delta = 10
      eprec = 0
      i = 0
   10 continue
      e = (1 + 1./n) ** n
      write(*, 90) int(n), delta, e
      n = n + 1000
      i = i + 1
      delta = dabs(e - eprec)
      eprec = e
      if (delta .gt. eps) go to 10
      write(*, 91) i, delta

   90 format(I8, 1X, D10.4, 1X, F14.12)
   91 format(/'Il valore accettabile è stato raggiunto con', I6,
     &       ' iterazioni,' / 'con delta = ', D10.4 /)

   	  end

Un vero fortranista avrebbe scritto tutto maiuscolo, io non ce l’ho fatta, perdono. Compilo

Eseguendo si ha (visualizzo solo la parte finale)

Mmmh! mica tanto buono, ecco

Sì lo so ci sono delle librerie, yadda yadda yaddaaa. Ma per oggi basta così.

Sulla ricorsività: perché sì e perché no

Il post Ancora sul sort di ‘Compilati vs. Interpretati’ – 2 ha fatto pensare a qualcuno che il Fortran non sia adatto per certe cose e/o che l’impossibilità di scrivere funzioni ricorsive sia un limite. Ebbene chi parla male del Fortran se la deve vedere con me: #sapevatelo 8)

Di seguito presenterò due script per il calcolo del fattoriale e del numero di Fibonacci, scritti in Python, sia in forma ricorsiva che no.

La definizione del fattoriale e della sequenza di Fibonacci si trovano al solito su Wikipedia. (Nota: lo sapete che ci sono anche scritte come si dovrebbe, qui il fattoriale). 😉

Ecco la versione ricorsiva del fattoriale

#!/usr/bin/ env python
# -*- coding: utf-8 -*-

import sys

def fact(n):
	print 'fact', n, ' : ',
	if n == 0 or n == 1:
		print
		return 1
	else:
		return n * fact(n - 1)

n = int(sys.argv[1])
print fact(n)

È da notare quante volte la funzione fact() venga chiamata. In forma non ricorsiva lo script diventa

#!/usr/bin/ env python
# -*- coding: utf-8 -*-

import sys

def fact(n):
	f = 1
	for c in range( 1, n + 1):
		f = f * c
	return f

n = int(sys.argv[1])
print fact(n)

praticamente un solo ciclo for.

La sequenza di Fibonacci, ricorsiva, è questa

#!/usr/bin/ env python
# -*- coding: utf-8 -*-

import sys

def fib(n):
	print 'fib', n, ' : ',
	if n == 0 or n == 1:
		print
		return 1
	else:
		return fib(n - 1) + fib(n - 2)

n = int(sys.argv[1])
print fib(n)

Per valori di n superiori il numero delle chiamate cresce in modo esplosivo, provate 8)

Ecco la versione newLISP

#!/usr/bin/newlisp

(define (fib n)
	(if (< n 2) 1
		(+ (fib (- n 1)) (fib (- n 2))
		)
	)
)

(set 'n (int (main-args 2)))
(println (fib n))

(exit)

Per il calcolo del fattoriale con newLISP utilizzo una scorciatoia. Dal manuale si ha

(exp (gammaln 6)) → 120
The example uses the equality of n! = gamma(n + 1) to calculate the factorial value of 5.

per cui

#!/usr/bin/newlisp

(define (fact n)
	(exp (gammaln (++ n))
	)
)

(set 'n (int (main-args 2)))
(println (fact n))

(exit)

A raccontarla tutta il fattoriale in Python è definito nel modulo math, per cui si potrebbe usare questa versione

#!/usr/bin/ env python
# -*- coding: utf-8 -*-

import sys, math

n = int(sys.argv[1])
print math.factorial(n)

Quindi benché la ricorsività sia buona e bella non sempre è utile o necessaria. Però piace ai matematici, ecco un blog che una volta seguivo (adesso mi manca il tempo), ricorsivo fin dal nome.

Ancora sul sort di ‘Compilati vs. Interpretati’ – 2

Panico!
Ebbene sì! Ci vuole coraggio per postare dopo l’intervento di ieri di glipari: ha reso Haskell sensato e –diciamola tutta– comprensibile 8)

Ma faccio finta di non accorgermi di niente e continuo con le mie robe 8)
Poi c’è un altro motivo per cui è stato difficile fare queste cose semplici-semplici: difficili da spiegare ai ggiovani –beati loro. Ma se leggete fino in fondo ci sarà la soluzione.
OK, dopo queste doverose premesse, vado e continuo con le mie osservazioni al post di dikiyvolk Compilati VS interpretati passo 1.

Nel post Verifica sul sort di ‘Compilati vs. Interpretati’ – 1 avevo provato a verificare le affermazioni di dikiyvolk (Walter ma in russo non ci sono parole più normali per un nickname sexy?) relative a Python 2.x e 3.x, confermandole in toto.

Ma non possiamo dimenticare newLISP. Perché mi è simpatico, è un mio amico (chissà con Haskell incombente come andrà a finire?) e allora ecco una versione modificata dello script di dikiyvolk, con relativi risultati. Riporto un solo caso rappresentativo, sempre sul solito PC athena.

#!/usr/bin/newlisp

(println "Lettura: " (time (set 'data (parse (read-file "LCbig.txt") "\n" 0))) " ms")
(println "Sort: " (time (sort data)) " ms")
(println "Scrittura: "(time (write-file "LCbig-out.txt" (join data "\n"))) " ms")

(exit)

OK, niente sorprese, newLISP fa il suo lavoro, onestamente. Il piccolino non delude, mai 8)

Ma, e qui vengo al dunque, per tanto tempo, quando ero giovane io, il linguaggio di programmazione era un altro, solo uno, quello.

Sì il Fortran. Si può fare in Fortran? Quasi-quasi ci provo!
Io, per quanto riguarda il Fortran non sono molto aggiornato. Nell’ambito professionale non sono mai andato oltre alla versione del ’78, il Fortran 77. Poi dal ’90 in poi ci sono state nuove versioni ma, almeno fino al 2000 si è continuato con quella. Anche perché i fortranisti macinano numeri e le librerie erano intoccabili (un esempio tra poco).

Allora, primo ostacolo: il Fortran 77 non c’è predefinita una subroutine per il sort; te la devi scrivere o, meglio ancora, trovarne una efficiente scritta da qualcuno più bravo di te. È quello che ho fatto io: googlata per fortran sort e –presto!– ecco:

• http://www.fortran.com/lcpsort.f95
• http://gd.tuwien.ac.at/languages/fortran/freemarket/quick_sort1.f
• http://gd.tuwien.ac.at/languages/fortran/freemarket/quick_sort2.f

Uh! vediamo 8)
lcpsort.f95: sembra bella, recente (2009), viene dall’università, ha riferimenti a tesi di dottorato, intrigante ma ha un difetto (per me): mi toccherebbe aggiornarmi sulle specifiche del Fortran moderno. No dai, vediamo se c’è qualcos’altro. Però è davvero scritta bene, mi piacerebbe lavorare con gente che scrive così; mi sento un po’ in colpa di non dare nemmeno un occhiata al PDF indicato;

quick_sort1.f: qui siamo sempre in ambito universitario, qualche anno prima (1994). Si tratta della trascrizione in Fortran di una procedura Algol (antenato del Pascal). Siccome la versione originale era in Fortran IV (o 66 come indicato nei commenti, i due termini sono quasi sinonimi) c’è un uso notevole di GOTO. Ha inoltre un’altra pecca: sorta interi e aggiustarla per trattare stringhe non è immediato (ci ho provato). Da notare, verso la fine, il commento fuorviante relativo alla ricorsione: non è vero, in Fortran non c’era e non c’è nella subroutine. Scritta come si usava una volta è difficile farsi un’idea di come funzioni, davvero.

quick_sort2.f: di Leonard J. Moss del SLAC. Fa la stessa cosa del precedente ma con uno stile decisamente migliore (imho). Scritta nell”86 partendo dal Knuth (ottima scelta!). Notare l’uso ragionevole dei GOTO (in Fortran 77 è difficile non usarne mai) e la scrittura ordinata. È immediato modificarla per ordinare un array di stringhe, come risulta dal codice riportato sotto.

Qui sono stato preso da una serie di dubbi; provo a elencarli.

In Fortran (ricordate che sto sempre parlando delle versioni che conosco) non c’è un istruzione per leggere un file in blocco, si deve leggere una riga o un record per volta. In realtà non risulterà un handicap come temevo.

Le stringhe in Fortran sono di lunghezza fissa, la definisci tu ma una volta dichiarata è quella: se leggi un numero minore di caratteri viene riempita con spazi alla fine; se leggi troppi caratteri perdi quelli in esubero alla capacità della stringa.

Problema relativo alla dimensione dei dati da trattare: duecento mega abbondanti! Questo per me è stato uno shock difficile da accettare e raccontare: il mio primo computer aveva mezzo mega di RAM (sì 512 KB di RAM) e un hard disk di 80 MB, ne ho parlato qui, c’è anche la foto. Sigh 8)

In ogni caso i parametri da passare alla subroutine devono essere modificati. Metto in due common i vettori (ai miei tempi gli array si chiamavano così, da noi) dei dati e dei puntatori. Sì, lo so, non sono puntatori, qui non ci sono ma si comportano esattamente come se lo fossero (sto parlando da vecchio fortrainer, sapevatelo).

Ecco con pochissime correzioni (le istruzioni originali sono state lasciare commentandole) il codice

c234567
      Program wsort

      parameter (numlinee = 5 700 000)
      common /datistr/ lines
      common /ordine/ index
      character*102 lines(numlinee)
      dimension index(numlinee)

      open(1, file = 'LCbig.txt')
*      open(1, file = 'T0.txt')
      read(1, 90, err = 10, end = 20) (lines(i), i = 1, 5 700 000)
      close(1)
   10 print*, i
      stop
   20 numlines = i-1
      write(*, 91) numlines

	  call SORTRX(numlines)

	  open(2, file = 'LCbig-o1.txt')
	  do 30 i = 1, numlines
	  	write(2, 90) lines(index(i))
   30 continue
	  close(2)

   90 format(A)
   91 format('Lette ', I10, ' linee')
   92 format(A)

      end

C From Leonard J. Moss of SLAC:

C Here's a hybrid QuickSort I wrote a number of years ago.  It's
C based on suggestions in Knuth, Volume 3, and performs much better
C than a pure QuickSort on short or partially ordered input arrays.

      SUBROUTINE SORTRX(N)
*      ,DATA,INDEX)
C===================================================================
C
C     SORTRX -- SORT, Real input, indeX output
C
C
C     Input:  N     INTEGER
C             DATA  REAL da me trasformato in string
C
C     Output: INDEX INTEGER (DIMENSION N)
C
C This routine performs an in-memory sort of the first N elements of
C array DATA, returning into array INDEX the indices of elements of
C DATA arranged in ascending order.  Thus,
C
C    DATA(INDEX(1)) will be the smallest number in array DATA;
C    DATA(INDEX(N)) will be the largest number in DATA.
C
C The original data is not physically rearranged.  The original order
C of equal input values is not necessarily preserved.
C
C===================================================================
C
C SORTRX uses a hybrid QuickSort algorithm, based on several
C suggestions in Knuth, Volume 3, Section 5.2.2.  In particular, the
C "pivot key" [my term] for dividing each subsequence is chosen to be
C the median of the first, last, and middle values of the subsequence;
C and the QuickSort is cut off when a subsequence has 9 or fewer
C elements, and a straight insertion sort of the entire array is done
C at the end.  The result is comparable to a pure insertion sort for
C very short arrays, and very fast for very large arrays (of order 12
C micro-sec/element on the 3081K for arrays of 10K elements).  It is
C also not subject to the poor performance of the pure QuickSort on
C partially ordered data.
C
C Created:  15 Jul 1986  Len Moss
C
C===================================================================

      parameter (numlinee = 5 700 000)
      common /datistr/ data
      common /ordine/ index
      character*102 data(numlinee)
      dimension index(numlinee)

*      INTEGER   N,INDEX(N)
*      REAL      DATA(N)

      INTEGER   LSTK(31),RSTK(31),ISTK
      INTEGER   L,R,I,J,P,INDEXP,INDEXT
*      REAL      DATAP
      character*102 DATAP

C     QuickSort Cutoff
C
C     Quit QuickSort-ing when a subsequence contains M or fewer
C     elements and finish off at end with straight insertion sort.
C     According to Knuth, V.3, the optimum value of M is around 9.

      INTEGER   M
      PARAMETER (M=9)

C===================================================================
C
C     Make initial guess for INDEX
      DO 50 I=1,N
         INDEX(I)=I
   50    CONTINUE

C     If array is short, skip QuickSort and go directly to
C     the straight insertion sort.

      IF (N.LE.M) GOTO 900

C===================================================================
C
C     QuickSort
C
C     The "Qn:"s correspond roughly to steps in Algorithm Q,
C     Knuth, V.3, PP.116-117, modified to select the median
C     of the first, last, and middle elements as the "pivot
C     key" (in Knuth's notation, "K").  Also modified to leave
C     data in place and produce an INDEX array.  To simplify
C     comments, let DATA[I]=DATA(INDEX(I)).

C Q1: Initialize
      ISTK=0
      L=1
      R=N

  200 CONTINUE

C Q2: Sort the subsequence DATA[L]..DATA[R].
C
C     At this point, DATA[l] <= DATA[m] <= DATA[r] for all l < L,
C     r > R, and L <= m <= R.  (First time through, there is no
C     DATA for l < L or r > R.)

      I=L
      J=R

C Q2.5: Select pivot key
C
C     Let the pivot, P, be the midpoint of this subsequence,
C     P=(L+R)/2; then rearrange INDEX(L), INDEX(P), and INDEX(R)
C     so the corresponding DATA values are in increasing order.
C     The pivot key, DATAP, is then DATA[P].

      P=(L+R)/2
      INDEXP=INDEX(P)
      DATAP=DATA(INDEXP)

      IF (DATA(INDEX(L)) .GT. DATAP) THEN
         INDEX(P)=INDEX(L)
         INDEX(L)=INDEXP
         INDEXP=INDEX(P)
         DATAP=DATA(INDEXP)
      ENDIF

      IF (DATAP .GT. DATA(INDEX(R))) THEN
         IF (DATA(INDEX(L)) .GT. DATA(INDEX(R))) THEN
            INDEX(P)=INDEX(L)
            INDEX(L)=INDEX(R)
         ELSE
            INDEX(P)=INDEX(R)
         ENDIF
         INDEX(R)=INDEXP
         INDEXP=INDEX(P)
         DATAP=DATA(INDEXP)
      ENDIF

C     Now we swap values between the right and left sides and/or
C     move DATAP until all smaller values are on the left and all
C     larger values are on the right.  Neither the left or right
C     side will be internally ordered yet; however, DATAP will be
C     in its final position.

  300 CONTINUE

C Q3: Search for datum on left >= DATAP
C
C     At this point, DATA[L] <= DATAP.  We can therefore start scanning
C     up from L, looking for a value >= DATAP (this scan is guaranteed
C     to terminate since we initially placed DATAP near the middle of
C     the subsequence).

         I=I+1
         IF (DATA(INDEX(I)).LT.DATAP) GOTO 300

  400 CONTINUE

C Q4: Search for datum on right <= DATAP
C
C     At this point, DATA[R] >= DATAP.  We can therefore start scanning
C     down from R, looking for a value <= DATAP (this scan is guaranteed
C     to terminate since we initially placed DATAP near the middle of
C     the subsequence).

         J=J-1
         IF (DATA(INDEX(J)).GT.DATAP) GOTO 400

C Q5: Have the two scans collided?

      IF (I.LT.J) THEN

C Q6: No, interchange DATA[I] <--> DATA[J] and continue

         INDEXT=INDEX(I)
         INDEX(I)=INDEX(J)
         INDEX(J)=INDEXT
         GOTO 300
      ELSE

C Q7: Yes, select next subsequence to sort
C
C     At this point, I >= J and DATA[l] <= DATA[I] == DATAP <= DATA[r],
C     for all L <= l < I and J < r <= R.  If both subsequences are
C     more than M elements long, push the longer one on the stack and
C     go back to QuickSort the shorter; if only one is more than M
C     elements long, go back and QuickSort it; otherwise, pop a
C     subsequence off the stack and QuickSort it.

         IF (R-J .GE. I-L .AND. I-L .GT. M) THEN
            ISTK=ISTK+1
            LSTK(ISTK)=J+1
            RSTK(ISTK)=R
            R=I-1
         ELSE IF (I-L .GT. R-J .AND. R-J .GT. M) THEN
            ISTK=ISTK+1
            LSTK(ISTK)=L
            RSTK(ISTK)=I-1
            L=J+1
         ELSE IF (R-J .GT. M) THEN
            L=J+1
         ELSE IF (I-L .GT. M) THEN
            R=I-1
         ELSE
C Q8: Pop the stack, or terminate QuickSort if empty
            IF (ISTK.LT.1) GOTO 900
            L=LSTK(ISTK)
            R=RSTK(ISTK)
            ISTK=ISTK-1
         ENDIF
         GOTO 200
      ENDIF

  900 CONTINUE

C===================================================================
C
C Q9: Straight Insertion sort

      DO 950 I=2,N
         IF (DATA(INDEX(I-1)) .GT. DATA(INDEX(I))) THEN
            INDEXP=INDEX(I)
            DATAP=DATA(INDEXP)
            P=I-1
  920       CONTINUE
               INDEX(P+1) = INDEX(P)
               P=P-1
               IF (P.GT.0) THEN
                  IF (DATA(INDEX(P)).GT.DATAP) GOTO 920
               ENDIF
            INDEX(P+1) = INDEXP
         ENDIF
  950    CONTINUE

C===================================================================
C
C     All done

      END

La lunghezza massima delle stringhe (102) è stata determinata con questo script Python

#!/usr/bin/ env python
# -*- coding: utf-8 -*-

import time

try:
  # lettura del file
  t0 = time.time()
  print "Tempo inizio: ", time.strftime("%H:%M:%S")
  file_input = open("LCbig-out.txt", "ru")
  lines = file_input.readlines()
  file_input.close()
  t1 = time.time()
  print "Tempo fine: ", time.strftime("%H:%M:%S")
  print "Tempo di caricamento dei dati: ", t1 - t0, "secondi"

  # controllo
  lung = 0
  nlines = len(lines)
  print "controllo", nlines

  t0 = time.time()
  print "Tempo inizio: ", time.strftime("%H:%M:%S")
  for c in range(0, nlines):
	  if len(lines[c]) > lung:
		  lung = len(lines[c])
		  print c, lung
  t1 = time.time()
  print "Tempo fine: ", time.strftime("%H:%M:%S")
  print "Tempo di controllo dei dati: ", t1 - t0, "secondi"
  print "linea più lunga", lung
except(IOError):
  print "Errore nell'elaborazione del file LCbig-out.txt"
  quit()

Notare inoltre come sia possibile scrivere bene i numeri grossi: siccome il compilatore elimina tutti gli spazi (tranne quelli all’interno di stringhe letterali) il numero 5700000 può essere scritto come 5 700 000. Il compilatore converte in maiuscolo tutto il codice (con l’eccezione ricordata prima) per cui io continuo a usare i caratteri minuscoli.
La versatilità del common mi consente di chiamare le linee da ordinare lines nel main e data nella subroutine.
La riga 12 meriterebbe un lungo discorso, gestisce l’errore, l’endfile, contiene quello che in Fortran si chiama DO implicito (DO è quello che tutti chiamano for).

Ecco cosa succede

Il file di output ha queste dimensioni

Troppo grosso, sì perche le stringhe sono paddate con blank fino alla lunghezza di 102 caratteri. Ecco uno script per sblankarle

#!/usr/bin/ env python
# -*- coding: utf-8 -*-

import time

try:
  # lettura del file
  file_input = open("LCbig-o1.txt", "ru")
  lines = file_input.readlines()
  file_input.close()
  t1 = time.time()
  nlines = len(lines)

  # sblank

  print "Tempo inizio: ", time.strftime("%H:%M:%S")
  t0 = time.time()

  for c in range(0, nlines):
	lines[c] = lines[c].strip() + '\n'

  t1 = time.time()
  print "Tempo fine: ", time.strftime("%H:%M:%S")
  print "Tempo di sblank dei dati: ", t1 - t0, "secondi"

  # scrivo
  file_output = open("LCbig-o3.txt", "w")
  file_output.writelines(lines)
  file_output.close()

except(IOError):
  print "Errore nell'elaborazione del file dei dati"
  quit()

Oh! anche il file originale aveva degli extra blank a destra! Modifico lo script e eseguo, ecco il risultato

OK. E questo è quanto. Non avrei mai immaginato di fare una cosa così in Fortran, non credevo fosse possibile, le dimensioni se ci penso mi fanno provare, ancora oggi, una vertigine difficile da descrivere.
In ogni caso certe cose non si devono fare in Fortran. Molto più semplice con Python. E stupendamente semplice e fenomenale newLISP.
Vero? 8) Non so se ci crederete ma scrivere questo post è stato per me difficile, quasi come una seduta psicoanalitica.

Ticol-tikì e l’interfaccia per il Fortran

Occhei, forse non si dice tikì ma tok. E poi si scrive tcl/tk. E quello di cui voglio parlare è wish.
Ma c’è chi l’ha già fatto, meglio di me, li introduco e lascio a loro la parola.
Tcl/Tk Cookbook di Lakshmi Sastry e Venkat VSS Sastry è una raccolta di esempi commentati su tcl/tk, come appunto dice il nome. Non è recentissimo, anzi per il web abbastanza stagionato ma funziona.

Il cookbook merita di essere letto tutto ma, continuando il discorso iniziato nei post precedenti considero in particolare il cap. 7 Tcl/Tk and FORTRAN in cui si sviluppa un piccolo programma per la soluzione delle equazioni di secondo grado in Fortran 77 e gli costruisce un interfaccia in tcl/tk per eseguirlo in finestra con i soliti pulsanti, campi di inserimento e visualizzazione cui siamo abituati.

Funziona tutto, con 2 piccolissime modifiche: l’invocazione del compilatore Fortran dev’essere cambiata da f77 a gfortran e nel file quads.tcl ho modificato

proc invokeQuads { } {
     set f [open |quads r+]

con

proc invokeQuads { } {
     set f [open |"quads" r+]

Inoltre l’opzione -f nella shebang di quads.tcl non è necessaria, l’evoluzione del linguaggio 😉

Naturalmente per lanciarlo il comando sarà

wish quads.tcl

a meno che si sia reso eseguibile con chmod +x quads.tcl. È anche possibile creare un lanciatore, per i più macisti: non c’è limite al peggio!

Allora: ok, funziona ma i programmi Fortran (quelli che conosco) servono a macinare numeri. Hanno bisogno di ingenti quantità di dati di input e producono outputs consistenti. Per questo la norma è di avere tutto su files. Per l’output spesso è utile visualizzarne una parte in forma grafica, cosa che conto di accennare in futuro. Inoltre se l’esecuzione dura più di qualche secondo è bene che il programma rassicuri l’operatore che tutto procede con messaggi su cosa sta facendo, a che punto è e anticipando qualche gossip. Quindi quando il gioco si fa duro cosa c’è di meglio del terminale e della riga di comando?

Tornando a tcl/tk sul web c’è parecchia roba, il linguaggio è vivo ma per quanto mi riguarda non è tra i miei preferiti. L’avevo usato agli inizi degli anni ’90 e mi ricordo parecchie cose ma non è sexy, preferisco usare altro, per esempio Python, pervasivo su Linux e che gode del supporto di entità come Google.

Ancora Fortran

Questa è la continuazione dal post precedente. Poi basta, c’è parecchia altra roba più interessante in giro.

Nel Fortran 90/95 abbiamo l’istruzione select case, proprio come nel Basic, il do, quello che tutti gli altri chiamano for, è molto generalizzato e fa anche quello che altrove fa while; c’è exit e cycle (che sarebbe la continue degli altri linguaggi, continue in Fortran c’è ma fa un’altra cosa, since 1954).
Finalmente sono stati introdotti gli array dinamici e i puntatori per gestirli; questi sono diversi da quelli del C, cioè si usano in modo diverso, anzi nemmeno questo, bisogna provare…, ci sono anche operatori globali per gli array che semplificano il codice e velocizzano l’esecuzione.
Dove i cambiamenti sono maggiori è nella parte non eseguibile (dichiarativa) del linguaggio. Sono stati introdotti i moduli, le procedure possono essere interne o esterne; questo influisce sulla vita e visibilità delle variabili.
Molto ampliata la fase di dichiarazione, con possibilità, tra l’altro, di definire record (struct in C). Come nel C abbiamo le interfacce (i file .h del C). Queste permettono (finalmente) di distinguere i parametri passati a funzioni e procedure come parametri in lettura, scrittura e lettura-scrittura.
Le procedure possono essere ricorsive, a volte serve, credo anche ai fortranisti.
Da sempre in Fortran esistono le funzioni generiche (es. la MAX() è in grado di sapere se gli argomenti sono interi o reali e restituire il valore appropriato in base al tipo di variabile a sinistra del segno =); adesso è possibile costruire le proprie funzioni generiche, ecco servito il polimorfismo. E c’è anche l’operator overloading.
Inoltre ci si è liberati del formato fisso, probabilmente nessuno usa più le schede perforate; le righe possono essere più lunghe, non bisogna più fare attenzione a colonna 73.
Nel 2005 sono uscite altre specifiche (Fortran 2003), e si continua tuttora (Fortran 2008), ci sono tentativi di standardizzare il calcolo parallelo.
Ma in ogni caso restano anche tutti i costrutti antichi: io qualche mese fa ho portato su gfortran un programma di calcolo vecchio più di 20 anni, le modifiche richieste sono state davvero minime. Certo un programmatore oggi lo scriverebbe in modo diverso: sono diventato vecchio!
Scrivendo queste note mi si è installato nella mente un virus quesito cui non so rispondere e del quale non riesco a liberarmi: se Micro$oft et al. non avessero diffuso la leggenda che il Basic è più semplice e avessero prodotto compilatori Fortran come si deve o se Linus Torvalds fosse nato 10-15 anni prima sono sicuro che il mondo informatico di oggi sarebbe diverso. O no? Ma qualcuno diceva così anche prima, riferendosi a Algol, poi diventato Pascal (che a me piaceva parecchio).

Operazione nostalgia

* : ~/bin $ whereis gfortran
gfortran: /usr/bin/gfortran /usr/share/man/man1/gfortran.1.gz
* : ~/bin $ ln -s /usr/bin/gfortran-4.4 ftn

d’ora in poi posso di nuovo invocare il compilatore Fortran, inteso come gfortran con ftn 😉

Argomenti sulla riga di comando

Clive Page’s list of Fortran Resources ha notizie interessanti, tra le quali le librerie grafiche e come leggere gli argomenti passati sulla linea di comando:

• COMMAND_ARGUMENT_COUNT() – This is an integer function returning the number of argument strings on the command-line (excluding the initial command-name if there is one).
• CALL GET_COMMAND_ARGUMENT(number, value, length, status) – This is a subroutine returning the string value of the Nth argument (argument 0 being the command-name). The length (optional) returns the length of the value returned, and status (also optional) returns 0 normally, or non-zero if there is an error such as a non-existing argument).

program arghhh!!! come direbbe il pirata pastafariano che è in me

implicit none
integer :: argnum, cont, arglung, error
character(len=40) :: arg

argnum = command_argument_count()
print *, argnum, 'argomenti sulla riga di comando'

do cont = 1, argnum
    call get_command_argument(cont, arg, arglung, error)
    print *, cont, arglung, error
    print *, arg
end do

end program arghhh

Compilazione ed esecuzione:

* : /media/algol/lab/OKp/fortran $ ftn -o arghhh -ffree-form arghhh.f
* : /media/algol/lab/OKp/fortran $ ./arghhh alpha beta "gamma delta" 'epsilon e tutti gli altri'
           4 argomenti sulla riga di comando
           1           5           0
 alpha
           2           4           0
 beta
           3          11           0
 gamma delta
           4          25           0
 epsilon e tutti gli altri
* : /media/algol/lab/OKp/fortran $

Introducendo i format si può controllare l’output

program arghhh!!! come direbbe il pirata pastafariano che è in me

implicit none
integer :: argnum, cont, arglung, error
character(len=40) :: arg

argnum = command_argument_count()
print 10, argnum, 'argomenti sulla riga di comando'

do cont = 1, argnum
    call get_command_argument(cont, arg, arglung, error)
    print 11, cont, arglung, error
    print 12, arg
end do

!formats
10 format(i4,1x,a)
11 format(3(i3))
12 format(a)

end program arghhh

* : /media/algol/lab/OKp/fortran $ ftn -o arghhh-f -ffree-form arghhh-f.f
* : /media/algol/lab/OKp/fortran $ ./arghhh-f alpha beta "gamma delta" 'epsilon e tutti gli altri' "e una riga troppo lunga che viene troncata"
5 argomenti sulla riga di comando
  1  5  0
alpha
  2  4  0
beta
  3 11  0
gamma delta
  4 25  0
epsilon e tutti gli altri
  5 40 -1
e una riga troppo lunga che viene tronca
* : /media/algol/lab/OKp/fortran $

Operazione nostalgia, nuova ricaduta

Questo è il minicomputer che usavo io, il Prime 550. Il mio aveva un solo hard disk (80 MB) e un’unità nastro sul pannello dove si intravede la scritta “PR1ME 550”. A proposito di dischi: le dimensioni si possono immaginare dai loro contenitori; se si toglievano dall’unità bisognava proteggerli dalla polvere perché non erano sottovuoto e (diceva il manuale) bastava una particella di fumo per mandarlo in crash. Vietato fumare! La cosa più grossa della stanza era comunque il condizionatore perché il ‘puter scaldava peggio di una stufa.
L’immagine è presa da qui. Se scorrete la pagina nell’immagine successiva vedete in primo piano anche la stampante: sigh! bei tempi, allora 😀

Il Fortran, dal IV al 77

Oggi Fortran. Fermi lì, niente panico, è semplice, davvero!

Quando io ho cominciato a lavorare, alla fine degli anni ’70 era usuale avere sul computer (ma allora si chiamava calcolatore o macchina) il compilatore Fortran e basta. A meno che uno fosse particolarmente ricco e allora c’era PL/I o Algol. Naturalmente i gestionali usavano il Cobol, ma è tutta un’altra cosa e noi li snobbavamo.

Con Google si scopre che ci sono millanta manuali di Fortran sul Web, sembra che lo usino ancora i fisici e gli ingegneri. Di sicuro ci sono state delle codifiche recenti, anche se in questo post racconterò delle origini; poi in futuro, chissà, capace che mi aggiorno.
Il Fortran ha origine sull IBM 704, come IBM Mathematical FORmula TRANslation System, John Backus, 1954.

Seguono poi versioni successive fino alla standardizzazione come FORTRAN IV del 1966, che qualcuno chiama Fortran 66. E il Fortran IV, o meglio varianti con estensioni rimase lo standard di fatto ben oltre l’avvento della versione successiva, Fortran 77, 1978.

Nel mondo Linux esiste, ovviamente, un ottimo compilatore fortran, installabile (almeno per Ubuntu) con il comando

sudo apt-get install gfortran

Un esempio pratico, cucinato secondo le diverse scuole di pensiero
Nota: purtroppo WordPress non supporta il fortran al pari di C/C++ e Python 😦
Seguendo la tradizione inizieremo con il classico “Hello world”.

print*, "Hello World!"
end

che compiliamo

* : /media/algol/lab/OKp/fortran $ gfortran -o hw --free-form hw.f
* : /media/algol/lab/OKp/fortran $

e eseguiamo

* : /media/algol/lab/OKp/fortran $ ./hw
 Hello World!
* : /media/algol/lab/OKp/fortran $

Ma mi sembra troppo poco. Propongo di verificare per alcuni casi la congettura di Collatz, dettagliatamente descritta qui. Nella versione italiana della stessa pagina è presente il codice Python che faccio mio con piccolissime modifiche

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import sys

def collatz(n):
	print n,
	if n > 1:
		if n % 2 == 0:
			n /= 2
		else:
			n = 3*n+1
		collatz(n)

# main
num = int(sys.argv[1])
collatz(num)
print

Provandolo si ottiene:

* : /media/algol/lab/OKp/fortran $ python collatz.py 12
12 6 3 10 5 16 8 4 2 1
* : /media/algol/lab/OKp/fortran $

Ma si diceva del Fortran. Quando io iniziai a lavorare avrei prodotto qualcosa come

C23456--ok fino a 72 compresa--
10    PRINT*,'NUMERO'
      READ(*,*,ERR=10) N
20    PRINT*, N
      IF(N.EQ.1) GOTO 999
      I=IDISP(N)
      IF(I.EQ.0) GOTO 30
C     CASO DISPARI
      N=3*N+1
      GOTO 20
30    CONTINUE
C     CASO PARI
      N=N/2
      GOTO 20
999   CONTINUE
C     FINITO
      END

      FUNCTION IDISP(K)
C     RITORNA 1 SE K E' DISPARI, 0 SE PARI
      J=K/2
      L=K-2*J
      IDISP=L
      END

dove:

• la prima riga è un commento (inizia con C in colonna 1);
• la riga 5 controlla se la variabile N è uguale a 1 (.EQ. significa equal), se sì salta all’istruzione con l’etichetta (label) 999;
• la riga 6 setta la variabile I chiamando la funzione IDISP(), definita sotto alle righe 19..24;
• se I vale 0 si salta all’istruzione con etichetta 30 altrimenti si esegue l’istruzione 9 e si salta all’istruzione con etichetta 20, etc, etc…;
• la riga 19 definisce la funzione IDISP() il cui corpo è rappresentato dalle righe 21 e 22;
• la funzione ritorna alla riga 23.

😦 sì, lo so, mancano i numeri di riga 😦

Notare che, siccome trattiamo sempre valori interi, tutte le variabili devono iniziare con una lettera compresa tra I e N (iniziali di integer). Tutte le variabili e le istruzioni devono essere MAIUSCOLE. In 24 righe di programma ci sono 4 GOTO: non male!
Inoltre dobbiamo rispettare la convenzione del formato fisso delle colonne:
i commenti devono essere su una riga a se stante in cui ci sia C in colonna 1;
le etichette (label) devono essere numeriche e inserite nelle colonne 1..5;
le istruzioni devono stare tra le colonne 7..72;
la colonna 6 è riservata al flag di continuazione: se un’istruzione è troppo lunga può essere interrotta e ripresa sulla riga successiva, inserendo un carattere a piacere in colonna 6. 0 (zero) non vale, viene considerato come spazio.
Secondo le convenzioni del tempo non ho inserito spazi inutili.
compilando ed eseguendo si ottiene:

* : /media/algol/lab/OKp/fortran $ gfortran -o CIV CIV.f
* : /media/algol/lab/OKp/fortran $ ./CIV
 NUMERO
12
          12
           6
           3
          10
           5
          16
           8
           4
           2
           1
* : /media/algol/lab/OKp/fortran $

Ma un vecchio avrebbe scritto qualcosa di diverso. Già allora c’erano delle istruzioni sconsigliate ma c’erano le abitudini vecchie anche di 15 anni.
Inoltre la mia versione ha una stranezza: usa una funzione. In Fortran esistono da sempre le funzioni che ritornano un valore, come in C o Python e le subroutines, sottoprogrammi, che non ritornano esplicitamente un valore ma solo come side effect. Qui entra in gioco una caratteristica del Fortran: tutti gli argomenti passati a funzioni e subroutines sono passati per indirizzo; se vengono modificati quando si torna al modulo chiamante i loro valori sono cambiati (come quando in pascal si usa var). Per qualche ragione che non ho mai capito il programmatore fortran IV serio (PF4S) non amava le funzioni; alle volte capitava che il capo mi facesse cambiare il mio codice per depurarlo di quelle schifezze.

La prima versione del PF4S usa il GOTO calcolato, riga 7. Il GOTO calcolato va all’istruzione con la prima label se la variabile di controllo vale 1, all’istruzione con la seconda etichetta se la variabile di controllo vale 2 e così via.

C23456--ok fino a 72 compresa--
10    PRINT*,'NUMERO'
      READ(*,*,ERR=10) N
20    PRINT*, N
      IF(N.EQ.1) GO TO 999
      CALL DISP(N,K)
      GOTO(30,50) K+1
C     CASO DISPARI
50    N=3*N+1
      GOTO 20
30    CONTINUE
C     CASO PARI
      N=N/2
      GOTO 20
999   CONTINUE
C     FINITO
      END

      SUBROUTINE DISP(K,L)
C     RITORNA L=1 SE K E' DISPARI, L=0 SE PARI
      J=K/2
      L=K-2*J
      END

La seconda versione del PF4S usa l’IF aritmetico. Il programma salta all’istruzione con la prima, seconda o terza label in funzione del valore della variabile di controllo, rispettivamente minore, uguale o maggiore di zero.

C23456--ok fino a 72 compresa--
10    PRINT*,'NUMERO'
      READ(*,*,ERR=10) N
20    PRINT*, N
      CALL DISP(N,K)
      IF(K) 999,30,50
C     CASO DISPARI
50    N=3*N+1
      GOTO 20
30    CONTINUE
C     CASO PARI
      N=N/2
      GOTO 20
999   CONTINUE
C     FINITO
      END

      SUBROUTINE DISP(K,L)
C     RITORNA L=1 SE K E' DISPARI, L=0 SE PARI
C     RITORNA L=-1 SE K=1
      L=-1
      IF(K.EQ.1) RETURN
      J=K/2
      L=K-2*J
      END

Poi alla fine del ’78 viene codificata la nuova versione del Fortran, il mitico Fortran 77. C’è da dire che parecchie novità erano false novità, già implementate dai singoli produttori, spesso in forme non compatibili con il nuovo standard. Ma c’erano in ogni caso parecchie cose nuove: l’IF non più limitato a una sola riga, ma nella forma IF THEN / [[ELSE IF] / [ELSE]] / END IF, la gestione delle stringhe, delle variabili logiche e dei files, il ciclo DO (quello che da tutte le altre parti si chiama for) migliorato e altro ancora.
Naturalmente parecchi PF4S si rifiutarono di usarlo. Una volta uno mi ha ricordato che per andare sulla Luna, e tornarci anche parecchie volte, il Fortran IV si era dimostrato più che sufficiente.
In ogni caso io adesso il programma l’avrei scritto così:

      program collatz
*23456--ok fino a 72 compresa--
10    print*,'numero'
      read(*,*,err=10) n
20    print*, n
      if(n.eq.1) go to 999
      if (mod(n, 2).eq.1) then
c       caso dispari
        n=3 * n + 1
      else
c       caso pari
        n = n / 2
      end if
      goto 20
999   continue
c     finito
      end

Questo standard non sarebbe stato modificato (pesantemente) fino al ’90.

Per chi vuole saperne di più

La prima immagine è presa da qui dove c’è un mucchio di cose interessanti. La seconda è pura nostalgia (il primo amore non si scorda mai).

Un manuale online (lo sto leggendo e mi sembra OK) in italiano qui.

Un altro, devo ancora vederlo qui, su questo sito trovate tanta di quella roba che poi vi perdete.

In futuro, forse, continua…