对Goldschmidt部门进行良好的初步估计

我正在使用Goldschmidt部门在Q22.10中计算固定点倒数,用于ARM上的软件光栅化器。

这是通过将分子设置为1来完成的,即分子成为第一次迭代的标量。 说实话,我在这里盲目地遵循维基百科算法。 文章说,如果分母在半开范围内缩放(0.5,1.0),那么良好的初步估计可以仅基于分母:设F为估计标量,D为分母,则F = 2 – D.

但是当这样做时,我会失去很多精确度。 如果我想找到512.00002f的倒数。 为了缩小数字,我在分数部分失去了10位精度,它被移出。 所以,我的问题是:

  • 有没有办法选择一个不需要标准化的更好估计? 为什么? 为什么不? 为什么这是或不可能的数学certificate将是伟大的。
  • 此外,是否可以预先计算第一个估计值,以便系列收敛得更快? 现在,它平均在第4次迭代后收敛。 在ARM上,这是大约50个周期的最坏情况,并且没有考虑到clz / bsr的仿真,也没有考虑内存查找。 如果可能的话,我想知道这样做是否会增加错误,以及增加多少错误。

这是我的测试用例。 注意:第13行的clz的软件实现来自我的post。 如果需要,可以用内在替换它。 clz应返回前导零的数量,并返回值32的32。

 #include  #include  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; else D <<= (-bitpos); /* (FNi / FDi) == (FN(i+1) / FD(i+1)) */ /* F = 2 - D */ F = (2ULL<>BASE; while(1){ DPREV = D; F = (2<>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; 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; } 

我无法忍受花一个小时解决你的问题……

该算法在Jean-Michel Muller的“Arithmetique des ordinateurs”第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的值,直到您有反向。 您只需要在范围0x80000000到0xFFFFFFFF范围内归一化的32位数X,并计算Y = 2 ^ 64 / X(Y最多为2 ^ 33)的近似值。

可以为您的Q22.10表示实现此简化算法,如下所示:

 // Fixed point inversion // EB Apr 2010 #include  #include  // 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<>(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>(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或负数,则结果不可表示(输入值太小)。

编辑。 作为我的评论的后续内容,这是第二个版本,在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 } 

有几个想法,虽然没有一个能直接解决你的问题。

  1. 为什么这个算法要分裂? 我在ARM中看到的大多数分歧使用了一些变量
     adcs hi, den, hi, lsl #1 subcc hi, hi, den adcs lo, lo, lo 

重复n位时间,通过clz的二进制搜索来确定从哪里开始。 这很快就结束了。

  1. 如果精度是一个大问题,那么您的固定点表示不限于32/64位。 它会慢一点,但你可以添加/ adc或sub / sbc来跨寄存器移动值。 mul / mla也是为这种工作而设计的。

同样,不是直接回答你的问题,但可能会有一些想法向前推进。 看到实际的ARM代码可能对我有所帮助。

Mads,你根本没有失去任何精确度。 当您将512.00002f除以2 ^ 10时,您只需将浮点数的指数减少10.尾数保持不变。 当然,除非指数达到其最小值,但这不应该发生,因为你缩放到(0.5,1)。

编辑:好的,所以你使用固定的小数点。 在这种情况下,您应该在算法中允许分母的不同表示。 D的值不仅在开始时(0.5,1)而且在整个计算过程中(很容易certificatex *(2-x)<1表示x <1)。所以你应该用十进制表示分母指向base = 32.这样,您将始终拥有32位精度。

编辑:要实现此function,您必须更改以下代码行:

  //bitpos = 31 - clz(val) - BASE; bitpos = 31 - clz(val) - 31; ... //F = (2ULL<>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,而是一些不同的值,我现在懒得弄清楚:)。