Martyr2專案實現——Number部分的問題求解 (1) Find Pi to Nth Digit
阿新 • • 發佈：2020-10-19
## 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開平方](https://blog.csdn.net/RickyIT/article/details/78051334)
**大數除法的精度問題**
在進行大數運算時，對於大數除法`BigDeciaml.divide()`，需要設定響應的計算精度和舍入方法(如何截斷數值)
這裡我們需要使用到`java.Math.MathContext`類，這個類描述了數字運算子的某些規則
我們可以使用預設的規則（比如`MathContext.DECIMAL128`)
也可以指定精度和舍入模式，定義自己的MathContext物件，構造方法為
`MathContext(int setPrecision, RoundingMode setRoundingMode)`
具體用法參考部落格：[java_math_MathContext](https://www.cnblogs.com/zjushuiping/archive/2012/05/31/2528212.html)
為了能夠實現我們的計算要求（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 |
#### 程式實現：
##### 主程式
```java
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);
}
}
```
##### 計算平方根程式：
```java
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(i