🥧 💜 ✌🏼 Cálculos rápidos y precisos para números de punto flotante utilizando el ejemplo de la función seno. Introducción y parte 1 🧙🏽 👩🏽‍🏭 😢

Lea atentamente muy buenos artículos de ArtemKaravaeven la adición de números de coma flotante. El tema es muy interesante y me gustaría continuarlo y mostrar con ejemplos cómo trabajar con números de coma flotante en la práctica. Tomemos la biblioteca GNU glibc (libm) como referencia. Y para que el artículo no sea demasiado aburrido, agregaremos un componente competitivo: intentaremos no solo repetir, sino también mejorar el código de la biblioteca, haciéndolo más rápido / preciso.

He elegido la función seno trigonométrica como ejemplo. Se trata de una función generalizada cuyas matemáticas son bien conocidas en la escuela y la universidad. Al mismo tiempo, con su implementación, habrá muchos ejemplos vívidos de trabajo "correcto" con números. Usaré el doble como número de punto flotante.

En esta serie de artículos, se planea mucho, desde matemáticas hasta códigos de máquina y opciones de compilación. El lenguaje para escribir un artículo es C ++, pero sin "adornos". A diferencia del lenguaje C, los ejemplos prácticos serán más legibles incluso para personas que no estén familiarizadas con el lenguaje y ocupen menos líneas.

Los artículos se redactarán mediante el método de inmersión. Se discutirán las subtareas, que luego se unirán en una única solución al problema.

Descomposición del seno en una serie de Taylor.

La función seno se expande en una serie infinita de Taylor.

\sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - \frac{x^{7}}{7!} + \frac{x^{9}}{9!} - \dots

$\sin(x)=x-{\frac {x^{3}}{3!}}+{\frac {x^{5}}{5!}}-{\frac {x^{7}}{7!}}+{\frac {x^{9}}{9!}}-\cdots$

Está claro que no podemos calcular una serie infinita, excepto en los casos en los que existe una fórmula analítica para una suma infinita. Pero este no es nuestro caso))) Supongamos que queremos calcular el seno en el intervalo

[0, \frac{π}{2}]

$[0, \frac{\pi}{2}]$ . Discutiremos el trabajo con intervalos con más detalle en la Parte 3. Sabiendo que

\sin (\frac{π}{2}) = 1

$\sin(\frac{\pi}{2})=1$ estimación, encontramos el primer término que se puede descartar en función de la condición de que

\frac{(π / 2)^{n}}{n!} < e

$\frac{(\pi/2)^n}{n!}<e$ , donde

e

$e$ es la diferencia entre el número 1 y el número más pequeño que es mayor que 1. Hablando en términos generales, este es el último bit de la mantisa (wiki). Es más fácil resolver esta ecuación por fuerza bruta. por

e \approx 2.22 \times 10^{- 16}

$e \approx 2.22\times10^{-16}$ . Me las arreglé

n = 23

$n=23$ ya se puede eliminar. La elección correcta del número de términos se discutirá en una de las siguientes partes, por lo que por hoy "reaseguraremos" y llevaremos los términos hasta

n = 25

$n=25$ inclusive.

El último término es aproximadamente 10,000 veces menor que

e

$e$ .

La solución más sencilla

Las manos ya están picando, escribimos:

Texto completo del programa para probar

#include <iostream>
#include <iomanip>
#include <cmath>
#include <array>
#include <bitset>
#include <quadmath.h>
//      clang
//#include "/usr/lib/gcc/x86_64-linux-gnu/10/include/quadmath.h"
#include <numeric>
#include <limits>
#include <vector>

#include <boost/timer/timer.hpp>
#include <boost/math/special_functions/factorials.hpp>

namespace bm = boost::math;

using namespace std;

typedef union { uint32_t i[2]; double x; } mynumber;

array<double, 25> fc;

