重磅消息–一种新的任意精度对数算法研究成功

经过许多个月的艰辛努力,我自己独立开发的高精度算术库(BCMP)终于可以告一段落了。另外,由我自己独创的对数log(x)算法也已经完成。测试报告显示,这个算法的性能是非常令人振奋的,在100-10000位精度下,可比现在最快软件Mathematica[注9]快2.8到8倍。

BCMP和GMP[注1]类似,利用BCMP可以计算任意精度整数和浮点数的算术运算,包括加,减,乘,除,开平方等。和GMP类似,我开发这套大数库的主要目标为追求极致的性能。现在,这一目标已基本达到。在中等精度(1万位十进制)以下,BCMP在大多数情况下,显著领先于GMP。当精度更高时,bcmp则慢于GMP,这是因为,大数运算库的关键是乘法,而不同大小的操作数需要使用不同的算法,没有一种算法可以适应于任意情况。BCMP中实现了basecase乘法[注2],karatsuba乘法[注3],Toom-Cook[注4]系列乘法,而尚未实现适应于操作数很大的SSA[注5]算法,或者FFT[注6]算法和NTT[注7]算法。和GMP不同,BCMP提供了2套接口,第一套接口支持2进制运算,内部表示采用\(2^{30}\)进制,另一套接口支持10进制运算,内部表示采用\(10^9\)进制。BCMP计划同时支持32位和64位方式编译,也计划同时支持在windows下编译和Linux下编译。现在已实现32位模式下的用Visual Studio编译,今后将逐步完善其他功能。

BCMP主要使用C语言和汇编语言开发,为了实现极高的性能,使用了大量的汇编代码,对关键部分做了精心的优化,除了普通的汇编指令外,还大量使用SSE2,SSE4和少数AVX2指令。

新的任意精度对数函数算法是我独创的,不同于任何一个公开发表的算法。新算法基于质数分解和比特突发算法,大大加速了计算速度,从算法复杂度来说,它和目前最快的算法AGM算法[注8]有相同的复杂度,都是O(log(P)M(P)),p为精度,其中M(P)代表两个P位整数相乘所需要的基本运算的次数。但是我的算法有更小的常数因子,速度更快。下表列出了几个最著名数学软件和大数库与BCMP在计算log(x)的性能对比结果。从这个表中可以看出,当精度为100到10000位时,我的新算法是MPFR[注11]的3.7到8.6倍

测试平台:
处理器: Intel(R) Core(TM) i7-2600K CPU @3.40GH
操作系统: Windows7 32-bit
Maple[注10] 版本:17.0
Mathematica版本: 10.3.0 32-bit
Pari[注12]版本: 2.9.3
mpfr版本:3.1.6

 

位数 时间,单位为毫秒
Maple Mathematica Pari MPFR BCMP_log1 BCMP_log2
20 0.029941 0.00247424 0.003114 0.004833 0.008654 0.000762
30 0.040354 0.00305579 0.003744 0.005667 0.008218 0.000975
40 0.043768 0.003999 0.004687 0.006456 0.010237 0.001043
50 0.045405 0.00464357 0.005267 0.006782 0.010515 0.001217
60 0.061336 0.00551004 0.007068 0.008939 0.012052 0.002432
70 0.063722 0.00660306 0.008129 0.009373 0.012467 0.002622
80 0.06523 0.00748533 0.008918 0.010451 0.014011 0.002847
90 0.067488 0.00873332 0.009401 0.010859 0.013989 0.002967
100 0.069835 0.00976343 0.011404 0.012145 0.017113 0.003194
125 0.072167 0.0121835 0.013434 0.014376 0.020085 0.003843
158 0.079578 0.0186991 0.020295 0.018568 0.026015 0.004825
199 0.084953 0.0257783 0.027464 0.023096 0.032686 0.005604
251 0.101058 0.0337484 0.039437 0.030076 0.041313 0.006338
316 0.111391 0.0496319 0.057234 0.03874 0.057398 0.008075
398 0.125592 0.0655544 0.090215 0.054713 0.076384 0.009971
501 0.140205 0.0873334 0.126219 0.076801 0.100851 0.012947
630 0.169863 0.143491 0.173471 0.1227 0.126775 0.025564
794 0.201257 0.213486 0.262401 0.192741 0.175562 0.022624
1000 0.26965 0.331134 0.378275 0.29729 0.260477 0.041692
1258 0.344599 0.4512 0.545905 0.425357 0.321757 0.056777
1584 0.437988 0.615793 0.836217 0.671066 0.504561 0.077545
1995 0.570668 0.871282 1.219222 1.042976 0.671209 0.133183
2511 0.807129 1.20226 1.723848 1.448282 0.983084 0.190615
3162 1.07266 1.68033 2.561828 2.257336 1.489376 0.273305
3981 1.47754 2.42274 3.66907 2.959093 2.155975 0.510987
5011 2.15234 3.27768 5.208427 3.788049 3.193957 0.727249
6309 2.96615 4.7273 7.723928 5.709551 4.507839 1.019573
7943 4.0625 6.69787 11.071721 7.894751 6.753858 1.748821
10000 6.33854 9.14068 15.950279 10.198301 9.649934 2.538205
17782 15.1146 21.4501 38.916218 23.862088 23.569465 6.929607
31622 42.25 56.5504 94.352581 68.980759 62.309453 18.861327
56234 681 117.001 207.879794 131.337302 146.628323 45.943156
100000 2012.67 260.002 587.256462 296.982964 364.771994 117.998387

性能对比折线图,以mpfr为基准,令mpfr的运行时间为1,值越小越好。

性能对比折线图2,去掉Maple的测试结果,仍以mpfr为基准,值越小越好。

说明:BCMP用两种方法实现了两个版本的log函数,第一个使用AGM算法,在上表中标记为BCMP_log1, 第二个使用新算法,标记为BCMP_log2

测试代码下载:
mathematica,maple,pari,mpfr,bcmp,time

测试结果下载:
test_result
更多的任意精度浮点软件性能比较,请参阅Comparison of multiple-precision floating-point software

