En varios contextos, por ejemplo, para la reducción de argumentos para funciones matemáticas, uno necesita calcular(a - K) / (a + K), dóndea es un argumento variable positivo yK Es una constante. En muchos casos,K es una potencia de dos, que es el caso de uso relevante para mi trabajo. Estoy buscando formas eficientes de calcular este cociente con mayor precisión de lo que se puede lograr con la división directa. Se puede suponer el soporte de hardware para la fusión múltiple agregada (FMA), ya que esta operación es proporcionada por todas las arquitecturas principales de CPU y GPU en este momento, y está disponible en C / C ++ a través de las funcionesfma() yfmaf().

Para facilitar la exploración, estoy experimentando confloat aritmética. Como planeo portar el enfoque adouble aritmética también, no se pueden utilizar operaciones que utilicen una precisión mayor que la nativa tanto del argumento como del resultado. Mi mejor solución hasta ahora es:

 /* 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 en el intervalo[K/2, 4.23*K], el código anterior calcula el cociente redondeado casi correctamente para todas las entradas (el error máximo es extremadamente cercano a 0.5 ulps), siempre queK es una potencia de 2, y no hay desbordamiento o subflujo en los resultados intermedios. porK No es una potencia de dos, este código es aún más preciso que el ingenuo algoritmo basado en la división. En términos de rendimiento, este código puede serMás rápido que el enfoque ingenuo en plataformas donde el recíproco de punto flotante se puede calcular más rápido que la división de punto flotante.

Hago la siguiente observación cuandoK = 2n: Cuando el límite superior del intervalo de trabajo aumenta a8*K, 16*K, ... el error máximo aumenta gradualmente y comienza a aproximarse lentamente al error máximo del cálculo ingenuo desde abajo. Desafortunadamente, lo mismo no parece ser cierto para el límite inferior del intervalo. Si el límite inferior cae a0.25*K, el error máximo del método mejorado anterior es igual al error máximo del método ingenuo.

¿Existe un método para calcular q = (a - K) / (a + K) que pueda lograr un error máximo menor (medido enulp frente al resultado matemático) en comparación con el método ingenuo y la secuencia de código anterior, en un intervalo más amplio,en particular para intervalos cuyo límite inferior es menor que0.5*K? La eficiencia es importante, pero es probable que se toleren algunas operaciones más de las que se usan en el código anterior.

En una respuesta a continuación, se señaló que podría mejorar la precisión al devolver el cociente como una suma no evaluada de dos operandos, es decir, como un par cabeza-colaq:qlo, es decir, similar al conocido doble-float y dobledouble formatos. En mi código anterior, esto significaría cambiar la última línea aqlo = r * e.

Este enfoque es ciertamente útil, y ya había contemplado su uso para un logaritmo de precisión extendida para usar enpow(). Pero no ayuda fundamentalmente con la ampliación deseada del intervalo en el que la computación mejorada proporciona cocientes más precisos. En un caso particular que estoy viendo, me gustaría usarK=2 (para precisión simple) oK=4 (para doble precisión) para mantener el intervalo de aproximación primario estrecho y el intervalo paraa es aproximadamente [0,28]. El problema práctico al que me enfrento es que para argumentos <0.25 * K, la precisión de la división mejorada no es sustancialmente mejor que con el método ingenuo.