紧急求助:如何只用"加减乘除"基础运算,求任意数的任意次方

[withpointer]

问题同上,如何只用"加减乘除"来求任意数(整数,小数...)的任意次方(整数,小数...)!不使用任何库函数

[yaos]

用牛顿叠代吧，很快的

：）

具体的算法要找书的

只用+-*/的

[liangbch]

http://numbers.computation.free.fr/Constants/constants.html 给出ln(2)的许多种计算方法，其最后的AGM迭代算算法可计算任意数的自然对数。

[withpointer]

to liangbch(宝宝):
实在抱歉,前几天才从北京回来,一直没有时间能仔细的看以下您的算法.
今天看了一下,一点收获都没有(见笑了,小弟的数学功底很差),如果您有空,可否将您的算法用浅显的语言或者程序算法表达出来???有些数学计算名词(比如:牛顿叠代等我都不太明白),
小弟先在这里谢谢了. 如果有代码可以发到我的邮箱里 [e-mail]

[gxqcn]

如果不要求速度，以下思路应不失为一种较好的选择，
实现起来也比较容易（仅针对实数集）。


目标：求 x 的 r 次方。


分析：

为简化讨论，不妨规定 x、r ≥ 0(若x<0，可能会无意义；当r<0时，将结果取倒即可)

(1)当 r 为非负整数时，很容易解决：
将 r 用二进制表示：(r[0],r[1],...,r[n])，（r[i] 皆取 0 或 1）其中r[0]表示“个位”，r[1]表示“十位”，。。。
则 x^r 可表示成一组数的乘积，它们来自于数列{ x, x^2, x^4, ..., x^(2^n) }，参与计算取舍决定于对应的 r[i] 是否为“1”；
而上述数列后项是前一项的平方，非常容易设计；

(2)当 r 不为整数时，令r = r1 + r0 ( r1 为非负整数，0＜r0＜1 )， r2 = 1/r0 > 1，
则 x^r = x^r1 * x^r0 = x^r1 *{[x^(-r0)]^(-1)} = x^r1 / (x^r2)；
上面的推导比较有技巧，我把小学生的大、中、小括号全用上了，以便阅读。
至此，我们把 x^r 表示成了两个幂的除法关系，
其中分子已在(1)中解决；分母可看作新的x^r项目，经过适当次数的迭代即可。

其实，更便于编程实现的是将 r 直接用二进制表示，(r[n], r[n-1], ...,r[0], r[-1], r[-2], ..., r[-m])，其中 r[0]表示“个位”，r[1]表示“十位”，r[-1]表示“十分位”，。。。
则 x^r 可表示成两组数的乘积之商，它们分别来自于数列{ x^(2^n), x^(2^(n-)), ..., x }及{ x, x^2, x^4, ..., x^(2^m) }，参与计算取舍决定于对应的bit是否为“1”。

为充分有效的共享数据，可计算max(m,n)，只构造一个最长的数列，且在构造上述两数列的同时即可计算连积。

以上算法，仅用到乘法和除法，但愿对你有参考价值；仓促行文，也许有失误之处。

[gxqcn]

我给的算法对于小数指数也适用啊。
仅用到平方、乘法、除法运算，未使用任何库函数，应符合你的要求。

比如说你要算的 9^1.6,

将 1.6 表示成二进制为 (1.[1001])，其中[]内的数据为循环节（无理数不存在循环节），
即 (1.6)10 = (1.1001 1001 ...)2，

则 9^1.6 = 9 ^ ( 1 + 1/2 + 1/16 + 1/32 + 1/256 + ... )
= 9 * ( 9 ^ 1/2 ) * ( 9 ^ 1/16 ) * ( 9 ^ 1/32 ) * ( 9 ^ 1/256 ) * ...
= 9 / [ ( 9^2 ) * ( 9^16 ) * ( 9^32 ) * ( 9^256 ) * ... ]
而 9^16 = ((( 9^2 )^2)^2)^2，9^32 = ( 9^16 )^2，...

不知这样说你可闹明白了没？

[gxqcn]

对不起，完全错了！因为 x^(1/r) ≠ 1/(x^r)，彻底失败！
如果对你造成误导，那就是罪过罪过了。。。

[withpointer]

to gxqcn(GxQ):

我用了您提供的(2)方法中的x^r1 / (x^r2)；但是计算结果不正确啊,我的代码在下面:
代码是VB的
Private Sub Command1_Click()
Dim X As Double, R As Double, R0 As Double, R1 As Double, _
R2 As Double
'9 ^ 1.6
X = 9: R = 1.6

R1 = 1: R0 = 0.6
R2 = 1 / R0

