Me sorprendió descubrir que es difícil encontrar información sobre este problema en ruso, como si a pocas personas les importara que las bibliotecas matemáticas utilizadas en los compiladores modernos a veces no den un resultado correctamente redondeado. Me preocupa esta situación, ya que solo estoy trabajando en el desarrollo de tales bibliotecas matemáticas. Este problema está bien cubierto en la literatura extranjera, por lo que decidí presentarlo en ruso en forma de ciencia popular, confiando en fuentes occidentales y aún con un poco de experiencia personal.

---

Amigos, para su comodidad, el artículo también está disponible como una presentación en video (casi 34 minutos), este formato es más adecuado para aquellos lectores que tienen dificultades para construir las imágenes matemáticas necesarias en su cabeza, ya que hay mucho material ilustrativo en la presentación. La información del video es completamente idéntica al contenido del artículo. Actúe según su conveniencia.

---

Repito que este no es un artículo científico, sino de divulgación científica, después de leerlo, lo conocerás brevemente.

• Las funciones elementales trascendentales (exp, sin, log, cosh y otras) que funcionan con aritmética de coma flotante se redondean incorrectamente, en ocasiones cometen un error en el último bit.
• La razón de los errores no siempre radica en la pereza o la baja calificación de los desarrolladores, sino en una circunstancia fundamental, que la ciencia moderna aún no ha podido superar.
• «», - .
• , , , , exp2(x) pow(2.0, x).

Para comprender el contenido de este artículo, debe estar familiarizado con el formato de punto flotante IEEE-754. Será suficiente si al menos entiendes que, por ejemplo, esto es: 0x400921FB54442D18 - número pi en formato de doble precisión (binary64, o double), es decir, simplemente entiendes lo que quiero decir con este registro; No necesito poder hacer tales transformaciones sobre la marcha. Y les recordaré los modos de redondeo en este artículo, esta es una parte importante de la historia. También es deseable saber inglés de "programador", porque habrá términos y citas de la literatura occidental, pero puede arreglárselas con un traductor en línea.



Primero, ejemplos, para que comprenda inmediatamente cuál es el tema de la conversación. Ahora daré el código en C ++, pero si este no es su idioma, estoy seguro de que comprenderá fácilmente lo que está escrito. Mire este código:

#include <stdio.h>
#include <cmath>

int main() {
  float x = 0.00296957581304013729095458984375f;  // ,  .
  float z;
  z = exp2f(x);  // z = 2**x  .
  printf ("%.8f\n", z);  //      8   .
  z = powf(2.0f, x);  // z = 2**x  
  printf ("%.8f\n", z);  //   .
  return 0;
}

El número x se escribe intencionalmente con tal número de dígitos significativos para que sea exactamente representable en el tipo flotante, es decir, para que el compilador lo convierta en un código binario sin redondeo. Después de todo, sabes muy bien que algunos compiladores no pueden redondear sin errores (si no lo sabes, indícalo en los comentarios, escribiré un artículo aparte con ejemplos). Más adelante en el programa, necesitamos calcular 2 ^x , pero hagámoslo de dos maneras: la función exp2f (x) y la exponenciación explícita de dos, powf (2.0f, x). El resultado, por supuesto, será diferente, porque dije anteriormente que las funciones elementales no pueden funcionar correctamente en todos los casos, y seleccioné especialmente un ejemplo para mostrar esto. Aquí está el resultado:

1.00206053
1.00206041

Cuatro compiladores me dieron estos valores: Microsoft C ++ (19.00.23026), Intel C ++ 15.0, GCC (6.3.0) y Clang (3.7.0). Se diferencian en un bit menos significativo. Aquí está el código hexadecimal para estos números:

0x3F804385  // 
0x3F804384  // 

Recuerde este ejemplo, en él veremos la esencia del problema un poco más tarde, pero por ahora, para que tenga una impresión más clara, consulte los ejemplos para el tipo de datos de doble precisión (doble, binary64) con algunas otras funciones elementales. Presento los resultados en la tabla. Las respuestas correctas (cuando estén disponibles) tienen * al final.







