负指数的平方

问题描述 投票：0回答：1

我不确定平方幂是否可以处理负指数。我实现了以下代码，该代码仅适用于正数。

    #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://en.wikipedia.org/wiki/Exponentiation_by_squaring 没有帮助，因为以下代码似乎是错误的。

    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).

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

    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 方法）：

  1. 假设你有x^0.5，那么你可以通过这种方法轻松计算平方根：从0开始到x/2，并不断检查x^2是否等于二分查找法的结果。

  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 次根的问题。

  3. 但是此时您会意识到此方法仅在您可以将根转换为 1/x 格式的情况下才有效。那么，如果我们无法转换，我们就会陷入困境吗？不，只要我们愿意，我们还是可以继续的。

  4. 将浮点转换为定点，然后计算 pow(a, b)。假设数字是 0.6，将其转换为 (24, 8) 浮点得到 Floor(0.6*1<<8) = 153(10011001). As you know 153 represents the fractional part so in fixed point this (10011001) represent (2^-1, 0, 0, 2^-3, 2^-4, 0, 0, 2^7).So we can again calculate the pow(a, 0.6) by calculating the 2,3, 4 and 7 root of x in fixed point. After calculating we again need to get the result in floating point by dividing with 1<<8.

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

还有一种基于日志的方法：

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

c algorithm math recursion
1个回答
6
投票

整数示例适用于 32 位 

int
算术，
DWORD
是 32 位 
unsigned int

  1. 浮动

    pow(x,y)
    =xy

    通常是这样评价的：

    因此可以计算小数指数：

    pow(x,y) = exp2(y*log2(x))
    。这也可以在定点上完成：

  2. 整数 

    pow(a,b)
    =ab 其中 
    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;
     }

  3. 整数 

    pow(a,b)
    =ab 其中 
    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;
    }

  4. 整数

    pow(a,b)
    =ab

    所以，如果 

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

  5. 整数 

    pow(a,b)
    =ab 其中 
    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)
    更改为 
    if (pow(a,bb)>x)
    ，其中 
    bb=1/b
    ...所以 
    b
    是您要查找的小数指数，
    bb
    是整数。另外 
    m
    是结果的位数，因此将 
    m=(bits(x)>>1);
    更改为 
    m=(bits(x)/bb);

[编辑1]定点 sqrt 示例

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

so this is unsigned fixed point. High `16` bits are integer and low `16` bits are fractional part.

- this is fp -> fx conversion:  `DWORD(float(x)*float(_fx32_one))`
- this is fp <- fx conversion:  `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 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
    是数学结果
  • c1
    是 fx32_pow 结果
  • d
    有区别

我希望没有忘记一些琐碎的事情，但看起来它工作正常。不要忘记定点的精度非常有限，因此结果会略有不同......

附注看看这个：

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.