Cálculos rápidos y precisos para números de punto flotante utilizando el ejemplo de la función seno. Parte 2: libm

Continúo la serie de artículos sobre cómo trabajar con punto flotante. En el primer artículo, hice una breve introducción matemática y mostré la forma más simple y obvia de calcular el seno con ejemplos de programas con diferentes "trampas". El artículo de hoy tendrá un estilo ligeramente diferente. No habrá práctica aquí, pero profundizaremos en las matemáticas y entraremos en el lugar santísimo: el código de la biblioteca estándar. También daré una respuesta a la pregunta al final del primer artículo. Entonces vamos.

Algo de matemáticas de nuevo

Obviamente, cuanto más rápido disminuye la serie de Taylor en valor absoluto, se necesitan menos términos para lograr la precisión requerida. Y así parece que el resultado será más preciso (a continuación se comentará con más detalle). Para comparar, por ejemplo, tome un término del séptimo grado de la serie de Taylor ($$$\frac{x^7}{7!}$) a $$$x=1.0$ y $$$x=0.1$... Los valores de la expresión serán$$$1.198\cdot 10^{-4}$ y $$$1.198\cdot 10^{-11}$respectivamente. Una gran diferencia, ¿no? Así que intentemos encontrar una manera de reducir el límite superior del intervalo de cálculo de la función seno.

Expansión de la serie alrededor de valores dados

Para comprender este método, es necesario volver al primer año del instituto y recordar las definiciones de la serie de Taylor ( wiki ). En pocas palabras: conociendo la función y sus derivadas en algún momento, puede encontrar los valores de la función en las proximidades de este punto expandiéndolo en una serie de Taylor. Para la función seno, esto significa lo siguiente

$$

$$\sin(x_0+\Delta x) \approx \sin(x_0)+\sin'(x_0)\frac{\Delta x}{1!}+\sin''(x_0)\frac{\Delta x^2}{2!}+\sin'''(x_0)\frac{\Delta x^3}{3!}+\cdots$$

¿Qué nos aporta este enfoque desde un punto de vista práctico? Imagina que tenemos un intervalo de$$$0$ antes de $$$\pi/2$... Elijamos 10 puntos distribuidos linealmente en este intervalo (la elección no es óptima):$$$x_0=0$, $$$x_1=\pi/20$, $$$x_2=2\pi/20$, $$$x_i=i\cdot \pi/20$... Para cada punto, calcule la placa con el seno y sus derivadas en este punto. Ahora puede modificar la función para que al obtener el valor$$$x$ la función toma el valor más cercano $$$x_i$ y se coloca en fila alrededor del punto $$$x_i$, no alrededor de cero$$$\Delta x = x - x_i$).

Usar transformaciones trigonométricas

Si retrocedemos aún más, a las clases superiores de la escuela, entonces podemos recordar una fórmula muy importante:

$$

$$\sin(x_0+\Delta x)=\sin(x_0)\cos(\Delta x)+\cos(x_0)\sin(\Delta x)$$

Y luego todo es igual que en el párrafo anterior. Seleccionamos puntos dentro del intervalo, calculamos el seno y el coseno para ellos, y al llamar a la función seno, buscamos el más cercano y, usando la fórmula anterior, calculamos el seno usando un valor pequeño$$$\Delta x$...

Piense cuál de estos dos métodos es mejor elegir, pero por ahora pasaremos de las matemáticas a los cálculos prácticos.

Propiedad de distribución de la multiplicación en el mundo de coma flotante

Tuve que pedir consejo a Internet, ¿cómo se llama? $$$a\cdot(b+c)=a\cdot b + a\cdot c$... Resulta ser una propiedad distributiva. Volvamos a la pregunta que hice al final de la primera parte. Es decir, por qué expresiones matemáticamente equivalentes$$$a\cdot(b+c)$ y $$$a\cdot b + a\cdot c$puede dar diferentes resultados en cálculos de coma flotante? La forma más sencilla de ilustrar esto es con un ejemplo. Tomemos un sistema hipotético que funciona con números de coma flotante en formato decimal con 4 dígitos de precisión. Pretendamos que$$$b=1.000\text{E}0$, $$$c=1.234\text{E-}2$y $$$a=5.678\text{E-}2$... Primero, tomemos una expresión entre corchetes y calculémosla paso a paso, recordando redondear en cada paso:

1)$$$b+c=1.000\text{E}0 + 1.234\text{E-}2=1.012\text{E}0$

2) $$$a\cdot(b+c)=5.678\text{E-}2\times 1.012\text{E}0=5.746\text{E-}2$

Respuesta recibida $$$y=5.746\text{E-}2$

Ahora calculemos la segunda expresión de la misma manera paso a paso:

$$$a\cdot b + a\cdot c=5.678\text{E-}2\times 1.000\text{E}0 + 5.678\text{E-}2\times 1.234\text{E-}2=5.678\text{E-}2 + 7.007\text{E-}4=5.748\text{E-}2$

Respuesta recibida $$$y=5.748\text{E-}2$

