I conti come li fa il computer

Zar, conosciuto anche come Roberto Zanasi è l’autore di un libricino piccolo come dimensioni ma, come dice qui in Padagna

cit ma cativ

(cit = piccolo).
Lo consiglio caldamente a tutti, trovate gli estremi sul suo blog.
Già perché ha anche un blog. Che dovete seguire. OK, genera dipendenza e a me non mi lascia dormire la notte mi ha costretto a una deviazione dal programma dei post che avevo in programma. Cioè questo post è conseguenza di questo. Poi ha continuato sul tema: essendo un matematico era prevedibile che a quell’1 sarebbe seguito un 2, eccolo. Comunque di questo non parlerò: mi basta il primo.

Allora la prima cosa che mi è venuta spontanea leggendo lo Zar è stata quella di aprire un terminale e invocare il fido Pythonin modo interattivo, con questi risultati (lo so mi ripeto, sono quelli pubblicati nel commento del suddetto post ma poi cambio)

* : _da_proc_1_ $ python
Python 2.6.6 (r266:84292, Sep 15 2010, 15:52:39) 
[GCC 4.4.5] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> un_terzo = 1./3
>>> un_terzo
0.33333333333333331
>>> uno = un_terzo * 3
>>> uno
1.0
>>> import math
>>> rad_2 = math.sqrt(2)
>>> rad_2
1.4142135623730951
>>> due = rad_2 * rad_2
>>> due
2.0000000000000004
>>> dueq = rad_2 ** 2
>>> dueq
2.0000000000000004
>>> print due == 2.
False
>>> due_q = rad_2 ** 2
>>> due_q
2.0000000000000004
>>> print due == due_q
True
>>> delta = abs(2 - due)
>>> delta
4.4408920985006262e-16
>>> quit()

C’è un arrotondamento! Per un’operazione semplice come calcolare la radice di 2 e poi farne il quadrato otteniamo un errore dell’ordine di del milionesimo di miliardo (10-15). Se poi si usassero questi numeri per fare altri calcoli in cascata l’errore, presumibilmente, aumenta.

Python è quello che viene in mente subito, come una volta si sarebbe usato il Basic (c’è chi lo usa ancora, non so se ci sia su Linux) ma io ho un altro interprete phigo, che lovvo, assaje: newLISP.
newLISP mi piace perché è piccolissimo, potente, multipiattaforma (su Windoze non c’è nemmeno bisogno di installarlo). Avete presente Visual Studio? Tutto l’opposto.
Allora vediamo cosa succede con newLISP, sempre in modalità interattiva


* : _da_proc_1_ $ newlisp
newLISP v.10.2.8 on Linux IPv4 UTF-8, execute 'newlisp -h' for more info.

