Cálculos precisos e rápidos para números de ponto flutuante usando o exemplo da função seno. Parte 3: ponto fixo

Continuamos o ciclo de palestras ( parte 1 e parte 2 ). Na parte 2, vimos o que está dentro da biblioteca libm e, neste trabalho, tentaremos alterar ligeiramente a função do_sin para aumentar sua precisão e velocidade. Vou citar esta função novamente do_sin ):



imagem



Conforme mostrado no artigo anterior, parte 132-145. Executado para x dentro do intervalo [0,126, 0,855469]. Bem. Vamos tentar escrever uma função que, dentro dos limites dados, seja mais precisa e, possivelmente, mais rápida.



A maneira como o usaremos é bastante óbvia. É necessário expandir a precisão dos cálculos para incluir mais casas decimais. A solução óbvia seria escolher o tipo duplo longo, contar nele e então converter de volta. Em termos de precisão, a solução deve ser boa, mas em termos de desempenho, pode haver problemas. Ainda assim, long double é um tipo de dado bastante exótico e seu suporte em processadores modernos não é uma prioridade. Em x86_64, as instruções SSE / AVX com este tipo de dados não funcionam. O coprocessador matemático será "explodido".



O que então você deve escolher? Vamos examinar mais de perto os limites de argumentos e funções.



Eles estão na região 1.0. Essa. na verdade, não precisamos de um ponto flutuante. Vamos usar um número inteiro de 64 bits ao calcular a função. Isso nos dará 10-11 bits adicionais à precisão original. Vamos descobrir como trabalhar com esses números. Um número neste formato é representado como a / d , onde a é um inteiro ed é um divisor que escolhemos constante para todas as variáveis ​​e armazenamos "em nossa memória", e não na memória do computador. Abaixo estão algumas operações para esses números:

cd=umad±bd=uma±bdcd=umadbd=umabd2cd=umadx=umaxd



Como você pode ver, não há nada de complicado nisso. A última fórmula mostra a multiplicação por qualquer inteiro. Observe também uma coisa bastante óbvia que o resultado da multiplicação de duas variáveis ​​inteiras sem sinal de tamanho N é mais frequentemente um número de tamanho até 2 * N inclusive. A adição pode causar um estouro de até 1 bit extra.



Vamos tentar escolher o divisor d . Obviamente, no mundo binário, o melhor é escolhê-lo como potência de dois, para não dividir, mas apenas movimentar o registrador, por exemplo. Que potência de dois você deve escolher? Encontre a dica nas instruções da máquina de multiplicação. Por exemplo, a instrução MUL padrão no sistema x86 multiplica 2 registradores e escreve o resultado em 2 registradores também, onde 1 dos registradores é a "parte superior" do resultado e o segundo é a parte inferior.



Por exemplo, se tivermos 2 números de 64 bits, o resultado será um número de 128 bits escrito em dois registradores de 64 bits. Vamos chamar RH - "maiúsculas" e RL - "minúsculas" 1 . Então, matematicamente, o resultado pode ser escrito comoR=RH264+Reu . Agora usamos as fórmulas acima e escrevemos a multiplicação parad=2-64

cd=uma264b264=umab2128=RH264+Reu2128=RH+Reu2-64264



E acontece que o resultado da multiplicação desses dois números de ponto fixo é o registrador R=RH .



É ainda mais fácil para o sistema Aarch64. A instrução "UMULH" multiplica dois registradores e escreve a parte "superior" da multiplicação no terceiro registrador.



Bem então. Especificamos um número de ponto fixo, mas ainda há um problema. Números negativos. Na série de Taylor, a expansão ocorre com um sinal variável. Para lidar com este problema, transformamos a fórmula de cálculo do polinômio pelo método Goner, para a seguinte forma:

pecado(x)x(1-x2(1/3!-x2(1/cinco!-x2(1/7!-x21/nove!))))



Verifique se é matematicamente exatamente igual à fórmula original. Mas em cada parêntese há uma série da forma1/(2n+1)!-x2() sempre positivo. Essa. esta conversão permite que a expressão seja avaliada como inteiros sem sinal.



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


