我忍不住花一个小时来解决你的问题......
Jean-Michel Muller 的“Arithmetique des ordinators”(法文)第 5.5.2 节描述了该算法。它实际上是以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或负数,则结果无法表示(输入值太小)。
编辑。作为我评论的后续内容,这是第二个版本,Q 上有隐式第 32 位。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
}