21 de abril de 2015

If 1+x is 1, how much is 1-x?

(The title should be "if (= (+ 1 x) 1), how much is (- 1 x)?", but ...)

Toy floating point

Mathematically, if we had infinite precision, x is 0, 1-x is 1, and we have finished. But if we do the calculations with floating point numbers everything is more complicated. First, if x is not 0 but very small, the result of 1+x and 1-x are also 1.

Almost all modern computers use the IEEE standard, but it is very complicated. To understand what is happening will use a toy version without denormalized numbers, infinities nor NANs. But we will try to keep the part of the behavior that interests us.

The toy version uses:
S|EXP|MANTISSA


  • S is for the sign: it uses a single bit, 0 indicates positive and 1 negative, it doesn't matter too much in what follows, so let's write it as + or - directly.
  • EXP stores the exponent: The idea is that the number is SOMETHING*2EXP. As all of this is in binary, so we use SOMETHING*2EXP instead of SOMETHING*10EXP. If EXP is positive then we have a large number, and if EXP is negative is a fraction. Since it is a toy version we will use magic and assume that EXP can be any integer (In the real version, EXP only takes values ​​in a certain range and also serves to identify the infinite and NAN. [*])
  • MANTISSA: The toy version will assume that we have 4 bits ABCD. It is not very important that it has exactly 4 bits, but it is easier do the calculations with 4 bits instead of the to the twenties of the real version.

If we normalize a number, we multiply or divide it by 2 until we get a number between 1 (included) and 2 (excluded). When we write it in base 2 the integer part is 1 and the rest are the decimal part. Then we can avoid to write the first 1 because it's always there, and then we interpret the mantissa ABCD as the number 1.ABCD2.

Examples

Putting all of this together, the number S|EXP|ABCD is
X = S 2EXP 1.ABCD2

For example +|-3|1001 represents the number
X = + 2-3 1.10012 = + 1/8 (1 + 1 1/2 + 0 1/4 + 0 1/8 + 1 1/16) = 1/8 * 25/16 = 25/128 = 0.1953125

Another example -|2|0011 represents the number
X = - 22 1.00112 = - 4 (1 + 0 + 0 1/2 1 1/8 1/4 + 1/16 + 1) = - 4 * 19/16 = - 19/4 = - 4.75

The astute reader should have noticed already that if we have normalized the number to be something between 1 and 2, the number can never be 0. Or in another words, when we added 1 to ABCD, there is no way to write the 0. In the toy version, there is no 0! (In the actual version this is also stored in EXP, and we have to add this case to [*].) For the calculations we will do, it's not very bad that there is no 0, so we just move on.

The numbers that are closer to 1

With this notation we can't write all the real numbers, only some fractions with a denominator that is a power of 2 (this includes integers).

