求助函数lnx 以及指数函数用汇编实现的算法

liangbch

应该说，我对计算高精度对数函数值还是很有研究的，目前仍然在研究高精度对数的快速算法，目标是在中等精度（大约几百位到几万位）的情况下，作出一个世界上最快的对数计算器。 在csdn 也曾经遇到过类似的问题，我也做了一些回答，下面给出一些链接，不知对你是否有用。 如何在汇编里面编写求对数的子程序？ 高分求助！（http://topic.csdn.net/t/20050908/10/4256707.html） 如何写一个程序算出lg(i)的值?（http://topic.csdn.net/t/20050113/16/3723330.html）

liangbch

楼上的程序只能精确到2位有效数字，明天给出一个可精确到12为有效数字的版本。

liangbch

下面给出两个计算 double数自然对数的函数，精度可达到12位10进制数的版本。第一个函数是Ln_2，第二个函数是Ln_3（ 实际上我写了三个函数，另外的一个是Ln_1，这个函数效率较低，没有给出）。对于任意双精度浮点数，Ln_3只需要做11次浮点乘除法，就可以达到12位以上的精度。 主要算法： 若$ 1

1. typedef unsigned short WORD;
2. typedef unsigned long DWORD;
4. #define LN2 0.69314718055994530941723212145818 //log(2)
5. #define DIG_NUM 12 //精确到12位有效数字
6. //#define DIG_NUM 15 //精确到15位有效数字
7. static double pow2[]=
8. {
9. 1.3407807929942597e+154,
10. 1.1579208923731620e+077,
11. 3.4028236692093846e+038,
12. 1.8446744073709552e+019,
13. 4.2949672960000000e+009,
14. 6.5536000000000000e+004,
15. 2.5600000000000000e+002,
16. 1.6000000000000000e+001,
17. 4.0000000000000000,
18. 2.0000000000000000,
19. 1.4142135623730951,
20. 1.1892071150027210,
21. 1.0905077326652577,
22. 1.0442737824274138,
23. 1.0218971486541166,
24. // 测试表明，当查表时查到 1.0218971486541166 这一项时，在其后
25. // 阶段，当精度达到10^-12时，应用泰勒公式，只需要计算3项，
26. // 即x, x^3/3, x^5/5,
27. // 从而使得整个计算过程达到最优，至多只需要11次浮点乘除法,具体算法见
28. // 函数Ln_3
31. 1.0108892860517005,
32. 1.0054299011128027,
33. 1.0027112750502025,
34. 1.0013547198921082,
35. 1.0006771306930664
36. };
38. static double pow2_exp[]=
39. {
40. 5.1200000000000000e+002,
41. 2.5600000000000000e+002,
42. 1.2800000000000000e+002,
43. 6.4000000000000000e+001,
44. 3.2000000000000000e+001,
45. 1.6000000000000000e+001,
46. 8.0000000000000000e+000,
47. 4.0000000000000000e+000,
48. 2.0000000000000000e+000,
49. 1.0000000000000000e+000,
50. 5.0000000000000000e-001,
51. 2.5000000000000000e-001,
52. 1.2500000000000000e-001,
53. 6.2500000000000000e-002,
54. 3.1250000000000000e-002,
55. 1.5625000000000000e-002,
56. 7.8125000000000000e-003,
57. 3.9062500000000000e-003,
58. 1.9531250000000000e-003,
59. 9.7656250000000000e-004
60. };
62. short GetExpBase2(double a) // 获得 a 的阶码
63. {
64. // 按照IEEE 754浮点数格式，取得阶码，仅仅适用于Intel 系列 cpu
65. WORD *pWord=(WORD *)(&a)+3;
66. WORD rank = ( (*pWord & 0x7fff) >>4 );
67. return (short)(rank-0x3ff);
68. }
70. double GetMantissa(double a) // 获得 a 的 尾数
71. {
72. // 按照IEEE 754浮点数格式，取得尾数，仅仅适用于Intel 系列 cpu
73. WORD *pWord=(WORD *)(&a)+3;
74. *pWord &= 0x800f; //清除阶码
75. *pWord |= 0x3ff0; //重置阶码
76. return a;
77. }
79. // 对数计算公式：
80. // 若x>1,y=(x-1)/(x+1),则有
81. // log(x)=2*(y+ y^3/3 + y^5/5 + y^7/7 + y^9/9 + ... y^(2k+1)/2k)(k为非负整数)
83. // 计算log(x)的第2个版本，
84. // x>0
85. double Ln_2(double x)
86. {
87. double r0,r1;
88. double item;
89. double s;
90. int i;
91. if (x<1)
92. return -(Ln_2(1.0/x));
94. for (r0=0,i=0;i<=9+5;i++)
95. {
96. if ( x>= pow2[i])
97. {
98. r0+= pow2_exp[i];
99. x/= pow2[i];
100. }
101. }
102. //此时x必然小于1.0218971486541
104. r0 *= LN2;
105. x=(x-1)/(x+1);
106. // 此时x小于0.010830001253374,2*x^7/7小于4.99265e-015,因此计算泰勒公式的前3项足够
108. s=x*x;
109. r1=item=x; //计算泰勒公式第1项
110. #if DIG_NUM ==12
111. for (i=3;i<=5;i+=2) //精确到12位有效数字
112. #else
113. for (i=3;i<=7;i+=2) //精确到15位有效数字
114. #endif
115. {
116. item *= s;
117. r1 += item/i; //泰勒公式其他3项
118. }
119. return r0+r1*2;
120. }
122. // 计算log(x)的第3个版本，速度更快，
123. // 精确到10^12,至多只需要11次浮点乘除法
124. // 要求：x必须是 IEEE754标准浮点格式，字节顺序是intel格式，在某些CPU可能不能正确工作
125. double Ln_3(double x)
126. {
127. double r0,r1;
128. double item;
129. double s;
130. int i;
131. if (x<1)
132. return -(Ln_3(1.0/x));
134. r0= GetExpBase2(x);
135. x= GetMantissa(x);
137. for (i=9+1;i<=9+5;i++)
138. {
139. if ( x>= pow2[i])
140. {
141. r0+= pow2_exp[i]; x/= pow2[i];
142. }
143. }
145. r0 *= LN2;
146. x=(x-1)/(x+1);
147. s=x*x;
149. r1=item=x; //计算泰勒公式第1项
150. #if DIG_NUM ==12
151. for (i=3;i<=5;i+=2) //精确到12位有效数字
152. #else
153. for (i=3;i<=7;i+=2) //精确到15位有效数字
154. #endif
155. {
156. item *= s;
157. r1 += item/i; //泰勒公式其他2项
158. }
159. return r0+r1*2;
160. }

