割圆法求高精度圆周率

由 三硝基豆腐 发布

前几天看几何原本的时候看到了对祖冲之的介绍,他利用割圆法求出圆周率后七位实实在在是让我膜拜,毕竟南北朝时期的计算工具还比较落后,我甚至无法想象祖冲之是怎么手动开根的,而且精确到第七位需要从最开始计算保留精度,这实实在在需要强大的计算能力。

为了重(xue)现(xi)昔(J)日(A)辉(V)煌(A),我写了一份并不算高效的割圆法求$\pi$的程序,用 JAVA 实现的,可以在 23 秒的时间内计算圆周率后 500 位。

好吧确实很慢,但我懒得优化了,因为主要还是学习 JAVA。

首先高精度类用的是 java 的 BigDecimal 类,这个类提供了基本的四则运算、比较以及绝对值和设置精度等方法,然而却没提供开根。

(真的是莫名其妙开根这么重要的操作都没有)

开根的话基本上有两种方法,一种是牛顿迭代法,另一种便是妇孺皆知(大雾)的二分法,它们复杂度都是$\log$级别的。

牛顿迭代法代码短好实现,并且听起来高级,然而需要做除法,这会较大地拖慢速度。

不过后面经过实验,即使二分法不需要做除法,二分法似乎仍然要慢一些,当然可能是我的二分写得太丑了。

牛顿迭代法其实很简单,设$a$的平方根为$x$,则方程为$x^2 = a, x^2 - a = 0$

首先随便把$x$的初始值$x_1$设为某个数,例如$a$

根据牛顿迭代公式:

$$x_{n + 1} = x_n - \frac {f(x_n)} {f^{'}(x_n)}$$

得到开方的迭代关系式:

$$x_{n + 1} = \frac {x_n + \frac {a} {x_n}} {2}$$

然后就是割圆法。

割圆法很简单,就是把圆切成多边形后计算周长,再根据$d = 2\pi r$计算$\pi$。

首先从正六边形开始,半径为$1$的话周长便为$6$,每边长设为$L_6 = 1$。

根据几何关系,不难求出:

$$L_{2n} = \sqrt {2 - \sqrt {4 - L_n^2}}$$

所以十二边形的边长$L_{12}$约为$0.517638$

于是就能迭代求出$\pi$值了。

以下代码中的 cnt 代表精确的位数(最末一两位并不保证精度,因为有四舍五入)

代码写得真的不行,毕竟刚开始学 JAVA。。。

import java.math.*;

class PiDemo
{
    static BigDecimal _05 = new BigDecimal("0.5");
    static BigDecimal _2 = new BigDecimal("2");
    static BigDecimal _4 = new BigDecimal("4");
    static int cnt = 500;
    static MathContext mc = new MathContext(cnt, RoundingMode.HALF_UP);
    static MathContext mc2 = new MathContext(cnt << 1, RoundingMode.HALF_UP);

    static private BigDecimal bigSqrt(BigDecimal a)
    {
        BigDecimal ans, lst = a;
        while (true)
        {
            ans = _05.multiply(lst.add(a.divide(lst, mc2), mc2), mc2);
            if (ans.compareTo(lst) == 0) return ans;
            lst = ans;
        }
    }

    public static void main(String args[])
    {
        BigDecimal l = new BigDecimal("1");
        BigDecimal ans, lst = new BigDecimal("0");
        BigDecimal s = new BigDecimal("6");
        while (true)
        {
            l = bigSqrt(_2.subtract(bigSqrt(_4.subtract(l.multiply(l, mc2), mc2)), mc2));
            s = _2.multiply(s, mc);
            ans = _05.multiply(l.multiply(s, mc), mc);
            if (ans.compareTo(lst) == 0) break;
            lst = ans;
        }
        System.out.println(ans.toString());
    }
};

仅有一条评论

  1. 雷明
    雷明 · 2020-06-18 13:13

    厉害了

发表评论