floating-point - factorizar - (x-1)^2 resolver

¿Por qué se prefiere(1-x)(1+x) a(1-x ^ 2)?

Question

Estaba buscando una implementación de la biblioteca en tiempo de ejecución de arcsin que se implementó calculando:

ArcTan(X, Sqrt(1 - X*X))

Sin embargo, el código que calculó 1 - X*X realmente evaluó (1-X)*(1+X) . ¿Hay una buena razón para preferir lo último? Sospecho que este último reduce la inexactitud de redondeo para X cerca de cero, pero no puedo razonar por qué sería así.

Answer 1

Escribí el programa a continuación para obtener algunos resultados empíricos para la precisión simple.

#include <float.h>
#include <math.h>
#include <stdio.h>
long double d1, d2, rel1, rel2;
float i1, i2;
int main() {
float f;
for (f = nextafterf(0, 2); f <= 1; f = nextafterf(f, 2))
{
long double o = 1.0L - ((long double)f * f);
float r1 = (1 - f) * (1 + f);
float r2 = 1 - f * f;
long double c1 = fabsl(o - r1);
long double c2 = fabsl(o - r2);
if (c1 > d1) d1 = c1;
if (c2 > d2) d2 = c2;
if (c1 / o > rel1) rel1 = c1 / o, i1 = f;
if (c2 / o > rel2) rel2 = c2 / o, i2 = f;
}
printf("(1-x)(1+x) abs:%Le relative:%Le/n", d1, rel1);
printf("1-x*x abs:%Le relative:%Le/n/n", d2, rel2);
printf("input1: %a 1-x:%a 1+x:%a (1-x)(1+x):%a o:%a/n", i1, 1-i1, 1+i1, (1-i1)*(1+i1), (double)(1 - ((long double)i1 * i1)));
printf("input2: %a x*x:%a 1-x*x:%a o:%a/n", i2, i2*i2, 1 - i2*i2, (double)(1 - ((long double)i2 * i2)));
}

Algunas observaciones:

• Usé long double 80 bits para calcular los metadatos. Esto no es suficiente para representar con precisión el error cometido para todos los valores de x , pero me temo que el programa se volverá prohibitivamente lento con mayor precisión.
• El valor de referencia o se calcula como 1.0L - ((long double)f * f) . Este es siempre el número long double más cercano al resultado real, porque (long double)f * f es exacto (ver, parece que la forma 1 - x*x veces puede ser mejor :)).

Obtuve los resultados a continuación:

(1-x)(1+x) abs:8.940394e-08 relative:9.447410e-08
1-x*x abs:4.470348e-08 relative:8.631498e-05
input1: 0x1.6a046ep-1 1-x:0x1.2bf724p-2 1+x:0x1.b50238p+0 (1-x)(1+x):0x1.0007bep-1 o:0x1.0007bc6a305ep-1
input2: 0x1.ffe96p-1 x*x:0x1.ffd2cp-1 1-x*x:0x1.6ap-12 o:0x1.69f8007p-12

De acuerdo con estos resultados, 1 - x*x tiene una precisión absoluta mejor y (1-x)*(1+x) tiene una precisión relativa mucho mejor. El punto flotante tiene que ver con la precisión relativa (todo el sistema está diseñado para permitir una representación relativamente precisa de los valores pequeños y grandes), por lo que se prefiere la última forma.

EDITAR: Computar el error final tiene más sentido, como se ilustra en la respuesta de Eric. Una subexpresión en una expresión como ArcTan(X, Sqrt(1 - X*X)) podría haber sido elegida no por su mejor precisión en general, sino porque era precisa donde más importaba. Agregar las líneas a continuación al cuerpo del bucle:

 long double a = atan2l(f, sqrtl(o));
float a1 = atan2f(f, sqrtf(r1));
float a2 = atan2f(f, sqrtf(r2));
long double e1 = fabsl(a - a1);
long double e2 = fabsl(a - a2);
if (e1 / a > ae1) ae1 = e1 / a, i1 = f;
if (e2 / a > ae2) ae2 = e2 / a, i2 = f;

Puede tener tanto sentido usar atan2l(f, sqrtf(r1)) porque no tengo la misma función de ArcTan que su sistema. De todos modos, con estas advertencias, para la expresión completa, el error relativo máximo en el intervalo [-1 ... 1] es 1.4e-07 para la versión (1-x) (1 + x) y 5.5e-7 para el 1 -x ² versión.

Answer 2