Espero que no tenga la impresión de que tomé deliberadamente algunas pruebas completamente únicas que difícilmente puede encontrar. Si es así, cocinemos de rodillas una enumeración completa de todos los posibles argumentos fraccionarios para la función 2 ^x para el tipo de datos flotante. Está claro que solo nos interesan los valores de x entre 0 y 1, porque otros argumentos producirán un resultado que solo difiere en el valor en el campo del exponente y no son de interés. Tú mismo entiendes:

$$

$$2^x = 2^{[x]}\cdot2^{\{x\}}.$$

Habiendo escrito dicho programa (el texto oculto estará debajo), verifiqué la función exp2f y cuántos valores erróneos produce en el intervalo x de 0 a 1.







Se desprende del programa siguiente que el número de argumentos x probados fue 197612997. Resulta que, por ejemplo, Microsoft C ++ calcula incorrectamente la función 2 ^x para casi el uno por ciento de ellos. No se regocijen, queridos fans de GCC y Clang, es solo que esta función está implementada correctamente en estos compiladores, pero llena de errores en otros.





No aburriré al lector con estos ejemplos, lo principal aquí fue mostrar que las implementaciones modernas de funciones trascendentales pueden redondear el último bit incorrectamente, y diferentes compiladores cometen errores en diferentes lugares, pero ninguno de ellos funcionará correctamente. Por cierto, el estándar IEEE-754 permite este error en el último bit (del que hablaré a continuación), pero todavía me parece extraño: ok doble, este es un tipo de datos grande, ¡pero el flotador se puede verificar por fuerza bruta! ¿Fue tan difícil de hacer? No es nada difícil y ya he mostrado un ejemplo.



Nuestro código de enumeración contiene una función "autoescrita" de cálculo correcto 2 ^xusando un polinomio aproximado del décimo grado, y fue escrito en unos pocos minutos, porque tales polinomios se derivan automáticamente, por ejemplo, en el sistema de álgebra computacional de Maple. Es suficiente establecer una condición para que el polinomio proporcione 54 bits de precisión (para esta función, 2 ^x ). ¿Por qué 54? Pero pronto lo descubrirás, justo después de que te cuente la esencia del problema y te diga por qué, en principio, ahora es imposible crear funciones trascendentales rápidas y correctas para el tipo de datos de precisión cuádruple (binary128), aunque ya hay intentos de atacar este problema en teoría.

Redondeo predeterminado y el problema con él

Si no está inmerso en el desarrollo de bibliotecas matemáticas, no hay nada de malo en olvidar la regla de redondeo predeterminada para números de coma flotante de acuerdo con el estándar IEEE-754. Por lo tanto, te lo recordaré. Si recuerda todo bien, eche un vistazo al menos al final de esta sección de todos modos, le espera una sorpresa: le mostraré una situación en la que redondear un número puede ser muy difícil.



Puede recordar fácilmente qué es "redondear hacia arriba" (a más infinito), "redondear hacia abajo" (a menos infinito) o "redondear a cero" por el nombre (en todo caso, existe Wikipedia). Las principales dificultades para los programadores surgen con el redondeo "al más cercano, pero en el caso de la misma distancia del más cercano, al que tiene el último dígito par". Sí, así se traduce este modo de redondeo, que la literatura occidental llama en definitiva: "Redondear los lazos más cercanos al par".



Este modo de redondeo se utiliza de forma predeterminada y funciona de la siguiente manera. Si, como resultado de los cálculos, la longitud de la mantisa resultó ser mayor que la que puede acomodar el tipo de datos resultante, el redondeo se realiza al más cercano de dos valores posibles. Sin embargo, puede surgir una situación en la que el número original resultó estar exactamente en el medio entre los dos más cercanos, luego se elige el resultado para el cual el último bit (después del redondeo) resulta ser par, es decir, igual a cero. Considere cuatro ejemplos en los que necesita redondear a dos bits después del punto decimal binario:

