求助函数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<x<2, y=(x-1)/(x+1) $,则有
  $log(x)=2*(y+ y^3/3 + y^5/5 + y^7/7 + y^9/9 + ... y^(2k+1)/(2k+1))(k为非负整数)$

1. typedef unsigned short WORD;
   
2. typedef unsigned long DWORD;
   
3. 
   
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
    
29. 
    
30. 
    
31. 1.0108892860517005,
    
32. 1.0054299011128027,
    
33. 1.0027112750502025,
    
34. 1.0013547198921082,
    
35. 1.0006771306930664
    
36. };
    
37. 
    
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. };
    
61. 
    
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. }
    
69. 
    
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. }
    
78. 
    
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为非负整数)
    
82. 
    
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));
    
93. 
    
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
     
103. 
     
104. r0 *= LN2;
     
105. x=(x-1)/(x+1);
     
106. // 此时x小于0.010830001253374,2*x^7/7小于4.99265e-015,因此计算泰勒公式的前3项足够
     
107. 
     
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. }
     
121. 
     
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));
     
133. 
     
134. r0= GetExpBase2(x);
     
135. x= GetMantissa(x);
     
136. 
     
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. }
     
144. 
     
145. r0 *= LN2;
     
146. x=(x-1)/(x+1);
     
147. s=x*x;
     
148. 
     
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

不错，学习了！