MsgBox X ^ R1 / (X ^ R2)
End Sub

[zzwu]

由于 R0 = 0.6,

故 R2 = 1/R0 = 1/0.6 = 1.66666...

这样
X ^ R1 / (X ^ R2)

分母中仍出现非整数次的乘方了!

[gxqcn]

对不起，(2)的推导是错误的。:(
但其后面的方法是正确的，且效率较高。

[withpointer]

to gxqcn:
还有没有类似(2)的方法呢?能不能推导出只做任意数的任意整数次方的公式呢?

[gxqcn]

“任意数的任意整数次方”用(1)即可，效率也较高。

[withpointer]

我的意思是 x ^ (1 / n) 这类算式给转换成 没有小数次方的公式,就是起到您的公式(2)的作用

[withpointer]

to gxqcn(gxq):
"9 / [ ( 9^2 ) * ( 9^16 ) * ( 9^32 ) * ( 9^256 ) * ... ]" 这个算式是怎么得出来的

9 ^ x的x的变换规律是什么?????

[gxqcn]

我休假完，重返岗位。

其实，求任意数的任意次方，仅需通过乘法、平方、开方实现（当然，开方的实现可能需要+-*/），仍举例说明，计算 9^1.6：

将 1.6 表示成二进制为 (1.[1001])，其中[]内的数据为循环节（无理数不存在循环节），
即 (1.6)10 = (1.1001 1001 ...)2，

则 9^1.6 = 9 ^ ( 1 + 1/2 + 1/16 + 1/32 + 1/256 + ... )
= 9 * ( 9 ^ 1/2 ) * ( 9 ^ 1/16 ) * ( 9 ^ 1/32 ) * ( 9 ^ 1/256 ) * ...
而 9^1/16 = ((( 9^1/2 )^1/2)^1/2)^1/2（连续进行4次开方运算），
9^1/32 在 9^1/16 基础上再开方，。。。

上次因为老婆要生儿子了，比较激动和兴奋，所以犯了低级错误，
在回家途中以找到问题的症结。

总结一下，整个算法的秘诀为（效率可能不会太高，但实现容易）：
将指数用二进制表示，小数点左边的连续平方，小数点右边的连续开方，将对应位为1的数连乘即可。

[withpointer]

to gxq:
你的算法，还要开平方，如何开平方啊(只用+-*/)???

[withpointer]

忘了恭喜gxq大哥了，祝baby健康成长!

[liangbch]

今天中午，写一个通用开方程序，可以对double型数开方，该程序效率很高，除了少数函数外，该程序的算法也适用于其他数据类型，甚至是自己定义的高精度大数。该程序的另一个特点是不使用除法，因为对大数运算来说，除法算法比较复杂，所以可以使用这里的算法轻松的实现大数的开方。

#include "stdafx.h"
#include "math.h"

//一个使用 加、减、乘来,对 double数 开方的算法，作者: [e-mail]
// 除 GetExpBase2，reciprocal_0 这两个函数仅仅适用于double型变量外，
// 其他函数的算法也适用于其他类型的数据

typedef unsigned short WORD;

short GetExpBase2(double a)
{
// 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu
WORD *pWord=(WORD *)(&a)+3;
WORD rank = ( (*pWord & 0x7fff) >>4 );
return (short)(rank-0x3ff);
}

//当n较大时，这个函数效率较低
//当n较大，可以使用平方算法优化
double pow(double a,int n)
{
if (n<=0)
return 0;
else if (n==0)
return 1;
else if (n==1)
return n;
//----------------------------
double r=a;

for (int i=1;i<n;i++)
r *=a;
return r;
}

//计算2 ^ n
double pow2(int n)
{
double r=1.00;
int i;
if ( n >0)
{
for (i=0;i<n;i++) r *=2 ;
}
else
{
for (i=0;i<-n;i++) r /=2 ;
}
return r;
}


//计算 a的绝对值
double myAbs(double a)
{
if (a <0) return 0-a;
else return a;
}

// 返回 1/a 的初值(近似值)
double reciprocal_0(double a)
{
// 按照IEEE 754浮点数格式，求其倒数，仅仅适用于Intel 系列 cpu
double x0=a;

WORD dwMantissa; //尾数高4位
WORD rank;
WORD *pWord=(WORD *)(&x0)+3;

dwMantissa= *pWord & 0x000f;

rank= (~(*pWord) + 0x7ff0);

*pWord = ((dwMantissa) | ( rank & 0xfff0 ));

switch ( dwMantissa)
{
case 0: break;
case 1: x0 *=0.8858;break;
case 2: x0 *=0.7901;break;
case 3: x0 *=0.7091;break;
case 4: x0 *=0.6400;break;
case 5: x0 *=0.5805;break;
case 6: x0 *=0.5289;break;
case 7: x0 *=0.4839;break;
case 8: x0 *=0.4444;break;
case 9: x0 *=0.4096;break;
case 10: x0 *=0.3787;break;
case 11: x0 *=0.3512;break;
case 12: x0 *=0.3265;break;
case 13: x0 *=0.3044;break;
case 14: x0 *=0.2844;break;
case 15: x0 *=0.2664;break;
}

return x0;
}

// 返回 1/a 的初值, 采用4步迭代，比reciprocal_0 精度更高
double reciprocal_1(double a)
{
double x0,x1;
x0=reciprocal_0(a);
for (int i=0;i<4;i++) // 4 步迭代
{
x1= x0 * (2-a*x0);
x0= x1;
}
return x0;
}

// 返回 1/a ，在reciprocal_0 的基础上，采用6步迭代算法，
double reciprocal(double a)
{
double x0,x1;
x0=reciprocal_0(a);
for (int i=0;i<6;i++)
{
x1= x0 * (2-a*x0);
x0= x1;
}
return x0;
}

// 计算 a 开 n 次方的初值
double sqrt_0(double a, int n)
{
double low,high,mid;
double r;
if (a <=0)
return 0;
//---------------
int e=GetExpBase2( a);

if (e >=0)
{ low=pow2(e/n); high=low*2; }
else
{ high=pow2(e/n); low=high/2; }

while (true) //使用2分法，求得更加精确的初值
{
mid=(low + high)/2;
r=pow(mid,n);

if (r > a ) high=mid;
else low=mid;

if ( r >a && r< 1.1*a || r <a && r >0.9*a )
break;
}
return mid;
}

double sqrt_n(double a, int n)
{
double x0,x1,x2,d;

if (a <=0)
return 0;

x0=sqrt_0(a,n); // a ^ (1/n)的初值
x1=reciprocal_1(x0); // x0 的倒数的近似值
d= x1 * 1e-16;

while (true)
{
double t=n+1 - a * pow(x1,n);
x2= x1 * t / n;
if ( myAbs(x2 - x1) < d)
break;
x1=x2;
}
return reciprocal(x1);
}

void test1()
{
double a;
double r1,r2;
int n;
printf("calc a ^ (1/n)\n");
printf("a=?"); scanf("%lf",&a);

printf("n=?(n>0)"); scanf("%d",&n);

r1=sqrt_n(a,n);
r2=exp(log(a)/n); //使用数学函数计算开方

printf("result1=%.15e,result2=%.15e\n", r1,r2);
}

int main(int argc, char* argv[])
{
test1(); return 0;
}

[liangbch]

to withpointer(NewOne) :
程序都给你写出来了，将给出的代码编程成程序，仔细调试研究去吧！

[gxqcn]

可见 http://community.csdn.net/Expert/FAQ/FAQ_Index.asp?id=182841

除了里面介绍的方法，也可以考虑二分法试算：
当被开方数 y 大于1时：选初始 max = y, min = 1.0，try = (max + min)/2；
当被开方数 y 小于1时：选初始 max = 1.0, min = y，try = (max + min)/2；
比较 try^2 与 y 的大小，修正 max 或 min 的值进行下一次迭代，直到达到预期精度。

[liangbch]

重写了计算n次方根的初值函数，将2分法换成速度更快的算法。

#include "stdafx.h"
#include "math.h"

// 一个使用 加、减、乘来,对 double数 开方的算法，作者: [e-mail]
// 除 GetExpBase2，reciprocal_0,GetMantissa 这3个函数仅仅适用于double型变量外，
// 其他函数的算法也适用于其他类型的数据

typedef unsigned short WORD;

const double A_2_pow[]=
{
2, //2 ^1
1.4142135623730950488016887242097, //2 ^0.5
1.1892071150027210667174999705605, //2 ^0.25
1.0905077326652576592070106557607, //2 ^0.125
1.0442737824274138403219664787399, //2 ^0.0625
1.0218971486541166782344801347833 //2 ^0.03125
};

/*
double 数 n 可表示为 m * 2^e ( 1<=m<2, e为整数)
m称为n的尾数，e称为n的阶码
*/

short GetExpBase2(double a) // 获得 a 的阶码
{
// 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu
WORD *pWord=(WORD *)(&a)+3;
WORD rank = ( (*pWord & 0x7fff) >>4 );
return (short)(rank-0x3ff);
}

double GetMantissa(double a) // 获得 a 的 尾数
{
// 按照IEEE 754浮点数格式，取得尾数，仅仅适用于Intel 系列 cpu
double r=a;
WORD *pWord=(WORD *)(&r)+3;

*pWord &= 0x800f; //清除阶码
*pWord |= 0x3ff0; //重置阶码
return r;
}

double myAbs(double a) //计算 a的绝对值
{
if (a <0) return 0-a;
else return a;
}

//当n较大时，这个函数效率较低
//当n较大，可以使用平方算法优化
double pow(double a,int n)
{
if (n<=0) return 0;
else if (n==0) return 1;
else if (n==1) return n;
//----------------------------
double r=a;
for (int i=1;i<n;i++)
r *=a;
return r;
}

double pow2(int n) //计算2 ^ n
{
double r=1.00;
int i;
if ( n >0)
for (i=0;i<n;i++) r *=2 ;
else
for (i=0;i<-n;i++) r /=2 ;
return r;
}

double pow2(int n1,int n2) //计算2 ^ (n1/n2),近似值
{
int i;
double exp =(double)n1 / (double)n2;
double r;
double baseExp;
int n= sizeof(A_2_pow)/sizeof(double); // 数组元素的个数

for (r=1,baseExp=1.0,i=0;i<n;)
{
if (exp>baseExp)
{
exp -= baseExp;
r *= A_2_pow[i];
}
else
{
i++;
baseExp /=2;
}
}
return r;
}

double reciprocal_0(double a) // 返回 1/a 的初值(近似值)
{
// 按照IEEE 754浮点数格式，求其倒数，仅仅适用于Intel 系列 cpu
double x0=a;

WORD dwMantissa; //尾数高4位
WORD rank;
WORD *pWord=(WORD *)(&x0)+3;

dwMantissa= *pWord & 0x000f;
rank= (~(*pWord) + 0x7ff0);
*pWord = ((dwMantissa) | ( rank & 0xfff0 ));

switch ( dwMantissa)
{
case 0: break;
case 1: x0 *=0.8858;break;
case 2: x0 *=0.7901;break;
case 3: x0 *=0.7091;break;
case 4: x0 *=0.6400;break;
case 5: x0 *=0.5805;break;
case 6: x0 *=0.5289;break;
case 7: x0 *=0.4839;break;
case 8: x0 *=0.4444;break;
case 9: x0 *=0.4096;break;
case 10: x0 *=0.3787;break;
case 11: x0 *=0.3512;break;
case 12: x0 *=0.3265;break;
case 13: x0 *=0.3044;break;
case 14: x0 *=0.2844;break;
case 15: x0 *=0.2664;break;
}

return x0;
}

// 返回 1/a 的初值, 采用4步迭代，比reciprocal_0 精度更高
double reciprocal_1(double a)
{
double x0,x1;
x0=reciprocal_0(a);
for (int i=0;i<4;i++) // 4 步迭代
{
x1= x0 * (2-a*x0);
x0= x1;
}
return x0;
}

// 返回 1/a ，在reciprocal_0 的基础上，采用6步迭代算法，
double reciprocal(double a)
{
double x0,x1;
x0=reciprocal_0(a);
for (int i=0;i<6;i++)
{
x1= x0 * (2-a*x0);
x0= x1;
}
return x0;
}

double sqrt_0(double a, int n) // 计算 a 开 n 次方的初值
{
double r;
int mod;

if (a <=0)
return 0;
//---------------
int e=GetExpBase2( a);
double Mantissa=GetMantissa(a);
//-------------------------
r=1.00 + (Mantissa-1)/n;

mod=e % n;
if (mod<0)
mod +=n;

r *= pow2((e-mod)/n);
r *= pow2(mod,n);
return r;
}

double sqrt_n(double a, int n) /// 计算 a 开 n 次方
{
double x0,x1,x2,d;
if (a <=0)
return 0;

x0=sqrt_0(a,n); // a ^ (1/n)的初值
x1=reciprocal_1(x0); // x0 的倒数的近似值
d= x1 * 1e-16;

while (true)
{
double t=n+1 - a * pow(x1,n);
x2= x1 * t / n;
if ( myAbs(x2 - x1) < d)
break;
x1=x2;
}
return reciprocal(x2);
}

void test1()
{
double a;
double r1,r2;
int n;
printf("calc a ^ (1/n)\n");
printf("a=?"); scanf("%lf",&a);

printf("n=?(n>0)"); scanf("%d",&n);

r1=sqrt_n(a,n);
r2=exp(log(a)/n); //使用数学函数计算开方

printf("result1=%.15e,result2=%.15e\n", r1,r2);
}

int main(int argc, char* argv[])
{
test1(); return 0;
}