Recogiendo buenas primeras estimaciones para la división Goldschmidt

Question

Estoy calculando FixedPoint recíprocos en Q22.10 con Goldschmidt división para su uso en mi rasterizer software en ARM.

Esto se hace simplemente establecer el numerador a 1, es decir el numerador se convierte en el escalar en la primera iteración. Para ser honesto, estoy un poco siguiendo el algoritmo de Wikipedia ciegamente aquí. El artículo dice que si el denominador se escala en el rango medio-abierto (0,5, 1,0], una buena primera estimación puede basarse en el denominador solo: Sea F el escalar estimado y D el denominador, entonces F = 2 - D.

Sin embargo, al hacer esto, pierde mucha precisión. Decir si quiero encontrar el recíproco de 512.00002f. Con el fin de ampliar el número hacia abajo, pierdo 10 bits de precisión en la parte de fracción, que se desplaza hacia fuera. Por lo tanto, mis preguntas son:

• ¿Hay una manera de recoger una mejor estimación que no requiere la normalización? ¿Por qué? Por qué no? Una prueba matemática de por qué esto es o no es posible sería grande.
• También, es posible calcular previamente las primeras estimaciones por lo que las serie converge más rápido? En este momento, converge después de la cuarta iteración en promedio. En este ARM es de unos ~ 50 ciclos peor de los casos, y que no está tomando la emulación de clz / BSR en cuenta, ni búsquedas de memoria. Si es posible, me gustaría saber si al hacerlo aumenta el error, y en qué medida.

Esta es mi caso de prueba. Nota: La aplicación de software de clz en la línea 13 es de mi puesto aquí . Se puede reemplazar con una intrínseca si lo desea. 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

No pudo resistirse a pasar una hora en su problema ...

Este algoritmo se describe en la sección 5.5.2 de "arithmétique des ordinateurs" por Jean-Michel Muller (en francés). En realidad, es un caso especial de iteraciones Newton con 1 como punto de partida. El libro da una formulación sencilla del algoritmo para calcular N / D, con D normalizó 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 duplica en cada iteración. En el caso de 32 bits, 4 iteraciones serán suficientes. También puede iterar hasta e se vuelve demasiado pequeña para modificar Q.

La normalización se utiliza, ya que proporciona el número máximo de bits significativos en el resultado. También es más fácil de calcular el error y el número de iteraciones necesarias cuando las entradas están en una gama conocida.

Una vez que se normaliza el valor de entrada, no es necesario preocuparse por el valor de la base hasta que tenga la inversa. Usted simplemente tiene un número X de 32 bits normalizado en el rango de 0x80000000 a 0xFFFFFFFF, y calcular una aproximación de Y = 2 ^ 64 / X (Y es como máximo de 2 ^ 33).

Este algoritmo simplificado puede ser implementado 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 señaló en el código, las multiplicaciones no están llenos 32x32-> 64 bits. E llegará a ser cada vez más pequeños y se ajusta inicialmente en 32 bits. Q será siempre en 34 bits. Tomamos sólo los 32 bits altos de los productos.

La derivación de 64-2*BASE-shl se deja como ejercicio para el lector :-). Si llega a ser 0 o negativo, el resultado no es representable (el valor de entrada es demasiado pequeño).

Editar. Como seguimiento a mi comentario, aquí hay una segunda versión con una implícita de 32 º poco en P. Tanto E y Q ahora se almacenan 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 2

Un par de ideas para usted, aunque ninguno que resolver su problema directamente como se ha dicho.

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

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

repetidas veces con n bits de una búsqueda binaria fuera de la clz para determinar dónde empezar. Eso es bastante Dang rápida.

2. Si la precisión es un gran problema, que no se limitan a 32/64 bits para su representación de punto fijo. Va a ser un poco más lento, pero se puede hacer de añadir / ADC o sub / SBC a los valores moverse a través de los registros. mul / mla también están diseñados para este tipo de trabajo.

Una vez más, no respuestas directas para ti, pero posiblemente algunas ideas para ir hacia adelante esto. Al ver el código ARM real sería probablemente me ayude un poco también.

Answer 3

Mads, que no están perdiendo precisión en absoluto. Cuando se divide 512.00002f por 2 ^ 10, que simplemente disminuir el exponente de su número de punto flotante por 10. Mantisa sigue siendo el mismo. Por supuesto, a menos que el exponente alcanza su valor mínimo, pero que no debería ocurrir ya que estás a escala (0,5, 1].

EDIT: Ok, así que está utilizando un punto decimal fijo. En ese caso, debería permitir una representación diferente del denominador en su algoritmo. El valor de D es de (0,5, 1] ??no sólo al principio, pero durante todo el cálculo (es fácil de demostrar que x * (2-x) <1 para x <1). Por lo que debe representar el denominador con decimales punto de base = 32. de esta manera usted tendrá 32 bits de precisión todo el tiempo.

EDIT: Para implementar este tendrá que cambiar las siguientes líneas de 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;

Además, al final vas a tener que cambiar de N no por bitpos pero algún valor diferente, que soy demasiado perezoso para averiguar en este momento:.)