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…

About these ads