For example number 1, we can write it as +|0|0000 because
X = + 20 1.00002 = + 1 (1 + 0 1/2 + 0 1/4 + 0 1/8 + 0 1/16) = 1 * 16/16 = 1
(It's strange, but all zeros is not 0 but 1.)

The next number we can write going up as little as possible, is +|0|0001 which is
X = + 20 1.00012 = + 1 (1 + 0 + 0 1/2 0 1/8 1/4 + 1/16 + 1) = 1 * 17/16 = 1.0625 = 1 + 0.0625

We can not write any number between 1 and 1.0625, because there are leaps. (In the notation IEEE mantissa has more bits and the leaps are much smaller, but they are still there.)

Going downward is harder. The numbers of the form +|0|ABCD are between 1 (inclusive) and 2 (excluded). So to grab a smaller number we have to lower the exponent, and look at the numbers +|-1|ABCD. The largest of these is clearly +|-1|1111 representing the number
X = + 2-1 1.11112 = + 1/2 (1 + 1 1/2 + 1 1/4 + 1 1/8 + 1 1/16) = 1/2 * 31/16 = 31/32 = 0.96875 = 1 -  0.03125

In other words, from the 1, there is a leap of 0.0625 = 1/16 upwards a leap of 0.03125 = 1/32 downwards. The interval of numbers that can not be written depends on the range where we are working and it just change at 1 as we saw in this example.

If 1+x is 1, how much is x?

Actually, we are interested only in the case where x is positive. If x is negative the idea is equivalent.

To get a result of 1+x that is not 1, it has to be at least 1+1/16, but in this representation the addition is rounded to the nearest representable number. Halfway between 1 and 1+1/16 is the number 1+1/32. The results below 1+1/32 are rounded to 1 and the results greater than 1+1/32 are rounded to 1+1/16. If the result is exactly 1+1/32 the rounding is not clear, let's suppose that the criterion is to choice for the even side and that it will be rounded to 1. 
Therefore if 1+x equals 1, x has to be between 0 and almost 1/32. (Case 1/32 is doubtful, we should read all the standard to be sure.)

If 1-x is 1, how much is x?

To get a result of 1-x that is not 1, the result has to be most 1-1/32, after rounding. Halfway between 1 and 1-1/32 is the number 1-1/64. The results greater than 1-1/64 are rounded to 1 and the smaller than 1-1/64 are rounded to 1-1/32. If the result is exactly  1-1/64, then it depends again on the rounding details, let's suppose that the criterion is to chooce for the even side and that this will be rounded to 1 too. Therefore if 1-x equals 1, x has to be between 0 and almost 1/64. (Case 1/64 is doubtful, if we had already read all standard we would know what happens.)

If 1+x is 1, how much is 1-x?

What matters is that the range of x (positive) where 1+x is 1 is different from the range of x in which 1-x is 1, and that's going to cause unexpected results.


x
1+x 1-x
0<x<1/64 1 1
x=1/64 1 1?
1/64<x<1/32 1 1-1/32
x=1/32 1? 1-1/32
1/32<x 1+1/16 or more 1-1/32 or less

So the interesting cases are when x is between 1/64 and 1/32 (and perhaps the edges). With exponential notation, when x is between 1/26 and 1/25. With some optimism, we can write this as 1/24+2 and 1/24+1, remembering that the mantissa has 4 bits. Looking carefully at the calculations we did, we can transform this into a proof and rewrite it as 1/2M+2 and 1/2M+1, where M is the number of bits we wrote explicitly in the mantissa.

If 1 + x is 1, how much is 1-x? (IEEE)

Using the IEEE standard rather than the toy version, the most important difference is that instead of having only 4 bits of mantissa we have 24 bits for single (32 bits counting those used to store the sign and exponent) and 52 bits for double (64 bits counting those used to store the sign and exponent). The exponent is no longer magical, and we have to be careful with the NAN, infinite and zero, but luckily these changes don't affect what we want to see. Assuming that the reader had done the proof, following the previous steps, to get a result of 1+x that is x but the result of 1-x is not 1, then x must be in the following range


Standar Total Mantisa Xmin Xmax
Toy ∞ :) 4 1/26 = 1/4 1/24 = 0.0156 1/25 = 1/2 1/24 = 0,0312
Single 32 23 1/225 = 1/4 1/223 = 2.98*10-8 1/224 = 1/2 1/223 = 5.96*10-08
Double 64 52 1/254 = 1/4 1/252 = 5.55*10-17 1/253 = 1/2 1/252 = 1.11*10-16

Experimentally

Let's see the results experimentally. Let's try the following program in Racket. First we make a generic function that takes values ​​of x between 0 and 1/2bits, spliced in 2extra steps. After that, we test it with double precision numbers (like 1.0) and single precision numbers (like 1.f0).

#lang racket