La respuesta verdadera es 0.0574806652.



Como puede ver, la respuesta obtenida en el segundo caso es mucho más cercana a la verdadera que en el primero. Si explicamos esto con los dedos, imagina que cuando en el primer caso sumamos a 1.0 el número$$$c=1.234\text{E-}2$simplemente descartamos los dos últimos dígitos. Ya no existen. En el segundo caso, el descarte se produce al final, después de la multiplicación. Aquellos. en el segundo caso, las operaciones de multiplicación son más precisas.



Parece que puedes terminar con esto, pero mira más de cerca el primer método y dime cuál será el resultado del cálculo.$$$b+c-b$... Y ... ¡tenemos una forma de redondear números de punto flotante! No se pierda este ejemplo. Date tiempo para resolverlo. El redondeo de números se utilizará de forma muy intensiva en este artículo y en los siguientes.



Notemos una característica más de esta expresión. Imagina que la precisión de 4 dígitos de la variable no es suficiente para nosotros. ¿Qué hacer? Y aquí ya tenemos la respuesta: representar el número en la forma$$$b+c$y almacenarlo en la memoria como la suma de dos dígitos. Y, en consecuencia, realice operaciones (por ejemplo, multiplicación) por separado para ambos términos. Esta técnica se describe con más detalle en el artículo Cómo agregar dos números de coma flotante sin pérdida de precisión .



En el artículo anterior, también escribí que el método$$$a\cdot(b+c)$hay una característica desagradable. Y es como sigue. Número$$$c$ siempre truncado en el último dígito significativo de un número $$$b$... Esto significa que independientemente del número$$$c$, si un $$$b+c \neq b$, entonces siempre es posible un error en el último signo, incluso para pequeños $$$c$... Esto no está permitido en el enfoque del capítulo siguiente.

Cómo funciona usando la biblioteca GNU como ejemplo

¿Cómo es? ¿Ha elegido cuál de los dos métodos descritos al principio del artículo eligió para el cálculo preciso del seno? Cualquiera que sea el método que elija, ambos son correctos. Además, son absolutamente idénticos. Créeme, compruébalo. A continuación, usaré fórmulas escolares. Son más fáciles de explicar.



Con los conocimientos adquiridos en el artículo anterior y en este artículo, puede comprender fácilmente el código de la biblioteca estándar. Abramos el archivo s_sin.c y busquemos la función __sin allí :



su código es bastante simple. Es fácil de entender que llama a un conjunto diferente de funciones dependiendo de los límites de la variable de entrada. En este artículo, discutiremos la sección de código 218-224 para ángulos 2 ^ -26 <| x | <0.855469. Puede ver que en esta sección del código, se llama a la función do_sin (x, 0). Nos detendremos en esta función con más detalle:

1. , dx=0 .
2. 129-130 , abs(x)<0.126, .. x , . , , , .
3. 136-137. , . x 2 . u x. , 0.345678. u=0.34, 0.005678.
4. 140-142. ( s ) ( c ) x . , cos(x)=1-c, 1.0, (. ), .
5. 143. u. , u=0.34 34. sin(u)=sn+ssn, cos(u)=cs+ccs. sn cs — «» u, ssn ccs — .
6. 144-145. sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. , , 144-145. — .

De hecho, he descrito solo la parte más simple de calcular el seno de esta manera. Quedan muchas matemáticas atrás. Por ejemplo, ¿cómo se calcula el tamaño de una tabla y los elementos que contiene? ¿De dónde vienen los números mágicos 0.126 y 0.855469? ¿Cuándo cortar el cálculo por el número de Taylor? Correcciones a los coeficientes de la serie de Taylor para refinar el resultado.



Todo esto, por supuesto, es interesante, pero, objetivamente, el método presentado tiene muchas desventajas: es necesario calcular el seno (s) y el coseno (c) simultáneamente, lo que requiere el doble de cálculos de la serie de Taylor ¹ . La multiplicación por valores tabulares, como podemos ver, tampoco es gratuita. Además, almacenar una tabla de 3520 bytes en RAM, por supuesto, no es un problema, pero acceder a ella (incluso en la caché) puede resultar caro.



Por lo tanto, en la siguiente parte intentaremos deshacernos de la placa y calcular el seno en el intervalo [0.126, 0.855469] directamente, pero con mayor precisión que en el primer capítulo.



Antes de terminar, una cuestión de ingenio rápido. El número grande en este ejemplo es 52776558133248 = 3 * 2 ⁴⁴ . ¿De dónde vino ese número, no, por ejemplo, 2 ⁴⁵ ? Formularé la pregunta con mayor precisión. ¿Por qué el número 3 * 2 ^{N es} óptimo al redondear números y no, por ejemplo, 2 ^{N + 1} ? Otra pregunta, ¿qué N debería elegir para redondear un número a un número entero?



¹ Cabe señalar que una ventaja significativa de este enfoque puede aparecer cuando el seno y el coseno se calculan simultáneamente desde el mismo ángulo. La segunda función se puede calcular de forma casi gratuita.