La guía paso a paso o código que encontrarás en este artículo es la resolución más rápida y válida que encontramos a tus dudas o dilema.
Solución:
Esta pregunta pide una forma de identificar los valores de la constante Y que hacen que sea seguro transformar x / Y en un cálculo más barato utilizando FMA para todos los valores posibles de x. Otro enfoque es utilizar static análisis para determinar una sobreaproximación de los valores x puede tomar, de modo que la transformación generalmente errónea se pueda aplicar sabiendo que los valores para los cuales el código transformado difiere de la división original no suceden.
Utilizando representaciones de conjuntos de valores de punto flotante que se adaptan bien a los problemas de los cálculos de punto flotante, incluso un análisis hacia adelante que comienza desde el principio de la función puede producir información útil. Por ejemplo:
float f(float z)
float x = 1.0f + z;
float r = x / Y;
return r;
Asumiendo el modo predeterminado de redondeo al más cercano x en la función anteriorsolo puede ser NaN (si la entrada es NaN), + 0.0f, o un número mayor que 2 -24en magnitud, pero no -0,0fo nada más cercano a cero que 2-24 Y. Esto justifica la transformación en una de las dos formas mostradas en la pregunta para muchos valores de la constante
. #pragma STDC FENV_ACCESS ON
static suposición sin la cual muchas optimizaciones son imposibles y que los compiladores de C ya hacen a menos que el programa utilice explícitamente x A delanteros
- análisis que predice la información para
trueanterior puede basarse en una representación de conjuntos de valores de punto flotante que una expresión puede tomar como una tupla de:falseuna representación para los conjuntos de posibles valores de NaN (dado que los comportamientos de NaN están subespecificados, una opción es usar solo un booleano, con - lo que significa que algunos NaN pueden estar presentes, y
- indicando que no hay NaN presente.),
- cuatro banderas booleanas que indican respectivamente la presencia de + inf, -inf, +0.0, -0.0,
un intervalo inclusivo de valores de coma flotante finitos negativos, y static un intervalo inclusivo de valores de coma flotante finitos positivos. + Para seguir este enfoque, todas las operaciones de punto flotante que pueden ocurrir en un programa C deben ser entendidas por el
- analizador. Para ilustrar, la suma entre conjuntos de valores U y V, que se utilizará para manejar
- en el código analizado, se puede implementar como:
- Si NaN está presente en uno de los operandos, o si los operandos pueden ser infinitos de signos opuestos, NaN está presente en el resultado.
- Si 0 no puede ser el resultado de la suma de un valor de U y un valor de V, use aritmética de intervalo estándar. El límite superior del resultado se obtiene para la suma redondeada al más cercano del valor más grande en U y el valor más grande en V, por lo que estos límites deben calcularse con redondeo al más cercano.
- Si 0 puede ser el resultado de la suma de un valor positivo de U y un valor negativo de V, entonces sea M el valor positivo más pequeño en U tal que -M esté presente en V.
- si succ (M) está presente en U, entonces este par de valores contribuye succ (M) – M a los valores positivos del resultado.
- si -succ (M) está presente en V, entonces este par de valores contribuye con el valor negativo M – succ (M) a los valores negativos del resultado.
- si pred (M) está presente en U, entonces este par de valores contribuye con el valor negativo pred (M) – M a los valores negativos del resultado.
si -pred (M) está presente en V, entonces este par de valores contribuye con el valor M – pred (M) a los valores positivos del resultado.
Haga el mismo trabajo si 0 puede ser el resultado de la suma de un valor negativo de U y un valor positivo de V. f Reconocimiento: lo anterior toma prestadas ideas de “Mejora de las restricciones de suma y resta de punto flotante”, Bruno Marre y Claude Michel
float f(float z, float t)
float x = 1.0f + z;
if (x + t == 0.0f)
float r = x / 6.0f;
return r;
return 0.0f;
Ejemplo: compilación de la función f debajo: x El enfoque de la pregunta se niega a transformar la división en función +0.0f en una forma alternativa, porque 6 no es uno de los valores para los que la división puede transformarse incondicionalmente. En cambio, lo que estoy sugiriendo es aplicar un análisis de valor simple comenzando desde el comienzo de la función que, en este caso, determina quees un flotador finito o al menos 2 x * C2 -24
en magnitud, y utilizar esta información para aplicar la transformación de Brisebarre et al, confiando en el conocimiento de que
- no se desborda.
YPara ser explícito, sugiero usar un algoritmo como el siguiente para decidir si transformar o no la división en algo más simple: - Es
- ¿Uno de los valores que se pueden transformar usando el método de Brisebarre et al según su algoritmo?
x¿C1 y C2 de su método tienen el mismo signo, o es posible excluir la posibilidad de que el dividendo sea infinito?x¿C1 y C2 de su método tienen el mismo signo, o puedenxtomar solo una de las dos representaciones de 0? Si en el caso donde C1 y C2 tienen signos diferentes y - solo puede ser una representación de cero, recuerde jugar (**) con los signos del cálculo basado en FMA para que produzca el cero correcto cuando
x * C2es cero.
¿Se puede garantizar que la magnitud del dividendo sea lo suficientemente grande como para excluir la posibilidad de que static subdesbordamientos?
Si la respuesta a las cuatro preguntas es “sí”, entonces la división se puede transformar en una multiplicación y un FMA en el contexto de la función que se está compilando. los
El análisis descrito anteriormente sirve para responder a las preguntas 2., 3. y 4.
q = x / y
(**) “jugar con los signos” significa usar -FMA (-C1, x, (-C2) * x) en lugar de FMA (C1, x, C2 * x) cuando sea necesario para que el resultado salga correctamente cuando x solo puede ser uno de los dos ceros con signo y Déjame reiniciar por tercera vez. Estamos tratando de acelerar qdónde xes una constante entera, y y , fmaf(a,b,c) , y a * b + c son todos valores de punto flotante IEEE 754-2008 binary32. Debajo,
indica una suma múltiple fusionada
C = 1.0f / y
utilizando valores binary32.
q = x * C
El algoritmo ingenuo es a través de un recíproco precalculado,
zh = 1.0f / y
zl = -fmaf(zh, y, -1.0f) / y
para que en tiempo de ejecución sea suficiente una multiplicación (mucho más rápida):
q = fmaf(x, zh, x * zl)
La aceleración Brisebarre-Muller-Raina utiliza dos constantes precalculadas,
C1 = 1.0f / y
C2 = -y
de modo que en tiempo de ejecución, una multiplicación y una multiplicación-suma fusionada son suficientes:
t1 = x * C1
t2 = fmaf(C1, t1, x)
q = fmaf(C2, t2, t1)
El algoritmo de Markstein combina el enfoque ingenuo con dos adiciones múltiples fusionadas que producen el resultado correcto si el enfoque ingenuo arroja un resultado dentro de 1 unidad en el lugar menos significativo, mediante el cálculo previo ypara que la división se pueda aproximar usando xEl enfoque ingenuo funciona para todos los poderes de dos
, pero por lo demás es bastante malo. Por ejemplo, para los divisores 7, 14, 15, 28 y 30, produce un resultado incorrecto para más de la mitad de todos los posibles y. x El enfoque Brisebarre-Muller-Raina falla de manera similar para casi todos los que no tienen poder de dos x, pero mucho menos yproducir el resultado incorrecto (menos del medio por ciento de todos los posibles
, varía dependiendo de
). yEl artículo de Brisebarre-Muller-Raina muestra que el error máximo en el enfoque ingenuo es ± 1,5 ULP. yEl enfoque de Markstein produce resultados correctos para potencias de dos
, y también para números enteros impares
. (No he encontrado un divisor entero impar fallido para el enfoque de Markstein). x Para el enfoque de Markstein, he analizado los divisores 1 – 19700 (datos sin procesar aquí).
donde falla el enfoque de Markstein para dicho divisor), podemos ver que ocurre un patrón simple:
Casos de falla de Markstein
(fuente: nominal-animal.net)
Si cambiamos el eje x al bit inverso (dígitos binarios en orden inverso, es decir, 0b11101101 → 0b10110111, datos) de los divisores, tenemos un patrón muy claro:
Casos de falla de Markstein, divisor inverso de bits 4194304/x(fuente: nominal-animal.net)
8388608/x Si dibujamos una línea recta a través del centro de los conjuntos de puntos, obtenemos una curva 2097152/x . (Recuerde, la trama considera solo la mitad de los posibles flotadores, así que cuando considere todos los posibles flotadores, duplíquelos).
y rev(y) entre paréntesis por completo el patrón de error completo. yPor lo tanto, si usamos 8388608/rev(y) para calcular el bit inverso del divisor y, luego 16777216/rev(x) es una buena aproximación de primer orden del número de casos (de todos los posibles flotantes) donde el enfoque de Markstein produce un resultado incorrecto para un divisor par, sin potencia de dos
. (O,
function markstein_failure_estimate(divisor):
if (divisor is zero)
return no estimate
if (divisor is not an integer)
return no estimate
if (divisor is negative)
negate divisor
# Consider, for avoiding underflow cases,
if (divisor is very large, say 1e+30 or larger)
return no estimate - do as division
while (divisor > 16777216)
divisor = divisor / 2
if (divisor is a power of two)
return 0
if (divisor is odd)
return 0
while (divisor is not odd)
divisor = divisor / 2
# Use return (1 + 83833608 / divisor) / 2
# if only nonnegative finite float divisors are counted!
return 1 + 8388608 / divisor
para el límite superior.) false Agregado 2016-02-28: encontré una aproximación para el número de casos de error usando el enfoque de Markstein, dado cualquier divisor entero (binary32). Aquí está como pseudocódigo:
Esto produce una estimación de error correcta dentro de ± 1 en los casos de falla de Markstein que he probado (pero aún no he probado adecuadamente divisores mayores que 8388608). La división final debe ser tal que no informe
ceros, pero no puedo garantizarlo (todavía). No tiene en cuenta divisores muy grandes (digamos 0x1p100, o 1e + 30, y mayores en magnitud) que tienen problemas de subdesbordamiento; definitivamente excluiría tales divisores de la aceleración de todos modos.
En las pruebas preliminares, la estimación parece increíblemente precisa. No dibujé una gráfica comparando las estimaciones y los errores reales para los divisores 1 a 20000, porque todos los puntos coinciden exactamente en las gráficas. (Dentro de este rango, la estimación es exacta o demasiado grande). Esencialmente, las estimaciones reproducen exactamente la primera gráfica de esta respuesta.
El patrón de fallas del enfoque de Markstein es regular y muy interesante. El enfoque funciona para todas las potencias de dos divisores y todos los divisores enteros impares.
Para divisores mayores que 16777216, veo consistentemente los mismos errores que para un divisor que se divide por la potencia más pequeña de dos para producir un valor menor que 16777216. Por ejemplo, 0x1.3cdfa4p + 23 y 0x1.3cdfa4p + 41, 0x1. d8874p + 23 y 0x1.d8874p + 32, 0x1.cf84f8p + 23 y 0x1.cf84f8p + 34, 0x1.e4a7fp + 23 y 0x1.e4a7fp + 37. (Dentro de cada par, la mantisa es la misma y solo varía el poder de dos).
Asumiendo que mi banco de pruebas no es en error, esto significa que el enfoque de Markstein también trabaja divisores mayores que 16777216 en magnitud (pero menores que, digamos, 1e + 30), si el divisor es tal que cuando se divide por la potencia más pequeña de dos, se obtiene un cociente de menos de 16777216 de magnitud y el cociente es impar.