「Java案例」求PI的值

案例解析

莱布尼茨级数法

π的近似值可以通过莱布尼茨级数公式计算：
π/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - …
这个公式虽然简单，但收敛速度较慢，先用它来理解基本原理。

# 源文件保存为“CalculatePI.java”
public class CalculatePI {public static void main(String[] args) {int terms = 1000000; // 计算的项数double pi = 0.0;for (int i = 0; i < terms; i++) {double term = 1.0 / (2 * i + 1);if (i % 2 == 0) {pi += term;} else {pi -= term;}}pi *= 4; // 因为计算的是π/4，所以需要乘以4System.out.println("计算得到的π值: " + pi);System.out.println("Math.PI实际值: " + Math.PI);System.out.println("误差: " + (Math.PI - pi));}
}

运行结果

计算得到的π值: 3.1415916535897743
Math.PI实际值: 3.141592653589793
误差: 1.0000000187915248E-6

代码解析：

• 设置计算项数为100万项，项数越多结果越精确
• 使用for循环逐项计算
• 根据项数的奇偶性决定加减
• 循环结束后乘以4得到π的近似值
• 最后输出计算结果并与Java内置的Math.PI比较

运行这个程序，可以发现虽然计算了100万项，但精度仍然只有小数点后5位左右，说明这个方法收敛确实很慢。

蒙特卡洛方法

现在尝试另一种有趣的方法——蒙特卡洛模拟。这个方法通过随机撒点来估算π值。

import java.util.Random;public class MonteCarloPI {public static void main(String[] args) {Random random = new Random();int totalPoints = 1000000; // 总点数int insideCircle = 0; // 落在圆内的点数for (int i = 0; i < totalPoints; i++) {double x = random.nextDouble(); // 0-1之间的随机x坐标double y = random.nextDouble(); // 0-1之间的随机y坐标// 判断点是否在单位圆内if (x * x + y * y <= 1) {insideCircle++;}}// 计算π值：面积比 = π/4double pi = 4.0 * insideCircle / totalPoints;System.out.println("蒙特卡洛法计算的π值: " + pi);}
}

运行结果

蒙特卡洛法计算的π值: 3.143272

代码解析：

• 在边长为1的正方形内随机撒点
• 计算落在半径为1的四分之一圆内的点数
• 根据面积比(圆面积/正方形面积=π/4)估算π值
• 点数越多，结果越精确

这个方法虽然直观有趣，但精度同样有限，适合理解概念而非高精度计算。

马青公式

为了提高计算效率，现在来看一个收敛更快的算法——马青公式。

# 源文件保存为“MachinPI.java”
import java.math.BigDecimal;
import java.math.RoundingMode;public class MachinPI {public static void main(String[] args) {int decimalPlaces = 15; // 小数位数BigDecimal pi = calculatePI(decimalPlaces);System.out.println("计算到小数点后" + decimalPlaces + "位的π值: " + pi);}private static BigDecimal calculatePI(int decimalPlaces) {BigDecimal sixteen = new BigDecimal(16);BigDecimal four = new BigDecimal(4);// 计算arctan(1/5)的级数展开BigDecimal arcTan1_5 = arctan(5, decimalPlaces);// 计算arctan(1/239)的级数展开BigDecimal arcTan1_239 = arctan(239, decimalPlaces);// 应用马青公式: π/4 = 4arctan(1/5) - arctan(1/239)BigDecimal pi = sixteen.multiply(arcTan1_5).subtract(four.multiply(arcTan1_239));return pi.setScale(decimalPlaces, RoundingMode.HALF_UP);}private static BigDecimal arctan(int denominator, int decimalPlaces) {BigDecimal result = BigDecimal.ZERO;BigDecimal term;int n = 0;// 设置计算精度int precision = decimalPlaces + 5;// 计算arctan(1/x)的泰勒级数展开BigDecimal x = new BigDecimal(denominator);BigDecimal xSquared = x.multiply(x);BigDecimal divisor;do {divisor = new BigDecimal(2 * n + 1);term = BigDecimal.ONE.divide(x.pow(2 * n + 1).multiply(divisor), precision, RoundingMode.HALF_UP);if (n % 2 == 0) {result = result.add(term);} else {result = result.subtract(term);}n++;} while (term.compareTo(new BigDecimal("1E-" + precision)) > 0);return result;}
}

运行结果

计算到小数点后15位的π值: 3.141592653589793

代码解析：

• 使用BigDecimal保证高精度计算
• 马青公式基于arctan函数的泰勒展开
• 通过设置小数位数控制计算精度
• 当新增项小于指定精度时停止计算
• 采用高精度计算模式避免舍入误差

这个方法收敛速度快得多，计算15位小数只需很少的迭代次数。

实践练习题

改进莱布尼茨级数

观察莱布尼茨级数的收敛情况，修改程序使其在达到指定精度时自动停止计算。

参考代码：

# 源文件保存为“LeibnizImproved.java”
public class LeibnizImproved {public static void main(String[] args) {double precision = 1e-6; // 期望精度double pi = 0.0;int i = 0;double term;do {term = 1.0 / (2 * i + 1);if (i % 2 == 0) {pi += term;} else {pi -= term;}i++;} while (Math.abs(term) > precision);pi *= 4;System.out.println("计算到精度" + precision + "的π值: " + pi);System.out.println("实际迭代次数: " + i);}
}

运行结果

计算到精度1.0E-6的π值: 3.1415946535856922
实际迭代次数: 500001

比较不同算法的效率

编写程序比较上述三种方法在相同精度下的计算效率，统计各自的迭代次数和耗时。

参考代码：

# 源文件保存为“CompareMethods.java”
import java.util.Random;
import java.math.BigDecimal;
import java.math.RoundingMode;class LeibnizImproved {public static void leibniz(double targetPrecision) {double precision = 1e-6; // 期望精度double pi = 0.0;int i = 0;double term;do {term = 1.0 / (2 * i + 1);if (i % 2 == 0) {pi += term;} else {pi -= term;}i++;} while (Math.abs(term) > precision);pi *= 4;System.out.println("计算到精度" + precision + "的π值: " + pi);System.out.println("实际迭代次数: " + i);}
}class MonteCarloPI {public static void monteCarlo(int num) {Random random = new Random();int totalPoints = 1000000; // 总点数int insideCircle = 0; // 落在圆内的点数for (int i = 0; i < totalPoints; i++) {double x = random.nextDouble(); // 0-1之间的随机x坐标double y = random.nextDouble(); // 0-1之间的随机y坐标// 判断点是否在单位圆内if (x * x + y * y <= 1) {insideCircle++;}}// 计算π值：面积比 = π/4double pi = 4.0 * insideCircle / totalPoints;System.out.println("蒙特卡洛法计算的π值: " + pi);}
}class MachinPI {public static BigDecimal calculatePI(int decimalPlaces) {BigDecimal sixteen = new BigDecimal(16);BigDecimal four = new BigDecimal(4);// 计算arctan(1/5)的级数展开BigDecimal arcTan1_5 = arctan(5, decimalPlaces);// 计算arctan(1/239)的级数展开BigDecimal arcTan1_239 = arctan(239, decimalPlaces);// 应用马青公式: π/4 = 4arctan(1/5) - arctan(1/239)BigDecimal pi = sixteen.multiply(arcTan1_5).subtract(four.multiply(arcTan1_239));return pi.setScale(decimalPlaces, RoundingMode.HALF_UP);}private static BigDecimal arctan(int denominator, int decimalPlaces) {BigDecimal result = BigDecimal.ZERO;BigDecimal term;int n = 0;// 设置计算精度int precision = decimalPlaces + 5;// 计算arctan(1/x)的泰勒级数展开BigDecimal x = new BigDecimal(denominator);BigDecimal xSquared = x.multiply(x);BigDecimal divisor;do {divisor = new BigDecimal(2 * n + 1);term = BigDecimal.ONE.divide(x.pow(2 * n + 1).multiply(divisor),precision,RoundingMode.HALF_UP);if (n % 2 == 0) {result = result.add(term);} else {result = result.subtract(term);}n++;} while (term.compareTo(new BigDecimal("1E-" + precision)) > 0);return result;}
}public class CompareMethods {public static void main(String[] args) {double targetPrecision = 1e-6;// 测试莱布尼茨级数法long start = System.nanoTime();LeibnizImproved.leibniz(targetPrecision);long end = System.nanoTime();System.out.println("莱布尼茨法耗时: " + (end - start) / 1e6 + "ms");// 测试蒙特卡洛法start = System.nanoTime();MonteCarloPI.monteCarlo((int) 1e6);end = System.nanoTime();System.out.println("蒙特卡洛法耗时: " + (end - start) / 1e6 + "ms");// 测试马青公式start = System.nanoTime();MachinPI.calculatePI(6);end = System.nanoTime();System.out.println("马青公式耗时: " + (end - start) / 1e6 + "ms");}
}

运行结果

计算到精度1.0E-6的π值: 3.1415946535856922
实际迭代次数: 500001
莱布尼茨法耗时: 32.9133ms
蒙特卡洛法计算的π值: 3.14054
蒙特卡洛法耗时: 64.2519ms
马青公式耗时: 1.8079ms