La derivada de ArcTan(X, Sqrt(1-X*X)) con respecto a X es 1/Sqrt(1-X*X) . Esto va al infinito como | X | va a 1. Por lo tanto, cuando X está cerca de 1 o -1, cualquier error en la evaluación tiene un gran efecto en el resultado. Por lo tanto, es fundamental que la evaluación minimice el error en estos casos.

Cuando X está cerca de 1, la evaluación de 1-X no tiene ningún error (en IEEE 754 o cualquier buen sistema de coma flotante, porque la escala del resultado es tal que su bit menos significativo es al menos tan bajo como el bit menos significativo en 1 o X, por lo que el resultado matemático exacto no tiene bits fuera de los bits significativos disponibles). Como 1-X es exacto, considere el efecto del error en 1+X considerando la derivada de ArcTan(X, Sqrt((1-X)*(1+X+e)) con respecto a e, donde e es el error introducido en la operación 1+X La derivada es, cuando X está cerca de 1 y e es pequeña, aproximadamente -1/10. (Tomando la derivada con Maple y sustituyendo 1 por x rinde -1/(sqrt(4+2e)*(5+2e) . Luego sustituyendo 0 por e produce -1/10.) Por lo tanto, el error en 1+X no es crítico.

Por lo tanto, evaluar la expresión como ArcTan(X, Sqrt((1-X)*(1+X)) es una buena forma de evaluarla.

La situación es simétrica para X cerca de -1. ( 1+X no tiene ningún error, y 1-X no es crítico).

Por el contrario, si consideramos el error en X*X , la derivada de ArcTan(X, Sqrt(1-X*X+e)) con respecto a e es, cuando X está cerca de 1, aproximadamente -1 / (2 * sqrt (e) * (1 + e)), entonces es grande cuando e es pequeño. Entonces, un pequeño error al evaluar X*X causará un gran error en el resultado, cuando X está cerca de 1.

Pregunte a Pascal Cuoq si, al evaluar una función f (x), generalmente estamos interesados ​​en minimizar el error relativo en el resultado final. Y, como he señalado, los errores que ocurren durante el cálculo son generalmente errores relativos en los resultados intermedios debido al redondeo en coma flotante. Pude ignorar esto en lo anterior porque estaba considerando la función cuando X está cerca de 1, así que tanto los valores intermedios considerados (1 + X y X * X) como el valor final tenían magnitudes cercanas a 1, por lo que se dividieron los valores por esas magnitudes no cambiaría nada significativamente.

Sin embargo, para completar, examiné la situación más de cerca. En Maple, escribí g := arctan(x, sqrt((1-x*x*(1+e0))*(1+e1))*(1+e2)) , permitiendo así los errores relativos e0, e1 y e2 en los cálculos de x*x , 1-x*x y sqrt , respectivamente, y escribí h:= arctan(x, sqrt((1-x)*(1+x)*(1+e0))*(1+e2)) para la alternativa. Tenga en cuenta que e0 en este caso combina los tres errores en 1-x , 1+x , y la multiplicación de ellos; el término de error completo podría ser (1+ea)*(1+eb)*(1+ec) , pero esto es efectivamente 1+e0 con un rango posible más grande para e0.

Luego examiné las derivadas de estas funciones con respecto a (una a la vez) e0, e1 y e2 dividido por abs (f (x)), donde f era la función ideal, arctan(x, sqrt(1-x*x)) . Por ejemplo, en Maple, examiné diff(g, e0) / abs(f(x)) . No realicé una evaluación analítica completa de estos; Examiné los valores de algunos valores de x cerca de 0 y cerca de 1 y para los valores de e0, e1 y e2 en uno de sus límites, -2 ^-54 .

Para x cerca de 0, los valores fueron todos de magnitud aproximadamente 1 o menos. Es decir, cualquier error relativo en el cálculo dio como resultado un error relativo similar en el resultado, o menos.

Para x cerca de 1, los valores con las derivadas de e1 y e2 eran pequeños, alrededor de 10 ^-8 o menos. Sin embargo, los valores con las derivadas de e0 fueron muy diferentes para los dos métodos. Para el método 1-x*x , el valor fue aproximadamente 2 • 10 ⁷ (usando x = 1-2 ^-53 ). Para el método (1-x)*(1+x) , el valor fue aproximadamente 5 • 10 ^-9 .

En resumen, los dos métodos no difieren mucho cerca de x = 0, pero el último método es significativamente mejor cerca de x = 1.