电源的平方为负指数 [英] Power by squaring for negative exponents
问题描述
我不知道,如果功率平方需要负指数的照顾。我采取了以下code这适用于只有正数。
I am not sure if power by squaring takes care of negative exponent. I implemented the following code which works for only positive numbers.
#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 不帮助如下code似乎是错误的。
Looking at https://en.wikipedia.org/wiki/Exponentiation_by_squaring doesn't help as the following code seems wrong.
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).
编辑: 由于阿米特该解决方案同时适用于正数和负数:
Thanks to amit this solution works for both negative and positive numbers:
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法):
For fractional exponent we can do below (Spektre method):
-
假设你有X ^ 0.5,那么您可以轻松地用这种方法计算平方根:从0到x / 2开始,并保持检查X ^ 2等于结果与否的binary搜索方法。
所以,如果你有X ^(1/3),则必须更换如果*中旬中期和LT; = N
到如果中期* *中旬中期&LT; = N
键,你会得到x.Same事情立方根适用于X ^(1/4),X ^(1/5)等。在x的^(2/5)的情况下,我们可以做(的x ^(1/5))^ 2,并再次降低发现x的第五根的问题。
So, in case you have x^(1/3) you have to replace if mid*mid <= n
to if mid*mid*mid <= n
and you will get the cube root of x.Same thing applies for x^(1/4), x^(1/5) and so on. In the case of x^(2/5) we can do (x^(1/5))^2 and again reduce the problem of finding the 5th root of x.
不过这个时候,你就会意识到,这种方法只有在工作的情况下,当你的可以的转换根,以1 / X格式。那么,我们坚持,如果我们不能转换?不,我们仍然可以继续,我们有这个意愿。
However by this time you would have realized that this method works only in the case when you can convert the root to 1/x format. So are we stuck if we can't convert? No, we can still go ahead as we have the will.
您浮点转换为固定点,然后计算POW(A,B)。假设数字是0.6,这个转换为(24,8)浮点产量楼(0.6 * 1&LT;&LT; 8)= 153(10011001)。大家知道153重新presents小数部分,从而在固定点今(10011001)重新present(2 ^ -1,0,0,2 ^ -3,2 ^ -4,0,0,2 ^ 7)。所以,我们可以通过计算2,3,4和7根x的固定点再次计算POW(一个,0.6)。计算后,我们再次需要得到的结果与1&LT除以浮点;&LT; 8。
Convert your floating point to fixed point and then calculate pow(a, b). Suppose the number is 0.6, converting this to (24, 8)floating point yields 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.
code以上方法可以接受的答案被找到。
Code for above method can be found in the accepted answer.
还有一个基于日志的方法
<一个href="http://stackoverflow.com/questions/19012100/how-math-pow-and-so-on-actualy-works/19072451#19072451">x^y = EXP2(Y * LOG2(X))
推荐答案
整数的例子是32位 INT
算术, DWORD
为32 无符号整型
Integer examples are for 32 bit int
arithmetics, DWORD
is 32bit unsigned int
-
浮动
POW(X,Y)= X ^ÿ
floating
pow(x,y)=x^y
- 在通常评估是这样的:如何Math.Pow(等)actualy工作
- 让分指数可以计算:
POW(X,Y)= EXP2(Y * LOG2(X))
- 在此也可以做到在定点定点BIGNUM战俘
- is usually evaluated like this: How Math.Pow (and so on) actualy works
- so the fractional exponent can be evaluated:
pow(x,y) = exp2(y*log2(x))
- this can be done also on fixed point fixed point bignum pow
整数 POW(A,B)=一^ B
,其中 A&GT; = 0,B&GT; = 0
integer pow(a,b)=a^b
where a>=0 , b>=0
-
这是很容易(你已经有一个)的平方做
this is easy (you already have that) done by squaring
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)=一^ B
,其中 B&GT; = 0
integer pow(a,b)=a^b
where b>=0
-
只需添加几IFS处理负一
just add few ifs to handle the negative 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)=一^ B
integer pow(a,b)=a^b
- 因此,如果
B&LT; 0
那么就意味着1 / powiu(A,-B)
- 正如你所看到的结果不是整数,在所有
- 所以要么忽略这种情况
- 或返回浮点值
- 或添加事半功倍的变量(这样你就可以评估对纯整数算术PI方程)
-
这是float结果:
- so if
b<0
then it means1/powiu(a,-b)
- as you can see the result is not integer at all
- so either ignore this case
- or return floating value
- or add a multiplier variable (so you can evaluate PI equations on pure Integer arithmetics)
this is float result:
float powfii(int a,int b)
{
if (b<0) return 1.0/float(powiu(a,-b));
else return powiu(a,b);
}
整数 POW(A,B)=一^ B
,其中 B
是小数
integer pow(a,b)=a^b
where b
is fractional
- 您可以做这样的事情
A ^(1 / BB)
,其中BB
是整数 - 在现实中,这是生根
- 这样您就可以使用二进制搜索来评价
-
A ^(1/2)
是平方根(一)
-
A ^(1 / BB)
是bb_root(一)
- 在这样做的二进制搜索c从MSB到LSB
- 和评估,如果
POW(C,BB)LT = A
然后离开位的是其他人清楚它 -
这是开方例如:
- you can do something like this
a^(1/bb)
wherebb
is integer - in reality this is rooting
- so you can use binary search to evaluate
a^(1/2)
issquare root(a)
a^(1/bb)
isbb_root(a)
- so do a binary search for c from MSB to LSB
- and evaluate if
pow(c,bb)<=a
then leave the bit as is else clear it this is sqrt example:
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;
}
所以现在只是改变了如果(A * A&GT; X)
与如果(POW(A,BB)X的催化剂)
定点开方例如
//---------------------------------------------------------------------------
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 - >外汇转换:
DWORD(浮点(X)*浮动(_fx32_one))
- 这是FP&LT; - 外汇转换:
浮动(DWORD(X))/浮点(_fx32_one))
-
fx32_mul(X,Y)
是X * Y
它采用80386+ 32位架构汇编(你可以改写为Karatsuba的或其他任何能独立于平台) -
fx32_sqrt(X)
是的sqrt(x)
- 在固定点你应该avare小数位移
- 乘法:
(A&LT;&LT; 16)*(B&LT;&LT; 16)=(A * B&LT;&LT; 32)
- 您需要通过
&GT移回;&GT; 16
来获得结果(A * B&LT;&LT; 16)
- 同样的结果可能溢出32位,因此我用64位的结果汇编
- so this is unsigned fixed point
- high 16 bits are integer
- 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)
isx*y
it uses assembler of 80386+ 32bit architecture (you can rewrite it to karatsuba or whatever else to be platform independent)fx32_sqrt(x)
issqrt(x)
- in fixed point you should be avare of the fractional bit shift
- for multiplication:
(a<<16)*(b<<16)=(a*b<<32)
- you need to shift back by
>>16
to get result(a*b<<16)
- also the result can overflow 32 bit therefore I use 64 bit result in assembly
32位有符号定点战俘C ++的例子
当你把所有的previous步骤在一起,你应该有这样的事情:
When you put all the previous steps together you should have something like this:
//---------------------------------------------------------------------------
//--- 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;
}
//---------------------------------------------------------------------------
我已经测试它是这样的:
I have tested it like this:
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
最接近定点重$ P $的A psentations,B
-
C0
是数学战俘结果 -
C1
是fx32_pow结果 -
D
的区别 - 的希望并没有忘记一些小事,但现在看来似乎正常工作
- 请不要忘记,固定点具有非常有限的precision
- 这样的结果会有所不同有点...
a,b
are floating pointx,y
are closest fixed point representations ofa,b
c0
is math pow resultc1
is fx32_pow resultd
is difference- hope did not forget something trivial but it seems like it works properly
- do not forget that fixed point has very limited precision
- so the results will differ a bit ...
这篇关于电源的平方为负指数的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!