1.GMP是一个用于任意精度算术的免费库,用于对有符号整数,有理数和浮点数进行操作。GMP的主要目标应用是密码应用和研究,互联网安全应用和计算机代数系统。GMP是一个非常优秀的大数库,已有26年的历史。GMP的作者曾宣称,GMP是地球上最快的大数库,原文“the fastest bignum library on the planet!”。因为GMP非常快,所以很多应用软件使用GMP,如计算机代数系统Mathematica和Maple中的整数运算就是调用GMP来实现的。
关于GMP更多的信息,请参阅https://gmplib.orghttps://en.wikipedia.org/wiki/GNU_Multiple_Precision_Arithmetic_Library

2.basecase乘法,也叫长乘法(Long multiplication),请参阅 https://en.wikipedia.org/wiki/Multiplication_algorithm

3.Karatsuba乘法,也称Karatsuba算法,这个算法由Anatoly Karatsuba在1960年发现的, 是一种比长乘法更快的算法,请参阅https://en.wikipedia.org/wiki/Karatsuba_algorithm

4.Toom–Cook乘法,狭义的Toom-Cook算法指Toom-3算法,可以看做是Karatsuba乘法的推广。广义的Toom-Cook乘法包括Toom4,Toom5等一大类算法。这个算法最早由俄罗斯数学家Andrei Toom提出,而Stephen Cook(是美国和加拿大计算机科学和数学家,图灵奖的获得者)给出更清晰的描述,请参阅维基百科https://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication

5.SSA算法指Schönhage–Strassen algorithm ,他是德国数学家和计算机科学家Arnold Schönhage和德国数学家Volker Strassen在1971开发出来的一个算法。关于这个算法的更多信息,请参阅维基百科https://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm 和百度百科https://baike.baidu.com/item/SSA/18691142,原始论文请参阅《A GMP-based Implementation of Schönhage-Strassen’s Large Integer Multiplication Algorithm》

6.FFT指快速傅里叶变换(fast Fourier transform),是基于复数的一种线性变换,这个算法在电子技术领域得到极其广泛的应用,是20世纪10大算法之一。另外,这个算法可用来计算大数乘法。关于FFT的更多信息,请参阅 https://en.wikipedia.org/wiki/Fast_Fourier_transform。关于20世纪10大算法,请参阅http://science.dataguru.cn/article-9643-1.html

7.NTT指数论变换(Number-theoretic transform),它是FFT推广到有限域上的一种变换算法,和FFT不同,他不使用复数域上的加,减和乘法,而是使用模算术。apfloat中的大数乘法用的就是NTT算法。感兴趣的读者可以参阅蒋增荣的《数论变换》或者孙琦的《快速数论变换》等书。

8.AGM指算法几何平均数,关于这个算法的描述请参阅澳大利亚数学家和计算机科学家Richard P. Brent《Fast Algorithms for High-Precision Computation of Elementary Functions》和的Borwein兄弟的《The arithmetic-geometricmean and fast computation of elementary functions》和。Borwein兄弟都是数学家,他们出生在苏格兰。

9.Mathematica 是一款功能强大的科学计算软件,和MATLAB、Maple 并称为三大数学软件。请参阅百度百科
https://baike.baidu.com/item/Mathematica

10.Maple和Mathematica类似,也是一款功能强大的科学计算软件。请参阅百度百科https://baike.baidu.com/item/maple/2306572

11.MPFR是基于GMP的一套2进制浮点库,请参阅维基百科https://en.wikipedia.org/wiki/GNU_MPFR

12.PARI/GP是一个计算机代数系统,请参阅维基百科https://en.wikipedia.org/wiki/PARI/GP

最短的计算大数乘法的c程序

说明:
1.这个程序接收2个从键盘输入的整数,计算他们的乘积,并输出结果。输入的两个整数的总长度不能大于99.
2.这个程序没什么大用,只是用来玩玩儿而已。
3.这个程序的主要目标是,使用尽可能短的代码来实现大数乘法。上面的代码
可在VC下编译并运行. 在GCC下编译,可省略#include语句和void关键字,
去除回车和不必要的空格,总长度仅仅194个字节。
另外,程序刻意避免使用数组来存贮中间结果和最终结果。
为此,使用了递归函数,同时,递归的使用也简化了代码。
4.在实际工作中,千万不要写这样的程序,否则会被骂死。
5.不要用这个程序考你的学生和面试者,即使他宣称精通C语言。

此类最短程序的特点
1.经常使用全局变量,全局变量的优点是
1).自动初始化数组和单变量为0,可省去某些变量初始化语句。
2).数组初始化为0也使得逻辑更简单,可省去某些边界值的判断。
3).在子程序,直接使用全局变量可省去某些参数定义和参数传递语句。

2.在表达式,大量使用“++”或者“–”之类运算符,此类语句往往起到
一箭双雕的效果,可有效的缩短代码长度.但在工作中,我强烈反对使用
这类运算符。

3.在比较语句中,很少使用if(i>0)这类语句,而是使用“if(i)”这样的
写法,这种写法比”>0″少了2个字母。

AGM函数近似值的估计

AGM是 Arithmetic-Geometric Mean的缩写,意为算术几何平均数。其定义为,给定两个正实数a0和b0,我们定义一个迭代过程

\( a_{i+1}= \frac {a_i+b_i} {2}, \qquad b_{i+1}= \sqrt {a_i \cdot b_i } \)

前者求两个数的算术平均数,后者求2个数的几何平均数。这两个序列有同样的极限,其值表示为A(a0,b0),在AGM迭代过程中,序列收敛地非常快,一旦\(a_i\) 和\(b_i\)大小相近时,每迭代一次,有效数字加倍。下表显示了这两个序列的收敛速度[1]

n \(a_n\) \(b_n\)
0 24 6
1 15 12
2 13.5 13.41640786499873817845
3 13.45820393249936908922 13.45813903099098487720
4 13.45817148174517698321 13.45817148170605385831
5 13.45817148172561542076 13.458171481725615420761

高斯发现,AGM函数可以用来求椭圆积分的值

\( A(1,x)= \frac{\pi} {2F(x)} \)

F(x)为椭圆积分函数

\( \displaystyle \int_0^{\pi \over 2} \frac{d\theta} { \sqrt {1-(1-x^2 ) \sin^2 \theta)}} ) \)