1. Redondea 1,00 1 001. El tercer bit después del punto decimal es 1, pero luego hay otro sexto bit, que es 1, lo que significa que el redondeo será hacia arriba, porque el número original está más cerca de 1,01 que de 1,00.
2. 1,001000. , 1,00 1,01, .
3. 1,011000. 1,01 1,10. , .
4. 1,010111. , 1,01, 1,10.

A partir de estos ejemplos, puede parecer que todo es simple, pero no lo es. El hecho es que a veces no podemos decir con certeza si realmente estamos en el medio entre dos valores. Vea un ejemplo. De nuevo, queremos redondear a dos bits después del punto decimal:



1,00 1 0000000000000000000000000000000000001



Ahora es obvio para usted que el redondeo debe ser hacia arriba, es decir, al número 1,01. Sin embargo, está viendo un número con 40 bits después del punto decimal. ¿Qué pasa si su algoritmo no puede proporcionar 40 bits de precisión y solo logra 30 bits? Luego dará otro número:



1.00 1 000000000000000000000000000



Sin saber que en la posición 40 (que el algoritmo no puede calcular) habrá una preciada, redondea este número hacia abajo y obtiene 1,00, lo cual es incorrecto. Redondeó el último bit incorrectamente, ese es el tema de nuestra discusión. De lo anterior, resulta que para obtener solo el segundo bit correcto, ¡debe calcular la función hasta 40 bits! ¡Guauu! ¿Y si la "locomotora" de ceros resulta ser aún más larga? De eso hablaremos en la siguiente sección.



Por cierto, este es el error que cometen muchos compiladores al convertir la notación decimal de un número de punto flotante al formato binario resultante. Si el número decimal original en el código del programa está demasiado cerca del medio entre dos valores binarios representables con precisión, entonces no se redondeará correctamente. Pero este no es el tema de este artículo, sino una razón para una historia separada.

La esencia del problema de redondear el último bit significativo

El problema se manifiesta por dos razones. El primero es el rechazo deliberado de cálculos laboriosos en favor de la velocidad. En este caso, siempre que se observe la precisión especificada, y qué bits habrá en la respuesta es un asunto secundario. La segunda razón es el dilema de Table Maker, que es el tema principal de nuestra conversación. Consideremos ambas razones con más detalle.

Primera razón

Usted, por supuesto, comprende que el cálculo de funciones trascendentales se implementa mediante algunos métodos aproximados, por ejemplo, mediante el método de aproximación de polinomios o incluso (raramente) mediante expansión de series. Para que los cálculos se realicen lo más rápido posible, los desarrolladores acuerdan realizar la menor cantidad posible de iteraciones del método numérico (o tomar un polinomio del menor grado posible), siempre que el algoritmo permita un error que no exceda la mitad del valor del último bit de la mantisa. En la literatura, esto se escribe como 0.5ulp (ulp = unidad en el último lugar ).



Por ejemplo, si estamos hablando de un número x de tipo float en el intervalo (0.5; 1), el valor ulp = 2 ^-23 . En el intervalo (1; 2) ulp = 2 ^-22 . En otras palabras, si x está en el intervalo (0; 1) entonces 2 ^xestará en el intervalo (1,2), y para asegurar una precisión de 0.5ulp, necesita, en términos generales, elegir EPS = 2 ^-23 (ya que denotaremos la constante "épsilon", en la gente común llamada "error" o "precisión", a quien como quieras, no le encuentres falta, por favor).



Para los cálculos aplicados, esto es suficiente, pero el hecho de que los últimos bits puedan no coincidir con el resultado absoluto no es importante para casi el 100% de los programadores, porque no es importante para ellos cuáles serán los bits, sino cuál será la precisión.



