Я не удержался, потратив час на твою проблему ...
Этот алгоритм описан в разделе 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 или отрицательным, результат не может быть представлен (входное значение слишком мало).
EDIT. В продолжение моего комментария, вот вторая версия с неявным 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
}