Java program to generate PI digits using Gauss-Legendre iterative algorithm using double and BigDecimal.
After 3 iterations as expected it generates 18 correct decimal digits. After 25 iterations it should generate 45 million of correct digits (assuming that the scale is set to be big enough).
3.141592646213542838751209274050779640674591064453125 3.1415926462135422821493444319826957743144372233456027945595394848214347672207952646946434489179913066
Main.java
import java.math.BigDecimal; import java.math.MathContext; import java.math.RoundingMode; // ? 2019 TheFlyingKeyboard and released under MIT License // theflyingkeyboard.net public class Main { private static final int scale = 100; public static void main(String[] args) { long iterations = 3; System.out.println(new BigDecimal(gaussLegendre(iterations))); System.out.println(gaussLegendreBigDecimal(iterations)); } private static double gaussLegendre(long iterations) { double a = 1.0d; double b = 1.0d / Math.sqrt(2); double t = 0.25d; double p = 1.0d; double an1; double bn1; double tn1; double pn1; double pi = 0.0d; for (long i = 0; i < iterations; ++i) { an1 = (a + b) / 2.0d; bn1 = Math.sqrt(a * b); tn1 = t - p * Math.pow(a - an1, 2); pn1 = 2 * p; pi = Math.pow(a + b, 2) / (4 * t); a = an1; b = bn1; t = tn1; p = pn1; } return pi; } private static BigDecimal gaussLegendreBigDecimal(long iterations) { MathContext mc = new MathContext(scale); BigDecimal two = new BigDecimal(2); BigDecimal a = new BigDecimal(1.0d); BigDecimal b = new BigDecimal(1.0d ).divide(new BigDecimal(2).sqrt(mc), scale, RoundingMode.CEILING); BigDecimal t = new BigDecimal(0.25d); BigDecimal p = new BigDecimal(1.0d); BigDecimal an1; BigDecimal bn1; BigDecimal tn1; BigDecimal pn1; BigDecimal pi = new BigDecimal(0); for (long i = 0; i < iterations; ++i) { an1 = a.add(b, mc).divide(two, scale, RoundingMode.CEILING); bn1 = a.multiply(b, mc).sqrt(mc); tn1 = t.subtract(p.multiply(a.subtract(an1, mc).pow(2, mc), mc), mc); pn1 = p.multiply(two, mc); pi = a.add(b, mc).pow(2, mc).divide(t.multiply(new BigDecimal(4), mc), scale, RoundingMode.CEILING); a = an1; b = bn1; t = tn1; p = pn1; } return pi; } }