Ancora confronti sui linguaggi

Intanto questa dovrebbe essere l’ultima puntata del tormentone sviluppato negli ultimi post: comincio a detestarlo al pari del caldo afoso che perdura.


Poi devo confessare un mio errore, grave: nella versione Fortran del listato ho postato un file non aggiornato, quello del modulo m_conf.f08. Questa è la versione corretta:

module m_conf
implicit none

real(8) :: f(13)

contains
    real(8) function fsin(a) result(r)
        real(8) :: a
        r = a - a ** 3 / f(3) + a ** 5 / f(5) &
            & - a ** 7 / f(7) + a ** 9 / f(9) &
            & - a ** 11 / f(11) + a ** 13 / f(13)
    end function fsin

    real(8) function fcos(a) result(r)
        real(8) :: a
        r = 1 - a ** 2 / f(2) + a ** 4 / f(4) &
            & - a ** 6 / f(6) + a ** 8 / f(8) &
            & - a ** 10 / f(10) + a ** 12 / f(12)
    end function fcos

end module

Conseguenza della mancata dichiarazione di real(8) per le funzioni fsin() e fcos() queste restituivano un valore real(4). Ma attenzione: siccome tutte le variabili erano in doppia precisione (ahemm, real(8), una volta si diceva così) il tempo di esecuzione non cambia!

Allora l’errore –OK, il grave errore– è quasi solo metodologico. E a mia giustificazione potrei citare il fatto di essere costantemente interrotto da infinite incombenze. E dal caldo, non dimentichiamo il caldo. C’è poi un altro motivo, più grave: a lavorare da solo non hai la possibilità di farlo vedere a qualcuno, raccontarglielo e allora queste cose salterebbero fuori subitissimo.

In ogni caso per verificare se il Fortran è davvero così più performante rispetto a Python ho modificato i codici facendo ripetere i calcoli 10 volte e introducendo il calcolo dei logaritmi decimali di seno e coseno.

Io da piccolo invece del liceo ho frequentato l’istituto per geometri (sono figlio di contadini). Uh! anche Piergiorgio Odifreddi e Flavio Briatore hanno fatto i geometri, pensa te!
A quei tempi non esistevano le calcolatrici elettroniche e i conti si dovevano fare a mano. In particolare per il corso di topografia questi assumevano un’importanza preponderante e si usavano i logaritmi, quelli del Bruhns. Chi fosse Karl Christian Bruhns l’ho scoperto solo l’altro giorno, grazie al solito Wiki. Devo confessarvi che adesso il Bruhns mi risulta simpatico; ma solo da adesso 😉

OK, fine dell’intermezzo nostalgico; il calcolo viene fatto usando la serie

presa da qui.

Poi la trasformazione da logaritmo naturale a decimale richiede solo una divisione ma nello spirito del test si è resa come una funzione. Questo perché anticamente, quando i ‘puter erano grossi e lenti, sarebbe stato un errore grave: chiamare una funzione richiede(va) aggiungere istruzioni, in questo caso inutili. Ma basta divagazioni, ecco le versioni aggiornate.

Python

#!/usr/bin/env python

import math

v = 2.0
f = ([1.0, 1.0, v])
for i in range(3, 14):
    v *= i
    f.append(v)

def fsin(a):
    return a - pow(a, 3) / f[3] + pow(a, 5) / f[5] \
             - pow(a, 7) / f[7] + pow(a, 9) / f[9] \
             - pow(a, 11) / f[11] + pow(a, 13) / f[13]

def fcos(a):
    return 1 - pow(a, 2) / f[2] + pow(a, 4) / f[4] \
             - pow(a, 6) / f[6] + pow(a, 8) / f[8] \
             - pow(a, 10) / f[10] + pow(a, 12) / f[12]

def ln(x):
    if x == 0.0:
        return -1.0e100
    r = (x - 1.0) / (x + 1.0)
    return 2.0 * (r + pow(r, 3) / 3.0 + pow(r, 5) / 5.0
                    + pow(r, 7) / 7.0 + pow(r, 9) / 9.0)

def log(x):
    return x / ln10

