Martyr2專案實現——Number部分的問題求解 (1) Find Pi to Nth Digit

sunfulv發表於2020-10-18
Git

Martyr2專案實現——Number部分的問題求解 (1) Find Pi to Nth Digit

Find Pi to Nth Digit

問題描述：

Find PI to the Nth Digit – Enter a number and have the program generate PI up to that many decimal places. Keep a limit to how far the program will go.

翻譯：

給定一個整數N，讓程式生成精確到小數點後N為的圓周率\(\pi\)

要保證程式執行的時間在一定限度下

計算原理:

常用的圓周率的數值計算方法有級數法，迭代法，隨機演算法

級數法：使用圓周率\(\pi\)的級數表示來計算

  1. 高斯提出的用於平方倒數和公式

    \[\frac{\pi^2}{6}=\frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}...+\frac{1}{(n-1)^2}+\frac{1}{n^2}+... \]

  2. 萊布尼茲公式

    \[\frac{\pi}{4}=\frac{1}{1}-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+...+(-1)^{n-1}\frac{1}{2n-1}+... \]

    不過萊布尼茲公式的收斂速度很慢

  3. 拉馬努金提出的公式

    \[\frac{1}{\pi}=\frac{2\sqrt{2}}{9801}\sum_{k=0}^{+\infty}\frac{(4k)!(1103+26390k)}{k!^4(396^{4k})} \]

使用級數法計算圓周率的收斂速度還是太慢

迭代演算法：適合計算機程式實現的計算圓周率的方法

​ 迭代演算法的收斂速度要比無窮級數快很多

​ 比較出名的演算法是高斯-勒讓德演算法

​ 高斯-勒讓德演算法：

​ 引入四個數列 \(\{a_n\},\{b_n\},\{t_n\},\{p_n\}\)

​ 他們的初值為：

\[a_0=1\qquad b_0=\frac{1}{\sqrt{2}}\qquad t_0=\frac{1}{4}\qquad p_0=1 \]

​ 遞推公式為：

\[a_{n+1}=\frac{a_n+b_n}{2}\;,b_{n+1}=\sqrt{a_nb_n}\;,t_{n+1}=t_n-p_n(a_n-a_{n+1})^2\;,p_{n+1}=2p_n. \]

​ 計算圓周率\(\pi\)近似值的方法：

\[\pi \approx\frac{(a_{n+1}+b_{n+1})^2}{4t_{n+1}} \]

​ 該演算法每執行一次迭代，計算出的圓周率的正確位數就會增加一倍多。

具體的實現：

我們準備將圓周率計算到小數點後1000位（N<=1000)

開方運算

考慮到java的浮點數最高只支援64位double雙精度浮點數，為了能夠計算的更精確，考慮使用java的大數類·java.Math.BigDecimal來進行計算。

注意到在使用高斯勒讓德演算法計算圓周率時，需要用到開平方運算，BigDecimal並沒有實現對大數物件的開方運算，我們需要自己實現。這裡使用牛頓迭代法來計算大數的開平方。

具體的計算方法參考部落格：java BigDecimal開平方

大數除法的精度問題

在進行大數運算時，對於大數除法BigDeciaml.divide()，需要設定響應的計算精度和舍入方法(如何截斷數值)

這裡我們需要使用到java.Math.MathContext類，這個類描述了數字運算子的某些規則

我們可以使用預設的規則（比如MathContext.DECIMAL128)

也可以指定精度和舍入模式，定義自己的MathContext物件，構造方法為

MathContext(int setPrecision, RoundingMode setRoundingMode)

具體用法參考部落格：java_math_MathContext

為了能夠實現我們的計算要求（1000位的圓周率），我們設定大數除法的計算精度為1002位（有效數字，自定義舍入方法）

MathContext mc = new MathContext(1002, RoundingMode.HALF_EVEN);

對於開平方運算，我們設定它的計算精度為500位（精確到小數點後100位）

下表是我們計算的每次迭代可以到達的計算精度

對於給定的引數N（要求計算小數位數）,我們通過查表來確定迭代次數，然後對得到的數值進行截斷。

迭代次數 精度（小數點後精確到的位數）
0 0
1 2
2 7
3 18
4 40
5 83
6 170
7 344
8 693
9 1000

程式實現：

主程式
import java.math.BigDecimal;
import java.math.MathContext;

public class CalculatePi {
    private static int[] map_array = {0,2,7,18,40,83,170,344,693,1000};

    public static String getPiValue(int N){ //獲取精確到小數點後N位的圓周率近似值
        if(N<0||N>1000) return "error:給定引數超出範圍！(預設引數範圍為[1,1000])";
        int index = 0;
        for(int i=map_array.length-1;i>=1;i--){
            if(N>map_array[i]) {
                index = i+1;
                break; //給定的引數N位於map_array[i,i+1]之間
            }
        }
        String value = calculate(N,index);
        return value;
    }

    private static String calculate(int N,int index) {
        //利用高斯-勒讓德迭代演算法來計算圓周率的近似值，index為迭代的次數
        if (index == 0) return "3";
        //設定初值
        BigDecimal a0 = new BigDecimal(1);
        BigDecimal a1 = new BigDecimal(1);
        BigDecimal b = CalculateSqrt.sqrt(new BigDecimal("0.5"));
        BigDecimal t = new BigDecimal("0.25");
        BigDecimal p = new BigDecimal(1);
        BigDecimal pi = new BigDecimal(3);
        MathContext mc = CalculateSqrt.mc;
        //進行迭代
        for (int i = 0; i < index; i++) {
            a1 = a0.add(b);
            a1 = a1.divide(new BigDecimal(2), mc);
            b = b.multiply(a0);
            b = CalculateSqrt.sqrt(b);
            BigDecimal temp = new BigDecimal(1);
            temp = a0.subtract(a1);
            temp = temp.multiply(temp);
            temp = temp.multiply(p);
            t = t.subtract(temp);
            p = p.multiply(new BigDecimal(2));
            temp = a1.add(b);
            temp = temp.multiply(temp);
            temp = temp.divide(new BigDecimal(4), mc);
            pi = temp.divide(t, mc);
            a0 = a1;
        }
        return pi.toString().substring(0, N + 2);
    }



    public static void main(String[] args) {
        int N = 10;
        String pi = getPiValue(1001);
        System.out.println(pi);
    }
}
計算平方根程式：
import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;

public class CalculateSqrt {
    private static int N = 1002;
    public static MathContext mc = new MathContext(N, RoundingMode.HALF_EVEN);
    private static String eps = "0."+repeatString("0",N/2)+"1";
    public static void main(String[] args) {
        BigDecimal n = new BigDecimal("2");
        BigDecimal r = sqrt(n);
        System.out.println(r.toString());
    }

    public static BigDecimal sqrt(BigDecimal num) {
        if(num.compareTo(BigDecimal.ZERO) < 0) {
            return BigDecimal.ZERO;
        }
        BigDecimal x = num.divide(new BigDecimal("2"), mc);
        while(x.subtract(x = sqrtIteration(x, num)).abs().compareTo(new BigDecimal(eps)) > 0);
        return x;
    }

    private static BigDecimal sqrtIteration(BigDecimal x, BigDecimal n) {
        return x.add(n.divide(x, mc)).divide(new BigDecimal("2"), mc);
    }
    private static String repeatString(String str,int n){
        StringBuffer sb = new StringBuffer();
        for(int i=0;i<n;i++){
            sb.append(str);
        }
        return sb.substring(0,sb.length());
    }
}

​

相關文章