代码之家 › 专栏 › 技术社区 › Mads Elvheim

为goldschmidt部门挑选良好的初始估计

taylor-series math c

Mads Elvheim · 技术社区 · 15 年前

我正在计算q22.10中的定点倒数 Goldschmidt division 在我的软件中使用手臂上的光栅。

只需将分子设置为1,即分子在第一次迭代时成为标量即可。老实说,我有点盲目地遵循维基百科的算法。文章指出,如果分母标度在半开放范围(0.5,1.0)内,一个好的初值估计只能基于分母:假设f是估计的标量,d是分母,那么f=2-d。

但当我这样做的时候,我失去了很多精确性。假设我想找到512.00002F的倒数,为了缩小这个数字,我在小数部分失去了10位精度,小数部分被移了出来。所以,我的问题是:

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

这是我的测试用例。注:软件实现 clz 13号线在我的岗位上 here . 如果你想的话,你可以用一个内在的来代替它。 CLZ 应返回前导零的数目,值0应返回32。

#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;
}

3 回复 | 直到 15 年前

Eric Bainville 15 年前

我忍不住花了一个小时来解决你的问题…

Jean-Michel Muller(法语)在“坐标算术”第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的值操心了,除非你有相反的值。你只需要在0x8000000到0xffffffff的范围内规范化一个32位数字x,然后计算一个近似值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
}

Michael Dorgan 15 年前

给你几个主意,尽管没有一个能直接解决你的问题。

为什么要用这个算法除法?我在手臂上看到的大部分伤口都使用了


      adcs hi, den, hi, lsl #1
      subcc hi, hi, den
      adcs lo, lo, lo

重复n位次,从clz开始进行二进制搜索,以确定从何处开始。太快了。

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

再说一次,不是直接给你答案,但可能有一些想法来推进这个。看到实际的arm代码可能对我也有帮助。

michalburger1 15 年前

玛兹,你一点也不失精准。当你用512.00002F除以2^10时,你只需将浮点数的指数减少10。尾数不变。当然,除非指数达到其最小值,但这不应该发生,因为你正在缩放到(0.5,1)。

编辑:好的,所以你用的是固定的小数点。在这种情况下,应该允许算法中分母的不同表示。d的值不仅从(0.5,1]开始,而且在整个计算过程中(很容易证明x*(2-x)<1代表x<1)。所以你应该用小数点表示分母,基数是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,而不是按位,而是一些不同的值,我现在懒得弄明白:)。