Aproximación optimizada de baja precisión a `rootn (x, n)`

rootn (float_t x, int_t n) es una función que calcula eln-th raíz x1 / n y es compatible con algunos lenguajes de programación comoOpenCL. Cuando se utilizan números de punto flotante IEEE-754, aproximaciones de inicio eficientes de baja precisión para cualquiern puede generarse en base a una simple manipulación del patrón de bits subyacente, suponiendo que solo operandos normalizadosx Necesita ser procesado.

El exponente binario deroot (x, n) será 1 / n del exponente binario dex. El campo exponente de un número de coma flotante IEEE-754 está sesgado. En lugar de desviar el exponente, dividirlo y volver a sesgar el resultado, simplemente podemos dividir el exponente sesgado porn, luego aplique un desplazamiento para compensar el sesgo previamente descuidado. Además, en lugar de extraer, luego dividir, el campo exponente, simplemente podemos dividir todo el operandox, reinterpretado como un entero. El desplazamiento requerido es trivial de encontrar ya que un argumento de 1 devolverá un resultado de 1 para cualquiern.

Si tenemos dos funciones auxiliares a nuestra disposición,__int_as_float() que reinterpreta un IEEE-754binary32 comoint32y__float_as_int() que reinterpreta unint32 operando comobinary32, llegamos a la siguiente aproximación de baja precisión arootn (x, n) de manera directa:

rootn (x, n) ~= __int_as_float((int)(__float_as_int(1.0f)*(1.0-1.0/n)) + __float_as_int(x)/n)

La división entera__float_as_int (x) / n puede reducirse a un turno o multiplicación más turno poroptimizaciones bien conocidas de la división de enteros por divisor constante. Algunos ejemplos trabajados son:

rootn (x,  2) ~= __int_as_float (0x1fc00000 + __float_as_int (x) / 2)  // sqrt (x)
rootn (x,  3) ~= __int_as_float (0x2a555556 + __float_as_int (x) / 3)  // cbrt (x)
rootn (x, -1) ~= __int_as_float (0x7f000000 - __float_as_int (x) / 1)  // rcp (x)
rootn (x, -2) ~= __int_as_float (0x5f400000 - __float_as_int (x) / 2)  // rsqrt (x)
rootn (x, -3) ~= __int_as_float (0x54aaaaaa - __float_as_int (x) / 3)  // rcbrt (x)

Con todas estas aproximaciones, el resultado será exacto solo cuandox = 2Nuevo Méjico para enterom. De lo contrario, la aproximación proporcionará una sobreestimación en comparación con el verdadero resultado matemático. Podemos reducir aproximadamente a la mitaderror relativo máximo al reducir ligeramente el desplazamiento, lo que lleva a una mezcla equilibrada de subestimación y sobreestimación. Esto se logra fácilmente mediante una búsqueda binaria del desplazamiento óptimo que utiliza todos los números de coma flotante en el intervalo [1, 2n) como casos de prueba. Al hacerlo, encontramos:

rootn (x, 2) ~= __int_as_float (0x1fbb4f2e + __float_as_int(x)/2) // max rel err = 3.47474e-2
rootn (x, 3) ~= __int_as_float (0x2a51067f + __float_as_int(x)/3) // max rel err = 3.15547e-2
rootn (x,-1) ~= __int_as_float (0x7ef311c2 - __float_as_int(x)/1) // max rel err = 5.05103e-2
rootn (x,-2) ~= __int_as_float (0x5f37642f - __float_as_int(x)/2) // max rel err = 3.42128e-2
rootn (x,-3) ~= __int_as_float (0x54a232a3 - __float_as_int(x)/3) // max rel err = 3.42405e-2

Algunos pueden notar que el cálculo pararootn (x,-2) es básicamente la parte inicial deRaíz cuadrada inversa rápida del terremoto.

Basado en la observación de las diferencias entre el desplazamiento bruto original y el desplazamiento finaloptimizado para minimizar el error relativo máximo, Podría formular una heurística para la corrección secundaria y, por lo tanto, el valor de compensación final, optimizado.

Sin embargo, me pregunto si es posible determinar el desplazamiento óptimo mediante alguna fórmula de forma cerrada, de modo que el valor absoluto máximo del error relativo, max (| (aprox (x, n) - x1 / n) / X1 / n|), se minimiza para todosx en [1,2n) Para facilitar la exposición, podemos restringir abinary32 Números (IEEE-754 de precisión simple).

Soy consciente de que, en general, no existe una solución de forma cerrada para aproximaciones minimax, sin embargo, tengo la impresión de que existen soluciones de forma cerrada para el caso de aproximaciones polinómicas a funciones algebraicas comon-th raíz. En este caso tenemos una aproximación lineal (por partes).

Respuestas a la pregunta(1)

Su respuesta a la pregunta