error - Computación eficiente(a-K)/(a+K) con precisión mejorada

floating point error

Question

En varios contextos, por ejemplo para la reducción de argumentos para funciones matemáticas, se necesita calcular (a - K) / (a + K) , donde a es un argumento de variable positiva y K 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 la que se puede lograr con la división directa. Se puede asumir el soporte de hardware para la adición múltiple fusionada (FMA), ya que esta operación es proporcionada por las principales arquitecturas de CPU y GPU en este momento, y está disponible en C / C ++ a través de las funciones fma() y fmaf() .

Para facilitar la exploración, estoy experimentando con aritmética float . Ya que planeo portar el enfoque de double aritmética también, no se pueden usar operaciones que usen una precisión mayor que la nativa del argumento y el 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 los argumentos a en el intervalo [K/2, 4.23*K] , el código anterior calcula el cociente casi correctamente redondeado para todas las entradas (el error máximo está muy cerca de 0.5 ulps), siempre que K sea ​​una potencia de 2, y haya Sin desbordamiento ni subdesbordamiento en resultados intermedios. Para K no es una potencia de dos, este código es aún más preciso que el algoritmo ingenuo basado en la división. En términos de rendimiento, este código puede ser má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 cuando K = 2 ⁿ : cuando el límite superior del intervalo de trabajo aumenta a 8*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 a 0.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 más pequeño (medido en ulp vs. el 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 inferior a 0.5*K ? La eficiencia es importante, pero es probable que se puedan tolerar unas cuantas operaciones más de las que se utilizan en el código anterior.

En una de las respuestas a continuación, se señaló que podía mejorar la precisión devolviendo el cociente como una suma no evaluada de dos operandos, es decir, como un par cabeza-cola q:qlo , es decir, similar al doble float conocido y double formato. En mi código anterior, esto significaría cambiar la última línea a qlo = r * e .

Este enfoque es ciertamente útil, y ya había contemplado su uso para un logaritmo de precisión extendida para uso en pow() . Pero no ayuda fundamentalmente con la ampliación deseada del intervalo en el que el cálculo mejorado proporciona cocientes más precisos. En un caso particular que estoy observando, me gustaría usar K=2 (para precisión simple) o K=4 (para precisión doble) para mantener el intervalo de aproximación primario estrecho, y el intervalo para a 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.

Answer 1

Dado que mi objetivo es simplemente ampliar el intervalo en el que se logran resultados precisos, en lugar de encontrar una solución que funcione para todos los valores posibles de a , el uso de la aritmética de doble float para todos los cálculos intermedios parece demasiado costoso.

Pensando un poco más sobre el problema, está claro que el cálculo del resto de la división, e en el código de mi pregunta, es la parte crucial para lograr un resultado más preciso. Matemáticamente, el resto es (aK) - q * (a + K). En mi código, simplemente utilicé m para representar (aK) y representé (a + k) como m + 2*K , ya que esto brinda resultados numéricamente superiores a la representación directa.

Con un costo computacional adicional relativamente pequeño, (a + K) se puede representar como un doble float , es decir, un par cabeza-cola p:plo , que lleva a la siguiente versión modificada de mi código original:

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 2 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
mx = fmaxf (a, K);
mn = fminf (a, K);
plo = (mx - p) + mn;
t = fmaf (q, -p, m);
e = fmaf (q, -plo, t);
q = fmaf (r, e, q);

Las pruebas muestran que esto entrega resultados casi correctamente redondeados para a en [K / 2, 2 ²⁴ * K), lo que permite un aumento sustancial en el límite superior del intervalo en el que se logran resultados precisos.

La ampliación del intervalo en el extremo inferior requiere una representación más precisa de (aK). Podemos calcular esto como un par de cabeza-cola de doble float m:mlo , lo que lleva a la siguiente variante de código:

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 3 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
plo = (a < K) ? ((K - p) + a) : ((a - p) + K);
mlo = (a < K) ? (a - (K + m)) : ((a - m) - K);
t = fmaf (q, -p, m);
e = fmaf (q, -plo, t);
e = e + mlo;
q = fmaf (r, e, q);

Las pruebas exhaustivas muestran que esto proporciona resultados casi redondeados correctamente en el intervalo [K / 2 ²⁴ , K * 2 ²⁴ ). Desafortunadamente, esto tiene un costo de diez operaciones adicionales en comparación con el código de mi pregunta, que es un precio elevado que se debe pagar para obtener el error máximo de alrededor de 1.625 ulps con el cálculo ingenuo hasta cerca de 0.5 ulp.

Como en mi código original de la pregunta, uno puede expresar (a + K) en términos de (aK), eliminando así el cálculo de la cola de p , plo . Este enfoque da como resultado el siguiente código:

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 4 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
mlo = (a < K) ? (a - (K + m)) : ((a - m) - K);
t = fmaf (q, -2.0f*K, m);
t = fmaf (q, -m, t);
e = fmaf (q - 1.0f, -mlo, t);
q = fmaf (r, e, q);

Esto resulta ventajoso si el enfoque principal disminuye el límite inferior del intervalo, que es mi enfoque particular como se explica en la pregunta. Las pruebas exhaustivas del caso de precisión simple muestran que cuando K = 2 ⁿ se producen resultados casi redondeados para valores de a en el intervalo [K / 2 ²⁴ , 4.23 * K]. Con un total de 14 o 15 operaciones (dependiendo de si una arquitectura admite la predicción completa o solo movimientos condicionales), esto requiere de siete a ocho operaciones más que mi código original.

Por último, uno podría basar el cálculo residual directamente en la variable original a para evitar el error inherente en el cálculo de m y p . Esto conduce al siguiente código que, para K = 2 ⁿ , calcula resultados casi redondeados para a en el intervalo [K / 2 ²⁴ , K / 3):

/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 5 */
m = a - K;
p = a + K;
r = 1.0f / p;       
q = m * r;
t = fmaf (q + 1.0f, -K, a);
e = fmaf (q, -a, t);
q = fmaf (r, e, q);

Answer 2

El problema es la suma en (a + K) . Cualquier pérdida de precisión en (a + K) se magnifica por la división. El problema no es la división en sí misma.

Si los exponentes de a y K son iguales (casi) no se pierde precisión, y si la diferencia absoluta entre los exponentes es mayor que el tamaño significativo, entonces (a + K) == a (si a tiene una magnitud mayor) o (a + K) == K (si K tiene una magnitud mayor).

No hay manera de prevenir esto. Aumentar el tamaño significativo (por ejemplo, usar "doble extendido" de 80 bits en 80x86) solo ayuda a ampliar ligeramente el "rango de resultados precisos". Para entender por qué, considere smallest + largest (donde smallest es el número de punto flotante positivo más pequeño de 32 bits puede ser). En este caso (para flotadores de 32 bits) necesitaría un tamaño significativo de aproximadamente 260 bits para que el resultado evite la pérdida de precisión por completo. Haciendo (por ejemplo) temp = 1/(a + K); result = a * temp - K / temp; temp = 1/(a + K); result = a * temp - K / temp; tampoco ayudará mucho porque todavía tiene exactamente el mismo problema (a + K) (pero evitaría un problema similar en (a - K) ). Tampoco puedes hacer result = anything / p + anything_error/p_error porque la división no funciona así.

Solo se me ocurren 3 alternativas para acercarme a 0.5 ulps para todos los valores positivos posibles de a que puede caber en un punto flotante de 32 bits. Ninguno es probable que sea aceptable.

La primera alternativa consiste en calcular previamente una tabla de búsqueda (con cálculos matemáticos con "números reales grandes") para cada valor de a , que (con algunos trucos) termina siendo de aproximadamente 2 GiB para punto flotante de 32 bits (y completamente insano para 64- punto flotante poco). Por supuesto, si el rango de valores posibles de a es menor que "cualquier valor positivo que pueda caber en una flotación de 32 bits", el tamaño de la tabla de búsqueda se reduciría.

La segunda alternativa es usar otra cosa ("número real grande") para el cálculo en tiempo de ejecución (y convertir a / desde un punto flotante de 32 bits).

La tercera alternativa implica "algo" (no sé cómo se llama, pero es caro). Establezca el modo de redondeo en "redondear a infinito positivo" y calcule temp1 = (a + K); if(a < K) temp2 = (a - K); temp1 = (a + K); if(a < K) temp2 = (a - K); luego cambie a "redondear a infinito negativo" y calcule if(a >= K) temp2 = (a - K); lower_bound = temp2 / temp1; if(a >= K) temp2 = (a - K); lower_bound = temp2 / temp1; . Luego haga a_lower = a y disminuya a_lower en la menor cantidad posible y repita el cálculo de "lower_bound", y continúe haciendo eso hasta que obtenga un valor diferente para lower_bound , luego vuelva al valor anterior de a_lower . Después de eso, esencialmente haces lo mismo (pero los modos de redondeo opuestos, y el incremento no disminuye) para determinar upper_bound y a_upper (comenzando con el valor original de a ). Finalmente, interpolar, como a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range; a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range; . Tenga en cuenta que deseará calcular un límite superior e inferior inicial y omitir todo esto si son iguales. También tenga en cuenta que esto es todo "en teoría, completamente sin probar" y probablemente lo hice en algún lugar.

Principalmente, lo que digo es que (en mi opinión) debes rendirte y aceptar que no hay nada que puedas hacer para acercarte a 0.5 ulp. Lo siento.. :)