pi = math.pi
deg = 60 * 60
nsteps = 180 * deg
step = pi / nsteps
ln10 = ln(10.0)
fmt = "{:3d} {:10.8f} {: 10.8f} {: 10.8f} {: 10.8e}"
fml = "{:12s} {: 12.8e} {: 12.8e}"

ssum = 0.0
csum = 0.0
for i in range (0, nsteps + 1):
    a = i * step
    s = fsin(a)
    ls = log(ln(s))
    ssum += s
    c = fcos(a)
    lc = log(ln(c))
    csum += c
    if (i % (10 * deg)) == 0:
        if c != 0.0:
            t = s / c
        print fmt.format(int(i / deg), a, s, c, t)
        print fml.format('', ls, lc)
print ssum, csum

Fortran

program fc10l
use m_c10l
implicit none

real(8) :: v, pi, step, ssum, csum
real(8) :: a, s, c, t, ls, lc
integer(4) :: deg, nsteps, i, j
character(30) fmt_, fml_

f(1) = 1.0
do i = 2, 13
    f(i) = i * f(i-1)
end do

pi = 4 * atan(1.0)
deg = 60 * 60
nsteps = 180 * deg
step = pi / nsteps
ln10 = myln(dble(10.0))
fmt_ = "(i3, 3(1x, f11.8), 1x, e12.4)"
fml_ = "(12x, 2(e12.4))"

do j = 1, 10
    ssum = 0.0
    csum = 0.0

    do i = 0, nsteps
        a = i * step
        s = fsin(a)
        ls = mylog(myln(s))
        ssum = ssum + s
        c = fcos(a)
        lc = mylog(myln(c))
        csum = csum + c
        if (mod(i,10 * deg) == 0) then
            if (c /= 0.0) t = s / c
            print fmt_, int(i / deg), a, s, c, t
            print fml_, ls, lc
        end if
    end do

    print *, ssum, csum
end do
end program

 

module m_c10l
implicit none

real(8) :: f(13)
real(8) ln10

contains
    real(8) function fsin(a) result(r)
        real(8) :: a
        r = a - a ** 3 / f(3) + a ** 5 / f(5) &
            & - a ** 7 / f(7) + a ** 9 / f(9) &
            & - a ** 11 / f(11) + a ** 13 / f(13)
    end function fsin

    real(8) function fcos(a) result(r)
        real(8) :: a
        r = 1 - a ** 2 / f(2) + a ** 4 / f(4) &
            & - a ** 6 / f(6) + a ** 8 / f(8) &
            & - a ** 10 / f(10) + a ** 12 / f(12)
    end function fcos

    real(8) function myln(x) result(res)
        real(8) :: x, r
        if(x == 0.0) then
            res = -1.0e20
        else
            r = (x - 1.0) / (x + 1.0)
            res = 2.0 *(r ** 3 / 3.0 + r ** 5 / 5.0 &
                & + r ** 7 / 7.0 + r ** 9 / 9.0)
        end if
    end function

    real(8) function mylog(x) result(res)
        real(8) :: x
        res = x / ln10
    end function

end module

Java

import static java.lang.Math.pow;
import static java.lang.System.out;

class jc10l {
    static double f[];
    //static double ln10;
    static double ln10 = myln(10.0);

    public static double fsin(double a) {
        return a - pow(a, 3.0) / f[3] + pow(a, 5.0) / f[5]
                 - pow(a, 7.0) / f[7] + pow(a, 9.0) / f[9]
                 - pow(a, 11.0) / f[11] + pow(a, 13.0) / f[13];
    }

    public static double fcos(double a) {
        return 1.0 - pow(a, 2.0) / f[2] + pow(a, 4.0) / f[4]
                   - pow(a, 6.0) / f[6] + pow(a, 8.0) / f[8]
                   - pow(a, 10.0) / f[10] + pow(a, 12.0) / f[12];
    }

    public static double myln(double x) {
        if(x == 0.0) {
            return -1.0e20;
        } else {
            double r = (x - 1.0) / (x + 1.0);
            return 2.0 * (r + pow(r, 3.0) / 3.0
                            + pow(r, 5.0) / 5.0
                            + pow(r, 7.0) / 7.0
                            + pow(r, 9.0) / 9.0);
        }
    }