(define (almost-one bits extra one)
  (for ([i (in-range (add1 (expt 2 extra)))])
    (define x (/ i (expt 2 (+ bits extra))))
    (displayln (format "~a/2^~a ~a ~a ~a"
                       (~r (* x (expt 2 bits)) #:precision '(= 4))
                       bits
                       (= (+ one x) one)
                       (= (- one x) one)
                       (~r x #:notation 'exponential)))))

(almost-one 52 4 1.0)
(newline)
(almost-one 23 4 1.f0)




The results are


Double precision (32 bits)  Single precision (16 bits)
(almost-one 52 4 1.0) (almost-one 23 4 1.f0)
0.0000/2^52 #t #t 0e+00
0.0625/2^52 #t #t 1.387779e-17
0.1250/2^52 #t #t 2.775558e-17
0.1875/2^52 #t #t 4.163336e-17
0.2500/2^52 #t #t 5.551115e-17
0.3125/2^52 #t #f 6.938894e-17
0.3750/2^52 #t #f 8.326673e-17
0.4375/2^52 #t #f 9.714451e-17
0.5000/2^52 #t #f 1.110223e-16

0.5625/2^52 #f #f 1.249001e-16
0.6250/2^52 #f #f 1.387779e-16
0.6875/2^52 #f #f 1.526557e-16
0.7500/2^52 #f #f 1.665335e-16
0.8125/2^52 #f #f 1.804112e-16
0.8750/2^52 #f #f 1.94289e-16
0.9375/2^52 #f #f 2.081668e-16
1.0000/2^52 #f #f 2.220446e-16
0.0000/2^23 #t #t 0e+00
0.0625/2^23 #t #t 7.450581e-09
0.1250/2^23 #t #t 1.490116e-08
0.1875/2^23 #t #t 2.235174e-08
0.2500/2^23 #t #t 2.980232e-08
0.3125/2^23 #t #f 3.72529e-08
0.3750/2^23 #t #f 4.470348e-08
0.4375/2^23 #t #f 5.215406e-08
0.5000/2^23 #t #f 5.960464e-08

0.5625/2^23 #f #f 6.705523e-08
0.6250/2^23 #f #f 7.450581e-08
0.6875/2^23 #f #f 8.195639e-08
0.7500/2^23 #f #f 8.940697e-08
0.8125/2^23 #f #f 9.685755e-08
0.8750/2^23 #f #f 1.043081e-07
0.9375/2^23 #f #f 1.117587e-07
1.0000/2^23 #f #f 1.192093e-07

If we increase the extra value, we see that the maximum is always in 1/2 and the minimum gets closer and closer to 1/4, as we marked in the other table.

To think about

  • What about (almost-one 52 4 0.5)?
  • What about (almost-one 52 4 1.5)?
  • What about (almost-one 52 4 1)?
Edit: Thanks to Colin Wright for many comments and corrections on the first published version. 

Si 1+x vale 1, ¿cuánto vale 1-x?

(El título debería ser “si (= (+ 1 x) 1), ¿cuánto vale (- 1 x)?, pero …)

Punto flotante de juguete

Matemáticamente, si tuviéramos infinita precisión, x vale 0, 1-x da 1 y a otra cosa. Pero si hacemos las cuentas con números en punto flotante todo se complica. Para empezar si x es distinto de 0 pero muy chiquito, 1+x y 1-x también valen 1.

Casi todas las computadoras actuales unan la norma IEEE, pero es muy complicada. Para entender lo que está pasando vamos a usar una versión de juguete sin números desnormalizados, ni infinitos, ni NANs. Pero tratando que se comporte igual en la parte que nos interesa.

La versión de juguete usa:
S|EXP|MANTISA


  • S es para el signo: usamos un solo bit, un 0 indica positivo y 1 negativo, mucho no importa en lo que sigue, así que vamos a escribirlo como + ó – directamente.
  • EXP guarda el exponente: La idea es que el número es ALGO*2EXP. Como todo esto es en binario, usamos ALGO*2EXP en vez de ALGO*10EXP. Si EXP es positivo es un número grande, si EXP es negativo es una fracción. Ya que es una versión de juguete vamos a usar magia y suponer que EXP puede ser cualquier número entero (En la versión real, EXP sólo toma valores en cierto rango y también sirve para identificar los infinitos y NAN. [*])
  • MANTISA: En la versión de juguete vamos a suponer que tiene 4 bits ABCD. No es muy importante que sean exactamente 4 bits, pero es más fácil hacer las cuentas con 4 bits que con los veintipico de la versión real.

Si normalizamos a un número lo multiplicamos o dividimos por 2 hasta que esté entre 1 (inclusive) y 2 (excluido). Al escribirlo en base 2 nos queda que la parte entera es 1 y el resto son decimales. Entonces podemos no escribir el primer 1 que está siempre ahí y a la mantisa ABCD la interpretamos como el número 1,ABCD2.

Unos ejemplos

Juntando todo esto, el número S|EXP|ABCD es
X = S 2EXP 1,ABCD2

Por ejemplo +|-3|1001 representa el número
X = + 2-3 1,10012 = + 1/8 (1 + 1 1/2 + 0 1/4 + 0 1/8 + 1 1/16) = 1/8 * 25/16 = 25/128 = 0,1953125

Otro ejemplo -|2|0011 representa el número
X = - 22 1,00112 = - 4 (1 + 0 1/2 + 0 1/4 + 1 1/8 + 1 1/16) = - 4 * 19/16 = - 19/4 = - 4,75

El lector astuto habrá notado que al normalizar el número para que esté entre 1 y 2, el número no puede ser 0. Pensado de otra forma, ya que agregamos el 1 a ABCD, no hay forma de escribir el 0. ¡En la versión de juguete no hay 0! (En la versión real esto también se anota en EXP, y hay que agregar este caso en [*].) Para lo que vamos a hacer no es tan malo que no haya 0, así que sigamos adelante.

Los números más cercanos a 1

Con esta notación no podemos escribir todos los números reales, solo alguna fracciones que tienen en el denominador una potencia de 2 (esto incluye a los enteros).

Por ejemplo el número 1 lo podemos escribir como +|0|0000 porque
X = + 20 1,00002 = + 1 (1 + 0 1/2 + 0 1/4 + 0 1/8 + 0 1/16) = 1 * 16/16 = 1
(Es extraño, pero todos ceros no representa al 0 sino al 1.)

El siguiente número que podemos escribir, subiendo lo menos posible, es +|0|0001 que es
X = + 20 1,00012 = + 1 (1 + 0 1/2 + 0 1/4 + 0 1/8 + 1 1/16) = 1 * 17/16 = 1,0625 = 1 + 0,0625

No podemos escribir ningún número entre 1 y 1,0625, porque hay saltitos. (En la notación IEEE la mantisa tiene más bits y los saltitos son mucho más chiquitos, pero igual hay saltitos.)

Para ir hacia abajo es más difícil. Los números de la forma +|0|ABCD están entre 1(inclusive) y 2(excluido). Así que para agarrar un número más chico tenemos que bajar el exponente, y ver los números +|-1|ABCD. El más grande de ellos es claramente +|-1|1111 que representa el número
X = + 2-1 1,11112 = + 1/2 (1 + 1 1/2 + 1 1/4 + 1 1/8 + 1 1/16) = 1/2 * 31/16 = 31/32 = 0,96875 = 1 – 0,03125

O sea que del 1 para arriba hay un salto de 0,0625 = 1/16 y del 1 para abajo hay un salto de 0,03125 = 1/32. Los saltos de números que no se pueden escribir dependen del rango en que uno esté trabajando y por ejemplo cambian justo en 1.

Si 1+x vale 1, ¿cuánto vale x?

En realidad nos interesa sólo si x es positivo. Si es negativo la idea es equivalente.
Para que 1+x no sea 1, al menos tiene que ser 1+1/16, pero al sumar en esta representación se redondea al número más cercano representable. A la mitad entre 1 y 1+1/16 está 1+1/32. Los resultados menores a 1+1/32 se redondean a 1 y los mayores a 1+1/32 se redondean a 1+1/16. Si el resultado es justo 1+1/32, depende del redondeo, supongo que se utiliza el criterio de elegir para el lado par, así que se debe redondear a 1. Con esto queda para que 1+x valga 1, x tiene que estar entre 0 y casi 1/32. (El caso 1/32 es dudoso, habría que leer todo el estandar para estar seguros.)

Si 1-x vale 1, ¿cuánto vale x?

Para que 1-x no sea 1, el resultado tiene que ser como máximo 1-1/32, después de redondear. A la mitad entre 1 y 1-1/32 está 1-1/64. Los resultados mayores a 1-1/64 se redondean a 1 y los menores a 1-1/64 se redondean a 1-1/32. Si el resultado es justo 1-1/64, depende de los detalles del redondeo, supongo que se utiliza el criterio de elegir para el lado par, así que se debe redondear a 1. Con esto queda para que 1-x valga 1, x tiene que estar entre 0 y casi 1/64. (El caso 1/64 es dudoso, si ya hubiéramos leído todo el estándar sabríamos que pasa.)

Si 1+x vale 1, ¿cuánto vale 1-x?

Lo importante es que el rango de x (positivos) en que 1+x es 1 es distinto del rango de x en que 1-x es 1, y eso nos va a traer resultados inesperados.


x 1+x 1-x
0<x<1/64 1 1
x=1/64 1 ¿1?
1/64<x<1/32 1 1-1/32
x=1/32 ¿1? 1-1/32
1/32<x 1+1/16 o más 1-1/32 o menos

O sea que los casos interesantes son cuándo x está entre 1/64 y 1/32 (y quizás los bordes). Con notación exponencial, esto es cuándo x está entre 1/26 y 1/25. Con un poco de optimismo, podemos escribir esto como 1/24+2 y 1/24+1, pensando que la mantisa tiene 4 bits. Mirando con cuidado las cuentas que hicimos, podemos transformar esto en una demostración y reescribirlo como 1/2M+2 y 1/2M+1, en donde M es la cantidad de bits que escribimos explícitamente en la mantisa.

Si 1+x vale 1, ¿cuánto vale 1-x? (IEEE)

Usando el estandard IEEE en vez de la versión de juguete, la diferencia más importante es que en vez de tener sólo 4 bits de mantisa tenemos 24 bits para los single (32 bits contando los que se usan para guardar el signo y el exponente) y 52 bits para los double (64 bits contando los que se usan para guardar el signo y el exponente). 

El exponente ya no es mágico, y hay que ver bien que pasa con los NAN, infinitos y cero, pero por suerte estos cambios no afectan lo que queremos ver. Suponiendo que el lector habrá hecho la demostración, siguiendo los pasos anteriores tenemos que para que 1+x sea 1 y 1-x sea distinto de 1 x tiene que estar en el siguiente rango


Estándar Total Mantisa Xmin Xmax
Juguete ∞ :) 4 1/26 = 1/4 1/24 = 0.0156 1/25 = 1/2 1/24 = 0,0312
Single 32 23 1/225 = 1/4 1/223 = 2,98*10-8 1/224 = 1/2 1/223 = 5,96*10-08
Double 64 52 1/254 = 1/4 1/252 = 5,55*10-17 1/253 = 1/2 1/252 = 1,11*10-16

Experimentalmente

Veamos como anda esto experimentalmente. Probemos con el siguiente programa en Racket. Primero hacemos una función genérica que toma valores de x entre 0 y 1/2bits, partidos en 2extra pasitos. Después la probamos con números de doble precisión (como 1.0) y de simple precisión (como 1.f0).
#lang racket

(define (almost-one bits extra one)
  (for ([i (in-range (add1 (expt 2 extra)))])
    (define x (/ i (expt 2 (+ bits extra))))
    (displayln (format "~a/2^~a ~a ~a ~a"
                       (~r (* x (expt 2 bits)) #:precision '(= 4))
                       bits
                       (= (+ one x) one)
                       (= (- one x) one)
                       (~r x #:notation 'exponential)))))

(almost-one 52 4 1.0)
(newline)
(almost-one 23 4 1.f0)


Los resultados son


Double precision (32 bits)  Single precision (16 bits)
(almost-one 52 4 1.0) (almost-one 23 4 1.f0)
0.0000/2^52 #t #t 0e+00
0.0625/2^52 #t #t 1.387779e-17
0.1250/2^52 #t #t 2.775558e-17
0.1875/2^52 #t #t 4.163336e-17
0.2500/2^52 #t #t 5.551115e-17
0.3125/2^52 #t #f 6.938894e-17
0.3750/2^52 #t #f 8.326673e-17
0.4375/2^52 #t #f 9.714451e-17
0.5000/2^52 #t #f 1.110223e-16

0.5625/2^52 #f #f 1.249001e-16
0.6250/2^52 #f #f 1.387779e-16
0.6875/2^52 #f #f 1.526557e-16
0.7500/2^52 #f #f 1.665335e-16
0.8125/2^52 #f #f 1.804112e-16
0.8750/2^52 #f #f 1.94289e-16
0.9375/2^52 #f #f 2.081668e-16
1.0000/2^52 #f #f 2.220446e-16
0.0000/2^23 #t #t 0e+00
0.0625/2^23 #t #t 7.450581e-09
0.1250/2^23 #t #t 1.490116e-08
0.1875/2^23 #t #t 2.235174e-08
0.2500/2^23 #t #t 2.980232e-08
0.3125/2^23 #t #f 3.72529e-08
0.3750/2^23 #t #f 4.470348e-08
0.4375/2^23 #t #f 5.215406e-08
0.5000/2^23 #t #f 5.960464e-08

0.5625/2^23 #f #f 6.705523e-08
0.6250/2^23 #f #f 7.450581e-08
0.6875/2^23 #f #f 8.195639e-08
0.7500/2^23 #f #f 8.940697e-08
0.8125/2^23 #f #f 9.685755e-08
0.8750/2^23 #f #f 1.043081e-07
0.9375/2^23 #f #f 1.117587e-07
1.0000/2^23 #f #f 1.192093e-07

Si subimos el valor de extra, vemos que el máximo siempre está en 1/2 y el mínimo se acerca cada vez más a 1/4, como habíamos marcado en la otra tabla.

Para pensar

  • ¿Qué pasa con (almost-one 52 4 0.5)?
  • ¿Qué pasa con (almost-one 52 4 1.5)?
  • ¿Qué pasa con (almost-one 52 4 1)?