Valores Tsx [i]
constexpr array<uint64_t, 18> tsx = { // 2^64/i!
    0x0000000000000000LL,
    0x0000000000000000LL,
    0x8000000000000000LL,
    0x2aaaaaaaaaaaaaaaLL, // Change to 0x2aaaaaaaaaaaaaafLL and check.
    0x0aaaaaaaaaaaaaaaLL,
    0x0222222222222222LL,
    0x005b05b05b05b05bLL,
    0x000d00d00d00d00dLL,
    0x0001a01a01a01a01LL,
    0x00002e3bc74aad8eLL,
    0x0000049f93edde27LL,
    0x0000006b99159fd5LL,
    0x00000008f76c77fcLL,
    0x00000000b092309dLL,
    0x000000000c9cba54LL,
    0x0000000000d73f9fLL,
    0x00000000000d73f9LL,
    0x000000000000ca96LL
};




Onde tsx[Eu]=1/Eu!em formato de ponto fixo. Desta vez, por conveniência, postei todo o código nogithubfast_sine, livrei-me do quadmath para compatibilidade com clang e arm. E mudei um pouco o método de cálculo do erro.



A comparação da função seno padrão e da função de ponto fixo é fornecida nas duas tabelas abaixo. A primeira tabela mostra a precisão do cálculo (é completamente a mesma para x86_64 e ARM). A segunda tabela é uma comparação de desempenho.



Função Número de erros Valor máximo ULP Desvio médio
sin_e7 0,0822187% 0,504787 7.10578e-20
sin_e7a 0,0560688% 0,503336 2.0985e-20
std :: sin 0,234681% 0,515376 ---




Durante o teste, o valor seno "verdadeiro" foi calculado usando a biblioteca MPFR... O valor máximo de ULP foi considerado como o desvio máximo do valor "verdadeiro". Porcentagem de erros - o número de casos em que o valor calculado da função seno por nós ou pelo binário libm não coincidiu com o valor do seno arredondado para cima. O valor médio do desvio mostra a "direção" do erro de cálculo: superestimação ou subestimação do valor. Como você pode ver na tabela, nossa função tende a superestimar os valores do seno. Isso pode ser consertado! Quem disse que os valores de tsx devem ser exatamente iguais aos coeficientes da série de Taylor. Uma ideia bastante óbvia sugere variar os valores dos coeficientes a fim de melhorar a precisão da aproximação e remover o componente constante do erro. É muito difícil fazer corretamente essa variação. Mas podemos tentar. Vamos tomar por exemplo4º valor da matriz de coeficientes tsx (tsx [3]) e altere o último número a para f. Vamos reiniciar o programa e ver a precisão (sin_e7a). Olha, é um pouco, mas aumentou! Adicionamos esse método ao nosso cofrinho.



Agora vamos ver qual é o desempenho. Para teste, peguei o que estava em mãos i5 mobile e um quarto raspberry ligeiramente overclockado (Raspberry PI 4 8GB), GCC10 da distribuição Ubuntu 20.04 x64 para ambos os sistemas.



Função x86_64 tempo, s Tempo ARM, s
sin_e7 0,174371 0,469210
std :: sin 0,154805 0,447807




Não pretendo ser mais preciso nessas medições. Variações de várias dezenas de por cento são possíveis dependendo da carga do processador. A principal conclusão pode ser feita assim. Mudar para aritmética inteira não oferece ganho de desempenho nos processadores modernos 2 . O número inimaginável de transistores nos processadores modernos torna possível realizar cálculos complexos rapidamente. Mas, eu acho que em processadores como Intel Atom, bem como em controladores fracos, esta abordagem pode dar um ganho de desempenho significativo. O que você acha?



Embora essa abordagem tenha resultado em um ganho de precisão, esse ganho de precisão parece mais interessante do que útil. Em termos de desempenho, essa abordagem pode se encontrar, por exemplo, na IoT. Mas, para computação de alto desempenho, não é mais mainstream. No mundo de hoje, SSE / AVX / CUDA prefere usar cálculo de função paralela. E na aritmética de ponto flutuante. Não existem análogos paralelos da função MUL. A função em si é mais uma homenagem à tradição.



No próximo capítulo, descreverei como você pode usar o AVX com eficácia para cálculos. Novamente, vamos entrar no código libm e tentar melhorá-lo.



1 Não há registros com esses nomes em processadores conhecidos por mim. Os nomes foram escolhidos por exemplo.

2Deve-se observar aqui que meu ARM está equipado com a versão mais recente do coprocessador matemático. Se o processador emular cálculos de ponto flutuante, os resultados podem ser dramaticamente diferentes.



All Articles