AGM函数不但可以快速的计算椭圆积分,也可用来计算对数函数ln(x).
[2]给出一个椭圆积分F(x)和对数函数ln(s)的关系式

\(
F(\frac{4}{s})=\ln(s)+ \frac{4\ln(s)-4}{s^2}+\frac{36\ln(s)-42}{s^4} + \frac{1200 \ln(s)-1480}{3s^6} + O(\frac{1}{s^8}) \)

如果s是足够的大,F(4/s)是非常好的ln(s)的近似值。为了计算ln(s)到p位2进制数。s必须大于\(2^{p \over 2 } \)。故为了计算ln(x),我们首先需要找到一个实数s和整数m,使得

\( s=x \cdot2^m \gt 2^{p \over 2} \)

[2]给出如下求ln(x)的公式

\( \displaystyle ln(x)=\frac{\pi}{2A({4 \over s},1)} -m \ln (2) \qquad (1)\)

我们将这个公式稍作一下变形

\( \displaystyle \ln(x)=\frac{{1 \over 8}\pi s}{A({s \over 4},1)} -m \ln (2)= \frac{{1 \over 2}\pi \cdot x \cdot 2^{m-2}}{A(x \cdot 2^{m-2},1)} – m \ln(2)
\qquad  (2) \)

反过来,利用这个公式,我们可以求AGM(1,y)的近似值,这里要求y>>1。在公式(2)中,我们取x为1,这样不但可以消去方程中的x,也可以消去ln(x),这样公式(2)变为

\(  \displaystyle \frac {{1 \over 2} \pi \cdot 2^{m-2}}{A(2^{m-2},1)}=m \ln(2) \quad (3)\)

令 \(y=2^{m-2}\), 则

\( \displaystyle A(y,1)= \frac { {1 \over 2} \pi \cdot y} {(log2(y)+2) \ln 2} = y \frac{{1 \over 2} \pi / \ln(2)}{log2(y)+2} \approx y \frac{2.26618007091}{log2(y)+2} \quad    (4) \)

下表是y取某些值时,使用公式(4)得到的近似值和真实值的对比。我们可以发现,当y>10,这个近似公式可以精确到3位以上有效数字。

y A(y,1) 公式(4)
10 4.25040709493 4.25819370444
\(2^4\) 6.03865834042 6.04314685577
\(2^6\) 18.1285335082 18.1294405673
\(2^8\) 58.0140204356 58.0142098154
\(2^{16}\) 8250.90983997 8250.90984041
\(2^{24}\) 1462315.09787 1462315.09787
\(2^{32}\) 286269685.042 286269685.042

参考文献
[1] agm(24, 6) at WolframAlpha. http://www.wolframalpha.com/input/?i=agm%2824%2C+6%29
[2] Muller J M. Elementary Functions: Algorithms and Implementation[M]. Springer Science+Business Media New York, 2016,130-131

对数常数:log(2)

本文译自 Xavier Gourdon 和 Pascal Sebah 的网站,原文http://numbers.computation.free.fr/Constants/Log2/log2.html

 

log 2 = 0.69314718055994530941723212145817656807550013436025…

通过减少劳动量,使天文学家的寿命翻了一番
– 皮埃尔·西蒙·德拉普拉斯 (1749-1827),关于对数

你不知道,对数表计算有多少诗意!
-卡尔·弗里德里希·高斯(1777-1855),给他的学生.

1 引言

对数的故事开始于苏格兰人约翰·内皮尔(1550-1617)在1614年出版的一本着作中. 这篇论文是拉丁文的,题目是“  Mirifici logarithmorum canonis descriptio[7]. 不过需要指出的是,瑞士钟表匠约斯特·比尔吉(JostBürgi,1552-1632)独立发明了对数, 但直到1620年,他的作品一直没有发表。.

由于可以分别用加法和减法来代替令人痛苦的乘法和除法,这个发明在欧洲大陆得到了非同寻常的欢迎并迅速传播(特别是得益于像开普勒这样的天文学家的热情)。在内皮尔三百周年纪念展上,乔治·吉布森写道:“对数的发明标志着科学史上的一个时代。

很快,对数表就被发布了. 其中之一是由英国数学家亨利·布里格斯(Henry Briggs)(1561-1631)于1624年出版了他的作品“ 算术对数 ”(Arithmetica logarithmica) [4]。这份对数表精确到十四位数字,并且包含了1-20,000以及9万到101,000的常用对数。. 1628年, 荷兰人Adrian Vlacq(1600-1667)完成了剩余的空缺部分[12]

数学家更喜欢使用一个所谓自然对数或双曲对数(表示为log或ln,即以e = 2.7182818284 为底的 对数,下面的定义允许简单推导出 对数的主要性质。

定义1 令x> 0

\(log ~x = \int_1^x \frac{dt}{t}=\int_1^{x-1} \frac{dt}{1+t}\)

它可被解释为双曲线y = 1 / t下的面积 ,t从1到x. 这个几何解释由耶稣会圣文森(1584-1667))在1647年所展示。

因此,log 2将被定义

\(log~ 2 = \int_1^2 \frac{dt}{t}=\int_0^{1} \frac{dt}{1+t}= \int_0^{1/2} \frac{dt}{1-t}\)

定理2 log 2 是无理数和超越数。

2 基于积分表示的算法

从log 2的积分表示,我们可以推导出一些公式。第一个是从积分的矩形近似(这里h =(b – a)/ n) 推导出来的

\(\displaystyle \int_a^b f(t)dt =h \sum_{k=1}^n f(a+kh) + O(h)\)

对区间[0,1] 应用f(t)= 1 /(1 + t)给出

\(\displaystyle log~ 2 = _{n \to \infty}^{lim} ( \frac{1}{n+1} + \frac{1}{n+2}+ \cdots + \frac{1}{2n})\)

随着这个序列的不断增加。得到的值也越来越精确。

x10
= 0.6(687…),
x100
= 0.69(065…),
x1,000
= 0.69(289…).

这是一个非常缓慢的收敛,但如果我们使用梯形法,则可以改进为

