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

math floating-point

Question

rootn (float_t x, int_t n) es una función que calcula la raíz número n x ^{1 / n} y es compatible con algunos lenguajes de programación como OpenCL . Cuando se utilizan los números de punto flotante IEEE-754, se pueden generar aproximaciones de inicio eficientes y de baja precisión para cualquier n basándose en la simple manipulación del patrón de bits subyacente, asumiendo que solo se deben procesar los operandos normalizados x .

El exponente binario de la root (x, n) será 1 / n del exponente binario de x . El campo de exponente de un número de punto flotante IEEE-754 está sesgado. En lugar de desviar el exponente, dividirlo y volver a sesgar el resultado, simplemente podemos dividir el exponente sesgado por n , luego aplicar un desplazamiento para compensar el sesgo que anteriormente se había olvidado. Además, en lugar de extraer, luego dividir, el campo exponente, podemos simplemente dividir todo el operando x , reinterpretado como un entero. El desplazamiento requerido es trivial para encontrar, ya que un argumento de 1 devolverá un resultado de 1 para cualquier n .

Si tenemos dos funciones auxiliares a nuestra disposición, __int_as_float() que reinterpreta un __int_as_float() IEEE- int32 como int32 , y __float_as_int() que reinterpreta un operando int32 como binary32 , llegamos a la siguiente aproximación de baja precisión para rootn (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 desplazamiento o multiplicación más el desplazamiento mediante optimizaciones bien conocidas de la división entera 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 cuando x = 2 ^{n * m} para el entero m . De lo contrario, la aproximación proporcionará una sobreestimación en comparación con el verdadero resultado matemático. Podemos reducir aproximadamente la mitad del error relativo máximo reduciendo ligeramente el desplazamiento, lo que lleva a una combinación equilibrada de subestimación y sobreestimación. Esto se logra fácilmente mediante una búsqueda binaria del desplazamiento óptimo que usa todos los números de punto flotante en el intervalo [1, 2 ⁿ ) 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 para rootn (x,-2) es básicamente la porción inicial de la raíz cuadrada inversa rápida de Quake .

Basándome en observar las diferencias entre el desplazamiento original sin procesar y el desplazamiento final optimizado 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 desplazamiento 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 (| (approx (x, n) - x ^{1 / n} ) / x ^{1 / n} |), se minimiza para todas las x en [1,2 ⁿ ). Para facilitar la exposición, podemos restringirnos a números binary32 (IEEE-754 de precisión simple).

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

Answer 1

Aquí hay un código de Octava (MATLAB) que calcula las compensaciones para n positivo, asumiendo las conjeturas a continuación. Parece que funciona para 2 y 3, pero sospecho que una de las suposiciones se rompe cuando n es demasiado grande. No hay tiempo para investigar en este momento.

% finds the best offset for positive n
% assuming that the conjectures below are true
function oint = offset(n)
% find the worst upward error with no offset
efrac = (1 / log(2) - 1) * n;
evec = [floor(efrac) ceil(efrac)];
rootnvec = 2 .^ (evec / n);
[abserr i] = max((1 + evec / n) ./ rootnvec);
relerr = abserr - 1;
rootnx = rootnvec(i);
% conjecture: z such that approx(z, n) = 1
% should have the worst downward error
fun = @(o) relerr - o / rootnx + (1 / (1 + o * n) ^ (1 / n) - 1);
oreal = fzero(fun, [0 1]);
oint = round((127 * (1 - 1 / n) - oreal) * 2 ^ 23);

Respuesta parcial solo para n positivo: solo voy a mordisquear un poco al conjeturar la peor sobreestimación, asumiendo que no ajustamos el desplazamiento hacia abajo.

Definamos una versión idealizada de la aproximación para x in [1, 2^n) .

rootn-A(x, n) = 1 + floor(lg(x))/n + ((x/2^floor(lg(x))) - 1) / n
                    ^^^^^^^^^^^^^^   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                    contribution of       contribution of the
                     the exponent             significand

Queremos maximizar rootn-A(x, n) / x^(1/n) .

Parece experimentalmente que el máximo ocurre cuando x es una potencia de dos. En este caso, el término significativo es cero y floor(lg(x)) = lg(x) , por lo que podemos maximizar

(1 + lg(x)/n) / x^(1/n).

Sustituya y = lg(x)/n , y podemos maximizar (1 + y) / 2^y para y in [0, 1) modo que n*y sea ​​un número entero. Eliminando la condición de integralidad, es un ejercicio de cálculo para mostrar que el máximo de esta función cóncava está en y = 1/log(2) - 1 , aproximadamente 0.4426950408889634 . Se deduce que el máximo para x una potencia de dos está en x = 2^floor((1/log(2) - 1) * n) o x = 2^ceil((1/log(2) - 1) * n) . Conjeturo que uno de estos es, de hecho, el máximo global.

En el extremo de la subestimación, parece que queremos x modo que la salida de rootn(x, n) sea 1 , al menos para n pequeña. Más por venir más adelante con suerte.