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.