\(\displaystyle  \int_a^b f(t)dt=~\frac{h}{2}(f(a)+f(b)+2~\sum_{k=1}^{n-1}f(a+kh))+O(h^2),\)

可得到

\(\displaystyle log~2 = _{n \to \infty }^{lim} \Big ( \frac{3}{4n} + ( \frac{1}{n+1} + \frac{1}{n+2} + \cdots + \frac{1}{2n-1}) \Big ) \)

某些迭代得到

x10
= 0.693(771…),
x100
= 0.6931(534…),
x1,000
= 0.693147(243…)..

其收敛速度有所改善,但仍为对数级收敛。我们现在对积分的计算作出最后的改进,即辛普森法

\( \displaystyle \int_a^b f(t)dt=\frac{h}{6} \sum_{k=0}^n \Big( f(a+kh)+f(a+(k-1)h) + 4f(a+(k-\frac{1}{2})h) \Big)+O(h^4),\)

给出

\( \displaystyle log ~2 = _{n \to \infty} ^{lim}  \sum_{k=0}^n ~ \Big( \frac{1}{n+k} + \frac{1}{n+k-1} + \frac{4}{n+k-1/2} \Big) \)

新的迭代得到

x10
= 0.693147(374…),
x100
= 0.6931471805(794…),
x1,000
= 0.69314718055994(726…).

从这些例子中我们可以注意到,这样的公式只能用来计算log 2的少数几位数,因为收敛速度太慢,以致无法得到高精度的常数。

3. 一个简单的极限

有趣的是,我们能够发现,对数可以通过一些非常简单的公式计算出来,这些公式与几何面积没有直接关系,只涉及开平方。

我们知道常数e可以由着名的极限定义来得到

\(e ~= _{n \to \infty}^{lim} ( 1+ \frac{1}{n}) ^ n,\)

值得注意的是,存在一个计算对数常数的类似的公式:

\(log ~ 2 ~ = ~ _{n \to \infty}^{lim} n( 2^{1/n}-1) ,\)

这显然是从那个扩展中推演出来的

\(2 ^{1/n} ~ = ~ e^{ log ~2/n}=1+\frac{log~2}{n} + O(\frac{1}{n^2}).\)

取n = 1000,我们发现

x1,000 = 0.693(387…).

一个更好的公式给出

\(log ~ 2 = ~ _{n \to \infty } ^\lim \frac{n}{2}(2^{1/n}-2^{-1/n}),\)

再次,当n = 1,000

x1,000 = 0.693147(236…),

这个误差是\(O(\frac{1} { n ^2})\)。总结本节,我们给出了赛德尔在1871年发布的不寻常的无限乘积公式

\( \displaystyle \frac { log ~ x }{x-1}= \prod _{k=1} ^{k= \infty} \frac {2} {x^\frac{1}{2^k}+1} \)

我们将它应用于x = 2

\( log ~ 2 = \frac{2}{2^{1/2}+1} \frac{2}{2^{1/4}+1} \frac{2}{2^{1/8}+1} \cdots \)

4 泰勒展开

1668年,尼古拉斯·墨卡托(1620-1687)在他的对数函数[6]中发表 了着名的对数函数级数展开式:

定理3 (墨卡托)

\(\displaystyle log (1+x)= x – \frac{x^2}{2} + \frac{x^3}{3} + \cdots= \sum _{k \geq1} \frac{(-1)^{k-1}x^k}{k}. \quad -1<x \leq 1 \)

墨卡托级数与詹姆斯•格雷戈里(James Gregory,1638-1675)研究的另一个级数非常相似:

\( arctan ~ x = x – \frac{x^3}{3} + \frac{x^5}{5} – \cdots \)

第一个使用墨卡托级数的数学家之一是伊莎克·牛顿(Isaac Newton,1643-1727),他在1671年写了着名的“De Methodis Serierum et
Fluxionum”(De Methodis Serierum et Fluxionum),此书包含了各种对数计算,但直到1736年,John Colson才将这本书翻译成英文并发表 。.
他对较小的x值(如± 0.1,1± 0.2)进行了精确的计算,然后使用log 2 = 2log(1.2)- log(0.9)- log(0.8)之类的关系计算大数的对数。.

对x = 1, 应用墨卡托级数会导致着名且慢的收敛

\( log ~ 2 = 1 – \frac{1}{2} + \frac{1}{3} – \frac{1}{4} + \frac{1}{5}- \cdots \)

分别取10,100,1000和10000项得到:

x10
= 0.6(456…)
x100
= 0.6(881…)
x1,000
= 0.69(264…)
x10,000
= 0.693(097…).

一个好得多的想法是应用x = – 1/2 , 这产生了一个经常用于数值计算的公式:

\( \displaystyle log ~2 =  \sum _{k \geq1} \frac{1}{k2^k} \qquad  (1)  \)

这个公式的收敛速度是几何级的(每获得1位数字,需要p次迭代)。

x1
= 0.(5)
x10
= 0.693(064…)
x100
= 0.6931471805599453094172321214581(688…)

应用这个公式,取k = 1000时,log 2的数字超过300位。 获得d位十进制数字所需的迭代次数k由下式给出(在前面的公式中N = 2)

\( \frac{1}{10^d}=\frac{1}{kN^k} \)

这给出了k的上限值

\( k \approx \frac{d}{log10(N)} \)

如果我们使用平凡的分解2 = (4/3) (3/2), 并在等号两端取对数,我们找到一个新的公式

\( \displaystyle log ~ 2= \sum_{k \geq 1} \frac{1}{k}  \Big( \frac{1}{3^k}+\frac{1}{4^k}\Big) \)

当K达到1000时,这个级数的部分总和将会超过470位有效数字。

从等式log 2 = – 1 / 2log(1 – 1/4)+ arctanh(1/2)(参见下一段关于arctanh的定义),我们推导出公式

\(  \displaystyle log ~2 = \sum _{k \geq 0} \Big ( \frac{1}{8k+8} + \frac{1}{4k+2} \Big ) \frac{1}{4^k}, \)

\(  \displaystyle log ~2 = \frac{2}{3} + \sum _{k \geq 1} \Big ( \frac{1}{2k} + \frac{1}{4k+1} + \frac{1}{8k+4} + \frac{1}{16k+12}\Big ) \frac{1}{16^k}. \)

