Escolhendo a primeira boa estimativas para Goldschmidt divisão

Question

Eu sou o cálculo fixedpoint reciprocals em Q22.10 com Goldschmidt divisão para usar no meu software quadrícula no BRAÇO.

Isso é feito por apenas definir o numerador 1, i.e o numerador torna-se o escalar na primeira iteração.Para ser honesto, eu estou seguindo as taxas algoritmo cegamente aqui.O artigo diz que se o denominador é escalado no meio-aberto (faixa de 0.5, 1.0], uma primeira estimativa pode ser baseado no denominador sozinho:Seja F ser estimado escalar e D ser o denominador, então F = 2 - D.

Mas, ao fazer isso, perde um monte de precisão.Dizer que se eu quero encontrar o recíproco da 512.00002 f.A fim de dimensionar o número para baixo, eu perco 10 bits de precisão na fração parte, que é deslocado para fora.Então, minhas perguntas são:

• Existe uma maneira de escolher a melhor estimativa que não necessitam de normalização?Por quê?Por que não?Uma prova matemática de por que isso é ou não é possível, seria ótimo.
• Também, é possível pré-calcular as primeiras estimativas para a série converge mais rapidamente?Agora, converge após a 4ª iteração, em média.No BRAÇO esta é cerca de ~50 ciclos de pior caso, e que não está tendo de emulação de van/bsr em conta, nem a memória de pesquisas.Se possível, eu gostaria de saber se isso aumenta o erro, e por quanto.

Aqui é o meu testcase.Nota:A implementação do software de clz na linha 13 é do meu post aqui.Você pode substituí-lo com um intrínseca, se você desejar. clz deve retornar o número de zeros à esquerda, e 32 para o 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

Não pude resistir a gastar uma hora no seu problema ...

Este algoritmo é descrito na Seção 5.5.2 de "arithmetique des ordinateurs", de Jean-Michel Muller (em francês). Na verdade, é um caso especial de newton iterações com 1 como ponto de partida. O livro fornece uma formulação simples do algoritmo para calcular N/D, com D normalizado em alcance [1/2,1 [:

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

O número de bits corretos dobra em cada iteração. No caso de 32 bits, 4 iterações serão suficientes. Você também pode iterar até e torna -se pequeno demais para modificar Q.

A normalização é usada porque fornece o número máximo de bits significativos no resultado. Também é mais fácil calcular o erro e o número de iterações necessárias quando as entradas estão em um intervalo conhecido.

Depois que seu valor de entrada é normalizado, você não precisa se preocupar com o valor da base até ter o inverso. Você simplesmente tem um número X de 32 bits normalizado na faixa de 0x80000000 a 0xffffffff e calcula uma aproximação de y = 2^64/x (y é no máximo 2^33).

Este algoritmo simplificado pode ser implementado para sua representação Q22.10 da seguinte maneira:

// 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);
}

Conforme observado no código, as multiplicações não estão cheias de 32x32-> 64 bits. E se tornará cada vez menor e se encaixa inicialmente em 32 bits. Q sempre estará em 34 bits. Tomamos apenas 32 bits altos dos produtos.

A derivação de 64-2*BASE-shl é deixado como um exercício para o leitor :-). Se se tornar 0 ou negativo, o resultado não será representável (o valor de entrada é muito pequeno).

EDITAR. Como acompanhamento do meu comentário, aqui está uma segunda versão com um bit implícito de 32º em Q. E e Q agora estão armazenados em 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

Algumas idéias para você, embora nenhuma que resolva seu problema diretamente, conforme declarado.

1. Por que esse algo para divisão? A maioria dos divisões que eu já vi no braço usar algum variente de

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

repetidos n bits vezes com uma pesquisa binária fora do CLZ para determinar por onde começar. Isso é muito rápido.

2. Se a precisão for um grande problema, você não está limitado a 32/64 bits para sua representação de ponto fixo. Será um pouco mais lento, mas você pode adicionar/adc ou sub/sbc para mover valores entre os registros. O MUL/MLA também foi projetado para esse tipo de trabalho.

Novamente, não respostas diretas para você, mas possivelmente algumas idéias para avançar isso. Ver o código do braço real provavelmente me ajudaria um pouco também.

Answer 3

Mads, você não perde precisão em tudo.Quando você dividir 512.00002 f por 2^10, você simplesmente diminuir o expoente de seu número em ponto flutuante por 10.Mantissa permanece o mesmo.É claro que, a menos que o expoente atinge o seu valor mínimo, mas isso não deve acontecer, já que você está de escala a (0.5, 1].

EDITAR:Ok, então você estiver usando um ponto decimal fixo.Nesse caso, você deve permitir uma representação diferente do denominador em seu algoritmo.O valor de D é de (0.5, 1], não só no início, mas ao longo de todo o cálculo (é fácil provar que x * (2-x) < 1 para x < 1).Portanto, você deve representar o denominador com o ponto decimal em base = 32.Desta forma, você terá 32 bits de precisão o tempo todo.

EDITAR:Para implementar isso, você terá que alterar as seguintes linhas 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;

Também no final, você vai ter de deslocar N não por bitpos mas algum valor diferente do que eu sou muito preguiçoso para descobrir agora :).