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

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

Этот алгоритм описан в разделе 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
}