这些公式还可以用来直接计算 log 2 的第n个二进制数字,而不用计算以前的数字。

5 Machin类公式

使用墨卡托级数已经大大提高了我们计算log2至很多位数字的能力,但是对于数p来说,可以走的更远。我们记得反双曲正切的定义。我们记
得反双曲正切的定义.

定义4 令| x | < 1:

\(  \displaystyle arctanh~x = \frac{1}{2} log \Big(\frac{1+x}{1-x} \Big) = \sum _{k \geq 0} \frac {x^{2k+1}}{2k+1}. \)

对于x = 1/3,这很容易导出与公式 \( \pi = arctan(1) \) 类似的基本公式log 2 = 2 arctanh(1/3)。我们将其重写为:

\( \displaystyle log ~2 = \frac{2}{3} \sum _{ k \geq 0 } \frac{1}{(2k+1)9^k} \)

下面是两个迭代:

x1
= 0.6(666…),
x10
= 0.6931471805(498…).

将公式应用到k = 1000时,log 2的有效数字超过了950位。. 可以做得更好吗?有一个著名的\( \pi \)的 公式,这个公式被称为梅钦的公式,由著名的
关系式子给出:

\(  \frac{\pi}{4}= 4~ arctan(\frac{1}{5}) – arctan(\frac{1}{239}).\)

这个想法是寻找计算log2的关系式. 如果我们注意 \( (a_1,a_2,\cdots,a_n)\)这一组有理数和 \( (p_1,p_2,\cdots, p_n)\)这样一组不断增加的整数, 这个关系式是下面这个样子的。

\( \displaystyle  log~2 =  \sum _{i=1} ^n a_i ~ arctanh(\frac{1}{p_i})\)

我们也希望公式是高效的。一个好的关系式是一个少的项数(n必须小)和大的\(p_i\)的折中。按照莱默的做法,我们通过下面的公式定义效率E.
定义5 (莱默度量)。 令\(p_i > 0\),  是n个元素的数组,则:

\( \displaystyle  E =  \sum _{i=1} ^n \frac{1}{log10({p_i}^2)} \)

例如log 2 = 2 arctanh(1/3),这里(n=1,p= 3),效率度量是1.05。使用同样的效率定义,Machin公式的效率度量为0.926(n=2, p1=5,p2=239), 所以它比log 2的公式效率更高。

可以推出如下分解定理:

定理6 令p>1 :

\( arctanh ( \frac{1} {p}) = arctanh (\frac{1}{2p-1}) + arctanh(\frac{1}{2p+1}) \)
\( arctanh ( \frac{1} {p}) = 2arctanh (\frac{1}{2p-1}) – arctanh(\frac{1}{2p^2-1}) \)
\( arctanh ( \frac{1} {p}) = 2arctanh (\frac{1}{2p+11}) + arctanh(\frac{1}{2p^2-1}) \)
\( arctanh ( \frac{1} {p}) = 2arctanh (\frac{1}{2p}) + arctanh(\frac{1}{4p^3-3p}) \)

对p = 3 应用这个定理给出了n = 2的关系式

推论7

\( log~2 = 2 arctanh(\frac{1}{5}) + 2 arctanh(\frac{1}{7}) \) Euler 1748   (2)
\( log~2 = 4 arctanh(\frac{1}{5}) – 2 arctanh(\frac{1}{17}) \)                     (3)
\( log~2 = 4 arctanh(\frac{1}{7}) + 2 arctanh(\frac{1}{17}) \)                    (4)
\( log~2 = 4 arctanh(\frac{1}{6}) + 2 arctanh(\frac{1}{99}) \)                    (5)

以上3个公式的效率分别为E = 1.307,E = 1.121,E = 0.998,E = 0.893。 使用其他分解公式和计算机计算,可以找到许多其他的这种性质的关系式(在[8]中给出)。为了说明这一点,我们选择另外几个公式,1个三项的和2个四项的。

定理8 (1997)。常数log 2可以通过以下任何1个快速收敛的级数获得:

\( log~2 = 18 arctanh(\frac{1}{26}) – 2 arctanh(\frac{1}{4801}) + 8 arctanh(\frac{1}{8749})  \)       (6)
\( log~2 = 144 arctanh(\frac{1}{251}) + 54 arctanh(\frac{1}{449}) – 38 arctanh(\frac{1}{4801}) +  62 arctanh(\frac{1}{8749})  \)      (7)
\( log~2 = 72 arctanh(\frac{1}{127}) + 54 arctanh(\frac{1}{449}) + 34 arctanh(\frac{1}{4801}) -10 arctanh(\frac{1}{8749}) \)       (8)

以上3个公式的效率分别是(E = 0.616,E = 0.659,E = 0.689)。最后两个公式用于计算log 2的\(10^8\)位。一个用于计算,另一个用于验证,正如你所看到的,只有5个arctanh的计算是必要的,因为有2项arctanh被重复使用。由于关系(6),几年后的记录增加到了\(5 \times 10^8\)以上。

请注意,还可以分别执行每一项的计算,使这些关系式的计算可以并行化。

例如,我们给出第二个公式的第一次迭代(包含251的那个),每次迭代会增加精度4.8位:

x1
= 0.6(456…)
x2
= 0.6931471805(304…)
x3
= 0.69314718055994(497…)
x4
= 0.69314718055994530941(317…)

5.1有理数的公式

寻找与有理数有关的关系可能是方便的,其形式是同一的:

\( \displaystyle log~2 = \sum _{i=1}^n a_i ~ arctanh( \frac{m_i}{p_i}) \)

\(m_i\)是整数,而上面的公式是1,莱默度量必须考虑这一点。

定理9 (1997)

\( log ~ 2= 6 arctanh(\frac{1}{9}) + 2 arctanh(\frac{3}{253}),\)        (9)
\( log ~ 2= 10 arctanh(\frac{1}{17}) + 4 arctanh(\frac{13}{499}),\)        (10)

其效率分别为 E = 0.784和 E = 0.722。注意关系式(10)用于几个常数的高精度计算。