Answer 3

Realmente no tengo una respuesta (los análisis de error de punto flotante apropiados son muy tediosos) pero algunas observaciones:

• Las instrucciones rápidas recíprocas (como RCPSS ) no son tan precisas como la división, por lo que puede ver una reducción en la precisión si las utiliza.
• m se calcula exactamente si a ∈ [0.5 × K _b , 2 ^{1 + n} × K _b ), donde K _b es la potencia de 2 debajo de K (o K en sí si K es una potencia de 2), y n es el número de ceros finales en el significando de K (es decir, si K es una potencia de 2, entonces n = 23).
• Esto es similar a una forma simplificada del algoritmo div2 de Dekker (1971) : para expandir el rango (particularmente el límite inferior), es probable que tenga que incorporar más términos de corrección (es decir, almacenar m como la suma de 2 float s, o utilizar un double ).

Answer 4

Si a es grande en comparación con K, entonces (aK) / (a ​​+ K) = 1 - 2K / (a ​​+ K) dará una buena aproximación. Si a es pequeño en comparación con K, entonces 2a / (a ​​+ K) - 1 dará una buena aproximación. Si K / 2 ≤ a ≤ 2K, entonces aK es una operación exacta, así que hacer la división dará un resultado decente.

Answer 5

Si puede relajar la API para devolver otra variable que modela el error, entonces la solución se vuelve mucho más sencilla:

