计算⾃然对数的快速算法
引⾔
在1982年，Tateaki. Sasaki 和 Yasumasa Kanada 发表了⼀篇论⽂：。在这篇只有四页的论⽂中，他们介绍了⼀个计算⾃然对数的快速算法。
实现该算法的 C# 程序
我们知道， Framework Class Library 中的 System.Math.Log ⽅法只能对 dobule 数据类型计算⾃然对数。现在我们为 decimal 数据类型实现⼀个 Log 扩展⽅法，如下所⽰：
1using System;
2
3namespace Skyiv.Extensions
4 {
5static class DecimalExtensions
6 {
7static readonly decimal pi2 = 6.283185307179586476925286766559m;
8static readonly decimal ln10 = 2.30258509299404568401799145468m;
9static readonly decimal eps2 = 0.00000000000001m; // 1e-14
10static readonly decimal eps1 = eps2 * eps2; // 1e-28
11
12public static decimal Sqrt(this decimal x)
13 {
14if (x < 0) throw new ArgumentException("Must be nonnegative");
15if (x == 0) return0;
16var y = (decimal)Math.Sqrt((double)x);
17return (y + x / y) / 2;
18 }
19
20public static decimal Log10(this decimal x)
21 {
22return Log(x) / ln10;
23 }
24
25public static decimal Log(this decimal x)
26 {
27if (x <= 0) throw new ArgumentException("Must be positive");
刘安琪
天津蓟县旅游28if (x == 1) return0;
29var k = 0;
30for (; x > 0.1m; k++) x /= 10;
31for (; x <= 0.01m; k--) x *= 10;
32return k * ln10 - NegativeLog(x);
33 }
34奇多音字
35static decimal NegativeLog(decimal q)
36 { // q in ( 0.01, 0.1 ]
37decimal r = q, s = q, n = q, q2 = q * q, q1 = q2 * q;
38for (var p = true; (n *= q1) > eps1; s += n, q1 *= q2)
39 r += (p = !p) ? n : -n;
40decimal u = 1 - 2 * r, v = 1 + 2 * s, t = u / v;
41decimal a = 1, b = Sqrt(1 - t * t * t * t);
42for (; a - b > eps2; b = Sqrt(a * b), a = t) t = (a + b) / 2;
43return pi2 / (a + b) / v / v;
44 }
45 }
46 }
在上述程序中：
第 12 ⾄ 18 ⾏实现 Sqrt 扩展⽅法。使⽤⽜顿迭代法，仅迭代⼀次 (第 17 ⾏)。这是因为第 16 ⾏使⽤ Math.Sqrt (double) ⽅法给出的初值的精度已经达到 15，⽽ decimal 的精度为 28，迭代⼀⾜够了。
相聚今天
第 25 ⾄ 33 ⾏实现 Log 扩展⽅法。该⽅法先将参数 x 变换到 ( 0.01, 0.1 ] 区间中 ( 第 29 ⾄ 31 ⾏)，然后调⽤ NegativeLog ⽅法进⾏计算。
第 35 ⾄ 44 ⾏的 NegativeLog ⽅法就是我们算法的核⼼，使⽤ 椭圆θ函数 ( 第 37 ⾄ 40 ⾏ ) 和 算术⼏何平均法 ( 第 41 ⾄ 43⾏ ) 来快速计算⾃然对数。该算法的原理请参阅参考资料[1]。
第 7 ⾏是事先计算出来的 2π 的值，在第 43 ⾏中使⽤。
第 8 ⾏是事先计算出来的 ln10 的值，在第 32 ⾏和 22 ⾏中使⽤。
我们的程序中使⽤ ln10 的值，⽽参考资料[1]中使⽤ ln2 的值。这是因为 decimal 使⽤⼗进制⽐较好。⽽⼀般的 double 使⽤⼆进制⽐较好。
第 10 ⾏和第 9 ⾏的 eps1 = 10-28 和 eps2 = 10-14 使⽤⼗进制控制计算精度。⽽参考资料[1]中使⽤⼆进制控制计算精度，即 2-p 和 2-p/2。
算法的性能
在上述程序中加⼊⼀些调试语句，就可以在程序运⾏时报告⼀些重要变量的值，结果如下所⽰：
q : 0.1000000000000000000000000000 k: 2
r : 0.0999000009999999000000001000
s : 0.1001000010000001000000001000
u : 0.8001999980000001999999998000
v : 1.2002000020000002000000002000 loop1: 4
b0: 0.8957696680606997488130397041
a : 0.9471678269798758062616205879
b : 0.9471678269798660936897543331 loop2: 3
NL: 2.3025850929940456840179914544
Ln: 2.3025850929940456840179914550
q : 0.0100000000000000000000000001 k: 28
r : 0.0099999900000000010000000001
s : 0.0100000100000000010000000001
u : 0.9800000199999999979999999998
v : 1.0200000200000000020000000002 loop1: 2
b0: 0.3845443947629648387969403232
a : 0.6556979528724023734005530455
b : 0.6556979528723956496570849678 loop2: 4
NL: 4.6051701859880913680359828988
Ln: 59.867212417845187784467777833
上⾯是对 10 和 100000000000000000000000001 分别调⽤ Log ⽅法计算⾃然对数的调试结果：
k 是 Log ⽅法执⾏到第 32 ⾏时的值。其他变量都是第 35 ⾄ 44 ⾏中 NegativeLog ⽅法执⾏时的值。
loop1 是第 38 ⾄ 39 ⾏的 for 循环的执⾏次数。
loop2 是第 42 ⾏的 for 循环执⾏次数。
q, r, s, u, v, a, b 都是计算终了时的值。
b0 是 b 的初值 (第 41 ⾏)。a 的初值是 1 。
NL 是 NegativeLog ⽅法的返回值，即参数 q 的⾃然对数的相反数。
Ln 是 Log ⽅法的返回值，即参数 x 的⾃然对数。
从上述调试结果可以看出，这个算法是⾮常⾼效的。算法中的两个 for 循环执⾏次数⼀般都不超过 4。
测试程序
下⾯是调⽤ decimal 数据类型的 Sqrt、Log 和 Log10 扩展⽅法的测试程序：
1using System;
2using Skyiv.Extensions;
3
4class Tester
5 {
6static void Main()
7 {
8foreach (var x in new decimal[] { 4 / decimal.MaxValue,
90.0000001m, 0.1m, 1, 10, 100000000, decimal.MaxValue })
10 {
11 Console.WriteLine("x : " + x);
12 Console.WriteLine("sqrt: " + x.Sqrt());
13 Console.WriteLine("ln : " + x.Log());
14 Console.WriteLine("lg : " + x.Log10());
15 Console.WriteLine();
16 }
17 }
18 }
运⾏结果如下所⽰：
work$ dmcs Tester.cs DecimalExtensions.cs
work$
x : 0.0000000000000000000000000001
sqrt: 0.00000000000001
ln : -64.472382603833279152503760731
lg : -28.000000000000000000000000000
x : 0.0000001
sqrt: 0.0003162277660168379331998894
ln : -16.118095650958319788125940182
幽默的聊天技巧lg : -6.9999999999999999999999999996
主板诊断卡
x : 0.1
sqrt: 0.3162277660168379331998893545
ln : -2.3025850929940456840179914544
lg : -0.9999999999999999999999999999
x : 1
sqrt: 1
ln : 0
lg : 0写给党的一封信
x : 10
sqrt: 3.1622776601683793319988935445
ln : 2.3025850929940456840179914550
lg : 1.0000000000000000000000000001
x : 100000000
sqrt: 10000
ln : 18.420680743952365472143931638
lg : 8.000000000000000000000000000
x : 79228162514264337593543950335
sqrt: 281474976710656.00000000000000
ln : 66.542129333754749704054283660
lg : 28.898879583742194740518933893
从中可以看出，这个算法的运算精度能够达到 27，只有最后⼀位可能有误差。参考资料
嘴唇薄的男人性格特点1.
2.
3.
4.
5.