Para aquellos que no entiendan, les daré un ejemplo en el sistema numérico decimal. Aquí hay dos números: 1.999999 y 2.0. Digamos que el primero es lo que recibió el programador y el segundo es el estándar de lo que debería haber sucedido si tuviéramos posibilidades ilimitadas. La diferencia entre ellos es solo una millonésima, es decir, la respuesta se calcula con un error de EPS = 10 ^-6 . Sin embargo, no hay un solo número correcto en esta respuesta. ¿Es mala? No, desde el punto de vista del programa de aplicación, esto es violeta, el programador redondeará la respuesta, digamos, a dos decimales y recibirá 2,00 (por ejemplo, se trataba de moneda, $ 2,00), no necesita más, pero el hecho de que puse EPS = 10 ^-6 en mi programa , luego bien hecho, tomé un margen para el error de los cálculos intermedios y resolví el problema correctamente.



En otras palabras, no se confunda: la precisión y el número de bits (o dígitos) correctos son dos cosas diferentes. Aquellos que necesitan precisión (esto es casi el 100% de los programadores), el problema discutido no les concierne en absoluto. Cualquiera que necesite una secuencia de bits para hacer coincidir una referencia correctamente redondeada está muy preocupado por este problema, por ejemplo, los desarrolladores de bibliotecas de funciones elementales. Sin embargo, es útil que todos sepan esto para el desarrollo general.



Permítanme recordarles que esta fue la primera dirección del problema: las últimas partes de la respuesta pueden estar equivocadas porque es una solución intencional. Lo principal es mantener la precisión de 0.5ulp (o más). Por lo tanto, el algoritmo numérico se selecciona solo a partir de esta condición, si solo funciona extremadamente rápido. En este caso, el Estándar permite la implementación de funciones elementales sin redondear correctamente el último bit. Estoy citando [1, sección 12.1]:

La versión de 1985 del estándar IEEE 754 para aritmética de coma flotante no especificó nada sobre la función elemental. Esto se debe a que durante años se ha creído que las funciones correctamente redondeadas serían demasiado lentas al menos para algunos argumentos de entrada. La situación cambió desde entonces y la versión 2008 de la norma recomienda (aunque no exige) que algunas funciones se redondeen correctamente.

Las siguientes son funciones recomendadas pero que no necesitan redondearse correctamente:

La segunda razón

Finalmente llegamos al tema de conversación: Table Maker's Dilemma (abreviado como TMD). No pude traducir adecuadamente su nombre al ruso, fue presentado por William Kahan (padre fundador de IEEE-754) en el artículo [2]. Quizás si lees el artículo, entenderás por qué el nombre es exactamente ese. En resumen, la esencia del dilema es que necesitamos obtener un redondeo absolutamente exacto de la función z = f (x), como si tuviéramos un registro de bits infinitos del resultado perfectamente calculado z a nuestra disposición. Pero está claro para todos que no podemos obtener una secuencia infinita. ¿Cuántos bits tomar entonces? Arriba, mostré un ejemplo en el que necesitamos ver 40 bits del resultado para obtener al menos 2 bits correctos después del redondeo. Y la esencia del problema TMD es que no sabemos de antemano, hasta cuántos bits calcular el valor de z para obtener tantos bits correctos después del redondeo como necesitemos. ¿Y si hay cien o mil? ¡No lo sabemos de antemano!



Por ejemplo, como dije, para la función 2 ^x , para el tipo de datos float, donde la parte fraccionaria de la mantisa tiene solo 23 bits, necesitamos realizar el cálculo con una precisión de ^2-54 para que el redondeo se produzca correctamente para todos los argumentos x posibles sin excepción. No es difícil obtener esta estimación mediante una búsqueda exhaustiva, pero para la mayoría de las otras funciones, especialmente para los tipos double o long double (ponga "clase" si sabe cuál es), tales estimaciones son desconocidas .



Entendamos ya por qué sucede esto. Deliberadamente di el primer ejemplo en este artículo con el tipo de datos flotante y le pedí que lo recordara, porque en este tipo solo hay 32 bits y será más fácil de ver, en otros tipos de datos la situación es similar.



Comenzamos con el número x = 0.00296957581304013729095458984375, este es un número exactamente representable en el tipo de datos flotante, es decir, está escrito de manera que se pueda convertir al sistema flotante binario sin redondeo. Calculamos 2 ^x , y si tuviéramos una calculadora con precisión infinita, entonces deberíamos obtener (Para que puedas comprobarme, el cálculo se hace en el sistema en línea WolframAlpha ):