    public static double mylog(double x) {
        return x / ln10;
    }

    public static void main(String args[]) {
        int i, j;
        f = new double[14];
        f[0] = 1.0;
        for (i = 1; i < 14; i++) {
            f[i] = i * f[i - 1];
        }
        int deg = 60 * 60;
        int nsteps = 180 * deg;
        double step = Math.PI / nsteps;
        double ssum;
        double csum;
        double a, s, c, t = 0.0;
        double ls, lc;

        for (j = 1; j < 11; j++) {
            ssum = 0.0;
            csum = 0.0;
            for (i = 0; i                 a = i * step;
                s = fsin(a);
                ls = mylog(myln(s));
                ssum += s;
                c = fcos(a);
                lc = mylog(myln(c));
                csum += c;
                if ((i % (10 * deg)) == 0) {
                    if (c != 0.0) {
                        t = s / c;
                    }
                    out.printf("%3d %11.8f %11.8f %11.8f %15.8e\n",
                          (i / deg), a, s, c, t);
                    out.printf("              %15.8e %15.8e\n",
                               ls, lc);
                }
            }
            out.println(ssum);
            out.println(csum);
        }
    }

}

newLISP

#!/usr/bin/newlisp

(set-locale "C") ;;; sono abituato così, perso un'ora!

