serie - Escogiendo buenas primeras estimaciones para la división de Goldschmidt

taylor swift

Question

Estoy calculando recíprocos de punto fijo en Q22.10 con la división Goldschmidt para usar en mi software rasterizador en ARM.

Esto se hace simplemente estableciendo el numerador en 1, es decir, el numerador se convierte en el escalar en la primera iteración. Para ser honesto, estoy siguiendo ciegamente el algoritmo de wikipedia aquí. El artículo dice que si el denominador se escala en el rango semiabierto (0.5, 1.0], una buena primera estimación puede basarse solo en el denominador: Sea F el escalar estimado y D el denominador, entonces F = 2 - RE.

Pero al hacer esto, pierdo mucha precisión. Di si quiero encontrar el recíproco de 512.00002f. Para reducir el número, pierdo 10 bits de precisión en la parte de la fracción, que se desplaza. Entonces, mis preguntas son:

• ¿Hay alguna manera de elegir una mejor estimación que no requiera normalización? ¿Por qué? Por qué no? Una prueba matemática de por qué esto es o no es posible sería genial.
• Además, ¿es posible calcular previamente las primeras estimaciones para que la serie converja más rápido? En este momento, converge después de la 4ta iteración en promedio. En ARM esto es alrededor de 50 ciclos en el peor de los casos, y eso no tiene en cuenta la emulación de clz / bsr, ni las búsquedas de memoria. Si es posible, me gustaría saber si al hacerlo aumenta el error y en qué medida.

Aquí está mi testcase. Nota: la implementación del software de clz en la línea 13 es de mi publicación here . Puedes reemplazarlo con un intrínseco si lo deseas. clz debe devolver el número de ceros a la izquierda y 32 para el valor 0.

#include <stdio.h>
#include <stdint.h>
const unsigned int BASE = 22ULL;
static unsigned int divfp(unsigned int val, int* iter)
{
/* Numerator, denominator, estimate scalar and previous denominator */
unsigned long long N,D,F, DPREV;
int bitpos;
*iter = 1;
D = val;
/* Get the shift amount + is right-shift, - is left-shift. */
bitpos = 31 - clz(val) - BASE;
/* Normalize into the half-range (0.5, 1.0] */
if(0 < bitpos)
D >>= bitpos;
else
D <<= (-bitpos);
/* (FNi / FDi) == (FN(i+1) / FD(i+1)) */
/* F = 2 - D */
F = (2ULL<<BASE) - D;
/* N = F for the first iteration, because the numerator is simply 1.
So don''t waste a 64-bit UMULL on a multiply with 1 */
N = F;
D = ((unsigned long long)D*F)>>BASE;
while(1){
DPREV = D;
F = (2<<(BASE)) - D;
D = ((unsigned long long)D*F)>>BASE;
/* Bail when we get the same value for two denominators in a row.
This means that the error is too small to make any further progress. */
if(D == DPREV)
break;
N = ((unsigned long long)N*F)>>BASE;
*iter = *iter + 1;
}
if(0 < bitpos)
N >>= bitpos;
else
N <<= (-bitpos);
return N;
}
int main(int argc, char* argv[])
{
double fv, fa;
int iter;
unsigned int D, result;
sscanf(argv[1], "%lf", &fv);
D = fv*(double)(1<<BASE);
result = divfp(D, &iter);
fa = (double)result / (double)(1UL << BASE);
printf("Value: %8.8lf 1/value: %8.8lf FP value: 0x%.8X/n", fv, fa, result);
printf("iteration: %d/n",iter);
return 0;
}

Answer 1

Mads, no estás perdiendo precisión alguna. Cuando divides 512.00002f por 2 ^ 10, simplemente disminuyes el exponente de tu número de punto flotante entre 10. Mantissa permanece igual. Por supuesto, a menos que el exponente alcance su valor mínimo, pero eso no debería suceder ya que está escalando a (0.5, 1].

EDIT: Ok, entonces estás usando un punto decimal fijo. En ese caso, debe permitir una representación diferente del denominador en su algoritmo. El valor de D es de (0.5, 1] ​​no solo al principio sino a lo largo de todo el cálculo (es fácil probar que x * (2-x) <1 para x <1). Por lo tanto, debe representar el denominador con decimal punto en la base = 32. De esta manera tendrá 32 bits de precisión todo el tiempo.

EDITAR: Para implementar esto, tendrá que cambiar las siguientes líneas de su código:

 //bitpos = 31 - clz(val) - BASE;
bitpos = 31 - clz(val) - 31;
...
//F = (2ULL<<BASE) - D;
//N = F;
//D = ((unsigned long long)D*F)>>BASE;
F = -D;
N = F >> (31 - BASE);
D = ((unsigned long long)D*F)>>31;
...
//F = (2<<(BASE)) - D;
//D = ((unsigned long long)D*F)>>BASE;
F = -D;
D = ((unsigned long long)D*F)>>31;
...
//N = ((unsigned long long)N*F)>>BASE;
N = ((unsigned long long)N*F)>>31;

También al final tendrás que cambiar N, no por bitpos, sino por un valor diferente que soy demasiado perezoso para entender ahora mismo :).

Answer 2

No pude resistir pasar una hora en tu problema ...

