Computando com eficiência (a - K) / (a + K) com precisão aprimorada
Em vários contextos, por exemplo, para a redução de argumentos para funções matemáticas, é necessário calcular(a - K) / (a + K)
, Ondea
é um argumento variável positivo eK
é uma constante. Em muitos casos,K
é um poder de dois, que é o caso de uso relevante para o meu trabalho. Estou procurando maneiras eficientes de calcular esse quociente com mais precisão do que é possível com a divisão direta. O suporte de hardware para FMA (Multiply Add) fundido pode ser assumido, pois esta operação é fornecida por todas as principais arquiteturas de CPU e GPU no momento e está disponível em C / C ++ através das funçõesfma()
efmaf()
.
Para facilitar a exploração, estou experimentandofloat
aritmética. Como pretendo portar a abordagem paradouble
aritmética, também, nenhuma operação usando maior que a precisão nativa do argumento e do resultado pode ser usada. Minha melhor solução até agora é:
/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 1 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
t = fmaf (q, -2.0f*K, m);
e = fmaf (q, -m, t);
q = fmaf (r, e, q);
Para argumentosa
no intervalo[K/2, 4.23*K]
, o código acima calcula o quociente arredondado quase corretamente para todas as entradas (o erro máximo é extremamente próximo de 0,5 ulps), desde queK
é uma potência de 2 e não há excedente ou insuficiente em resultados intermediários. ParaK
não um poder de dois, esse código ainda é mais preciso que o ingênuo algoritmo baseado na divisão. Em termos de desempenho, esse código pode serMais rápido do que a abordagem ingênua em plataformas nas quais o ponto recíproco de ponto flutuante pode ser calculado mais rapidamente que a divisão de ponto flutuante.
Eu faço a seguinte observação quandoK
= 2n: Quando o limite superior do intervalo de trabalho aumenta para8*K
, 16*K
, ... o erro máximo aumenta gradualmente e começa a aproximar lentamente o erro máximo da computação ingênua a partir de baixo. Infelizmente, o mesmo não parece ser verdadeiro para o limite inferior do intervalo. Se o limite inferior cair para0.25*K
, o erro máximo do método aprimorado acima é igual ao erro máximo do método ingênuo.
Existe um método para calcular q = (a - K) / (a + K) que pode obter um erro máximo menor (medido emulp vs o resultado matemático) em comparação com o método ingênuo e a sequência de códigos acima, em um intervalo maior,em particular para intervalos cujo limite inferior é inferior a0.5*K
? A eficiência é importante, mas é provável que tolerem mais algumas operações do que as usadas no código acima.
Em uma resposta abaixo, foi apontado que eu poderia melhorar a precisão retornando o quociente como uma soma não avaliada de dois operandos, ou seja, como um par de cauda-cabeçaq:qlo
, ou seja, semelhante à bem conhecida duplafloat
e duplodouble
formatos. No meu código acima, isso significaria alterar a última linha paraqlo = r * e
.
Essa abordagem é certamente útil, e eu já havia contemplado seu uso em um logaritmo de precisão estendida para uso empow()
. Mas isso não ajuda fundamentalmente com a ampliação desejada do intervalo no qual a computação aprimorada fornece quocientes mais precisos. Em um caso específico que estou olhando, gostaria de usarK=2
(para precisão única) ouK=4
(para precisão dupla) para manter o intervalo de aproximação primário estreito e o intervalo paraa
é aproximadamente [0,28]. O problema prático que estou enfrentando é que, para argumentos <0,25 * K, a precisão da divisão aprimorada não é substancialmente melhor do que com o método ingênuo.