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

Continuamos el ciclo de conferencias ( parte 1 y parte 2 ). En la parte 2, analizamos lo que hay dentro de la biblioteca libm y en este trabajo intentaremos alterar ligeramente la función do_sin para aumentar su precisión y velocidad. Citaré esta función nuevamente do_sin ):







Como se muestra en el artículo anterior, parte 132-145. Ejecutado para x dentro del rango [0,126, 0,855469]. Bien. Intentemos escribir una función que, dentro de los límites dados, sea más precisa y, posiblemente, más rápida.



La forma en que lo usaremos es bastante obvia. Es necesario ampliar la precisión de los cálculos para incluir más lugares decimales. La solución obvia sería elegir el tipo doble largo, contar en él y luego volver a convertir. En términos de precisión, la solución debería ser buena, pero en términos de rendimiento, puede haber problemas. Aún así, long double es un tipo de datos bastante exótico y su soporte en procesadores modernos no es una prioridad. En x86_64 SSE / AVX, las instrucciones con este tipo de datos no funcionan. El coprocesador matemático quedará "impresionado".



Entonces, ¿qué debería elegir? Echemos un vistazo más de cerca a los límites de los argumentos y las funciones.



Están en la región 1.0. Aquellos. de hecho, no necesitamos un punto flotante. Usemos un entero de 64 bits al calcular la función. Esto nos dará 10-11 bits adicionales a la precisión original. Averigüemos cómo trabajar con estos números. Un número en este formato se representa como a / d , donde a es un número entero yd es un divisor que elegimos constante para todas las variables y almacenamos "en nuestra memoria", y no en la memoria de la computadora. A continuación se muestran algunas operaciones para dichos números:

$$

$$\frac{c}{d}=\frac{a}{d}\pm\frac{b}{d}=\frac{a\pm b}{d}\\ \frac{c}{d}=\frac{a}{d}\cdot\frac{b}{d}=\frac{a \cdot b}{d^2}\\ \frac{c}{d}=\frac{a}{d}\cdot x=\frac{a \cdot x}{d}$$

Como puede ver, no tiene nada de complicado. La última fórmula muestra la multiplicación por cualquier número entero. Tenga en cuenta también una cosa bastante obvia que el resultado de multiplicar dos variables enteras sin signo de tamaño N es más a menudo un número de tamaño hasta 2 * N inclusive. La adición puede provocar un desbordamiento de hasta 1 bit adicional.



Intentemos elegir el divisor d . Obviamente, en el mundo binario, lo mejor es elegirlo como una potencia de dos, para no dividir, sino simplemente mover el registro, por ejemplo. ¿Qué potencia de dos deberías elegir? Encuentra la pista en las instrucciones de la máquina de multiplicar. Por ejemplo, la instrucción MUL estándar en el sistema x86 multiplica 2 registros y también escribe el resultado en 2 registros, donde 1 de los registros es la "parte superior" del resultado y el segundo es la parte inferior.



Por ejemplo, si tenemos 2 números de 64 bits, el resultado será un número de 128 bits escrito en dos registros de 64 bits. Llamemos RH - mayúscula "y RL -" minúscula " ¹ . Entonces, matemáticamente, el resultado se puede escribir como$$$R=R_H \cdot 2^{64}+R_L$... Ahora usamos las fórmulas anteriores y escribimos la multiplicación para$$$d=2^{-64}$

$$

$$\frac{c}{d}=\frac{a}{2^{64}}\cdot\frac{b}{2^{64}}=\frac{a \cdot b}{2^{128}}=\frac{R_H \cdot 2^{64} + R_L}{2^{128}}=\frac{R_H + R_L \cdot 2^{-64}}{2^{64}}$$

Y resulta que el resultado de multiplicar estos dos números de coma fija es el registro $$$R=R_H$...



Es incluso más fácil para el sistema Aarch64. La instrucción "UMULH" multiplica dos registros y escribe la parte "superior" de la multiplicación en el 3er registro.



Bien entonces. Hemos especificado un número de punto fijo, pero todavía hay un problema. Números negativos. En la serie de Taylor, la expansión va con un signo variable. Para hacer frente a este problema, transformamos la fórmula para calcular el polinomio por el método de Goner, a la siguiente forma:

$$

$$\sin(x)\approx x(1-x^2(1/3!-x^2(1/5!-x^2(1/7!-x^2\cdot1/9!))))$$

