new*_*old 4 c algorithm math recursion
我不确定平方的幂是否可以处理负指数。我实现了以下仅适用于正数的代码。
#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).
编辑：感谢 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 得到浮点数的结果。
上述方法的代码可以在接受的答案中找到。
还有一个基于日志的方法：
整数示例用于 32 位int算术，DWORD是 32 位unsigned int
漂浮的 pow(x,y)=x^y
通常是这样评价的：
因此可以计算分数指数：pow(x,y) = exp2(y*log2(x))。这也可以在固定点上完成：
整数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
只需添加几个ifs 来处理负数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) 是 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);
[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是小数部分。
DWORD(float(x)*float(_fx32_one))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以获得 result (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,bc0 是数学战绩c1 是 fx32_pow 结果d 是区别
希望没有忘记一些琐碎的事情，但它似乎可以正常工作。不要忘记固定点的精度非常有限，因此结果会有所不同......
PS看看这个：
| 归档时间： |
|
| 查看次数： |
3766 次 |
| 最近记录： |