负指数的平方幂算法

Question

负指数的平方幂算法

6

我不确定指数为负数时幂运算的正确性。我写了下面这段代码，但仅适用于正数。

    #include <stdio.h>
    int powe(int x, int exp)
    {
         if (x == 0)
            return 1;
         if (x == 1)
            return x;
         if (x&1)
                return powe(x*x, exp/2);
         else
                return x*powe(x*x, (exp-1)/2);       
    }

查看https://zh.wikipedia.org/wiki/平方求幂算法并不能解决问题，因为以下代码似乎是错误的。

    Function exp-by-squaring(x, n ) 
      if n < 0  then return exp-by-squaring(1 / x, - n );
      else if n = 0  then return  1;
      else if n = 1  then return  x ; 
      else if n is even  then return exp-by-squaring(x * x,  n / 2);
      else if n is odd  then return x * exp-by-squaring(x * x, (n - 1) / 2).

编辑：感谢 Amit，这个解决方案适用于负数和正数：

    float powe(float x, int exp)
    {
            if (exp < 0)
                    return powe(1/x, -exp);
            if (exp == 0)
                    return 1;
            if (exp == 1)
                    return x;
            if (((int)exp)%2==0)
                    return powe(x*x, exp/2);
            else
                    return x*powe(x*x, (exp-1)/2);
    }

对于分数指数，我们可以使用以下方法（Spektre方法）：

假设你有x^0.5，那么你可以通过以下方法轻松计算平方根：从0开始到x/2并在二分查找方法中不断检查x^2是否等于结果。
所以，在情况下x^(1/3)，你必须将if mid*mid <= n替换为if mid*mid*mid <= n，然后您将得到x的立方根。对于x^(1/4)，x^(1/5)等也适用相同的方法。在x^(2/5)的情况下，我们可以进行(x^(1/5))^2，并再次减少寻找x的5次方根的问题。
但是这时候你会意识到，这种方法仅适用于当您可以将根转换为1/x格式的情况。如果我们无法转换，那么我们就卡住了吗？不，我们仍然可以继续前进，因为我们有意愿。
将您的浮点数转换为定点数，然后计算pow(a，b)。假设数字是0.6，将其转换为（24，8）浮点数得到Floor（0.6*1<<8）= 153（10011001）。正如您所知，153表示小数部分，因此在定点中，此（10011001）表示（2^-1,0,0,2^-3,2^-4,0,0,2^7）。因此，我们可以再次通过计算定点中的x的2、3、4和7根来计算pow(a，0.6)。计算后，我们还需要通过除以1<<8来以浮点格式获取结果。

上述方法的代码可在已接受的答案中找到。

还有一种基于对数的方法：

x^y = exp2(y*log2(x))

- newbie_old

1

你为什么认为这是错的？第二段代码实际上解决了负指数问题，因为 x^n = (x^-1)^(-n)。 - amit

1

第一段代码确实没有考虑到负指数的情况，所以表现出那样的行为是可以预期的。但是第二段代码是正确的，尽管你说它“看起来不对”。我的问题是：为什么你认为第二段代码看起来不对？并且我提到它能够正常工作是因为 x^n = (x^-1)^(-n)。 - amit

在你的问题中：上面那个是“第一段代码”，下面那个是“第二段代码”。 - amit

你的实现有误，请将 if (x&1) 替换为 if (exp % 2)。 - Chan Kha Vu

1

@FalconUA：是的，愚蠢的错误。我想要执行exp&1来检查奇数。 - newbie_old

显示剩余4条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Spektre · Accepted Answer

整数示例适用于32位的int算术运算，DWORD是32位的unsigned int

浮点数pow(x,y)=x^y

通常是这样计算的：

如何实际工作 Math.Pow(等等)

因此，分数指数可以计算：pow(x,y) = exp2(y*log2(x))。这也可以在定点上完成：

定点 bignum pow

整数 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;
            }

整数 pow(a,b)=a^b，其中b>=0

只需添加几个if来处理负数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;
         }

整数 pow(a,b)=a^b

如果b<0，那么意味着1/powiu(a,-b)。你可以看到结果根本不是整数，所以要么忽略这种情况，要么返回浮点值，或者添加一个乘数变量（这样你就可以在纯整数算术上评估PI方程）。这是浮点数结果：

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

整数 pow(a,b)=a^b，其中 b 是分数

您可以这样做：a^(1/bb)，其中 bb 是整数。实际上，这是求根运算，因此可以使用二分查找进行评估：

a^(1/2) 是 平方根(a)
a^(1/bb) 是 bb次方根(a)

因此，从最高位到最低位对c进行二分查找，并计算是否满足 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)更改为if (pow(a,bb)>x)，其中bb=1/b...因此b是您要查找的分数指数，而bb是整数。另外，m是结果的位数，所以将m=(bits(x)>>1);更改为m=(bits(x)/bb); [编辑1] 固定点平方根示例

//---------------------------------------------------------------------------
const int _fx32_fract=16;       // fractional bits count
const 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位为小数部分。

将fx转换为fp： 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位结果。

[编辑2] 32位有符号定点数pow C++示例

当您将所有先前的步骤结合起来时，应该得到类似于以下内容：

//---------------------------------------------------------------------------
//--- 32bit signed fixed point format (2os complement)
//---------------------------------------------------------------------------
// |MSB              LSB|
// |integer|.|fractional|
//---------------------------------------------------------------------------
const int _fx32_bits=32;                                // all bits count
const int _fx32_fract_bits=16;                          // fractional bits count
const 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 mask
const int _fx32_integ_mask=0xFFFFFFFF-_fx32_fract_mask; // integer bits mask
const int _fx32_sMSB_mask =1<<(_fx32_bits-1);           // max signed bit mask
const 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; /*Error*/ }     // underflow
    if (sx) { if (_y) a=0; /*Error*/ 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 是 math pow 的结果
c1 是 fx32_pow 的结果
d 是差异

希望没有忘记什么琐碎的东西，但看起来它似乎正常工作。请注意，定点具有非常有限的精度，因此结果会略有不同...

P.S. 看一下这个：