El algoritmo más rápido para identificar la x más pequeña y más grande que hace que la ecuación de doble precisión x + a == b sea verdadera

En el contexto del análisis estático, me interesa determinar los valores de x en la twig de entonces del condicional a continuación:

 double x; x = …; if (x + a == b) { … 

se puede suponer que a y b son constantes de doble precisión (la parte más fácil del problema es generalizar a expresiones arbitrarias), y se puede asumir que el comstackdor sigue estrictamente IEEE 754 ( FLT_EVAL_METHOD es 0). Se puede suponer que el modo de redondeo en tiempo de ejecución es el más cercano al par.

Si el cálculo con racionales fuera barato, sería simple: los valores para x serían los números de doble precisión contenidos en el intervalo racional (b – a – 0.5 * ulp1 (b) … b – a + 0.5 * ulp2 (b) ). Los límites deben incluirse si b es par, excluido si b es impar, y ulp1 y ulp2 son dos definiciones ligeramente diferentes de “ULP” que pueden tomarse idénticas si a uno no le importa perder una pequeña precisión en las potencias de dos.

Desafortunadamente, la computación con racionales puede ser costosa. Considere que otra posibilidad es obtener cada uno de los límites por dicotomía, en 64 adiciones de doble precisión (cada operación decidiendo un bit del resultado). Las 128 adiciones de punto flotante para obtener los límites inferior y superior pueden ser más rápidas que cualquier solución basada en matemáticas.

Me pregunto si hay una manera de mejorar sobre la idea de “128 adiciones de punto flotante”. En realidad, tengo mi propia solución que implica cambios de modo de redondeo y nextafter llamadas, pero no me gustaría limitar el estilo de nadie y hacer que pierdan una solución más elegante que la que tengo actualmente. Además, no estoy seguro de que cambiar el modo de redondeo dos veces sea realmente más barato que 64 adiciones de punto flotante.

Ya has dado una solución agradable y elegante en tu pregunta:

Si el cálculo con racionales fuera barato, sería simple: los valores para x serían los números de doble precisión contenidos en el intervalo racional (b – a – 0.5 * ulp1 (b) … b – a + 0.5 * ulp2 (b) ). Los límites deben incluirse si b es par, excluido si b es impar, y ulp1 y ulp2 son dos definiciones ligeramente diferentes de “ULP” que pueden tomarse idénticas si a uno no le importa perder una pequeña precisión en las potencias de dos.

Lo que sigue es un esbozo medio razonado de una solución parcial al problema basado en este párrafo. Espero tener la oportunidad de desarrollarlo pronto. Para obtener una solución real, tendrás que manejar subnormales, ceros, NaN y todas esas otras cosas divertidas. Voy a suponer que a y b son, digamos, tales que 1e-300 < |a| < 1e300 1e-300 < |a| < 1e300 y 1e-300 < |b| < 1e300 1e-300 < |b| < 1e300 para que no se produzcan locuras en ningún punto.

En ausencia de desbordamiento y subdesbordamiento, puede obtener ulp1(b) de b - nextafter(b, -1.0/0.0) . Puede obtener ulp2(b) de nextafter(b, 1.0/0.0) - b .

Si b/2 <= a <= 2b , entonces el teorema de Sterbenz le dice que b - a es exacto. Entonces (b - a) - ulp1 / 2 será el double más cercano al límite inferior y (b - a) + ulp2 / 2 será el double más cercano al límite superior. Pruebe estos valores y los valores inmediatamente antes y después, y elija el intervalo más amplio que funcione.

Si b > 2a , b - a > b/2 . El valor calculado de b - a está desactivado a lo sumo la mitad de una ulp. Una ulp1 es a lo sumo dos ulp, al igual que una ulp2 , por lo que el intervalo racional que dio es como máximo dos ulp de ancho. Averigua cuál de los cinco valores más cercanos a ba funciona.

Si a > 2b , una ulp de ba es al menos tan grande como una ulp de b ; Si algo funciona, apuesto a que tendrá que estar entre los tres valores más cercanos a ba . Imagino que el caso donde b tienen diferentes signos funciona de manera similar.

Escribí una pequeña stack de código C ++ implementando esta idea. No fallaron las pruebas aleatorias de fuzz (en diferentes rangos) antes de que me aburriera de esperar. Aquí está:

 void addeq_range(double a, double b, double &xlo, double &xhi) { if (a != a) return; // empty interval if (b != b) { if (aa != 0) { xlo = xhi = -a; return; } else return; // empty interval } if (bb != 0) { // TODO: handle me. } // b is now guaranteed to be finite. if (aa != 0) return; // empty interval if (b < 0) { addeq_range(-a, -b, xlo, xhi); xlo = -xlo; xhi = -xhi; return; } // b is now guaranteed to be zero or positive finite and a is finite. if (a >= b/2 && a <= 2*b) { double upulp = nextafter(b, 1.0/0.0) - b; double downulp = b - nextafter(b, -1.0/0.0); xlo = (ba) - downulp/2; xhi = (ba) + upulp/2; if (xlo + a == b) { xlo = nextafter(xlo, -1.0/0.0); if (xlo + a != b) xlo = nextafter(xlo, 1.0/0.0); } else xlo = nextafter(xlo, 1.0/0.0); if (xhi + a == b) { xhi = nextafter(xhi, 1.0/0.0); if (xhi + a != b) xhi = nextafter(xhi, -1.0/0.0); } else xhi = nextafter(xhi, -1.0/0.0); } else { double xmid = ba; if (xmid + a < b) { xhi = xlo = nextafter(xmid, 1.0/0.0); if (xhi + a != b) xhi = xmid; } else if (xmid + a == b) { xlo = nextafter(xmid, -1.0/0.0); xhi = nextafter(xmid, 1.0/0.0); if (xlo + a != b) xlo = xmid; if (xhi + a != b) xhi = xmid; } else { xlo = xhi = nextafter(xmid, -1.0/0.0); if (xlo + a != b) xlo = xmid; } } }