【ZT】算pi的几类算法引自: (转)为什么要算 PI? 计算机xxx的用途就是进行人类无法完成的复杂运算,算 PI 就是这样的运算之一。虽然算 PI 本身没有多大的实际意义,但是对于计算机爱好者来说作为一种编程的挑战,还是很有意思的。算 PI 看似简单,其实它还牵涉到一些有用的数学知识。 {dy}类算法:arctan 的级数展开
很容易想到,要得到超高精度的 PI 值,实数在计算机中必须以数组的形式进行存取,数组的大小跟所需的有效位数成正比。在这个算法中,PI 的有效位数 n 随 (2) 的求和项数线性增加。而为计算 (2) 中的每一项,需要进行超高精度实数除以小整数(52, 2392, 2k+1)的循环,循环所需次数也跟 n 成正比。所以,这个算法总的时间复杂度为 O(n2)。 这个算法的优点是简单,而且只需要进行整数运算。下面给出我写的算 PI 程序。在程序中,我采用了一些提高速度的措施:超高精度实数以数组的形式进行存取,数组元素的类型为 64 位整数(long long),每个元素储存 12 个十进制位;对 xk (x = 1/5, 1/239) 的头部和尾部的 0 的数量进行估计,只对非 0 的部分进行计算。 另外,还有许多跟 (1) 类似的式子,但不常用。例如:
第二类算法:与 1/PI 有关的级数
以上两个级数(还有其它类似形式的级数,但不常用)比起 arctan 的泰勒级数要复杂得多。虽然仍然是线性收敛,总的时间复杂度也仍然是 O(n2),但它们的收敛速度相当快, (Ramanujan) 每项可以增加 8 位有效数字, (Chudnovsky) 每项可以增加 14 位。 在这个算法中,除了要进行超高精度实数(数组形式)和小整数的运算外,还有一次超高精度实数的开方和倒数的运算,这需要用到 FFT(快速傅立叶变换),在下文叙述。 第三类算法:算术几何平均值和迭代法 算术几何平均值(Arithmetic-Geometric Mean, AGM) M(a, b) 定义如下:
然后,由椭圆积分的一系列理论(抱歉,过程我不懂)可以推导出如下公式:
根据这条公式可以制定适当的迭代算法。在迭代过程中,有效位数随迭代次数按 2 的指数增加,即每迭代一次有效位数乘 2。算法中的超高精度实数的乘、除、开方等运算需要使用 FFT,在下文叙述。综合考虑 FFT 的时间复杂度,整个算法的时间复杂度约为 O(n log(n)2)。 除了 (AGM) 以外,还有其它的迭代序列,它们具有同样的时间复杂度。例如下面的这个序列将按 4 的指数收敛到 1/PI:
FFT 如上所述,第二和第三类算法不可避免地要涉及超高精度实数(数组形式存取的多位数)的乘、除、开方等运算。多位数乘法如果按照常规方法来计算,逐位相乘然后相加,其时间复杂度将达到 O(n2)。使用 FFT 可大大减少计算量。 设有复数数组 a[k] 和 b[k] (k=0~n-1),正向和反向的离散傅立叶变换(DFT)定义如下: (i = sqrt(-1))
(3) 和 (4) 中的 (1/n) 可以放在任何一个式子中,也可以拆成 (1/sqrt(n)) 同时放在两个式子中,目的是保证正向和反向傅立叶变换以后不会相差一个因子。 当 n 的所有素因子均为小整数,尤其是当 n 为 2 的整数次幂的时候,使用适当的算法经过仔细的协调,可以避免多余的计算,使离散傅立叶变换 (3) 和 (4) 减少至 O(n log(n)) 的时间复杂度,即所谓的快速傅立叶变换(FFT)。具体的细节请查阅相关书籍。下面给出我写的一段 FFT 程序,仅供参考。另外也有已经开发的 FFT 函数库,例如 ,可以直接使用。
利用 FFT,要计算 n1 位和 n2 位的两个多位数乘法,可以这样进行:开辟两个长度为 n(n>=n1+n2,取 2m {zj0}) 的复数数组,将两个多位数从低位到高位分别填入,高位补 0。对两个数组分别进行正向傅立叶变换。将得到的两个变换后的数组的对应项相乘,然后进行反向傅立叶变换,{zh1}得到一个结果数组。由于傅立叶变换是在复数域中进行的,因此还要对结果数组进行取整和进位,才能得到最终的乘积。 值得留意的是傅立叶变换的精度问题。我们知道,在计算机中实数用单精度数或双精度数表示,它们会存在一定的误差。在计算多位数乘法时,n 往往是一个很大的数字,傅立叶变换过程中需要对数组的每一项进行求和,如何保证精度带来的误差不会因为求和而超出允许的范围?我的观点是必须使用双精度实数,而且由于统计特性,精度带来的误差在求和过程中不会很大,一般不会影响计算的正确性。如果需要保证计算的正确性,我想到两种检查方法。{dy}种是取模验算。例如,如果乘数和被乘数对 17 的模分别是 8 和 6,那么积对 17 的模就应该是 14。第二种是检查运算结果中浮点数偏离整数的{zd0}值。如果偏差只有比如 10-3 量级,我们可以认为这个尺度的乘法运算很安全;如果偏差达到 0.5,说明运算已经出错了;如果偏差达到 0.1 量级,那也比较危险,也许换个别的乘数和被乘数就溢出了。 多位数的倒数和开方可以通过牛顿迭代求根法转化为乘法运算。例如,要计算 x = 1/a ,根据牛顿迭代法令 f(x) = 1/x - a ,可以得到以下迭代序列:
要计算 x = sqrt(a) ,可以先计算 x = 1 / sqrt(a) ,令 f(x) = 1/x2 - a ,可以得到以下迭代序列:
(5) 和 (6) 均以 2 的指数收敛到所求结果。还存在其它更复杂一些的迭代序列,它们以更高的指数收敛,在此不提。不过需要提醒的是,跟 (AGM) 不同,这里 (5) 和 (6) 中的 x0 只是 1/a 和 1 / sqrt(a) 的约值,在前几次的迭代中不必进行满 n 位数的乘法运算,因而可以减少计算量。 Comments (3)
|