double sin_e1(double x) {
  double result = 0;
  int sign = 1;
  for(int i = 1; i < 25; i += 2) {
    result += sign * pow(x, i) / bm::unchecked_factorial<double>(i);
    sign = -sign;
  }
  return result;
}

double sin_e2(double x) {
  double result = 0;
  int sign = 1;
  double xx = x * x;
  double pw = x;
  double fti = 1.0;
  for(int i = 1; i < 25; i += 2) {
    fti /= i;
    result += sign * pw * fti;
    fti /= ( i + 1 );
    sign = -sign;
    pw  *= xx;
  }
  return result;
}

double sin_e3(double x) {
  double result = 0;
  for(int i = 25; i >= 1; i -= 2) {
    result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);
  }
  return result;
}

double sin_e4(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 1; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x * res;
}

double sin_e5(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 3; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x + x * xx * res;
}

inline
double fsin(double x) {
  double result;
  asm ("fsin" :"=t" (result) : "0" (x));
  return result;
}

#define SIN(a) fsin(a)
//#define SIN(a) sin(a)
//#define SIN(a) sin_e5(a)
// ^^     . ^^

int main() {

  __uint128_t ft = 1;
  fc[1] = 1.0; //3 * 5;
  for(int i = 2; i < fc.size(); i++) {
    ft *= i;
    // factorial with sign for Taylor series
    fc[i] = (1.0 / ft) * (( (i - 2) % 4 < 2) ? -1 : 1);
  }
  vector<double> xv;
  xv.resize(8 * 2000000);
  //      0  M_PI/2
  for (int i = 0; i < xv.size(); i++) {
    //      .
    xv[i] = (M_PI / 2) * i / double(xv.size());
  }

  double res = 0;
  {
    boost::timer::auto_cpu_timer at;
    for(int i = 0; i < xv.size(); i++) {
      res += SIN(xv[i]);
    }
  }

  int co = 0, cn = 0;
  //      .
  __float128 avg = 0.0, div = 0.0;
  for(int i = 0; i < xv.size(); i++) {
    mynumber dold, dnew;
    dold.x = sin(xv[i]);
    dnew.x = SIN(xv[i]);
    __float128 q = sinq(xv[i]); // <= sinq  .
    __float128 dd = __float128(dnew.x) - q;
    //     .
    div += dd * dd;
    avg += dd;
    //  ,         .
    //   ,         .
    if( dold.i[0] != dnew.i[0] || dold.i[1] != dnew.i[1] ) {
      if( fabsq(q - dold.x) <= fabsq(q - dnew.x) )
        co++;
      else
        cn++;
    }
  }
  avg /= xv.size();
  div /= xv.size();

  cout << res << endl;

  //  ,           .
  cout << "Better libm: " <<  co << " / " << xv.size() << "(" << 100.0 * co / xv.size() << "%)" << endl;

  //  ,  ""         .
  cout << "Better new: " <<  cn << " / " << xv.size() << "(" << 100.0 * cn / xv.size() << "%)" << endl;

  //         .
  cout << "  Avg / std new: " << double(avg) << " / " << double(sqrtq( div - avg * avg )) << endl;
  return 0;
}

double sin_e1(double x) {
  double result = 0;
  int sign = 1;
  for(int i = 1; i < 25; i += 2) {
    result += sign * pow(x, i) / bm::factorial<double>(i);
    sign = -sign;
  }
  return result;
}

Cómo acelerar el programa, creo que mucha gente se dio cuenta de inmediato. ¿Cuántas veces cree que sus cambios pueden acelerar el programa? La versión optimizada y la respuesta a la pregunta debajo del spoiler.

Versión optimizada del programa.

double sin_e2(double x) {
  double result = 0;
  int sign = 1;
  double xx = x * x;
  double pw = x;
  double fti = 1.0;
  for(int i = 1; i < 25; i += 2) {
    fti /= i;
    result += sign * pw * fti;
    fti /= ( i + 1 );
    sign = -sign;
    pw  *= xx;
  }
  return result;
}

