Выбор хороших первых оценок для подразделения Goldschmidt

StackOverflow https://stackoverflow.com/questions/2661541

  •  27-09-2019
  •  | 
  •  

Вопрос

Я вычисляю обратные связи с фиксированной точкой в Q22.10 с помощью Подразделение Гольдшмидта для использования в моем программном растеризаторе на ARM.

Это делается простым установлением числителя равным 1, т.е. числитель становится скаляром на первой итерации.Честно говоря, здесь я как бы слепо следую алгоритму википедии.В статье говорится, что если знаменатель масштабируется в полуоткрытом диапазоне (0,5, 1,0], хорошая первая оценка может быть основана только на знаменателе:Пусть F - предполагаемый скаляр, а D - знаменатель, тогда F = 2 - D.

Но, делая это, я сильно теряю точность.Скажите, хочу ли я найти величину, обратную 512.00002f.Чтобы уменьшить масштаб числа, я теряю 10 бит точности в дробной части, которая смещается в сторону уменьшения.Итак, мои вопросы таковы:

  • Есть ли способ выбрать лучшую оценку, которая не требует нормализации?Почему?Почему бы и нет?Математическое доказательство того, почему это возможно или не возможно, было бы отличным.
  • Кроме того, возможно ли предварительно рассчитать первые оценки, чтобы ряд сходился быстрее?Прямо сейчас он сходится в среднем после 4-й итерации.В ARM это примерно ~ 50 циклов в худшем случае, и это без учета эмуляции clz / bsr и поиска в памяти.Если это возможно, я хотел бы знать, увеличивает ли это ошибку и на сколько.

Вот мой тестовый пример.Примечание:Программная реализация clz в строке 13 - это из моего поста здесь.Вы можете заменить его встроенным, если хотите. clz должно возвращать количество начальных нулей и 32 для значения 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;
}
Это было полезно?

Решение

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

Этот алгоритм описан в разделе 5.5.2 «Арифметик-ординатуров» Джин-Мишель Мюллеру (по-французски). На самом деле это особый случай Newton Iterations с 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.

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

Как только ваше входное значение нормализуется, вам не нужно беспокоить значение базы, пока у вас не будет обратным. У вас просто 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);
}

Как отмечено в коде, умножения не заполнены 32х32-> 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
}

Другие советы

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

  1. Почему это алгос для разделения? Большинство деливок, которые я видел в руке, используйте какой-то вариант
    
      adcs hi, den, hi, lsl #1
      subcc hi, hi, den
      adcs lo, lo, lo

Повторяются N бита времени с двоичным поиском Of CLZ, чтобы определить, где начать. Это довольно чертовски быстро.

  1. Если точность является большой проблемой, вы не ограничены 32/64 битами для вашего представления фиксированного точка. Это будет немного медленнее, но вы можете сделать добавление / ADC или SUB / SBC для перемещения значений на регистры. MUL / MLA также предназначены для этой работы.

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

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

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

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

  //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 не на bitpos, а на какое-то другое значение, которое мне лень вычислять прямо сейчас :).

Лицензировано под: CC-BY-SA с атрибуция
Не связан с StackOverflow
scroll top