复制代码

无心人

他需要的是51系列的代码 可能不能用某些技巧的 因为硬件并不支持浮点操作 且是8位CPU

liangbch

我是知道这一点，所以给出2个版本的算法。 我不可能知道每一种CPU的汇编，只能给出算法和C语言代码。具体到某一种CPU,当然需要他自己用汇编语言去实现这个算法了，其中包括实现 类似C语言浮点数的 表示方法和运算方法。

无心人

呵呵 C51不知道具体如何运算的

liangbch

下面给出计算指数函数的简单算法 a,b为任意浮点数，求$a^x$. step1: 将x表示为2进制，整数部分存入数组 $a_i, 0<=i<=k$,小数部分存入数组$b_i,1<=i<=m$，如13.75可表示为 1101.11, 数组a[0..3]={1,0,1,1},数组b[1..2]={1,1} step2: 计算$a^(2^0),a^(2^1),a^(2^2),直到a^(2^k)$, 存入数组$e1_i,0<=i<=k$, 对于13.75，需要计算到$a^8$. step3: 计算$a^(2^-1),a^(2^-2),b^(2^-3)$, 直到$b^(2^-m)$,存入数组$e2_i,1<=i<=m$, 对于13.75，需要计算到$a^(1/4)$ step4: $r=1$ step5: 如 $a_i>0,0<=i<=k$,则计算$r = r ** e1_i$ step6: 如 $b_i>0,1<=i<=m$,则计算$r = r ** e2_i$

liangbch

楼上用到开平方算法，具体算法可在本站中查找。

silitex

不错，学习了！