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… 😉