1.0020604729652405753669743044108123031635398201893943954577320057 ...







Traduzcamos este número al sistema binario, digamos que 64 bits serán suficientes: 1.00000000100001110000100 1 000000000000000000000000000001101111101



El bit de redondeo (24 bits después del punto decimal) está subrayado. Pregunta: ¿dónde redondear? ¿Arriba o abajo? Claramente, usted sabe esto porque ve suficientes bits y puede tomar una decisión. Pero mira con cuidado ...



Después del bit de redondeo, tenemos 29 ceros. Esto significa que estamos muy, muy cerca del medio entre los dos números más cercanos y es suficiente con bajar un poco, ya que la dirección del redondeo cambiará. Pero la pregunta es: ¿dónde será este cambio? El algoritmo numérico puede secuencialmente, paso a paso, acercarse al valor exacto desde diferentes lados, y hasta que pasemos todos estos 29 ceros y alcancemos una precisión que exceda el valor del último cero en esta "locomotora", no sabremos la dirección de redondeo. ... ¿Qué pasa si, en realidad, la respuesta correcta debería ser:



1.00000000100001110000100 0 11111111111111111111111111111?



Entonces el redondeo disminuirá.



No lo sabemos hasta que nuestra precisión alcanza el 54º bit después del punto decimal. Cuando se conozca exactamente el bit 54, sabremos exactamente a cuál de los dos números más cercanos estamos más cerca. Estos números se denominan puntos más difíciles de redondear [1, sección 12.3] (puntos críticos para el redondeo), y el número 54 se denomina dureza a redondear y se indica con la letra m en el libro citado.



La complejidad del redondeo (m) es el número de bits que es el mínimo necesario para asegurar que para todos los argumentos de alguna función f (x) y para un rango preseleccionado, la función f (x) se redondea correctamente al último bit (para diferentes modos de redondeo, puede haber diferentes valor m). En otras palabras, para el tipo de datos float y para el argumento x desde el rango (0; 1) para el modo de redondeo al tiempo de redondeo "par más cercano" m = 54. Esto significa que para absolutamente todo x del intervalo (0; 1) podemos poner en el algoritmo la misma precisión ESP = ^2-54 , y todos los resultados se redondearán correctamente a 23 bits después del punto decimal binario.



De hecho, algunos algoritmos pueden proporcionar un resultado exacto y, en base a 53 e incluso 52 bits, la fuerza bruta lo muestra, pero en teoría, necesitas exactamente 54. Si no fuera por la posibilidad de generar la fuerza bruta, no podríamos "hacer trampa". y guardar un par de bits, como hice en el programa de fuerza bruta anterior. Tomé un polinomio con un grado más bajo de lo que debería, pero aún funciona, solo porque tuve suerte.



Entonces, independientemente del modo de redondeo, tenemos dos situaciones posibles: o surge una "locomotora de vapor" de ceros en el área de redondeo, o una "locomotora de vapor" de unos. La tarea del algoritmo correcto para calcular la función trascendental f (x) es refinar el valor de esta función hasta que la precisión exceda el valor del último bit de esta "locomotora", y hasta que quede claro que como resultado de las fluctuaciones posteriores del algoritmo numérico para calcular f (x) los ceros no se convertirán en unos, o viceversa. Tan pronto como todo se haya estabilizado y el algoritmo haya alcanzado una precisión tal que esté más allá de los límites de la "locomotora de vapor", entonces podemos redondear como si tuviéramos un número infinito de bits. Y este redondeo será con el último bit correcto. Pero, ¿cómo se puede lograr esto?

"Muletas"

Como se mencionó, el principal problema es conseguir que el algoritmo supere la locomotora de ceros o unos que viene inmediatamente después del bit de redondeo. Cuando se supera la locomotora y la vemos como un todo, entonces esto equivale al hecho de que estos ceros o unos ya se han calculado exactamente , y ya sabemos exactamente en qué dirección se producirá ahora el redondeo. Pero si no conocemos la longitud de la locomotora, ¿cómo podemos diseñar un algoritmo?

