平方负指数的权力

我不确定平方的权力是否负指数。 我实现了以下只适用于正数的代码。

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

编辑:感谢阿米特这个解决scheme适用于负数和正数:

  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),你必须replaceif mid*mid <= nif mid*mid*mid <= n ,那么你将得到x的立方根。同样的, (1/4),x ^(1/5)等等。 在x ^(2/5)的情况下,我们可以做(x ^(1/5))^ 2,再次减lessfindx的第5个根的问题。

  3. 然而到了这个时候,你会意识到这种方法只适用于可以将根转换为1 / x格式的情况。 所以我们卡住了,如果我们不能转换? 不,我们仍然可以按照我们的意愿继续前进。

  4. 将你的浮点数转换成固定点,然后计算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来得到浮点结果。

以上方法的代码可以在接受的答案中find。

还有一个基于日志的方法

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

整数例子是32位int运算, DWORD是32位unsigned int

  1. pow(x,y)=x^y

    通常评估如下:

    • Math.Pow(等等)实际上是如何工作的

    所以可以评估分数指数: pow(x,y) = exp2(y*log2(x)) 。 这也可以在固定点上完成:

    • 定点高粱粉
  2. 整数幂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; } 
  3. 整数pow(a,b)=a^b其中b>=0

    只要加上less数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)=a^b

    所以如果b<0那么它意味着1/powiu(a,-b)正如你所看到的,结果根本不是整数,所以要么忽略这种情况,要么返回浮点值或者添加一个乘数variables(这样就可以计算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)=a^b其中b是分数

    你可以做这样的事情a^(1/bb)其中bb是整数。 事实上,这是生根,所以你可以使用二进制search来评估:

    • a^(1/2)square root(a)
    • a^(1/bb)bb_root(a)

    所以做二进制searchcMSBLSB,并评估如果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);

[edit1]定点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; } //--------------------------------------------------------------------------- 

所以这是无符号的定点。 高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位结果。

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

我已经testing过,像这样:

 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,ya,b最接近的不动点表示
  • c0是math战果的结果
  • c1是fx32_pow的结果
  • d是不同的

    希望没有忘记一些微不足道的事情,但似乎是正常的。 不要忘记,定点精度非常有限,所以结果会有所不同。

PS看看这个:

  • 如何在一个时钟周期内获得32位input的平方根?
  • 定点log2,exp2
  • 整数C ++ pow(x,1 / y)实现