Este algoritmo se describe en la sección 5.5.2 de "Arithmetique des ordinateurs" por Jean-Michel Muller (en francés). En realidad es un caso especial de iteraciones de Newton con 1 como punto de partida. El libro ofrece una formulación simple del algoritmo para calcular N / D, con D normalizada en el rango [1 / 2,1 [:

e = 1 - D
Q = N
repeat K times:
Q = Q * (1+e)
e = e*e

El número de bits correctos se duplica en cada iteración. En el caso de 32 bits, 4 iteraciones serán suficientes. También puede iterar hasta que e sea ​​demasiado pequeño para modificar Q

Se utiliza la normalización porque proporciona el número máximo de bits significativos en el resultado. También es más fácil calcular el error y la cantidad de iteraciones necesarias cuando las entradas están en un rango conocido.

Una vez que su valor de entrada se normaliza, no necesita molestarse con el valor de BASE hasta que tenga el inverso. Simplemente tiene un número X de 32 bits normalizado en el rango 0x80000000 a 0xFFFFFFFF, y calcula una aproximación de Y = 2 ^ 64 / X (Y es como máximo 2 ^ 33).

Este algoritmo simplificado se puede implementar para su representación Q22.10 de la siguiente manera:

// Fixed point inversion
// EB Apr 2010
#include <math.h>
#include <stdio.h>
// Number X is represented by integer I: X = I/2^BASE.
// We have (32-BASE) bits in integral part, and BASE bits in fractional part
#define BASE 22
typedef unsigned int uint32;
typedef unsigned long long int uint64;
// Convert FP to/from double (debug)
double toDouble(uint32 fp) { return fp/(double)(1<<BASE); }
uint32 toFP(double x) { return (int)floor(0.5+x*(1<<BASE)); }
// Return inverse of FP
uint32 inverse(uint32 fp)
{
if (fp == 0) return (uint32)-1; // invalid
// Shift FP to have the most significant bit set
int shl = 0; // normalization shift
uint32 nfp = fp; // normalized FP
while ( (nfp & 0x80000000) == 0 ) { nfp <<= 1; shl++; } // use "clz" instead
uint64 q = 0x100000000ULL; // 2^32
uint64 e = 0x100000000ULL - (uint64)nfp; // 2^32-NFP
int i;
for (i=0;i<4;i++) // iterate
{
// Both multiplications are actually
// 32x32 bits truncated to the 32 high bits
q += (q*e)>>(uint64)32;
e = (e*e)>>(uint64)32;
printf("Q=0x%llx E=0x%llx/n",q,e);
}
// Here, (Q/2^32) is the inverse of (NFP/2^32).
// We have 2^31<=NFP<2^32 and 2^32<Q<=2^33
return (uint32)(q>>(64-2*BASE-shl));
}
int main()
{
double x = 1.234567;
uint32 xx = toFP(x);
uint32 yy = inverse(xx);
double y = toDouble(yy);
printf("X=%f Y=%f X*Y=%f/n",x,y,x*y);
printf("XX=0x%08x YY=0x%08x XX*YY=0x%016llx/n",xx,yy,(uint64)xx*(uint64)yy);
}

Como se indica en el código, las multiplicaciones no son completas 32x32-> 64 bits. E se hará cada vez más pequeño y encaja inicialmente en 32 bits. Q siempre estará en 34 bits. Tomamos solo los 32 bits altos de los productos.

La derivación de 64-2*BASE-shl se deja como un ejercicio para el lector :-). Si se convierte en 0 o negativo, el resultado no se puede representar (el valor de entrada es demasiado pequeño).

EDITAR. Como seguimiento a mi comentario, aquí hay una segunda versión con un bit implícito 32-th en Q. Tanto E como Q ahora están almacenados en 32 bits:

uint32 inverse2(uint32 fp)
{
if (fp == 0) return (uint32)-1; // invalid
// Shift FP to have the most significant bit set
int shl = 0; // normalization shift for FP
uint32 nfp = fp; // normalized FP
while ( (nfp & 0x80000000) == 0 ) { nfp <<= 1; shl++; } // use "clz" instead
int shr = 64-2*BASE-shl; // normalization shift for Q
if (shr <= 0) return (uint32)-1; // overflow
uint64 e = 1 + (0xFFFFFFFF ^ nfp); // 2^32-NFP, max value is 2^31
uint64 q = e; // 2^32 implicit bit, and implicit first iteration
int i;
for (i=0;i<3;i++) // iterate
{
e = (e*e)>>(uint64)32;
q += e + ((q*e)>>(uint64)32);
}
return (uint32)(q>>shr) + (1<<(32-shr)); // insert implicit bit
}

Answer 3

Un par de ideas para usted, aunque ninguna que resuelva su problema directamente como se indica.

1. ¿Por qué este algo para la división? La mayoría de las divisiones que he visto en ARM usan alguna variante de

   adcs hi, den, hi, lsl #1
   subcc hi, hi, den
   adcs lo, lo, lo

repite n bits veces con una búsqueda binaria fuera del clz para determinar dónde comenzar. Eso es bastante rápido.

2. Si la precisión es un gran problema, no está limitado a 32/64 bits para su representación de punto fijo. Será un poco más lento, pero puede hacer agregar / adc o sub / sbc para mover valores a través de los registros. mul / mla también están diseñados para este tipo de trabajo.

Nuevamente, no son respuestas directas para usted, sino posiblemente algunas ideas para avanzar en esto. Ver el código ARM real probablemente también me ayudaría un poco.