La primera "muleta"

Al lector le puede parecer que la respuesta es obvia: tome la aritmética con precisión infinita y coloque un número deliberadamente excesivo de bits, y si no es suficiente, coloque otro y vuelva a calcular. En general, es correcto. Esto se hace cuando la velocidad y los recursos de la computadora no juegan un papel especial. Este enfoque tiene un nombre: estrategia multinivel de Ziv [1, sección 12.3]. Su esencia es sumamente simple. El algoritmo debe admitir cálculos en varios niveles: un cálculo preliminar rápido (en la mayoría de los casos resulta ser definitivo), un cálculo más lento pero más preciso (guarda en la mayoría de los casos críticos), un cálculo aún más lento, pero aún más preciso (cuando es absolutamente "malo "Tenía que hacerlo) y así sucesivamente.



En la inmensa mayoría de los casos, basta con tomar la precisión un poco más de 0.5ulp, pero si aparece una "locomotora", la aumentamos. Mientras permanezca la "locomotora de vapor", aumentamos la precisión hasta que quede estrictamente claro que las fluctuaciones adicionales del método numérico no afectarán a esta "locomotora de vapor". Entonces, por ejemplo, en nuestro caso, si hemos llegado a ESP = ^2-54 , entonces en la posición 54 aparece una unidad que, por así decirlo, "protege" a la locomotora de ceros y garantiza que ya no habrá una resta de un valor mayor o igual a ^2-53 y los ceros no se convertirán en unos, arrastrando el bit de redondeo a cero con él.



Fue una presentación de divulgación científica, de todos modos con la prueba de redondeo de Ziv, donde se muestra qué tan rápido, en un solo paso, comprobar si hemos alcanzado la precisión deseada, se puede leer en [1, Capítulo 12], o en [3, Sección 10,5].



El problema con este enfoque es obvio. Es necesario diseñar un algoritmo para el cálculo de cada función trascendental f (x) para que en el transcurso de la pieza sea posible aumentar la precisión de los cálculos. Para una implementación de software, esto todavía no da tanto miedo, por ejemplo, el método de Newton permite, en términos generales, duplicar el número de bits exactos después del punto decimal en cada iteración. Puede doblar hasta que sea "suficiente", aunque este es un proceso que requiere bastante tiempo, debo admitir que el método de Newton no siempre está justificado, porque requiere calcular la función inversa f ^-1(x), que en algunos casos puede no ser más simple que calcular f (x). Para la implementación de hardware, la "estrategia Ziva" es completamente inadecuada. El algoritmo, cableado en el procesador, debe realizar una serie de acciones con el número de bits ya preestablecido, y esto es bastante problemático de implementar si no conocemos este número de antemano. ¿Hacer inventario? ¿Cuánto cuesta?



El enfoque probabilístico para resolver el problema [1, Sección 12.6] nos permite estimar el valor de m (recuerde, este es el número de bits, que es suficiente para un redondeo correcto). Resulta que la longitud de la "locomotora" en el sentido probabilístico es ligeramente mayor que la longitud de la mantisa del número. Así, en la mayoría de los casos será suficiente tomar m un poco más del doble del valor de la mantisa, y solo en casos muy raros será necesario tomar aún más. Cito a los autores de este trabajo: "deducimos que en la práctica, m debe ser ligeramente mayor que 2p" (tienen p - la longitud de la mantisa junto con la parte entera, es decir, p = 24 para float). Más adelante en el texto, muestran que la probabilidad de error con tal estrategia es cercana a cero, pero aún positiva, y esto se confirma mediante experimentos.



Sin embargo, todavía hay casos en los que el valor de m debe tomarse aún más y el peor de los casos no se conoce de antemano. Existen estimaciones teóricas del peor de los casos [1, sección 12.7.2], pero producen millones de bits impensables, lo que no es bueno. Aquí hay una tabla del trabajo citado (esto es para la función exp (x) en el intervalo de -ln (2) a ln (2)):

