Cálculos precisos e rápidos para números de ponto flutuante usando o exemplo da função seno. Introdução e parte 1

Leia com atenção artigos muito bons de ArtemKaravaevna adição de números de ponto flutuante. O tópico é muito interessante e eu gostaria de continuar e mostrar com exemplos como trabalhar com números de ponto flutuante na prática. Vamos tomar a biblioteca GNU glibc (libm) como referência. E para que o artigo não fique muito enfadonho, adicionaremos um componente competitivo: tentaremos não só repetir, mas também melhorar o código da biblioteca, tornando-o mais rápido / preciso.



Escolhi a função seno trigonométrica como exemplo. Esta é uma função amplamente difundida, cuja matemática é bem conhecida na escola e na universidade. Ao mesmo tempo, quando for implementado, haverá muitos exemplos vívidos de trabalho "correto" com números. Vou usar double como um número de ponto flutuante.



Nesta série de artigos, muito é planejado, desde matemática até códigos de máquina e opções de compilador. A linguagem de escrever um artigo é C ++, mas sem "enfeites". Em contraste com a linguagem C, os exemplos de trabalho serão mais legíveis mesmo para pessoas não familiarizadas com a linguagem e ocuparão menos linhas.



Os artigos serão escritos pelo método de imersão. Serão discutidas subtarefas, que então se reunirão em uma única solução para o problema.



Decomposição do seno em uma série de Taylor.



A função seno se expande em uma série infinita de Taylor.



sin(x)=xx33!+x55!x77!+x99!



É claro que não podemos calcular uma série infinita, exceto nos casos em que existe uma fórmula analítica para uma soma infinita. Mas este não é o nosso caso))) Suponha que queremos calcular o seno no intervalo[0,π2]... Discutiremos o trabalho com intervalos com mais detalhes na Parte 3. Sabendo quesin(π2)=1 estimativa, encontre o primeiro termo que pode ser descartado com base na condição de (π/2)nn!<eOnde eé a diferença entre o número 1 e o menor número que é maior do que 1. Grosso modo, este é o último bit da mantissa ( wiki ). É mais fácil resolver essa equação pela força bruta. Parae2.22×1016... eu conseguin=23já pode ser descartado. A escolha correta do número de termos será discutida em uma das próximas partes, portanto, por hoje iremos "ressegurar" e levar os termos atén=25inclusive.

O último termo é aproximadamente 10.000 vezes menor quee...



Solução mais simples



As mãos já estão coçando, escrevemos:



Texto completo do programa para teste
#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;
}


Como acelerar o programa, acho que muitas pessoas descobriram na hora. Quantas vezes você acha que suas mudanças podem acelerar o programa? A versão otimizada e a resposta à pergunta sob o spoiler.



Versão otimizada do 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;
}


10000 (GNU C++ v10; -O2)



Melhorando a precisão



Metodologia



A precisão do cálculo da função será determinada por 2 parâmetros padrão.



O desvio padrão do pecado verdadeiro (float128) e a média desse desvio. O último parâmetro pode fornecer informações importantes sobre como nossa função se comporta. Ela pode subestimar ou superestimar sistematicamente o resultado.



Além desses parâmetros, apresentaremos mais dois. Junto com nossa função, chamamos a função sin (double) incorporada à biblioteca. Se os resultados de duas funções: nossa e a integrada não corresponderem (bit a bit), então adicionamos às estatísticas qual das duas funções está mais longe do valor verdadeiro.



Ordem de soma



Vamos voltar ao exemplo original. Como você pode aumentar sua precisão "rapidamente"? Aqueles que leram o artigo com atenção É possível adicionar N números de tipo duplo com mais precisão? provavelmente dará uma resposta imediata. É necessário virar o ciclo na direção oposta. Para adicionar do menor módulo ao maior.



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;
}


Os resultados são apresentados na tabela.



Função Erro médio STD O nosso é melhor Melhor 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%


Pode parecer que usar os algoritmos de soma inteligente removerá o erro para quase 0, mas não é o caso. É claro que esses algoritmos irão aumentar a precisão, mas para se livrar completamente dos erros, algoritmos de multiplicação inteligentes também são necessários. Eles existem, mas são muito caros: há muitas operações desnecessárias. Seu uso não se justifica aqui. No entanto, voltaremos a eles mais tarde em um contexto diferente.



Resta muito pouco. Combine algoritmos rápidos e precisos. Para fazer isso, vamos voltar à série Taylor novamente. Vamos restringi-lo a 4 membros, por exemplo, e fazer a seguinte transformação.



sin(x)x(1+x2(1/3!+x2(1/5!+x2(1/7!+x21/9!))))





Você pode expandir os parênteses e verificar se a expressão original foi obtida. Essa visualização é muito fácil de encaixar em um loop.



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, mas perdeu precisão em comparação com e3. Novamente, o arredondamento é um problema. Vamos dar uma olhada na última etapa do loop e transformar um pouco a expressão original.

sin(x)x+xx2(1/3!+))





E o código correspondente.



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;
}


A precisão dobrou em comparação com libm. Se você consegue adivinhar por que a precisão aumentou, escreva nos comentários. Além disso, há mais uma coisa muito mais desagradável sobre sin_e4, que está faltando em sin_e5, relacionada à precisão. Tente adivinhar qual é o problema. Na próxima parte, com certeza falarei sobre isso em detalhes.



Se você gostou do artigo, no próximo irei lhe contar como a GNU libc calcula um seno com um ULP máximo de 0,548.



All Articles