6 各种超几何级数

高斯研究了以下级数家族

\(  F(a,b,c,z)= 1 + \frac {a \cdot b} {c} \frac{z}{1 ! } + \frac {a(a+1) \cdot b (b+1)}{ c(c+1)} \frac{z^2}{2!} + \cdots \)

其中a,b,c是实数。该级数收敛于 | z | <1.。在圆| z | = 1时,级数在c > a + b时收敛。这些函数被称称之为超几何函数 [1]。通过使用合适的a,b,c的值,许多通常的函数可以表示为超几何函数。

例10

\( ln(1+x)= xF(1,1;2;-x), \)
\( arctanh~x = xF(\frac{1}{2},1;\frac{3}{2};x^2), \)
\( arctan~x = xF(\frac{1}{2},1;\frac{3}{2};-x^2), \)
\( \frac{1}{1-x}= F(1,1;1;x). \)

1797年,高斯的朋友约翰·弗里德里希·普法夫(Johann Friederich Pfaff,1765-1825)给出了有趣的关系式

\( F(a,b,c,z)= (1-z)^{-b} F(c-a,b,c, \frac{z}{z-1} ). \)

如果我们将这个公式应用于arctanh函数

\( arctanh~x = xF(\frac{1}{2},1;\frac{3}{2};x^2), \)

我们发现,

\( arctanh~x =\frac{x}{1-x^2} F(1,1;\frac{3}{2};\frac{x^2}{x^2-1}), \)

对这个函数,我们给出一个新的级数

\( arctanh~x = \frac{x}{1-x^2} \Big ( 1- \frac{2}{3}(\frac{x^2}{1-x^2} ) + \frac{2 \cdot 4 } {3 \cdot 5} (\frac{x^2}{1-x^2})^2 – \frac{2 \cdot 4 \cdot 6}{3 \cdot 5 \cdot 7}  ( \frac{x^2}{1-x^2})^3 + \cdots \Big ) \)

使用 log2 = 2 arctanh(1/3); 经过一些简单的操作,该级数

推论11

\( log 2 = \frac{3}{4} \Big ( 1- \frac{1}{12} + \frac{1 \cdot 2}{12 \cdot 20}- \frac{1 \cdot 2 \cdot 3}{12 \cdot 20 \cdot 28} + \cdots \Big ). \)

这个级数的每个新项都乘以- k /(8 k +4),其中k = 1,2,…欧拉给出了同样类型的公式来计算\( \pi \)。这个推论的一个很好的应用就是编写一个很小的代码来计算log 2,就计算\( \pi \)一样。

7. log2和AGM

所有先前的级数的收敛率是对数或几何的。根据AGM(算术几何平均值)(参见[2]),可以找出一个计算log 2的公式,这个公式就像\( \pi \)一样,是二次收敛的。
基于AGM的二次收敛算法中最经典的log 2如下所示。从开始一个\(a_0\)和 \(b_0>0 \),我们考虑AGM迭代

我们使用符号

\( a_{n+1}=\frac{a_n+b_n}{2}, \qquad b_{n+1}= \sqrt {a_n b_n },  \)

我们使用迭代

\( \displaystyle  R(a_0,b_0)= \frac{1} {(1- \sum _{n \geq 0} 2^{n-1}( {a_n}^2-{b_n}^2) )}  \)

定理12 令 \( N \geq 3\) 和 \( x \leq 1\)

\(  \vert log ~ x – R(1,10^{-N}) + R(1,10^{-N}) x  \vert \leq \frac{N}{10 ^{2(N-2)}}\)