float foo(float a, float k, float *res)
{
    float ret=(a-k)/(a+k);
    *res = fmaf(-ret,a+k,a-k)/(a+k);
    return ret;
}

Esta solución solo maneja el error de corte de división, pero no maneja la pérdida de precisión de a+k y ak .

Para manejar esos errores, creo que necesito usar doble precisión, o bithack para usar un punto fijo.

El código de prueba se actualiza para generar artificialmente bits no menos significativos en la entrada

código de prueba

https://ideone.com/bHxAg8

Answer 6

Una posibilidad es rastrear el error de myp en m1 y p1 con Dekker / Schewchuk clásico:

m=a-k;
k0=a-m;
a0=k0+m;
k1=k0-k;
a1=a-a0;
m1=a1+k1;

p=a+k;
k0=p-a;
a0=p-k0;
k1=k-k0;
a1=a-a0;
p1=a1+k1;

Luego, corrige la división ingenua:

q=m/p;
r0=fmaf(p,-q,m);
r1=fmaf(p1,-q,m1);
r=r0+r1;
q1=r/p;
q=q+q1;

Eso te costará 2 divisiones, pero debería estar cerca de la mitad de la vida si no te equivoco.

Pero estas divisiones pueden reemplazarse por multiplicaciones con inverso de p sin ningún problema, ya que la primera división redondeada incorrectamente se compensará con el resto r, y la segunda división redondeada incorrectamente no importa (los últimos bits de corrección q1 no cambiarán nada ).