通过平方来求负指数的幂

通过平方来求负指数的幂

整数示例适用于32位

int

算术，

DWORD

为32位

unsigned int

1. 漂浮的

   pow(x,y)=x^y

通常这样评估：

* Math.Pow（等等）如何实际工作

因此可以计算分数指数：

pow(x,y) = exp2(y*log2(x))

。这也可以在固定点上完成：

* 定点大号战俘

1. 整数

   pow(a,b)=a^b

   其中

   a>=0 , b>=0

通过平方很容易（您已经拥有）：

        DWORD powuu(DWORD a,DWORD b)        {int i,bits=32;        DWORD d=1;        for (i=0;i<bits;i++) { d*=d; if (DWORd(b&0x80000000)) d*=a; b<<=1; }        return d;        }

1. 整数

   pow(a,b)=a^b

   其中

   b>=0

只需添加几个

if

s来处理负数

a

        int powiu(int a,DWORD b)     {     int sig=0,c;     if ((a<0)&&(DWORd(b&1)) { sig=1; a=-a; } // negative output only if a<0 and b is odd     c=powuu(a,b); if (sig) c=-c;     return c;     }

1. 整数

   pow(a,b)=a^b

因此，如果

b<0

这样就意味着

1/powiu(a,-b)

您可以看到结果根本不是整数，因此可以忽略这种情况或返回浮点值或添加一个乘数变量（以便您可以

PI

使用纯Integer算术对方程进行求值）。这是浮动结果：

        float powfii(int a,int b)     {     if (b<0) return 1.0/float(powiu(a,-b));     else return powiu(a,b);     }

1. 整数

   pow(a,b)=a^b

   其中

   b

   小数

您可以执行以下 *** 作

a^(1/bb)

，其中

bb

整数是。实际上，这是扎根的，因此您可以使用二进制搜索来评估：

* `a^(1/2)` 是 `square root(a)`* `a^(1/bb)` 是 `bb_root(a)`

因此，请

c

从 MSB 到 LSB 进行二进制搜索，

pow(c,bb)<=a

然后评估是否保留

bit

原样。这是

sqrt

示例：

        int bits(DWORD p) // count how many bits is p        {        DWORD m=0x80000000; int b=32;        for (;m;m>>=1,b--)         if (p>=m) break;        return b;        }    DWORD sqrt(const DWORD &x)        {        DWORD m,a;        m=(bits(x)>>1);        if (m) m=1<<m; else m=1;        for (a=0;m;m>>=1) { a|=m; if (a*a>x) a^=m; }        return a;        }

因此，现在只需更改

if (a*a>x)

with

if (pow(a,bb)>x)

，其中

bb=1/b

…

b

是您要寻找的小数指数，并且

bb

是整数。也是

m

结果的位数，因此更改

m=(bits(x)>>1);

为

m=(bits(x)/bb);

[edit1]定点sqrt示例

//---------------------------------------------------------------------------const int _fx32_fract=16;       // fractional bits countconst int _fx32_one  =1<<_fx32_fract;DWORD fx32_mul(const DWORD &x,const DWORD &y)   // unsigned fixed point mul    {    DWORD a=x,b=y;   // asm has access only to local variables    asm { // compute (a*b)>>_fx32_fract        mov eax,a    // eax=a        mov ebx,b    // ebx=b        mul eax,ebx  // (edx,eax)=eax*ebx        mov ebx,_fx32_one        div ebx      // eax=(edx,eax)>>_fx32_fract        mov a,eax;        }    return a;    }DWORD fx32_sqrt(const DWORD &x) // unsigned fixed point sqrt    {    DWORD m,a;    if (!x) return 0;    m=bits(x);       // integer bits    if (m>_fx32_fract) m-=_fx32_fract; else m=0;    m>>=1;// sqrt integer result is half of x integer bits    m=_fx32_one<<m;  // MSB of result mask    for (a=0;m;m>>=1)// test bits from MSB to 0        {        a|=m;        // bit set        if (fx32_mul(a,a)>x)    // if result is too big         a^=m;       // bit clear        }    return a;    }//---------------------------------------------------------------------------

所以这是无符号的定点。高位

16

是整数，低位

16

是小数部分。

• 这是fp-> fx转换：

  DWORd(float(x)*float(_fx32_one))

• 这是fp <-fx转换：

  float(DWORd(x))/float(_fx32_one))

• fx32_mul(x,y)

  是否

  x*y

  使用80386+ 32位体系结构的汇编程序（您可以将其重写为karatsuba或其他任何独立于平台的文件）

• fx32_sqrt(x)

  是

  sqrt(x)

在固定点上，您应该知道乘法的小数位移位：

(a<<16)*(b<<16)=(a*b<<32)

您需要向后移

>>16

以获得结果

(a*b<<16)

。结果也可能溢出

32

位，因此我

64

在汇编中使用位结果。

[edit2] 32位带符号定点Pow C ++示例

将前面的所有步骤放在一起时，应具有以下内容：

//---------------------------------------------------------------------------//--- 32bit signed fixed point format (2os complement)//---------------------------------------------------------------------------// |MSB   LSB|// |integer|.|fractional|//---------------------------------------------------------------------------const int _fx32_bits=32;          // all bits countconst int _fx32_fract_bits=16;    // fractional bits countconst int _fx32_integ_bits=_fx32_bits-_fx32_fract_bits; // integer bits count//---------------------------------------------------------------------------const int _fx32_one       =1<<_fx32_fract_bits;         // constant=1.0 (fixed point)const float _fx32_onef    =_fx32_one;        // constant=1.0 (floating point)const int _fx32_fract_mask=_fx32_one-1;      // fractional bits maskconst int _fx32_integ_mask=0xFFFFFFFF-_fx32_fract_mask; // integer bits maskconst int _fx32_sMSB_mask =1<<(_fx32_bits-1);// max signed bit maskconst int _fx32_uMSB_mask =1<<(_fx32_bits-2);// max unsigned bit mask//---------------------------------------------------------------------------float fx32_get(int   x) { return float(x)/_fx32_onef; }int   fx32_set(float x) { return int(float(x*_fx32_onef)); }//---------------------------------------------------------------------------int fx32_mul(const int &x,const int &y) // x*y    {    int a=x,b=y;     // asm has access only to local variables    asm { // compute (a*b)>>_fx32_fract        mov eax,a        mov ebx,b        mul eax,ebx  // (edx,eax)=a*b        mov ebx,_fx32_one        div ebx      // eax=(a*b)>>_fx32_fract        mov a,eax;        }    return a;    }//---------------------------------------------------------------------------int fx32_div(const int &x,const int &y) // x/y    {    int a=x,b=y;     // asm has access only to local variables    asm { // compute (a*b)>>_fx32_fract        mov eax,a        mov ebx,_fx32_one        mul eax,ebx  // (edx,eax)=a<<_fx32_fract        mov ebx,b        div ebx      // eax=(a<<_fx32_fract)/b        mov a,eax;        }    return a;    }//---------------------------------------------------------------------------int fx32_abs_sqrt(int x) // |x|^(0.5)    {    int m,a;    if (!x) return 0;    if (x<0) x=-x;    m=bits(x);       // integer bits    for (a=x,m=0;a;a>>=1,m++);  // count all bits    m-=_fx32_fract_bits;        // compute result integer bits (half of x integer bits)    if (m<0) m=0; m>>=1;    m=_fx32_one<<m;  // MSB of result mask    for (a=0;m;m>>=1)// test bits from MSB to 0        {        a|=m;        // bit set        if (fx32_mul(a,a)>x)    // if result is too big         a^=m;       // bit clear        }    return a;    }//---------------------------------------------------------------------------int fx32_pow(int x,int y)       // x^y    {    // handle special cases    if (!y) return _fx32_one;     // x^0 = 1    if (!x) return 0;  // 0^y = 0  if y!=0    if (y==-_fx32_one) return fx32_div(_fx32_one,x);    // x^-1 = 1/x    if (y==+_fx32_one) return x;  // x^+1 = x    int m,a,b,_y; int sx,sy;    // handle the signs    sx=0; if (x<0) { sx=1; x=-x; }    sy=0; if (y<0) { sy=1; y=-y; }    _y=y&_fx32_fract_mask;      // _y fractional part of exponent     y=y&_fx32_integ_mask;      //  y integer part of exponent    a=_fx32_one;     // ini result    // powering by squaring x^y    if (y)        {        for (m=_fx32_uMSB_mask;(m>_fx32_one)&&(m>y);m>>=1);     // find mask of highest bit of exponent        for (;m>=_fx32_one;m>>=1) { a=fx32_mul(a,a); if (int(y&m)) a=fx32_mul(a,x); }        }    // powering by rooting x^_y    if (_y)        {        for (b=x,m=_fx32_one>>1;m;m>>=1)      // use only fractional part { b=fx32_abs_sqrt(b); if (int(_y&m)) a=fx32_mul(a,b); }        }    // handle signs    if (sy) { if (a) a=fx32_div(_fx32_one,a); else a=0;  }     // underflow    if (sx) { if (_y) a=0;  else if(int(y&_fx32_one)) a=-a; }  // negative number ^ non integer exponent, here could add test if 1/_y is integer instead    return a;    }//---------------------------------------------------------------------------

我已经像这样测试过：

float a,b,c0,c1,d;int x,y;for (a=0.0,x=fx32_set(a);a<=10.0;a+=0.1,x=fx32_set(a)) for (b=-2.5,y=fx32_set(b);b<=2.5;b+=0.1,y=fx32_set(b))    {    if (!x) continue; // math pow has problems with this    if (!y) continue; // math pow has problems with this    c0=pow(a,b);    c1=fx32_get(fx32_pow(x,y));    d=0.0;    if (fabs(c1)<1e-3) d=c1-c0; else d=(c0/c1)-1.0;    if (fabs(d)>0.1)     d=d; // here add breakpoint to check inconsistencies with math pow    }

• a,b

  是浮点数

• x,y

  是最接近的定点表示形式

  a,b

• c0

  是数学战俘结果

• c1

  是fx32_pow结果

• d

  是差异

希望并没有忘记一些琐碎的事情，但是看起来它运作正常。不要忘记固定点的精度非常有限，因此结果会有所不同…