;         0 1 2  3  4 5 6  7  8 9 10  11  12 13
(set 'p '(1 1 -1 -1 1 1 -1 -1 1 1 -1  -1  1  1))
(set 'f '(1))
(for (i 1 13)
    (set 't (* (p i) i (abs (f (- i 1)))))
    (push t f i)
)

(define (fsin a)
    (add a (div (pow a 3) (f 3)) (div (pow a 5) (f 5))
           (div (pow a 7) (f 7)) (div (pow a 9) (f 9))
           (div (pow a 11) (f 11)) (div (pow a 13) (f 13))
    )
)

(define (fcos a)
    (add 1 (div (pow a 2) (f 2)) (div (pow a 4) (f 4))
           (div (pow a 6) (f 6)) (div (pow a 8) (f 8))
           (div (pow a 10) (f 10)) (div (pow a 12) (f 12))
    )
)

(define (myln x)
    (if (= x 0.0)
        (set 'r -1.0e100)
        (begin
            (set 'r (div (dec x 1.0) (add x 1.0)))
            (mul 2.0 (add r (div (pow r 3) 3.0)
                            (div (pow r 5) 5.0)
                            (div (pow r 7) 7.0)
                            (div (pow r 9) 9.0)
                      )
            )
        )
    )
)

(define (mylog x)
    (div x ln10)
)

(set 'pi (mul 4 (atan 1.0))
     'deg (* 60 60)
     'nsteps (* 180 deg)
     'step (div pi nsteps)
     'ln10 (myln 10)
)

(for (j 1 10)
    (set 'ssum 0
         'csum 0
    )

    (for (i 0 nsteps)
        (set 'a (mul i step)
             's (fsin a)
             'ls (mylog (myln s))
             'c (fcos a)
             'lc (mylog (myln c))
        )
        (inc ssum s)
        (inc csum c)
        (if (!= c 0)
            (set 't (div s c)))
        (if (= (mod i (* deg 10)) 0)
            (begin
            (println (format "%3d %11.8f %11.8f %11.8f %15.8e"
                     (/ i deg) a s c t))
            (println (format "                %15.8e %15.8e"
                     ls lc))
            )
        )
    )
    (println ssum " " csum)
)

(exit)

Con newLISP ho perso un’ora a inseguire un errore mysteryuoso assay! Continuava a dirmi che le parentesi non erano bilanciate anche se l’editor mi diceva che erano OK. Dimenticare qualche parentesi con il Lisp è la cosa peggiore che un niubbo possa fare, ma io tanto niubbo non sono e poi ormai una certa pratica ce l’ho. Altro errore è quando apri e chiudi dove non devi ma neanche questo mi sembrava il caso. Poi l’illuminazione, la virgola al posto del punto. Atz! dimenticato (set-locale "C"), sì io i numeri sono abituato a scriverli alla ‘mericana.
A dire il vero ho dimenticato Nmila punti-e-virgola nella versione Java ma lì il compilatore te lo dice.

Nell’esecuzione dei confronti è tornato mooolto utile il fatto che time scrive su stderr; quindi ridirigendo l’output (stdout) sul teminale si può avere un risultato leggibile

Il risultato finale conferma sostanzialmente quello della prima versione. I rapporti rispetto al Fortran dei tempi di esecuzione diminuiscono leggermente ma rimangono elevati.

Per motivi non solo informatici pare che si userà il Fortran per i calcoli e Python per cose complementari, come la costruzione dei grafici riassuntivi delle elaborazioni. Mi sa che ne riparlerò, forse… 😉

Posta un commento o usa questo indirizzo per il trackback.

Commenti

  • dikiyvolk  Il 10 luglio 2012 alle 15:36

    Ciao juhan,

    Anche tu a fare i test 🙂 comparativi. Io li trovo molto utili per evidenziare le differenze tra i vari linguaggi e il perché a volte si passa per retrogradi su certe scelte… sempre che chi ti ascolta sia disposto ad ascoltare 🙂
    Aspetto anche l’articolo sui grafici di python.

    Solo per sottolineare quanto tu sia fortunato io non solo muoio dal caldo a casa ma anche in ufficio (si è rotto il condizionamento).
    Appena riesco provo anche io a fare i tuoi test 🙂 come mai hai escluso il C?

    Qui servirebbe una bella birra gelata da buon programmatore 🙂

    Buona giornata

    • juhan  Il 10 luglio 2012 alle 15:45

      Huela! Intanto caldo anche qui, troppo.
      Il test doveva essere solo Fortran e Python, poi è saltato fuori Java (e C# che ho detto di no, non c’è per Linux (spero sia vero)). Allora ci ho aggiunto newLISP che si comporta bene. Il C ha la fama di essere difficilissimo e poi lì hanno un sacco di roba vecchia in Fortran (gli ing. del Poli lo usano tuttora).
      Adesso si dovrebbe partire con il Fortran (forse) e grafici con Python (forse). Tutto questo caldo permettendo. Sui grafici conto di postare qualcosa, prossimamente; prima finisco ancora con il Fortran.

  • dikiyvolk  Il 18 luglio 2012 alle 11:58

    Ciao juhan,

    Finalmente ho riscritto il tuo codice in C e i tempi sulla mia macchina virtuale sono pessimi 😦

    real 0m21.049s
    user 0m21.027s
    sys 0m0.005s

    #include <stdio.h>
    #include <math.h>

    static double f[14];

    static double fsin(double a) {
    return a – pow(a, 3.0) / f[3] + pow(a, 5.0) / f[5]
    – pow(a, 7.0) / f[7] + pow(a, 9.0) / f[9]
    – pow(a, 11.0) / f[11] + pow(a, 13.0) / f[13];
    }

    static double fcos(double a) {
    return 1.0 – pow(a, 2.0) / f[2] + pow(a, 4.0) / f[4]
    – pow(a, 6.0) / f[6] + pow(a, 8.0) / f[8]
    – pow(a, 10.0) / f[10] + pow(a, 12.0) / f[12];
    }

    static double myln(double x) {
    if (x == 0.0) {
    return -1.0e20;
    } else {
    double r = (x – 1.0) / (x + 1.0);
    return 2.0 * (r + pow(r, 3.0) / 3.0
    + pow(r, 5.0) / 5.0
    + pow(r, 7.0) / 7.0
    + pow(r, 9.0) / 9.0);
    }
    }

    static double ln10 = myln(10.0);

    static double mylog(double x) {
    return x / ln10;
    }

    int main(int argc, char **argv)
    {
    int i, j;
    f[0] = 1.0;
    for (i = 1; i < 14; i++)
    f[i] = i * f[i – 1];

    int deg = 60 * 60;
    int nsteps = 180 * deg;
    double step = M_PI / nsteps;
    double ssum;
    double csum;
    double a, s, c, t = 0.0;
    double ls, lc;

    for (j = 1; j < 11; j++)
    {
    ssum = 0.0;
    csum = 0.0;
    for (i = 0; i < nsteps; i++)
    {
    a = i * step;
    s = fsin(a);
    ls = mylog(myln(s));
    ssum += s;
    c = fcos(a);
    lc = mylog(myln(c));
    csum += c;

    if ((i % (10 * deg)) == 0)
    {
    if (c != 0.0)
    t = s / c;

    printf("%3d %11.8f %11.8f %11.8f %15.8e\n", (i / deg), a, s, c, t);
    printf(" %15.8e %15.8e\n", ls, lc);
    }
    }
    }
    printf("%d\n", ssum);
    printf("%d\n", csum);
    }

    Per compilarlo ho dato il comando:

    g++ main.cpp -o cc10l -O3

    Riesci a compilarlo tu e vedere i tempi? Magari ho sbagliato qualche cosa nel codice, l’ho tradotto al volo dal sorgente java prima della pausa caffe 🙂

    • glipari  Il 18 luglio 2012 alle 14:34

      Ho provato io sul mio Intel(R) Core(TM) i7-2760QM CPU @ 2.40GHzPentium, con Kubuntu, gcc 4.6.3, e compilato come g++ prova.c -o prova. Ecco i risultati

      real 0m9.976s
      user 0m9.953s
      sys 0m0.008s

      Come posso compilare il programma fortran?

      • juhan  Il 18 luglio 2012 alle 14:39

        Installare gfortran e poi gfortran -o fc10l m_c10l.f08 fc10l.f08
        (i nomi dei due file sorgenti).
        Appena posso provo anch’io.

    • glipari  Il 18 luglio 2012 alle 16:44

      dikiyvolk,

      cambia tutti le chiamate a pow() con la funzione mypow() qui sotto, e otterrai uno speed-up di circa 10 volte:

      inline double mypow(double x, int y)
      {
      double r = 1;
      for (int i=0; i<y; i++) r *= x;
      return r;
      }

      I tempi ottenuti sono adesso:

      real 0m0.931s
      user 0m0.928s
      sys 0m0.004s

      E con il Fortran:

      real 0m1.136s
      user 0m1.116s
      sys 0m0.016s

      Probabilmente si può ottimizzare ancora un pochino, sia il C che il Fortran…

      • dikiyvolk  Il 19 luglio 2012 alle 09:11

        Ciao,

        Ho fatto le modifiche che mi hai suggerito e i tempi sempre sulla macchina virtuale sono stati:

        real 0m1.285s
        user 0m1.282s
        sys 0m0.002s

        Quindi molto meglio. Però la mia domanda è:

        – Perché la pow standard è così lenta? Non ho trovato il codice, ma googlando in giro si parla di controlli che deve fare per esponenziali frazionari.

        Tu ne sai qualche cosa?

        Buona giornata

      • glipari  Il 19 luglio 2012 alle 12:02

        La pow() deve trattare il caso generico di esponenziale con la virgola. Ad esempio 2^0.5 equivale alla radice di 2. In generale, la funzione x^y non è una semplice moltiplicazione di x per y volte…. Juhan ne dovrebbe sapere qualcosa!!! 🙂 Però, se l’esponente è intero, possiamo semplicifare alla grande, ed ecco l’accellerazione!

      • dikiyvolk  Il 19 luglio 2012 alle 12:56

        Si avevo capito che il problema fosse quello… ma la domanda che mi viene in mente è come mai quella (la funzione pow) degli altri linguaggi tipo java è più veloce?

      • juhan  Il 19 luglio 2012 alle 12:59

        poi ci scrivo qualcosa, mi manca il tempo e ho troppo caldo 😦

  • dikiyvolk  Il 24 luglio 2012 alle 15:03

    Ciao Juhan,

    Ho scritto anche il codice in Ruby. Come linguaggio l’ho trovato un po’ ostico da capire, ma è sicuramente una limitazione del mio approccio alla programmazione un po’ troppo retrò.

    Comunque ho fatto ben 3 versioni del programma.
    Quella che più si avvicina a quello python è la seguente:

    include Math
    
      def fsin(a)
          return a - a**3 / $f[3] + a**5 / $f[5] \
                 - a**7 / $f[7] + a**9 / $f[9] \
                 - a**11 / $f[11] + a**13 / $f[13]
      end
    
      def fcos(a)
          return 1 - a**2 / $f[2] + a**4 / $f[4] \
                 - a**6 / $f[6] + a**8 / $f[8] \
                 - a**10 / $f[10] + a**12 / $f[12]
      end
    
      def myln(x)
        if x == 0.0
          return -1.0e20
        else
          r = (x - 1.0) / (x + 1.0);
          return 2.0 * (r + (r**3.0) / 3.0 \
                 + (r**5.0) / 5.0 \
                 + (r**7.0) / 7.0 \
                 + (r**9.0) / 9.0)
        end
    
      end
    
      def mylog(x)
          return x / Math.log(10)
      end
    
      $f = Array.new(14)
      $f[0] = 1.0
      for i in 1..$f.length
        $f[i] = i * $f[i - 1]
      end
    
      deg = 60 * 60
      nsteps = 180 * deg
      step = Math::PI / nsteps
      ssum = 0.0
      csum = 0.0
      t = 0.0
    
      for j in 1..(11 - 1)
        ssum = 0.0;
        csum = 0.0;
        for i in 0..(nsteps -1)
          a = i * step;
          s = fsin(a);
          ls = mylog(myln(s));
          ssum += s;
          c = fcos(a);
          lc = mylog(myln(c));
          csum += c;
          if (i % (10 * deg)) == 0
            if c != 0.0
              t = s / c;
            end
            printf("%3d %11.8f %11.8f %11.8f %15.8e\n", (i / deg), a, s, c, t)
            printf("              %15.8e %15.8e\n", ls, lc)
          end
        end
      end
      print(ssum)
      print(csum)
    

    Le performance non sono state buone… anzi sulla mia macchina virtuale 😦 i tempi sono stati:

    # time ruby r_conf_2.rb > rl2

    real 1m25.174s
    user 1m25.043s
    sys 0m0.033s

    Con java e i tuoi sorgenti invece i miei tempi sono stati:

    # time java -classpath . jc10l > jo

    real 0m37.901s
    user 0m37.984s
    sys 0m0.117s

    Con python 2.7.2:

    # time /tmp/pc10l.py > op

    real 0m8.464s
    user 0m8.377s
    sys 0m0.070s

    Su ruby ho anche provato a cambiare l’operatore ** con una funzione come avevamo fatto con il C, ma i tempi sono quadruplicati. Magari sono io che ho sbagliato a scrivere il codice in Ruby, ci vorrebbe un esperto.
    Tu con che versione di python hai fatto i test?

    • juhan  Il 24 luglio 2012 alle 15:10

      Ho postato stamattina un aggiornamento con le cose tue e del prof.
      Ruby è usato anche standalone o solo per il Web? E avrebbe vantaggi rispetto a Python che mi sembra più popolare?
      Inoltre qui si parla di Lua; qualcuno ne sa qualcosa?
      Poi basta vorrei/dovrei dedicarmi alla grafica 😉
      O fuggo dal blog e ne faccio uno tutto sul Lisp: (solo-lambda) 😀

      • dikiyvolk  Il 24 luglio 2012 alle 15:26

        Ciao… non vale mi aggiorni le pagine a tradimento 🙂 è che OKPanico è la mia pagina di default…e aprendola al mattino non l’aggiorno se non chiudo il browser.
        Ruby si può usare standalone ma è molto usato per applicativi WEB. All’estero vedo molte richieste di sviluppatori, molto ben pagati (più del C/C++) ma in Italia nulla. Ha avuto un discreto successo tanto che hanno creato JRuby che è un Ruby in Virtual Machine.
        Lua??? Mi manca… ma sono troppi!!! Altro che la carica dei 101 🙂 Scappo a programmare in COBOL su un’isola deserta e niente messaggio in bottiglia per farmi venire a salvare 🙂

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo di WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione /  Modifica )

Google photo

Stai commentando usando il tuo account Google. Chiudi sessione /  Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione /  Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione /  Modifica )

Connessione a %s...

Questo sito utilizza Akismet per ridurre lo spam. Scopri come vengono elaborati i dati derivati dai commenti.

%d blogger hanno fatto clic su Mi Piace per questo: