¡Descarga Punto Flotante y más Apuntes en PDF de Matemáticas solo en Docsity!
Cap´ıtulo 6
Aritm´etica de punto flotante.
Estudiaremos (muy superficialmente) la representaci´on de los n´umeros reales en los ordenadores y su ‘aritm´etica’. Las referencias pertinentes a la bibliograf´ıa ir´an apareciendo a lo largo de la exposici´on. En todo este cap´ıtulo, cuando no se especifica lo contrario, 〈〈^ n´umero〉〉^ significa n´umero real.
Aperitivo.
Lee atentamente estos enunciados:
(1) Seg´un tu calculadora, ¿qu´e vale
1 3
(2) Comprueba con tu calculadora que
a) a =
b) b =
a − 1 2 − a
de modo que a > b.
Responde ahora cuidadosamente a estas cuestiones:
(i) ¿C´omo se explica el resultado de (1)?
(ii) ¿Son correctas las operaciones de (2)? S´I ut NO ut NO SE´ ut
(iii) En (2), lo cierto es que: a > b ut a = b ut a < b ut.
¿Qu´e consecuencias sacas?
166 CAP´ITULO 6. ARITM ETICA DE PUNTO FLOTANTE.´
6.1. Representaci´on en punto flotante
En una reuni´on se plantea la siguiente pregunta: “¿Cu´anto vale π?” El ingeniero dice: “Es aproximadamente 3 y 1/7.” El f´ısico responde: “Es 3.141592.” El matem´atico piensa un momento, y replica: “Es igual a π.” (Del folklore)
Hemos mostrado la posibilidad te´orica de representar los n´umeros reales mediante sucesiones de d´ıgitos en distintas bases. Ahora vamos a enfrentarnos con su representaci´on pr´actica en un orde- nador, con las limitaciones efectivas que la realidad impone.
6.1.1. N´umeros de m´aquina
El conjunto de los n´umeros reales no puede representarse plenamente en un ordenador (¡ni siquiera el de los n´umeros enteros!). Si queremos guardar las cifras decimales (o m´as frecuentemente, binarias) de un n´umero en un ordenador, s´olo disponemos de un n´umero finito de ‘huecos’ para hacerlo. Podemos, simplemente, ignorar las cifras del desarrollo a partir de un cierto lugar; pero se aprovecha mejor el espacio evitando “ceros innecesarios” como en la ‘notaci´on cient´ıfica’, con una representaci´on en punto flotante:
n´umero = signo × mantisa × base exponente
Por ejemplo, para la constante de gravitaci´on universal G = 6. 673 × 10 −^11 , la mantisa es 6.673, la base es 10 y el exponente −11. Este y otros ejemplos de magnitudes muy grandes o muy peque˜nas, como
Masa del electr´on 9. 109 × 10 −^28 g. Masa del sol 1. 989 × 1030 kg. Constante de Plank 6. 626 × 10 −^34 J·s N´umero de Avogadro 6. 022 × 1023 mol−^1
ponen de manifiesto el ahorro de cifras que conlleva esta notaci´on. Pero cada n´umero admite infinitas representaciones de este tipo: − 111 /2 = − 55. 5 × 100 = − 5. 55 × 101 = − 0. 555 × 102 =.. ., o tambi´en − 111 /2 = − 555 × 10 −^1 = − 5550 × 10 −^2 =... Para tener unicidad hace falta normalizar la representaci´on. Habitualmente, como hemos hecho en los ejemplos de notaci´on cient´ıfica, el primer d´ıgito no nulo de la mantisa aparece a la izquierda del punto, la parte entera del valor absoluto consta de una sola cifra: − 111 /2 = − 5. 55 × 101. Esto excluye el n´umero 0, que queda como excepci´on, y ha de ser representado aparte. En matem´aticas es tambi´en frecuente otra normalizaci´on, colocando la primera cifra no nula inmediatamente despu´es del punto; el valor absoluto de la mantisa es siempre menor o igual que 1 y su parte entera es nula: − 111 /2 = − 0. 555 × 102. La diferencia entre ambas normalizaciones es inesencial: todo se reduce a sumar o restar 1 al exponente. Lo importante es que el punto decimal flota para eliminar los ceros ‘innecesarios’: 2 /111 = 0. 018018018... = 1. 8018018 × 10 −^2 , o bien 2/111 = 0. 18018018 × 10 −^1.
168 CAP´ITULO 6. ARITM ETICA DE PUNTO FLOTANTE.´
Por ser la notaci´on m´as c´omoda, ´esta es la forma de los valores que vamos a usar para definir los sistemas num´ericos en punto flotante.
Definici´on 6.1.1. Un sistema num´erico en punto flotante es un subconjunto finito de R, F = F (β, p, emin, emax), caracterizado por cuatro par´ametros enteros:
- β ∈ Z (la base o radix de F )
- p ∈ Z (la precisi´on de F )
- emin ∈ Z (el menor exponente de F )
- emax ∈ Z (el mayor exponente de F )
Dados valores espec´ıficos de β, p, emin, emax, se define
FN = {y = ±m × βe−p^ : m, e ∈ Z, βp−^1 ≤ m ≤ βp^ − 1 , emin ≤ e ≤ emax},
y finalmente
F = { 0 } ∪ FN
El mayor y el menor n´umero positivo de F son
f max := m´ax{y ∈ F : y > 0 } = βemax^ (1 − β−p) f minN := m´ın{y ∈ FN : y > 0 } = βemin−^1
Llamaremos rango de F al conjunto [−f max, −f minN ] ∪ { 0 } ∪ [f minN , f max] (que contiene a F estrictamente).
N´otese que ahora la mantisa m de un n´umero y = 0.d 1 d 2... dp × βe^ de FN es el entero d 1 d 2... dp (escritos en base β). Llamaremos a d 1 el d´ıgito m´as significativo y a dp el d´ıgito menos significativo de y.
En la norma LIA-1, los par´ametros β y p han de satisfacer
β ≥ 2 y p ≥ 2 ,
y β tiene que ser par. As´ı mismo, “para que haya bastantes elementos en F ”, los par´ametros emin, emax y p han de satisfacer
p − 2 ≤ −emin ≤ βp^ − 1 , p ≤ emax ≤ βp^ − 1.
En IEEE 754 y 854 se consideran cuatro formatos distintos: precisi´on simple, precisi´on doble, precisi´on simple extendida, precisi´on doble extendida, con valores para los par´ametros (ver [Hig], p. 41):
formato β p emin emax IEEE simple 2 24 − 125 128 IEEE doble 2 53 − 1021 1024 IEEE simple extendida 2 ≥ 32 ≤ − 1021 ≥ 1024 IEEE doble extendida 2 ≥ 64 ≤ − 16381 ≥ 16384
(reajustados a la normalizaci´on 〈〈^0 .d 1.. .〉〉).
6.1. REPRESENTACI ON EN PUNTO FLOTANTE´ 169
Ejemplos
1.- Si β = 2, p = 5 [= 2^2 + 1], y = 11/2 = 2^2 + 1 + 1/2 = 2^3 (2−^1 + 2−^3 + 2−^4 ), tendr´ıamos (con mantisa y exponente en base 2, escrito el cero como O y el uno como I para distinguir) y = O.IOIIO × 2 II^ = IOIIO × 2 II−IOI^ ; la mantisa m es IOIIO y el exponente e es II.
2.- ‘Minimodelo’ con β = 2, p = 3, emin = − 1 , emax = 3. Los n´umeros no negativos de F ser´ıan (en escritura decimal) el 0 y
fracci´on 1 / 2 1 /2 + 1/ 8 1 /2 + 1/ 4 1 /2 + 1/4 + 1/ 8 exponente
- 25 0. 3125 0. 3750 0. 4375 × 2 −^1
- 5 0. 625 0. 750 0. 875 × 20
- 0 1. 25 1. 50 1. 75 × 21
- 0 2. 5 3. 0 3. 5 × 22
- 0 5. 0 6. 0 7. 0 × 23
que se distribuyen gr´aficamente as´ı:
Junto con sus opuestos (los sim´etricos respecto del origen) forman el sistema F = F (2, 3 , − 1 , 3). Este ejemplo ilustra c´omo el espaciado de los n´umeros de F salta un factor 2 (= β) con cada potencia de 2 (= β).
6.1.2. Aproximaci´on por redondeo
Elegido un determinado sistema num´erico en punto flotante para implementar en el ordenador, el resto de los n´umeros reales (¡casi todos!) no existen para ´el. Esta situaci´on no es nueva en absoluto: cada vez que hemos usado un valor num´erico de π, hemos puesto π = 3.14, o π = 3.1416, o π = 3.141592654 (cuando disponemos de calculadora); no hemos podido manejar su ‘valor exacto’, que, como dec´ıa el matem´atico, ‘es π’. La realidad fuerza inevitablemente a tomar aproximaciones, por lo que, si no podemos evitar los errores que generan, tratemos al menos de minimizarlos. Eso hacemos cuando ponemos π = 3.1416 (aproximaci´on por exceso) y no π = 3.1415 (aproxi- maci´on por defecto o por corte): tomamos el valor m´as pr´oximo a π entre los que tienen cinco cifras decimales. Para simplificar(*), definiremos:
Definici´on 6.1.2 (Redondeo y desbordamiento.). Dado F = F (β, p, emin, emax), el redondeo ( respecto de F ) es la aplicaci´on f l : R → F ∪ {overflow, underflow} dada por
f l(x) = overflow si |x| > f max = m´ax F (hay un desbordamiento por exceso) f l(x) = underflow si 0 < |x| < f minN = m´ın{|y| : y ∈ F, y 6 = 0} (hay un desbordamiento por defecto).
f l(x) = elemento de F m´as pr´oximo a x si x est´a en el rango de F. Cuando haya dos n´umeros en F a la misma distancia de x, convendremos en tomar como f l(x) el que tenga el ´ultimo d´ıgito dp par.
(*)La norma IEEE 754 contempla distintos modos de redondeo: redondeo hacia +∞ (round down), redondeo hacia −∞ (round up), redondeo hacia 0 (round towards zero) y redondeo al m´as pr´oximo (round to nearest). Ver [Gold]
6.1. REPRESENTACI ON EN PUNTO FLOTANTE´ 171
Pero seg´un dice [Hig], ‘la cantidad m´as ´util asociada con F , que est´a por todas partes en el mundo del an´alisis de los errores de redondeo’, es la llamada unidad de redondeo, que no es otra cosa que la mitad del ´epsilon de la m´aquina.
Definici´on 6.1.5. Dado un sistema num´erico en punto flotante F = F (β, p, emin, emax), el n´umero
u =
β^1 −p^ =
M se llama unidad de redondeo del sistema F.
En el est´andar IEEE con simple precisi´on se tiene u = 2−^24 ≈ 5. 96 × 10 −^8 ; en doble precisi´on, u = 2−^53 ≈ 1. 11 × 10 −^16. A partir del lema anterior, se obtiene:
Proposici´on 6.1.6. Si x ∈ R est´a en el rango de F , entonces
| f l(x) − x| ≤ u |x| =
β^1 −p^ |x|;
precisando m´as,
f l(x) = x(1 + δ), |δ| < u =
β^1 −p.
Demostraci´on. Ver [Hig], p. 42, Th. 2.2.
Esta proposici´on nos da cotas superiores del error absoluto y del error relativo en la aproxima- ci´on de x por f l(x). Recordamos la definici´on de estos conceptos.
Definici´on 6.1.7. Sea x un n´umero real. Si ˜x es otro n´umero real, que tomamos como una aproxi- maci´on de x, se llama error absoluto (de la aproximaci´on) a la diferencia
Ea(x) = |x − x˜|,
y error relativo (cuando x 6 = 0) al cociente
Er(x) = Ea(x) |x|
|x − x˜| |x|
∣∣ 1 − ˜x x
∣∣x˜ x
As´ı ˜x = x + α con |α| = Ea(x); y para x 6 = 0, x˜ = x(1 + ρ) con |ρ| = Er(x).
La proposici´on 6.1.6 significa, por tanto, que el error relativo del redondeo es siempre estricta- mente menor que la unidad de redondeo.
Nota. En algunos textos, el error se expresa en t´erminos de la unidad en el ´ultimo lugar, que se define para un y ∈ F por
ulp(y) = ulp(± 0 .d 1 d 2... dp × βe) = 0. 00... 01 × βe^ = βe−p
(ulp: unity in the last place). As´ı, por ejemplo, si F corresponde a β = 10 y p = 3, y se redondea .0314159 a. 314 × 10 −^1 , entonces hay un error de .159 unidades en el ´ultimo lugar. Dada la distri- buci´on de los elementos de F , es preferible usar el error relativo, m´as preciso; mientras que para x = 0.9994999999 tendr´ıamos f l(x) = 0. 999 × 100 , con ulp(f l(x)) = 10−^3 y un error de 0. 4999999 ulp ≈ 0 .5 ulp, el error relativo es 0. 0005002502252 ≈ 0. 1 u, y sin embargo para x = 1.005 ser´ıa f l(x) = 0. 100 × 101 , con ulp(f l(x)) = 10−^2 , un error de 0.5 ulp igualmente, pero el error relativo correspondiente ser´ıa 0. 004975124378 ≈ u, aproximadamente diez veces mayor que el anterior (este es el fen´omeno llamado wobbling precision, ‘precisi´on bamboleante’ u ‘oscilante’).
172 CAP´ITULO 6. ARITM ETICA DE PUNTO FLOTANTE.´
6.2. Aritm´etica ‘aproximada’
Pr´acticamente en todas las operaciones que realiza un ordenador los datos llegan inevitablemen- te afectados de errores, o se generan errores de redondeo en el curso de los c´alculos. Al efectuar las operaciones m´as elementales, las aritm´eticas (sumas, restas, multiplicaciones, divisiones) los errores de los resultados ¿son del mismo orden que los iniciales o, por el contrario, se produce un aumento significativo en los mismos? ¿son las propias operaciones una nueva fuente de error, se realicen como se realicen? Responder con precisi´on a esta pregunta es una tarea ciertamente complicada, y es el objeto de estudios espec´ıficos fuera de lugar en esta asignatura. Pero podemos hacernos una idea de las dificultades y los problemas que se suscitan examinando qu´e sucede en nuestro modelo de sistemas de n´umeros en punto flotante. Para empezar, ¿c´omo se propagan los errores en una suma?
Lema 6.2.1. Sean x, y, αx, αy, εa(x), εa(y), n´umeros reales tales que |αx| ≤ εa(x), |αy| ≤ εa(y), y sean ˜x = x + αx y ˜ = y + αy
. Entonces Ea = |(x + y) − (x˜ + ˜y)| ≤ |αx| + |αy|,
pudiendo tomar Ea cualquier valor desde Ea = 0 hasta Ea = εa(x) + εa(y). An´alogamente, si x + y 6 = 0, y para ρx, ρy, es
x ˜ = x(1 + ρx) y ˜ = y(1 + ρy)
, entonces Er = |(x + y) − (x˜ + ˜y)| |x + y|
|x ρx + y ρy| |x + y|
y nada puede decirse en general: Er puede tomar cualquier valor en [0, +∞).
Demostraci´on. Por la desigualdad triangular
Ea = |(x + y) − (˜x + ˜y)| = | − (αx + αy)| ≤ |αx| + |αy|.
Como cualquier valor r ∈ [0, εa(x) + εa(y)] puede ponerse en la forma r = s + t con 0 ≤ s ≤ εa(x), 0 ≤ t ≤ εa(y) (¿por qu´e?), Ea puede alcanzar el valor r. En cuanto a Er, sustituyendo ˜x y y˜ por su expresi´on en t´erminos de x e y aparece la igualdad del enunciado. Que Er puede alcanzar cualquier valor no negativo es claro: tomemos, por ejemplo, x > 0, h > 0 cualesquiera, y = −x + h, ρx ≥ 0 arbitrario, ρy = 0; as´ı x + y = h y
Er = xρx h
que puede tomar tanto el valor Er = 0 para ρx = 0 como cualquier valor positivo R (basta que sea ρx > 0, h = xρx R
Vemos, pues, que aunque el error absoluto de los sumandos est´e acotado por un cierto εa, el de la suma puede llegar a ser 2εa, y al cabo de N sumas, N εa, y que el error relativo queda descontrolado. Como generalmente lo que importa es el error relativo, afinaremos el resultado anterior.
Proposici´on 6.2.2. Sean x, y, ρx, ρy, εr(x), εr(y), n´umeros reales tales que x + y 6 = 0 y |ρx| ≤ εr(x), |ρy| ≤ εr(y), y sean x ˜ = x(1 + ρx) ˜y = y(1 + ρy)
Si x e y tienen el mismo signo,
Er ≤ m´ax{|ρx|, |ρy|},
pudiendo tomar Er cualquier valor desde Er = 0 hasta Er = m´ax{εr(x), εr(y)}.
174 CAP´ITULO 6. ARITM ETICA DE PUNTO FLOTANTE.´
Veamos algunos ejemplos de los errores que se originan al realizar estas operaciones. Con base 10 y precisi´on p = 5, tomamos
y = 5/ 7 f l(y) = 0. 71428 × 100 u = 0. 714251 f l(u) = 0. 71425 × 100 v = 98765. 9 f l(v) = 0. 98765 × 105 w = 0. 111111 × 10 −^4 f l(w) = 0. 11111 × 10 −^4
Calculando, obtenemos:
Operaci´on Resultado Error absoluto Error relativo y u 0. 30000 × 10 −^4 0. 472 × 10 −^5 0. 136 (y u) w 0. 27000 × 101 0. 425 0. 136 (y u) w 0. 29629 × 101 0. 466 0. 136 u ⊕ v 0. 98766 × 105 0. 162 × 101 0. 164 × 10 −^4 y w 0. 64286 × 105 0. 779 0. 122 × 10 −^4
Se observa que el comportamiento es bastante desigual. Adem´as, con estas nuevas operaciones se pierden algunas de las buen´ısimas propiedades algebraicas y estructurales de R. Comprob´emoslo. Siguiendo con base 10 y precisi´on p = 5, prescindiendo del 0 previo al punto por comodidad de escritura, tomamos
x = 1000000 f l(x) =. 10000 × 107 y = − 1000000 f l(y) = −. 10000 × 107 z = 1 f l(z) =. 10000 × 101
Calculando, resulta:
x + y x ⊕ y y + z y ⊕ z x + y + z (x ⊕ y) ⊕ z x ⊕ (y ⊕ z)
- 00000 × 100 − 999999 −. 10000 × 107 1. 10000 × 101 (= 1 ). 00000 × 100 (= 0 )
¡No hay asociatividad para la ‘suma flotante’ ⊕! Ni funciona tampoco la propiedad cancelativa: en este mismo ejemplo,
y ⊕ z = y ⊕ 0 6 =⇒ z = 0.
Tambi´en podemos observar que para x = 3 × 107 , por ejemplo,
f l(1/x) =. 33333 × 10 −^7 , x (1/x) = f l((. 30000 × 108 )(. 33333 × 10 −^7 )) =. 99999 ,
mientras que si y =. 33334 × 10 −^7 , z =. 33335 × 10 −^7 ,
x y = y x = x z = z x =. 10000 × 101 = 1.
Esto implica que tampoco el producto es asociativo: para estos valores,
z (x y) = z 1 = z 6 = y = 1 y = (z x) y.
Adem´as, x tiene dos inversos, y y z, ninguno de los cuales coincide, por cierto, con f l(1/x).
6.2. ARITM ETICA ‘APROXIMADA’´ 175
Comparando, en general, con las propiedades fundamentales de R, encontramos:
PROPIEDADES VALIDEZ de cuerpo Suma cerrada NO Puede haber underflow/overflow Suma asociativa NO a veces (a ⊕ b) ⊕ c 6 = a ⊕ (b ⊕ c) Suma conmutativa SI siempre a ⊕ b = b ⊕ a Elemento neutro SI siempre 0 ∈ F y a ⊕ 0 = 0 ⊕ a = a para todo a Elemento opuesto SI siempre −a ∈ F y a ⊕ (−a) = (−a) ⊕ a = 0
Producto cerrado NO Puede haber underflow/overflow Producto asociativo NO a veces (a b) c 6 = a (b c) Producto conmutativo SI siempre a b = b a Elemento unidad SI siempre 1 ∈ F y 1 a = a 1 = a para todo a Elemento inverso ?? ???
Distributiva NO a veces a (b ⊕ c) 6 = (a b) ⊕ (a c) de cuerpo ordenado Reflexiva SI siempre a ≤ a Sim´etrica SI siempre a ≤ b y b ≤ a implica a = b (a, b ∈ F ) Transitiva SI siempre a ≤ b, b ≤ c implica a ≤ c Invariancia por traslaci´on SI siempre a ≤ b implica a ⊕ c ≤ b ⊕ c Invariancia por escalado SI siempre a ≤ b, 0 ≤ c implica a c ≤ b c Completitud SI trivial: F es finito Arquimediana SI trivial: F es acotado Densidad NO F es finito topol´ogicas Conexi´on NO todos los puntos son aislados
Revisemos ahora el funcionamiento de estas operaciones y el error cometido al tomarlas como aproximaciones de las operaciones exactas.
Errores en sumas y restas
Estudiaremos solamente x ⊕ y y x y cuando x y y son elementos positivos del rango de F , porque todos los dem´as casos puede reducirse a ´este teniendo en cuenta que f l(−x) = − f l(x), x y = x ⊕ (−y), etc.
Proposici´on 6.2.4. Dados x, y ≥ f minN , tales que x ⊕ y ∈ F , si
f l(x) = x(1 + ρx) (|ρx| < u), f l(y) = y(1 + ρy) (|ρy| < u),
se tiene
Er = |(x + y) − (x ⊕ y)| |x + y|
< 2 u + u^2 ,
y as´ı x ⊕ y = (x + y)(1 + ρ), |ρ| < 2 u + u^2 ≈ M.
Demostraci´on. Si z = f l(x) + f l(y) = x + y + xρx + yρy, x ⊕ y = f l(z) = z(1 + ρz ) = z + zρz (|ρz | < u), ser´a
x ⊕ y − (x + y) = x ⊕ y − z + z − (x + y) = zρz + xρx + yρy = (x + y + xρx + yρy)ρz + xρx + yρy = (x + y)ρz + (xρx + yρy)(ρz + 1),
6.2. ARITM ETICA ‘APROXIMADA’´ 177
Proposici´on 6.2.5. Dados n´umeros reales no nulos x e y, tales que x y ∈ F , si
f l(x) = x(1 + ρx) (|ρx| < u), f l(y) = y(1 + ρy) (|ρy| < u),
se tiene
Er = |(x y) − (x y)| |x y| < 3 u + 3u^2 + u^3 ,
o sea
x y = (x y)(1 + ρ), |ρ| < 3 u + 3u^2 + u^3 ≈ 3 u.
Demostraci´on. Poniendo z = f l(x) f l(y) = xy(1 + ρx)(1 + ρy), x y = f l(z) = z(1 + ρz ) (|ρz | < u), queda x y = xy(1 + ρx)(1 + ρy)(1 + ρz ), y as´ı
Er = |(x y) − (x y)| |x y|
= |(1 + ρx)(1 + ρy)(1 + ρz ) − 1 |
= |ρx + ρy + ρz + ρxρy + ρxρz + ρyρz + ρxρyρz | < 3 u + 3u^2 + u^3 ≈ 3 u,
suponiendo 3u^2 + u^3 ≈ 0.
An´alogamente, para la operaci´on :
Proposici´on 6.2.6. Dados n´umeros reales x e y, tales que x y ∈ F , si
f l(x) = x(1 + ρx) (|ρx| < u), f l(y) = y(1 + ρy) (|ρy| < u),
se tiene
Er = |(x/y) − (x y)| |x/y|
3 u + u^2 1 − u
o sea
x y = (x/y)(1 + ρ), |ρ| <
3 u + u^2 1 − u ≈ 3 u.
Demostraci´on. Poniendo z = f l(x)/ f l(y) = (x/y)
1 + ρx 1 + ρy , x y = f l(z) = z(1 + ρz ) (|ρz | < u),
queda x y = (x/y) (1 + ρx)(1 + ρz ) 1 + ρy , y as´ı, usando que |1 + ρy| ≥ 1 − |ρy| > 1 − u > 0,
Er =
|(x y) − (x/y)| |x/y|
∣∣^ (1 +^ ρx)(1 +^ ρz^ ) 1 + ρy
(1 + ρx)(1 + ρz ) − (1 + ρy) 1 + ρy
ρx + ρz − ρy + ρxρz 1 + ρy
3 u + u^2 1 − u
≈ 3 u
suponiendo 4u^2 /(1 − u) ≈ 0.
178 CAP´ITULO 6. ARITM ETICA DE PUNTO FLOTANTE.´
6.3. Estrategias de c´alculo para reducir los errores de redondeo
La computaci´on en punto flotante es inexacta por naturaleza, y no es dif´ıcil utilizarla mal de modo que las respuestas computadas consistan casi enteramente en ‘ruido’. Uno de los principales problemas del an´alisis num´erico es determinar c´omo ser´an de exactos los resultados de ciertos m´etodos num´ericos; aqu´ı est´a involu- crado un problema de ‘salto de credibilidad’: no sabemos cu´antas respuestas del ordenador creer. Los usuarios novatos del ordenador resuelven este problema confiando impl´ıci- tamente en el ordenador como autoridad infalible; tienden a creer que todos los d´ıgitos de una respuesta impresa son significativos. Los usuarios desilusionados del ordenador mantienen justamente la actitud opuesta, constantemente temen que sus respuestas carezcan casi de sentido. (D. E. Knuth) El an´alisis elemental de las operaciones que acabamos de efectuar muestra que es necesario un estudio previo a cualquier proceso de c´alculo con ordenador, para evitar una acumulaci´on exage- rada de errores, especialmente cuando haya que realizar operaciones consecutivas (a veces en gran n´umero) con el ordenador. Hay tarea para personas con ingenio hasta en los casos m´as simples, como se˜nala J. M. Muller en su art´ıculo Ordenadores en busca de aritm´etica [MC1, pp. 92–99], al- tamente recomendable: 〈〈^... Si el lector deduce que no hay que depositar jam´as una confianza ciega en los resultados facilitados por un ordenador y que la implantaci´on r´apida de una operaci´on tan elemental como la suma compete todav´ıa al campo de la investigaci´on, me dar´e por absolutamente satisfecho. 〉〉 Los ejemplos siguientes ilustran algunas de las dificultades y muestran caminos para paliar sus consecuencias.
6.3.1. Reformulaci´on del problema: Ejemplos
C´alculo de una diferencia de cuadrados
Consideremos la evaluaci´on de x^2 − y^2 cuando x = 1.21, y = 1.20 en un sistema con β = 10 y p = 3. Resulta
x y x^2 y^2 x^2 − y^2 Error relativo Valor exacto 1. 21 1. 20 1. 4641 1. 4400 0. 0241 Redondeo 0. 121 × 10 0. 120 × 10 0. 146 × 10 0. 144 × 10 0. 200 × 10 −^1 0 , 170... ≈ 34 u
ya que en este caso u =
10 −^2 = 0.005.
Pero x^2 − y^2 = (x − y)(x + y). ¿Qu´e sucede si rehacemos el c´alculo usando la segunda expresi´on?
x y x − y x + y (x − y)(x + y) Error relativo Valor exacto 1. 21 1. 20 0. 01 2. 41 0. 0241 Redondeo 0. 121 × 10 0. 120 × 10 0. 1 × 10 −^1 0. 241 × 10 0. 241 × 10 −^1
¡El resultado es exacto! Todas las cantidades x, y, x + y, x − y, (x + y)(x − y) son ‘n´umeros de m´aquina’ y no hemos introducido ning´un error, mientras que al calcular x x ha sido necesario el redondeo y x x − y y ha originado una cancelaci´on con los dos ´ultimos d´ıgitos falsos. Si bien el ejemplo es un poco ‘tramposo’, incluso cuando x e y deban redondearse la segunda expresi´on es preferible para el c´alculo, puesto que si f l(x) y f l(y) tienen varias cifras iniciales coincidentes, sus cuadrados tendr´an iguales doble n´umero de d´ıgitos, y la cancelaci´on es ‘a´un m´as catastr´ofica’.
180 CAP´ITULO 6. ARITM ETICA DE PUNTO FLOTANTE.´
Veamos una ilustraci´on clara (muy simplificada). Con precisi´on p = 3 y base 10, y sumandos 100, 0.9, 0.9, 0.9, podemos proceder as´ı:
(((100 ⊕ 0 .9) ⊕ 0 .9) ⊕ 0. 9 = (100 ⊕ 0 .9) ⊕ 0. 9 = 100 ⊕ 0. 9 = 100 (100 ⊕ 0 .9) ⊕ (0. 9 ⊕ 0 .9) = 100 ⊕ 1. 8 = 102 ((0. 9 ⊕ 0 .9) ⊕ 0 .9) ⊕ 100 = (1. 8 ⊕ 0 .9) ⊕ 100 = 2. 7 ⊕ 100 = 103 valor exacto: 100 + 0.9 + 0.9 + 0. 9 = 100 .9 + 0.9 + 0. 9 = 101 .8 + 0. 9 = 102. 7
En las evaluaciones primera y tercera hemos procedido recurrentemente: en cada etapa sumamos la suma parcial previa con el sumando siguiente. En la segunda hemos emparejado los sumandos, reduciendo el n´umero de operaciones. Observamos que el primer procedimiento da mayor error, el segundo es ‘menos costoso’, y que el tercero da el valor redondeado del resultado exacto, el mejor valor posible. La causa del ´exito final es aqu´ı muy aparente: al sumar n´umeros de tama˜no muy distinto, hay que tener precauci´on de que ‘el pez grande no se coma al chico’, como sucede en la evaluaci´on 100 ⊕ 0 .9 = 100. Parece conveniente ordenar los sumandos de menor a mayor para homogeneizar tama˜nos en la medida de lo posible, permitiendo as´ı que la adici´on de los sumandos peque˜nos lleve a un valor que no desaparezca cuando lleguen los sumandos ‘grandes’. (Esta es una peque˜na muestra de las dificultades que tienen que superar los programadores incluso en situaciones que no parecen a priori tan complicadas). Adem´as, los dos primeros sumandos est´a involucrados en n − 1 operaciones y el ´ultimo en una s´ola, con lo que el error absoluto final es menor. Por otra parte, cuando la rapidez de c´alculo es importante, la opci´on del emparejamiento parece ventajosa. Adem´as, al disminuir el n´umero de operaciones efecuadas, cabe suponer que disminuir´a el crecimiento de los errores (con datos de magnitudes similares, al menos). N´otese que en una suma de n = 2k^ sumandos, cada sumando estar´ıa ahora involucrado en k = log 2 (n) operaciones, lo que lleva a un error absoluto del orden de n log 2 (n) veces el m´aximo de los errores de los sumandos, ver [1] p. 15. Sobre las dificultades de la implementaci´on de algoritmos eficaces en la adici´on de muchos sumandos, ver m´as detalles en [3], Floating point summation, y [4], as´ı c´omo en el art´ıculo de J. P. Muller ya citado [MC1]. En [Brz], pp. 10–11, hay un ejemplo de ‘correcci´on de la aritm´etica’ para efectuar sumas.
Costo operativo y eficiencia: Evaluaci´on de polinomios
〈〈 (^) Un mismo problema suele poder resolverse por diversos m´etodos num´ericos. ¿C´omo decidirse por uno? El tama˜no de los errores es ciertamente un criterio ´util: podr´ıamos elegir el m´etodo que diese menores errores. No obstante hay un criterio mejor: nos quedaremos con el m´etodo que, dando errores dentro de unos l´ımites predeterminados, necesite el menor trabajo. Esta preocupaci´on por estudiar la eficiencia de los diversos m´etodos es un rasgo distintivo del an´alisis num´erico. La eficiencia depende tanto del tama˜no de los errores como del costo ope- rativo del m´etodo. Aquellas tareas para las que, en un momento hist´orico dado, no se dispone de un m´etodo aceptablemente eficiente no se pueden llevar a cabo. Ejemplo. M´etodo de Horner (1744–1834) para evaluar polinomios. Supongamos que desea- mos hallar el valor de un polinomio cu´artico a 0 + a 1 x + a 2 x^2 + a 3 x^3 + a 4 x^4 (los ai son n´umeros reales conocidos) para un valor dado de x 0 de la variable x. El m´etodo m´as ingenuo halla sucesivamente los productos x^20 = x 0 ∗ x 0 , x^30 = x^20 ∗ x 0 , x^40 = x^30 ∗ x 0 con un coste operativo de tres multiplicaciones, luego los productos ai ∗ xi 0 (otras cuatro multiplicaciones), y finalmente hace las cuatro sumas/restas. (Que sean efectivamente sumas o restas depende del signo con- creto que vayan teniendo los sumandos; para sumar 48 y -27 hay que restar.) En total cuatro sumas/restas y siete multiplicaciones. Sin embargo haciendo las operaciones en orden
a 0 + x 0 ∗
[
a 1 + x 0 ∗
[
a 2 + x 0 ∗ [a 3 + x 0 ∗ a 4 ]
]]
6.3. ESTRATEGIAS DE C ALCULO´ 181
s´olo son precisas cuatro sumas/restas y cuatro multiplicaciones. Para un polinomio de grado N el m´etodo ingenuo requiere N sumas/restas y 2N − 1 multiplicaciones, el m´etodo (1.6) llamado de Horner o de multiplicaci´on encajada usa N sumas/restas y N multiplicaciones. Como ninguno de los dos m´etodos genera errores —excepto los derivados de usar mantisa de longitud finita— el de Horner, que tiene menor costo operativo es m´as eficiente. Si s´olo se va a evaluar una vez un polinomio y se usa un ordenador no hay diferencia entre usar uno u otro m´etodo. Si, como a veces ocurre, dentro de un programa grande hay que evaluar millones de veces, usar uno u otro m´etodo puede significar tener que esperar el resultado dos horas o s´olo una. Si calculamos con papel y l´apiz en un examen usar un mal m´etodo puede implicar no tener el resultado a tiempo. 〉〉 ([S-S], pp. 13 y 14.)
6.3.2. Refinamientos en la representaci´on y en la precisi´on
Dos recursos que permiten mejorar la exactitud de los c´alculos son el uso de n´umeros desnorma- lizados y de precisi´on m´ultiple.
N´umeros desnormalizados
Para paliar los problemas que ocasiona el redondeo de n´umeros ‘cercanos al cero’, se ampl´ıa a veces el sistema F de n´umeros en punto flotante incluyendo los n´umeros desnormalizados (tambi´en llamados subnormales en ingl´es). Tanto los est´andares IEEE 754 y 854 como LIA-1 y unas cuantas implementaciones ‘no-IEEE’ incluyen n´umeros desnormalizados. Siguiendo con la notaci´on habitual, un n´umero en punto flotante desnormalizado es un n´umero real de la forma
y = ±m × βemin−p
donde m es un entero del intervalo [1, βp−^1 − 1]. La fracci´on correspondiente g queda en el intervalo [β−p, 1 /β −β−p]; su d´ıgito m´as significativo es 0, y su exponente emin. Los n´umeros desnormalizados rellenan parcialmente los “huecos de subdesbordamiento” que se producen entre ±βemin−^1 y 0. Tomados juntos, componen el conjunto que en LIA-1 se denota FD. La admisi´on o no de n´umeros desnormalizados se controla mediante el par´ametro booleano denorm, que toma el valor true o false seg´un se permitan n´umeros desnormalizados o no. En el primer caso, se redefine F = { 0 }∪FN ∪FD; en el segundo, se mantiene la definici´on previa. El m´ınimo n´umero positivo en punto flotante desnormalizado es βemin−p^ = f minN · εM. Los valores de FD est´an igualmente espaciados, con el mismo espacio que el de los n´umeros de F entre βemin−^1 y βemin^ , que es βemin−p^ = f minN · εM. Volviendo al ejemplo con β = 2, p = 3, emin = −1, emax = 3, la representaci´on gr´afica de F con los n´umeros desnormalizados incluidos (2−^4 = 0.0625, 2 × 2 −^4 = 0.125, 3 × 2 −^4 = 0.1875) ser´ıa
↓
Los valores de FD tienen un m´aximo absoluto de error de representaci´on de βemin−p. Sin embargo, como los n´umeros desnormalizados tienen menos de p d´ıgitos de precisi´on, el error relativo de representaci´on puede variar ampliamente. Este error relativo va desde M = β^1 −p^ en el extremo superior de FD hasta 1 en el extremo inferior de FD. Cerca de 0, el error relativo crece sin tope. Siempre que una adici´on o sustracci´on de n´umeros de F produce un resultado en FD, dicho resultado es exacto, el error relativo es cero. Un ejemplo trivial: en el modelo anterior, 1 0 .875 = 0.125 con n´umeros desnormalizados incluidos, 1 0 .875 = 0 en caso contrario.
6.3. ESTRATEGIAS DE C ALCULO´ 183
est´e manchada solamente con un error de afectaci´on. Demos un ejemplo para ilustrar nuestro prop´osito: calculemos a = b × c
memoria acumulador b 0. 842345 → 0. 8423450 c 0. 516287 · 10 −^2 → 0. 5162870 · 10 −^2 a 0. 434892 · 10 −^2 ← 0. 4348917 · 10 −^2
Para una suma, tras la transferencia de los operandos al acumulador, el ordenador comienza por modificar el exponente del operando menor (en valor absoluto) para que tenga el mismo valor que el de el m´as grande; de hecho, el ordenador alinea como nosotros las cifras de las potencias id´enticas de 10, unas sobre otras:
memoria acumulador b 0. 842345 → 0. 8423450 c 0. 516287 · 10 −^2 → 0. 0051629 a 0. 847508 ← 0. 8475079
Y lo mismo, evidentemente, para una sustracci´on. 〉〉
Extra˜namente, hay fabricantes (como cray) que construyen algunas de sus m´aquinas sin au- mentar los bits disponibles para los datos intermedios, limit´andose a recortar ´estos al tama˜no correspondiente a la precisi´on, lo que puede suponer un error relativo a˜nadido hasta de β − 1, es decir, del 100 % cuando β = 2 y del 900 % cuando β = 10. En efecto: [Gold] Con sumandos en F :
Proposici´on 6.3.1. Usando un formato de punto flotante con par´ametros β y p, y computando diferencias usando p d´ıgitos, el error relativo del resultado puede llegar a ser β − 1.
Demostraci´on. Para x =
β = 0. 10 · · · 0 × 100 , α = β − 1, y = α β^2
α βp+^ = 0.αα · · · α × 10 −^1
resulta x − y = β−p−^1 , x y = β−p, Er = β − 1.
En cambio:
Proposici´on 6.3.2. Si x e y son n´umeros en punto flotante en un formato con par´ametros β y p, y si la sustracci´on se hace con p + 1 d´ıgitos (i.e., con un d´ıgito ‘de seguridad’) entonces el error relativo de redondeo del resultado es menor que 2 , donde es el ´epsilon de la m´aquina.
(La adici´on est´a incluida en el teorema anterior ya que x e y pueden ser positivos o negativos.) Para m´aquinas que trabajan sin d´ıgitos extras hay muchos algoritmos que dejan de funcionar correctamente, como lo hacen en caso contrario. Existen adem´as resultados como el siguiente, v´alidos cuando se opera con digitos ‘de seguridad’ (v. [Hig], p. 50):
Teorema 6.3.3 (Sterbenz). Sean x e y n´umeros en punto flotante con y/ 2 ≤ x ≤ 2 y. Si la sustracci´on se realiza con un d´ıgito de seguridad, entonces x − y es computado exactamente (supo- niendo que |x − y| ≥ f minN .)
6.3.3. Buscar algoritmos estables
〈〈 (^) Se puede oir a veces el siguiente argumento ingenuo:
¿Por q´ue perder el tiempo con todos esos peque˜nos errores de punto flotante? Incluso IEEE en precisi´on simple nos da 7 d´ıgitos significativos por lo menos, y seguro que esto es suficiente para la mayor´ıa de las aplicaciones, y adem´as est´a la precisi´on doble, y seguramente esta tiene que bastar.
184 CAP´ITULO 6. ARITM ETICA DE PUNTO FLOTANTE.´
Un contraargumento es que los errores pueden propagarse y crecer r´apidamente, y como son inevitables “peque˜nos” errores de redondeo, existe siempre el peligro de que nuestras computaciones produzcan basura.
Los problemas que tienen una tendencia inherente a la “amplificaci´on de errores” se llaman mal condicionados; requieren usar paquetes de aritm´etica multiprecisi´on y un estrecho control de errores. Un ejemplo cl´asico de problema mal condicionado es cierto polinomio de grado 20; tras multi- plicar uno de sus coeficientes (el que multiplica a la potencia 19) por un factor 1.0000000005, la localizaci´on de la mayor parte de las ra´ıces cambia completamente. En el anterior ejemplo, simplemente entrar los datos (los coeficientes del polinomio) en el ordenador puede ser suficiente para estropear todo el c´alculo, por cuanto la conversi´on de la representaci´on decimal externa a la interna, binaria, puede introducir errores demasiado grandes. Un problema bien condicionado puede tener un algoritmo que amplifique los errores: tal algoritmo se dice que es inestable. En t´erminos sencillos, un algoritmo inestable transforma un “peque˜no” cambio en las cantidades de entrada en un cambio “grande” en las cantidades de salida. 〉〉 ([4], § 4-7.)
La estabilidad de un algoritmo es una de sus propiedades m´as importantes, pero no es sencilla de definir con precisi´on. En [1], p. 13, se describe as´ı: 〈〈 (^) Supongamos que g(x) es la soluci´on exacta de un problema con inputs x, y g∗(x) es el valor calculado por un algoritmo particular. La soluci´on calculada g∗(x) puede ser pensada a menudo como la soluci´on exacta de un problema perturbado con inputs x˜, esto es, g∗(x) = g(˜x). Si siempre hay un ˜x cercano a x tal que g∗(x) = g(x˜), entonces el algoritmo se dice estable. Si no, el algoritmo es inestable. 〉〉 N´otese que la estabilidad es una propiedad relativa al algoritmo empleado para resolver un pro- blema, mientras que el condicionamiento se refiere al propio problema: un problema de computaci´on se dice mal condicionado si peque˜nas perturbaciones en el problema producen perturbaciones grandes en la soluci´on exacta. Adem´as del ejemplo anteriormente se˜nalado sobre c´alculo de ra´ıces de un polinomio, otro problema t´ıpico mal condicionado es el de la resoluci´on de un sistema lineal cuando la matriz de los coeficientes es casi singular: peque˜nos cambios en dichos coeficientes o en los t´erminos independientes pueden llevar a soluciones exactas muy diferentes; como al almacenar los datos en la m´aquina se producen errores de redondeo, es sustancialmente imposible dar soluciones precisas a los problemas mal condicionados con cualquier algoritmo que se utilice. Por contra, a veces un algoritmo inestable puede modificarse convenientemente para lograr uno estable.
Estos son los consejos a los programadores de un experto ([Hig], pp. 30 y 31):
〈〈 (^) Dise˜no de algoritmos estables
No hay una receta simple para dise˜nar algoritmos num´ericamente estables. Mientras que esto ayuda a mantener a los analistas num´ericos en el negocio (¡incluso probando que los algoritmos de los dem´as son inestables!) no es una buena noticia para los cient´ıficos computacionales en general. El mejor consejo es ser consciente de la necesidad de la estabilidad num´erica cuando se dise˜na un algoritmo y no concentrarse exclusivamente en otros aspectos, tales como el costo computacional y la paralelizabilidad. Se pueden dar unas pocas pautas.
- Intentar evitar restar cantidades contaminadas de error (aunque tales restas pueden ser inevi- tables).