Aceleración de más de 10,000 veces (GNU C ++ v10; -O2)

Mejorando la precisión

Metodología

La precisión del cálculo de la función estará determinada por 2 parámetros estándar.

La desviación estándar del verdadero pecado (float128) y la media de esta desviación. El último parámetro puede dar información importante sobre cómo se comporta nuestra función. Puede subestimar o sobrestimar sistemáticamente el resultado.

Además de estos parámetros, introduciremos dos más. Junto con nuestra función, llamamos a la función sin (doble) integrada en la biblioteca. Si los resultados de dos funciones: la nuestra y la incorporada no coinciden (bit a bit), entonces agregamos a las estadísticas cuál de las dos funciones está más lejos del valor real.

Orden de suma

Volvamos al ejemplo original. ¿Cómo puede aumentar su precisión "rápidamente"? Aquellos que leen el artículo con atención ¿Es posible sumar N números de tipo doble con mayor precisión? lo más probable es que dé una respuesta de inmediato. Es necesario girar el ciclo en sentido contrario. Para agregar desde el módulo más pequeño al más grande.

double sin_e3(double x) {
  double result = 0;
  for(int i = 25; i >= 1; i -= 2) {
    result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);
  }
  return result;
}

Los resultados se muestran en la tabla.

Función	Error medio	ETS	Mejor nuestro	Mejor libm
sin_e1	-1.28562e-18	8.25717e-17	0.0588438%	53,5466%
sin_e3	-3.4074e-21	3.39727e-17	0,0423%	10.8049%
sin_e4	8.79046e-18	4.77326e-17	0,0686%	27,6594%
sin_e5	8.78307e-18	3.69995e-17	0.0477062%	13,5105%

Puede parecer que el uso de algoritmos de suma inteligente eliminará el error casi a 0, pero este no es el caso. Por supuesto, estos algoritmos aumentarán la precisión, pero para eliminar completamente los errores, también se requieren algoritmos de multiplicación inteligentes. Existen, pero son muy caras: hay muchas operaciones innecesarias. Su uso no está justificado aquí. Sin embargo, volveremos a ellos más adelante en un contexto diferente.

Queda muy poco. Combine algoritmos rápidos y precisos. Para hacer esto, volvamos a la serie de Taylor nuevamente. Limitémoslo a 4 miembros, por ejemplo, y hagamos la siguiente transformación.

\sin (x) \approx x (1 + x^{2} (- 1 / 3! + x^{2} (1 / 5! + x^{2} (- 1 / 7! + x^{2} \cdot 1 / 9!))))

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

Puede expandir los paréntesis y comprobar que se obtiene la expresión original. Esta vista es muy fácil de encajar en un bucle.

double sin_e4(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 1; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x * res;
}

Funciona rápido, pero pierde precisión en comparación con e3. Nuevamente, el redondeo es un problema. Echemos un vistazo al último paso del ciclo y transformemos ligeramente la expresión original.

\sin (x) \approx x + x \cdot x^{2} (- 1 / 3! + \dots))

$\sin(x)\approx x+x \cdot x^2(-1/3!+ \cdots))$

Y el código correspondiente.

double sin_e5(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 3; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x + x * xx * res;
}

La precisión se ha duplicado en comparación con libm. Si puede adivinar por qué ha aumentado la precisión, escriba los comentarios. Además, hay una cosa mucho más desagradable acerca de sin_e4, que falta en sin_e5, relacionada con la precisión. Intenta adivinar cuál es el problema. En la siguiente parte, definitivamente lo hablaré en detalle.

Si te gusta el artículo, en el siguiente te diré cómo la libc GNU calcula un seno con un ULP máximo de 0.548.

Cálculos rápidos y precisos para números de punto flotante utilizando el ejemplo de la función seno. Introducción y parte 1

Descomposición del seno en una serie de Taylor.

La solución más sencilla

Mejorando la precisión

Metodología

Orden de suma

More articles: