¿Hay alguna forma de conseguir un redondeo correcto con la instrucción i387 fsqrt?

Question

¿Hay alguna manera de obtener el redondeo correcto con la instrucción i387 fsqrt?...

...además de cambiar el modo de precisión en la palabra de control x87: sé que es posible, pero no es una solución razonable porque tiene problemas desagradables de tipo reentrada en los que el modo de precisión será incorrecto si se interrumpe la operación sqrt.

El problema que estoy tratando es el siguiente:el x87 fsqrt opcode realiza una operación de raíz cuadrada correctamente redondeada (según IEEE 754) en la precisión de los registros fpu, que asumiré que es precisión extendida (80 bits).Sin embargo, quiero usarlo para implementar funciones eficientes de raíz cuadrada de precisión simple y doble con los resultados correctamente redondeados (según el modo de redondeo actual).Dado que el resultado tiene un exceso de precisión, el segundo paso es convertir el resultado a redondeos de precisión simple o doble nuevamente, posiblemente dejando un resultado no redondeado correctamente.

Con algunas operaciones es posible solucionar este problema con sesgos.Por ejemplo, puedo evitar el exceso de precisión en los resultados de la suma agregando un sesgo en forma de potencia de dos que fuerce los 52 bits significativos de un valor de doble precisión a los últimos 52 bits de la mantisa de precisión extendida de 63 bits. .Pero no veo ninguna forma obvia de hacer ese truco con la raíz cuadrada.

¿Alguna idea inteligente?

(También etiquetado como C porque la aplicación prevista es la implementación de C sqrt y sqrtf funciones.)

Answer 1

Primero, dejemos de lado lo obvio:Deberías usar SSE en lugar de x87.La ESS sqrtss y sqrtsd Las instrucciones hacen exactamente lo que usted desea, son compatibles con todos los sistemas x86 modernos y también son significativamente más rápidas.

Ahora, si insistes en usar x87, comenzaré con las buenas noticias:No necesitas hacer nada para flotar.Necesitas 2p + 2 bits para calcular una raíz cuadrada correctamente redondeada en un formato de punto flotante de p bits.Porque 80 > 2*24 + 2, el redondeo adicional a precisión simple siempre redondeará correctamente y tendrá una raíz cuadrada correctamente redondeada.

Ahora las malas noticias: 80 < 2*53 + 2, así que no hay tanta suerte con la doble precisión.Puedo sugerir varias soluciones;aquí hay uno fácil y agradable que se me viene a la cabeza.

1. dejar y = round_to_double(x87_square_root(x));
2. utilizar un producto Dekker (cabeza-cola) para calcular a y b tal que y*y = a + b exactamente.
3. calcular el residual r = x - a - b.
4. if (r == 0) return y
5. if (r > 0), dejar y1 = y + 1 ulp, y calcular a1, b1 calle. y1*y1 = a1 + b1.Comparar r1 = x - a1 - b1 a r, y regresar y o y1, dependiendo de cuál tenga el residual más pequeño (o el que tenga cero bits de orden inferior, si los residuos son iguales en magnitud).
6. if (r < 0), haz lo mismo para y1 = y - 1 ulp.

Este procedimiento sólo maneja el modo de redondeo predeterminado;sin embargo, en los modos de redondeo dirigido, simplemente redondear al formato de destino es lo correcto.

Answer 2

Bien, creo que tengo una mejor solución:

1. Calcular y=sqrt(x) en precisión extendida (fsqrt).
2. Si los últimos 11 bits no son 0x400, simplemente convierta a doble precisión y regrese.
3. Agregar 0x100-(fpu_status_word&0x200) a la palabra baja de la representación de precisión extendida.
4. Convertir a doble precisión y regresar.

El paso 3 se basa en el hecho de que el bit C1 (0x200) de la palabra de estado es 1 si y sólo si fsqrtEl resultado fue redondeado.Esto es válido porque, debido a la prueba del paso 2, x no era un cuadrado perfecto;si fuera un cuadrado perfecto, y no tendría bits más allá de la doble precisión.

Puede ser más rápido realizar el paso 3 con un punto flotante condicional funcionando en lugar de trabajar en la representación de bits y recargar.

Aquí está el código (parece funcionar en todos los casos):

sqrt:
    fldl 4(%esp)
    fsqrt
    fstsw %ax
    sub $12,%esp
    fld %st(0)
    fstpt (%esp)
    mov (%esp),%ecx
    and $0x7ff,%ecx
    cmp $0x400,%ecx
    jnz 1f
    and $0x200,%eax
    sub $0x100,%eax
    sub %eax,(%esp)
    fstp %st(0)
    fldt (%esp)
1:  add $12,%esp
    fstpl 4(%esp)
    fldl 4(%esp)
    ret

Answer 3

Puede que no sea lo que quieres, ya que no aprovecha el 387. fsqrt instrucción, pero hay una sorprendentemente eficiente sqrtf(float) en glibc implementado con aritmética de enteros de 32 bits.Incluso maneja correctamente NaN, Infs y subnormales; podría ser posible eliminar algunas de estas comprobaciones con instrucciones x87 reales/indicadores de palabras de control FP.ver: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

El dbl-64/e_sqrt.c El código no es tan amigable.Es difícil decir qué suposiciones se hacen a simple vista.Curiosamente, el i386 de la biblioteca sqrt[f|l] implementaciones solo llamen fsqrt, pero carga el valor de manera diferente. flds para SP, fldl para DP.