Compruebe que matemáticamente sea exactamente igual a la fórmula original. Pero en cada paréntesis hay un número como$$$1/(2n + 1)! - x^2\cdot(\cdots)$siempre positivo. Aquellos. esta conversión permite que la expresión se evalúe como enteros sin signo.

constexpr mynumber toint    = {{0x00000000, 0x43F00000}};  /*  18446744073709551616 = 2^64     */
constexpr mynumber todouble = {{0x00000000, 0x3BF00000}};  /*  ~5.42101086242752217003726400434E-20 = 2^-64     */

double sin_e7(double xd) {
  uint64_t x = xd * toint.x;
  uint64_t xx = mul2(x, x);
  uint64_t res = tsx[19]; 
  for(int i = 17; i >= 3; i -= 2) {
    res = tsx[i] - mul2(res, xx);
  }
  res = mul2(res, xx);
  res = x - mul2(x, res);
  return res * todouble.x;
}

Dónde $$$tsx[i]=1/i!$en formato de punto fijo. Esta vez, por conveniencia, publiqué todo el código en el github fast_sine , me deshice de quadmath por compatibilidad con clang y arm. Y cambié un poco el método de cálculo del error.



La comparación de la función seno estándar y la función de punto fijo se da en las dos tablas siguientes. La primera tabla muestra la precisión del cálculo (es completamente igual para x86_64 y ARM). La segunda tabla es una comparación de desempeño.







Durante la prueba, el valor del seno "verdadero" se calculó utilizando la biblioteca MPFR... El valor ULP máximo se consideró como la desviación máxima del valor "verdadero". Porcentaje de errores: el número de casos en los que el valor calculado de la función seno por nosotros o por libm no coincidió con el valor seno redondeado al doble. El valor medio de la desviación muestra la "dirección" del error de cálculo: sobreestimación o subestimación del valor. Como puede ver en la tabla, nuestra función tiende a sobreestimar los valores del seno. ¡Esto se puede arreglar! Quién dijo que los valores de tsx deberían ser exactamente iguales a los coeficientes de la serie de Taylor. Se sugiere una idea bastante obvia para variar los valores de los coeficientes para mejorar la precisión de la aproximación y eliminar el componente constante del error. Es bastante difícil hacer correctamente tal variación. Pero podemos intentar. Tomemos por ejemplo4to valor de la matriz de coeficientes tsx (tsx [3]) y cambie el último número a por f. Reiniciemos el programa y veamos la precisión (sin_e7a). ¡Mira, es un poco, pero aumentado! Agregamos este método a nuestra alcancía.



Ahora veamos cuál es el rendimiento. Para las pruebas, tomé lo que tenía a mano i5 mobile y una cuarta frambuesa ligeramente overclockeada (Raspberry PI 4 8GB), GCC10 de la distribución Ubuntu 20.04 x64 para ambos sistemas.







No pretendo ser más preciso en estas medidas. Son posibles variaciones de varias decenas de por ciento según la carga del procesador. La conclusión principal se puede hacer así. El cambio a la aritmética de números enteros no mejora el rendimiento en los procesadores modernos ² . La increíble cantidad de transistores en los procesadores modernos hace posible realizar rápidamente cálculos complejos. Pero creo que en procesadores como Intel Atom, así como en controladores débiles, este enfoque puede proporcionar una ganancia de rendimiento significativa. ¿Qué piensas?



Si bien este enfoque ha dado como resultado una ganancia en precisión, esta ganancia en precisión parece más interesante que útil. En términos de rendimiento, este enfoque puede encontrarse, por ejemplo, en el IoT. Pero para la informática de alto rendimiento, ya no es la corriente principal. En el mundo actual, SSE / AVX / CUDA prefieren utilizar el cálculo de funciones paralelas. Y en aritmética de coma flotante. No hay análogos en paralelo de la función MUL. La función en sí es más bien un tributo a la tradición.



En el próximo capítulo, describiré cómo puede usar AVX de manera efectiva para los cálculos. De nuevo, vayamos al código libm e intentemos mejorarlo.



¹ No hay registros con tales nombres en procesadores que conozco. Los nombres se han elegido, por ejemplo.

²Cabe señalar aquí que mi ARM está equipado con la última versión del coprocesador matemático. Si el procesador emulara cálculos de coma flotante, los resultados podrían ser dramáticamente diferentes.