Segunda "muleta"

En la práctica, m no será tan grande. Y para determinar el peor de los casos, se aplica una segunda "muleta", que se denomina "cálculo previo exhaustivo". Para el tipo de datos float (32 bits), si la función f tiene un argumento (x), entonces podemos "ejecutar" fácilmente todos los valores posibles de x. El problema surgirá solo con funciones que tengan más de un argumento (entre ellas pow (x, y)), para las cuales no podríamos pensar en nada parecido. Después de verificar todos los valores posibles de x, calculamos nuestra constante m para cada función f (x) y para cada modo de redondeo. Luego, los algoritmos de cálculo que deben implementarse en el hardware están diseñados para proporcionar una precisión de 2 ^m . Entonces se garantiza que el redondeo de f (x) es correcto en todos los casos.



Para el tipo doble (64 bits), la enumeración simple es casi imposible. Sin embargo, ¡se están resolviendo! ¿Pero cómo? La respuesta se da en [4]. Te lo contaré muy brevemente.



El dominio de la función f (x) se divide en segmentos muy pequeños de modo que dentro de cada segmento es posible reemplazar f (x) con una función lineal de la forma b-ax (los coeficientes ayb son, por supuesto, diferentes para diferentes segmentos). El tamaño de estos segmentos se calcula analíticamente de modo que tal función lineal sería de hecho casi indistinguible de la original en cada segmento.



Luego, después de algunas operaciones de escalado y desplazamiento, llegamos al siguiente problema: ¿puede un eje b en línea recta "acercarse lo suficiente" a un punto entero?



Resulta que es relativamente fácil dar una respuesta de sí o no. Es decir, "sí" - si los puntos potencialmente peligrosos están cerca de una línea recta, y "no" - si no, tal punto, en principio, puede acercarse a la línea. La belleza de este método es que la respuesta “no” en la práctica se obtiene en la inmensa mayoría de los casos, y la respuesta “sí”, que rara vez se obtiene, obliga a recorrer el segmento con fuerza bruta para determinar qué puntos específicos fueron críticos.



Sin embargo, la iteración sobre los argumentos de f (x) se reduce muchas, muchas veces y permite detectar puntos de ruptura para números como double (binary64) y long double (¡80 bits!). Esto se hace en supercomputadoras y, por supuesto, tarjetas de video ... en tu tiempo libre de minería. Sin embargo, nadie sabe todavía qué hacer con el tipo de datos binary128. Permítame recordarle que la parte fraccionaria de la mantisa de tales números es de 112 bits . Por lo tanto, en la literatura extranjera sobre este tema hasta el momento, solo se pueden encontrar razonamientos semifilosóficos que comienzan con “esperamos ...” (“esperamos ...”).



Los detalles del método, que le permiten determinar rápidamente el paso de una línea cerca de puntos enteros, son inapropiados aquí. Para aquellos que deseen aprender el proceso, recomiendo mirar más detenidamente el problema de encontrar la distancia entre una línea recta y Z ² , por ejemplo, en el artículo [5]. Describe un algoritmo mejorado, que en el curso de la construcción se asemeja al famoso algoritmo de Euclides para encontrar el máximo común divisor. Daré la misma imagen de [4] y [5], que muestra la transformación posterior del problema:







hay tablas extensas que contienen los peores casos de redondeo a diferentes intervalos para cada función trascendental. Se encuentran en [1 sección 12.8.4] y en [3, sección 10.5.3.2], así como en artículos separados, por ejemplo, en [6].



Daré algunos ejemplos tomando filas aleatorias de tales tablas. Hago hincapié en que estos no son los peores casos para todo x, pero solo para algunos intervalos pequeños, consulte la fuente si está interesado.







¿Cómo leer la tabla? El valor x se especifica en notación doble de coma flotante hexadecimal. Primero, como se esperaba, hay uno inicial, luego 52 bits de la parte fraccionaria de la mantisa y la letra P. Esta letra significa "multiplica por dos a una potencia" seguida de un grado. Por ejemplo, P-23 significa que la mantisa especificada debe multiplicarse por 2 ^-23 .