> (set 'un_terzo (div 1,0 3,0))
0,3333333333
> (set 'uno (mul un_terzo 3))
1
> (set 'rad_2 (sqrt 2,0))
1,414213562
> (set 'due (mul rad_2 rad_2))
2
> (set 'dueq (pow rad_2 2))
2
> (set 'uguale (= due 2,0))
nil
> (set 'delta (abs (- due 2,0)))
0
> (set 'mag (> (- due 2,0) 0))
nil
> (set 'mn (> 0 (- due 2,0)))
nil
> * : _da_proc_1_ $ 

Notare come newLISP rilevi le impostazioni locali relative ai numeri per cui scrivo 1,0 e 2,0 e non 1.0 e 2.0 come i programmatori sono soliti fare. Questa opzione è peraltro modificabile con una semplice istruzione (che non vi dico, andate sul sito di Lutz).

Ma nel mondo di Linux esistono altri programmi per questo genere di problemi. Io sono affezionato a calc (dovete installarlo), ecco cosa succede


* : arrotondamento $ calc
C-style arbitrary precision calculator (version 2.12.3.3)
Calc is open software. For license details type:  help copyright
[Type "exit" to exit, or "help" for help.]

; rad_2 = sqrt(2.0)
; print rad_2
1.4142135623730950488
; due = rad_2 * rad_2
; print due
~2.00000000000000000000
; print due - 2
~-0.00000000000000000000
; uguale = (due == 2)
; uguale
	0
; mag = (due > 2)
; print mag
0
; mn = (due < 2)
; print mn
1
; delta = 2 - due
; print delta
~0.00000000000000000000
; print delta / 1.e-15
~0.00000477643336092562
; print delta / 1.e-20
0.477643336092561856
; quit

Notare come il delta qui sia dell’ordine di 10-20: bravo calc!

Ma fin dalle origini di Unix esiste bc, sento dire da geek, nerd e affini. OK eccolo


* : arrotondamento $ bc
bc 1.06.95
Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
scale = 20
rad_2 = sqrt(2)
print rad_2
1.41421356237309504880
due = rad_2 * rad_2
print due
1.99999999999999999999
scale = 100
rad_2c = sqrt(2)
print rad_2c
1.4\
14213562373095048801688724209698078569671875376948073176679737990732\
4784621070388503875343276415727
due_c = rad_2c * rad_2c
print due_c
1.99999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999
delta = 2 - due_c
print delta
.00\
00000000000000000000000000000000000000000000000000000000000000000000\
000000000000000000000000000001
scale = 1000
rad_2m = sqrt(2)
print rad_2m
1.414213562373095048801688724209698078569671875376948073176679737990\
73247846210703885038753432764157273501384623091229702492483605585073\
72126441214970999358314132226659275055927557999505011527820605714701\
09559971605970274534596862014728517418640889198609552329230484308714\
32145083976260362799525140798968725339654633180882964062061525835239\
50547457502877599617298355752203375318570113543746034084988471603868\
99970699004815030544027790316454247823068492936918621580578463111596\
66871301301561856898723723528850926486124949771542183342042856860601\
46824720771435854874155657069677653720226485447015858801620758474922\
65722600208558446652145839889394437092659180031138824646815708263010\
05948587040031864803421948972782906410450726368813137398552561173220\
40245091227700226941127573627280495738108967504018369868368450725799\
36472906076299694138047565482372899718032680247442062926912485905218\
10044598421505911202494413417285314781058036033710773091828693147101\
71111683916581726889419758716582152128229518488472
due_m = rad_2m * rad_2m
print due_m
1.999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999
print 2 - due_m
.00000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000001
print due_c - due_m
-.00000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999999999999999999999999999999999999999\
9999999999999
quit
* : arrotondamento $ 

Notare come cambiano le cose al cambiare della variabile scale. Io mi sono fermato a 1000 ma se qualcuno vuole…

Sia Python che newLISP vengono di solito usati creando uno script e eseguendolo. Certo, ecco, Python (arr.py)

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

import math

un_terzo = 1./3
uno = un_terzo * 3
print "un_terzo =", un_terzo
print "uno = ", uno
if uno == 1.:
	print "gli uni sono uguali"
else:
	print "gli uni sono diversi"

rad_2 = math.sqrt(2.0)
due = rad_2 * rad_2
print "rad_2 =", rad_2
print "due =", due
if due == 2.:
	print "i due sono uguali"
else:
	print "i due sono diversi"

delta = abs(2. - due)
print "delta =", delta

che produce


* : arrotondamento $ python arr.py
un_terzo = 0.333333333333
uno =  1.0
gli uni sono uguali
rad_2 = 1.41421356237
due = 2.0
i due sono diversi
delta = 4.4408920985e-16
* : arrotondamento $

e con newLISP (arr.nl)

;; verifica arrotondamento con newLISP
; notare che nl usa la virgola perché legge le impostazioni locali
; naturalmente si può cambiare ma lo lascio così.

(set 'rad_2 (sqrt 2,0))
(print "rad_2 = " rad_2 "\n")
(set 'due (mul rad_2 rad_2))
(print "due = " due "\n")

(if (= 2,0 due)
	(print "uguali\n")
	(print "diversi\n"))
(set 'delta (- 2,0 due))
(print "delta = " delta "\n")

(if (> delta 0)
	(print "delta è positivo\n"))
(if (< delta 0)
	(print "delta è negativo\n"))
(if (= delta 0)
	(print "delta è nullo\n"))

(exit)

si ottiene


* : arrotondamento $ newlisp arr.nl 
rad_2 = 1,414213562
due = 2
diversi
delta = 0
delta è nullo
* : arrotondamento $ 

Qui ci sarebbe da perplimersi sui test di delta: non è zero (altrimenti due sarebbe uguale a 2), non è positivo nè negativo ma nullo 😦

Sta di fatto che i conti giusti, per i matematti, finora non si sono visti. Ma non abbiamo provato con linguaggi compilati; chissà magari cambia qualcosa 8)

Allora una versione scolastica in C (arr.c)

#include
#include

void main(void) {
  	double un_terzo, uno, rad_2, due, delta;

	un_terzo = 1./3.;
	uno = un_terzo * 3.;
	printf("un_terzo = %f\nuno = %f\n", un_terzo, uno);
	if (uno == 1.) {
		printf("gli uni sono uguali\n");
	}else{
		printf("gli uni sono diversi\n");
	}
	rad_2 = sqrt(2.0);
	due = rad_2 * rad_2;
	printf("rad_2 = %f\ndue = %f\n", rad_2, due);
	if (due == 2.) {
		printf("i due sono uguali\n");
	}else{
		printf("i due sono diversi\n");
	}
	delta = abs(due - 2.0);
	printf("delta = %e\n", delta);
	printf("delta / 1e-15 = %e\n", delta / 1.e-15);
}

che produce


* : arrotondamento $ cc -lc arr.c
* : arrotondamento $ ./a.out
un_terzo = 0.333333
uno = 1.000000
gli uni sono uguali
rad_2 = 1.414214
due = 2.000000
i due sono diversi
delta = 0.000000e+00
delta / 1e-15 = 0.000000e+00

Mistero, anche qui, sul valore di delta. Io che sono vecchio ho passato buona parte della mia vita con il Fortran. Vediamo cosa capita con il Fortran 77. C’è poco da ridere! questa è la versione normale del Fortran, vero che poi ne sono uscite altre 3 o 4 ma per altri motivi che qui non c’interessano (e che non saprei usare senza consultare il manuale ;-)).
Ecco il codice (arr.f)

	  program arr
C234567
	  real*8 un_terzo, uno, rad_2, due, delta

	  un_terzo = 1./3.
	  uno = un_terzo * 3
	  print *, 'un_terzo =', un_terzo
	  print *, 'uno =', uno
	  if(uno.eq.1.) then
	  	print *, 'gli uni sono uguali'
	  else
	  	print *, 'gli uni sono diversi'
	  end if

	  rad_2 = sqrt(2.0)
	  due = rad_2 * rad_2
	  print *, 'rad_2 = ', rad_2
	  print *, 'due = ', due
	  if (due.eq.2.0) then
	  	print *, 'i due sono uguali'
	  else
	  	print *, 'i due sono diversi'
	  end if

	  delta = abs(due - 2.0)
	  print *, 'delta = ', delta

	  end

che produce


* : arrotondamento $ gfortran -o arr arr.f
* : arrotondamento $ ./arr
 un_terzo =  0.33333334326744080     
 uno =   1.0000000298023224     
 gli uni sono diversi
 rad_2 =    1.4142135381698608     
 due =    1.9999999315429164     
 i due sono diversi
 delta =   6.84570835574049852E-008
* : arrotondamento $ 

Ecco adesso bisognerebbe trarre le conclusioni, aprire il dibattito (No il dibattito no!, Nanni Moretti ai tempi in cui il Fortran era giovane).
Facciamo così: il post è già molto lungo, io altro materiale (non mio) ce l’ho; se l’argomento akkiappa ci sarà una continuazione, forse, anzi molto probabilmente 8)

Posta un commento o usa questo indirizzo per il trackback.

Commenti

  • zar  Il 1 giugno 2011 alle 16:19

    Ma che significa che zero non è zero ma è nullo?

    • juhan  Il 1 giugno 2011 alle 16:22

      Non lo so! Ma non dirlo a nessuno che devo fare quello saputo.
      Sei stato velocissimo 8)

  • glipari  Il 1 giugno 2011 alle 17:24

    vi lascio dei puntatori:

    http://en.wikipedia.org/wiki/Floating_point

    e

    http://en.wikipedia.org/wiki/IEEE_754-2008

    Tutti i maggiori processori usano questa rappresentazione, e i calcoli vengono fatti in hardware (una volta si diceva c’era il coprocessore, ora è tutto integrato). I vari linguaggi si comportano in modo leggermente diverso perché usano un “layer” software diverso. Alcuni linguaggi usano rappresentazioni alternative e routine software per rendere arbitrariamente precisa la rappresentazione.

    Ora non ho tempo, ma magari stasera (se non crollo dal sonno) metto qualche commento più puntuale.

  • glipari  Il 1 giugno 2011 alle 17:26

    Dimenticavo: una cosa è la rappresentazione interna, un’altra cosa quello che vedi a video. Per il programma C, devi mettere il format giusto nella printf, e ti verrà fuori lo stesso output del python.

    • juhan  Il 1 giugno 2011 alle 17:48

      Grazie 😀
      Sono andato un po’ di corsa, dopo aver letto il post di Zar non riuscivo a non pensarci e l’ho scritto quasi di getto e non ponderato come faccio di solito. E sono successe cose strane, per dirne una mi sono dimenticato che in Fortran si usa “end if”.
      Un’altra cosa: con gli arrotondamenti i programmatori fortran ci lottavano tutto il tempo; spero si sia capito.

      • glipari  Il 1 giugno 2011 alle 18:27

        è uno dei motivi per cui il fortran è ancora in giro: qualcuno ha scritto delle buone librerie matematiche in fortran in cui tutti i problemi di arrotondamento sono stati accuratamente tenuti sotto controllo. per cui i fisici preferiscono fidarsi e usare queste librerie invece di stare a riscriversele e a risbatterci la testa in un altro linguaggio.

        e chiamali fessi!

Trackback

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 )

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: