Выбор хороших первых оценок для дивизиона Гольдшмидта

Ответ 1

Я не мог удержаться от часа на вашу проблему...

Этот алгоритм описан в разделе 5.5.2 "Arithmetique des ordinateurs" Жан-Мишеля Мюллера (на французском). На самом деле это особый случай итераций Ньютона с 1 в качестве отправной точки. В книге дается простая формулировка алгоритма для вычисления N/D, причем D нормирована в диапазоне [1/2,1 [:

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

Количество правильных бит удваивается на каждой итерации. В случае 32 бит будет достаточно 4 итераций. Вы также можете выполнить итерацию до тех пор, пока e не станет слишком мал, чтобы изменить Q.

Нормализация используется, поскольку она обеспечивает максимальное количество значимых бит в результате. Также легче вычислить ошибку и количество итераций, необходимых, когда входы находятся в известном диапазоне.

Как только ваше входное значение нормализовано, вам не нужно беспокоиться о значении BASE, пока не получится обратное. Вы просто имеете 32-разрядное число X, нормированное в диапазоне от 0x80000000 до 0xFFFFFFFF, и вычисляете приближение Y = 2 ^ 64/X (Y не более 2 ^ 33).

Этот упрощенный алгоритм может быть реализован для вашего представления Q22.10 следующим образом:

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

Как отмечено в коде, умножения не заполнены 32x32- > 64 бит. E будет уменьшаться и уменьшаться и вначале помещается на 32 бита. Q всегда будет на 34 бита. Мы принимаем только 32 разряда продуктов.

Вывод 64-2*BASE-shl оставлен в качестве упражнения для читателя:-). Если он становится 0 или отрицательным, результат не представляется (входное значение слишком мало).

ИЗМЕНИТЬ. В качестве продолжения моего комментария здесь представлена ​​вторая версия с неявным 32-м битом в Q. И E, и Q теперь хранятся на 32 битах:

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
}

Ответ 2

Несколько идей для вас, хотя никто не решает вашу проблему напрямую, как указано.

• Почему этот алго для деления? Большинство разделов, которые я видел в ARM, используют некоторые переменные

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

повторяется n бит раз с двоичным поиском вне clz, чтобы определить, с чего начать. Это довольно быстро.

2. Если точность является большой проблемой, вы не ограничены 32/64 бит для вашего представления с фиксированной точкой. Это будет немного медленнее, но вы можете добавить /adc или sub/sbc, чтобы перемещать значения в регистры. mul/mla также предназначены для такого рода работ.

Опять же, не прямые ответы для вас, но, возможно, несколько идей, чтобы идти вперед. Видеть фактический код ARM, вероятно, тоже поможет мне.

Ответ 3

Безумие, вы совсем не теряете точности. Когда вы делите 512.00002f на 2 ^ 10, вы просто уменьшаете показатель вашего числа с плавающей запятой на 10. Mantissa остается прежним. Конечно, если экспонент не достигнет своего минимального значения, но этого не должно произойти, так как вы масштабируетесь до (0,5, 1).

EDIT: Хорошо, вы используете фиксированную десятичную точку. В этом случае вы должны разрешить другое представление знаменателя в вашем алгоритме. Величина D равна (0,5, 1) не только в начале, но и во всем вычислении (легко доказать, что x * (2-x) < 1 для x < 1). Таким образом, вы должны представить знаменатель с десятичной точкой в ​​базе = 32. Таким образом, вы будете иметь 32 бита точности все время.

EDIT: для этого вам придется изменить следующие строки вашего кода:

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

Кроме того, в конце вам придется сдвинуть N не по битпотам, а какое-то другое значение, которое мне слишком ленив, чтобы разобраться прямо сейчас:).