Además, imagine que la función f (x) se calcula con precisión infinita y los primeros 53 bits se cortan de ella (¡sin redondeo!). Son estos 53 bits (uno de ellos hasta la coma) los que se indican en la columna f (x). Los bits posteriores se indican en la última columna. El signo de "grado" de la secuencia de bits en la última columna significa el número de repeticiones de bits, es decir, por ejemplo, 10 ⁵³1011 significa que primero viene el bit igual a 1, luego 53 ceros y luego 1011. Luego la elipsis, lo que significa que, en general, no necesitamos el resto de los bits en absoluto.



Además, es una cuestión de tecnología: conocemos los peores casos para cada intervalo de una función tomada por separado y podemos elegir para este intervalo una aproximación tal que cubra el peor de los casos con su precisión. Con solo años de supercomputación, es posible crear implementaciones de hardware rápidas y precisas de funciones elementales. El asunto es pequeño: queda por enseñar al menos a los desarrolladores de compiladores a utilizar estas tablas.

¿Por qué es necesario?

¡Gran pregunta! Después de todo, he hablado en repetidas ocasiones anteriormente de que casi el 100% de los programadores no necesitan conocer una función elemental hasta el último bit correctamente redondeado (a menudo ni siquiera necesitan la mitad de los bits), ¿por qué los científicos manejan supercomputadoras y compilan tablas para resolver un problema "inútil"?



Primero, el desafío es fundamental. Es bastante interesante no obtener un redondeo exacto en aras de un redondeo exacto, pero en principio para comprender cómo se podría resolver este interesante problema, ¿qué secretos de la matemática computacional nos revelará su solución? ¿Cómo se podrían utilizar estos secretos en otras tareas? Ciencias fundamentales: son así, puedes hacer algún tipo de "tonterías" durante décadas, y luego cien años después, gracias a estas "tonterías", se produce un avance científico en alguna otra área.



En segundo lugar, la cuestión de la portabilidad del código. Si una función puede permitirse manejar los últimos bits del resultado de la forma que desee, significa que en diferentes plataformas y en diferentes compiladores, se pueden obtener resultados ligeramente diferentes (incluso si están dentro del error especificado). En algunos casos esto no es importante, pero en algunos puede ser significativo, especialmente cuando el programa tiene un error que aparece en una plataforma, pero no aparece en otra plataforma precisamente por los diferentes bits del resultado. Pero, ¿por qué le estoy describiendo el conocido dolor de cabeza asociado con diferentes comportamientos de programas? Sabes todo esto sin mí. Sería genial tener un sistema matemático que funcione exactamente igual en todas las plataformas, sin importar cuán compilado esté. Para esto necesitas correctamente Redondea el último bit.

Lista de fuentes

[1] Jean-Michel Muller, "Funciones elementales: algoritmos e implementación", 2016



[2] William Kahan, " Un logaritmo demasiado inteligente por la mitad ", 2004



[3] Jean-Michel Muller, "Manual de aritmética de punto flotante" , 2018



[4] Vincent Lefèvre, Jean-Michel Muller, “Hacia lo trascendental correctamente redondeado”, IEEE TRANSACTIONS ON COMPUTERS, VOL. 47, NO. 11, NOVIEMBRE de 1998. págs. 1235-1243



[5] Vincent Lefèvre. “Nuevos resultados sobre la distancia entre un segmento y Z ² ”. Aplicación al redondeo exacto. 17º Simposio de IEEE sobre aritmética informática - Arith'17, junio de 2005, Cape Cod, MA,

Estados Unidos. págs. 68-75



[6] Vincent Lefèvre, Jean-Michel Muller, “Peores casos para el redondeo correcto de las funciones elementales en doble precisión”, Rapport de recherche (INSTITUT NA TIONAL DE RECHERCHE EN INFORMA TIQUE ET EN AUTOMA TIQUE) n˚4044 - Novembre 2000 - 19 páginas.