通过选择x = 1/2,该算法给出log 2的约2N位十进制数​​。这个是二次收敛的,并且当 \(2^{n-1} ( {a_n}^2 – {b_n}^2 \)小于 \(10^{-2N} \) 就该停止,为 \( n \approx log_2(N) \)

请注意,这个算法不是专门用来计算log2的,对于任何实数X,都可以用来计算logX的值。. 使用FFT乘法,以获得 n位有效数字的x,它的复杂性是O(n log^2n)。

8 计算log p

从log 2的一个很好的近似开始,由于这个关系式的缘故,很容易找到所有整数的对数的表达式

\( \displaystyle 2 ~ arctanh( \frac{1}{2p+1}) = log \Big ( \frac{1+\frac{1}{2p+1}}{1-\frac{1}{2p+1} } \Big )= log( \frac{p+1}{p}), \)

因此我们推断:

定理13 令p> 0

\( log(p+1)= log ~ p + 2 \Big ( \frac{1}{2p+1} + \frac{1}{3} \frac{1}{ {(2p+1)}^3} + \frac{1}{5} \frac{1}{ {(2p+1)}^5} + \cdots \Big ) . \)

当你需要计算多个连续的整数的对数时,这可能是有用的。

实例14 用p=2,p=4 和p=7

\( log ~ 3 = log ~ 2 + \frac{2}{5} \Big ( 1+ \frac{1}{3} \frac{1}{25} + \frac{1}{5}\frac{1}{25^2} + \cdots \Big ),\)
\( log ~ 5 = 2log2 + \frac{2}{9} \Big ( 1+ \frac{1}{3} \frac{1}{81} + \frac{1}{5}\frac{1}{81^2} + \cdots \Big ), \)
\( log ~ 7 = 3log2 + \frac{2}{15} \Big ( 1+ \frac{1}{3} \frac{1}{225} + \frac{1}{5}\frac{1}{225^2} + \cdots \Big ). \)

9 log2公式的集合

log2的积分和级数公式可以在log2的公式集合中找到。

10 log 2的计算记录

 

位数 什么时候 备注
16 约1671 牛顿 他用log2 = 2log(1.2)-log(0.9)-log(0.8)和墨卡托级数对数。
25 1748 欧拉 关系式(2)被使用。欧拉还计算了整数3 到10的的自然对数,并精确到25位有效数字[5]
137 1853 Shanks 还计算了log3,log5和log10的值,并精确到137位数字。他的工作发表了一个p计算[9]
273 1878 亚当斯 亚当斯还以相同的精度计算了log3,log5和log7。
330 1940 H. S. Uhler 罗康瑞还使用相同的精度,计算出了log3,log5,log7和log17的值[11]
3 683 1962 D.W. Sweeney 在同一篇文章指出,计算log2是计算欧拉常数最多达3566位的必要条件)。 用于log2的公式是扩展(1)(见[10]
2 015 926 1997 P. Demiche 他的计算使用公式(2)这个经典方法。HP9000 / 871在160MHz下花了16个小时。
5 039 926 1997 P. Demiche 同样的技术被使用,在同一台计算机上做计算,花了130个小时
10 079 926 1997年12月 P. Demiche 相同的技术,在同一台计算机上做计算,400小时。
29 243 200 1997年12月 X. Gourdon 使用AGM公式(11),用到FFT技术。计算是在SGI R10000工作站上完成的,耗时20小时58分钟。
58 484 499 1997年12月 X. Gourdon 相同的AGM技术,在同一台计算机上83小时。
108 000 000 1998年9月 X. Gourdon 四项公式(7)与二分法一起使用。公式(8)被用于验证。在一台具有256M内存PII 450,这个算过程化
了47个小时的。
200 001 000 2001年9月 X. Gourdon 和 S. Kondo 使用公式(6)和(5)。 在一台具有256M内存的Athlon 1300,花了8个小
时。
240 000 000 2001年9月 X. Gourdon 和 S. Kondo 使用公式(6)和(10)。在一台具有1024M内存的PIV 1400上,花了14个小
时.
500 000 999 2001年9月 X. Gourdon 和 S. Kondo 使用公式(6)和(10)。在一台具有1024M内存的PIV 1400上,花了28个小
时.

参考文献
[1] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover, New York, (1964)
[2] J.M. Borwein and P.B.Borwein, “Pi and the AGM – A study in Analytic Number Theory and Computational Complexity”, A
Wiley-Interscience Publication, New York, (1987)
[3] R.P. Brent, Fast multiple-Precision evaluation of elementary functions, J. Assoc. Comput. Mach., (1976), vol. 23,
p.242-251
[4] H.Briggs, Arithmetica logarithmica sive logarithmorum Chiliades Triginta, London, (1624)
[5] L. Euler, Introduction à l analyse infinitésimale (french traduction by Labey), Barrois, aîné, Librairie, (original
1748, traduction 1796), vol. 1, p. 89-90
[6] N. Mercator, Logarithmotechnia, London, (1668)
[7] J. Napier, Mirifici logarithmorum canonis descriptio, Edimboug, (1614)
[8] P. Sebah, Machin like formulae for logarithm, Unpublished, (1997).
[9] W. Shanks, Contributions to Mathematics Comprising Chiefly the Rectification of the Circle to 607 Places of Decimals,
G. Bell, London, (1853)
[10] D.W. Sweeney, On the Computation of Euler’s Constant, Mathematics of Computation, (1963), p. 170-178
[11] H.S. Uhler, Recalculation and extension of the modulus and of the logarithms of 2, 3, 5, 7 and 17, Proc. Nat. Acad.
Sci., (1940), vol. 26, p. 205-212
[12] A. Vlacq, Arithmetica logarithmica, Gouda, (1628)

也谈圆周率计算

  1. 圆周率计算历史

[1] [2] [3] [4] 给出详尽的关于圆周率计算的历史。我们这里仅给出几个有代表性的事件。

1.1. 第一个将π计算到小数点后2位数字的是古希腊数学家阿基米德 Archimedes,阿基米德是有史以来最伟大的三个数学家之一,另外两个是牛顿和高斯。

1.2. 第一个将圆周率计算到小数点后7位以上者是中国南北朝时期的数学家祖冲之,他的记录保持了1000多年,直到1424年才被中亚数学家阿尔卡西(AI-Kashi)打破。

1.3. 第一个将圆周率计算到100位以上者,是英国数学家梅钦(John Machin) 他在1706将圆周率计算到100位。

1.4. 第一个突破500位,是威廉·香克斯( William Shanks),他在1874年,将圆周率计算至707位,但只有527位是正确的。

1.5. 圆周率的手工计算最高纪录由是英国的弗格森(DF Ferguson),1948他和美国的伦奇共同发表了π的808位小数值。

1.6. 第一个突破1000位的是DF Ferguson和John Wrench,他们使用台式计算器将pi计算至1120位,这也是计算机出现之前的最高纪录。

1.7. 第一个突破10000位,Francois Genuys, 他在1958年1月用一台IBM 704将圆周率计算到1万位。

1.8. 第一个突破1亿位的,是日本人金田康正(Yasumasa Kanada),他在1987年1月,他在一台NEC SX-2电脑上,将圆周率计算到134,214,700位。

1.9. 第一个突破1万亿位的,也是日本人金田康正,他使用的就是一种Machine公式。

2. 圆周率的算法

2.1.多边型法

在文艺复兴时代(或者说17世纪)之前,圆周率都是通过计算多边形的周长来近似得到的。这类算法很慢,祖冲之将圆周率至小数点后7位花了15年时间。

2.2 梅钦公式或类梅钦公式

这个公式是英国数学家梅钦(John Machin)提出,1706年梅钦计算π值突破100位小数大关,他利用了如下公式:

 

这个公式有2项,每项都是一个常数乘以arctan(1/p)形式的数。而arctan(1/p)的展开式是一个级数形式,可以使用普通的方法来计算,也可以使用Binary splitting方法使得计算效率更高。这个公式非常好用,甚至在2002年12月,金田康正首次将圆周率计算到1万亿位以上时,他使用的就是一种Machine公式,这个公式是挪威人Størmer提出来的,见 [3][7]

2.3 拉马努金公式和Chudnovsky公式

印度传奇数学家给出多个计算圆周率率的公式,见[10], 这里我们贴出至今仍被广为使用的那个公式。

 

楚德诺夫斯基(Chudnovsky)公式,楚德诺夫斯基兄弟是美国数学家和工程师,楚德诺夫斯基兄弟发明了一个新的公式,称之为Chudnovsky公式, 这个公式是Ramanujan公式改良的基础上改良而成的。1994年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。

Chudnovsky公式的另一个更方便于计算机编程的形式是:
2.3 AGM算法

在1976年,理查德·布兰特(Richard Peirce Brent)和尤金·萨拉明(Eugene Salamin)独立的发表了新的计算圆周率的算法。称之为布伦特-萨拉明(或萨拉明-布伦特)算法。这个算法是基于是一个迭代算法,每迭代一次,有效数字增加1倍。这个算法用到高斯-勒让德AGM迭代算法。其方法见下

初值:

重复计算:

最后计算:

3. 当代圆周率计算的几个重要人物

3.1. 金田康正(Kanada Yasumasa

是一位日本数学家,东京大学信息科学系教授,在过去30年中他的21次计算中,创造了11次世界纪录。他的名字出现在[2]中19次。

3.2. Chudnovsky兄弟

是美国数学家和工程师,他们开发了Chudnovsky算法,并多次创造了圆周率的计算纪录。在[2]中,他们的名字出现了5次。

3.3. 法布里斯•贝拉

是一位杰出的软件工程师,FFmpeg 和QEMU 软件的创建者。他用自己开发的软件在一台CPU为Core i7 CPU 2.93GHz台式机上运行了90多天,将圆周率算到小数点后27000亿位,创造当新的世界纪录。

3.4. 近藤茂(Shigeru Kondo)

是一名狂热的圆周率爱好者。他并非数学家,他总是使用别人写的程序来计算圆周率。在近几年,他使用余智恒(Alexander Yee)开发的y-cruncher 程序创创造了3次世界纪录,分别将圆周率计算到5万亿,10万亿,和12万亿位。他也曾用另一个程序PiFast创造了在PC上圆周率计算位数的世界纪录,但计算位数没有超过当时的巨型机,故在[2]中没有提及。在[3]中,他的名字出现了23次,比如在2003年9月,他在一台P4 3.2GHZ的电脑上,将圆周率250亿。

3.5.余智恒(Alexander Yee)

是一位年轻的计算机天才,在2010年,近藤茂用他的程序创造了新的世界记录时,他仅仅22岁。如今,他的程序世界上计算圆周率最快的程序。

4. 几种计算圆周率的软件

4.1.  SuperPi

是金田康正基于运行巨型机上的程序改写的,可运行在普通PC的版本。Super PI在超频社区非常流行,既可以作为测试这些系统性能的基准,也可以作为压力测试来检查它们是否仍然正常工作。该程序使用浮点FFT算法来做大数乘法,使用AGM算法来计算圆周率。相对于后来者,这个程序的速度并不是很快。我自己编写的仅仅使用Toom-4乘法和AGM算法的计算Pi的程序居然可和SuperPI打成平手。关于SuperPi,请参见[4]

4.2. PiFast

Xavier Gourdon的PiFast是2003年Microsoft Windows最快的程序。据其作者称,它可以在2.4 GHz Pentium 4上以3.5秒计算一百万位数。PiFast还可以计算其它无理数等e√2。它也可以在效率较低的情况下以很少的内存工作(低至几十兆以计算十亿位以上的数字),PiFast可使用拉马努金公式和和Chudnovsky公式来计算圆周率,比SuperPi要快很多,但他依然是一个单线程的程序,运行在32位的操作系统,不能重发发挥CPU的性能。

4.3. QuickPi

Steve Pagliarulo写的QuickPiye 也是同样优秀的计算Pi的程序,当计算精度在4亿位以下时,他的速度快于PiFast。像PiFast一样,QuickPi也可以计算其他的无理数,比如e,√2和√3。

4.4. y-cruncher

他是余智恒写的一个计算PI的程序,可运行在32位或者64位操作系统,可工作于多线程模式或者单线程模式,可充分发挥CPU的计算能力。当然,这个程序也是最复杂的。[8]提到,这个程序的源码竟包含47万行以上的代码。[9]中提到,这个程序使用了多至9种计算大数乘法的算法。

4.5. 其他计算圆周率的程序

除了上述几个程序外,还有其他一些程序同样可计算圆周率,这里就不再重点介绍了。他们有wPrime,IntelBurnTest,Prime95,Montecarlo superPI,OCCT等。

5.y-cruncher 和Pifast的性能比较

由于piFast是32位程序,仅仅运行单线程模式。所以我的测试仅Focus在单线程模式,测试计算2500万位圆周率的运行时间,时间单位为秒。测试结果1

运行环境 Win7 64bit, CPU: i7-4700HQ @ 2.40GHz, 8G RAM
程序 y-cruncher PiFast
命令

行或参数

y-cruncher custom pi -dec:25m -TD:1 -PF:none 0,0,25000000,2048,1
Final Multiply 0.213 1.51
Total time 8.039 45.55
Note: Total time不包括从2进制转化10进制的时间。以Final Multiply步骤为参照。y-cruncher的速度是PiFast的7倍

测试结果2

运行环境 Win7 32bit, CPU: i7-2600K @ 3.40GHz, 16G RAM
程序 y-cruncher PiFast
命令行或参数 y-cruncher custom pi -dec:25m -TD:1 -PF:none 0,0,25000000,2048,1
Final Multiply 0.704 1.29
Total time 22.513 40.34
Note: Total time不包括从2进制转化10进制的时间,Final Multiply步骤为参照y-cruncher的速度是PiFast的1.8倍

参考文献:

[1] https://en.wikipedia.org/wiki/Approximations_of_%CF%80 π的近似值

[2] https://en.wikipedia.org/wiki/Chronology_of_computation_of_%CF%80 π计算年表

[3] http://numbers.computation.free.fr/Constants/PiProgram/computations.html π and its computation through the ages

[4] https://baike.baidu.com/item/%E5%9C%86%E5%91%A8%E7%8E%87/139930 百度百科圆周率

[5] https://en.wikipedia.org/wiki/Super_PI

[6] https://en.wikipedia.org/wiki/Yasumasa_Kanada金田康正

[7] https://en.wikipedia.org/wiki/Carl_St%C3%B8rmer Carl stomer

[8] http://www.numberworld.org/y-cruncher/algorithms.html

[9] http://www.numberworld.org/y-cruncher/internals/multiplication.html乘法

[10] http://numbers.computation.free.fr/Constants/Pi/piSeries.html