In this article, I will try to calculate pi in Java using the following formula.
\frac{4}{\pi}=\sum_{n=0}^\infty \frac{(-1)^n(4n)!(1123+21460n)}{882^{2n+1}(4^nn!)^4}
This is the official version of Ramanujan (1887/12/22 --1920/4/26). Looking at it sternly, it makes me wonder how I came up with this, and how this really gives me pi.
(I don't know which formula the exchange is about) The following exchanges are also mysterious.
Professor Hardy "Please tell me how you came up with this formula." Ramanujan "Please believe me, Goddess Namagiri stood at the bedside and taught me."
I want to actually calculate it.
Basically, it is a program without any ingenuity.
--The numerator and denominator are calculated by BigInteger respectively. --Big Decimal calculates from the point of dividing a fraction (in this program, the precision is specified for the first time here). —— Factorial is too terrible to do the same calculation over and over, so we are “developing” to reuse the previous calculation results. --The last digit contains the error (the error is proportional to the "number of repetitions").
PamaPi.java
/**
*Calculate pi.
*/
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.RoundingMode;
public class RamaPi {
/**
*Calculate factorial(Reuse the previous calculation result)。
*/
class Facto {
/**Last argument*/
int n = 1;
/**Last value*/
BigInteger val = BigInteger.ONE;
BigInteger calc(final int arg) {
for ( ; n <= arg ; n++) {
val = val.multiply(BigInteger.valueOf(n));
}
return val;
}
}
/**Factorial calculation 1*/
Facto workFacto1 = new Facto();
/**Factorial calculation 2*/
Facto workFacto2 = new Facto();
/**Constant "882"*/
static final BigInteger VAL_882 = BigInteger.valueOf(882);
/**Constant "4"*/
static final BigInteger VAL_4 = BigInteger.valueOf(4);
/**
*Calculate the formula numerator of Ramanujan.
*/
BigInteger bunsi(final int n) {
final int sign = (n&1) == 0 ? 1 : -1;
final BigInteger a = workFacto1.calc(4 * n);
final BigInteger b = BigInteger.valueOf(sign * (1123 + 21460*n));
return a.multiply(b);
}
/**
*Calculate the official denominator of Ramanujan.
*/
BigInteger bunbo(final int n) {
final BigInteger a = VAL_882.pow(2*n+1);
final BigInteger b = VAL_4.pow(n);
final BigInteger c = workFacto2.calc(n);
final BigInteger d = b.multiply(c).pow(4);
return a.multiply(d);
}
/**
*Find the pi from Ramanujan's formula.
*/
void calc(final int digits) {
BigDecimal sum = BigDecimal.ZERO;
/**Number of digits to discontinue addition (less than this number of digits cannot be added)*/
final BigDecimal thr = BigDecimal.ONE.divide(BigDecimal.TEN.pow(digits));
/**Time to start calculation*/
long t0 = System.nanoTime();
int iter = 0;
for ( ;; iter++ ) {
final BigDecimal bunbo = new BigDecimal(bunbo(iter));
final BigDecimal bunsi = new BigDecimal(bunsi(iter));
final BigDecimal adder = bunsi.divide(bunbo, digits, RoundingMode.HALF_EVEN);
if (adder.abs().compareTo(thr) < 0) {
break;
}
/*Rewrite sum with the result of addition*/
sum = sum.add(adder);
/*If the number of digits is small, the state of convergence is displayed.*/
if (digits < 100) {
System.out.println( BigDecimal.valueOf(4).divide(sum, digits, RoundingMode.HALF_EVEN));
}
}
/**Time when the calculation was completed*/
long t1 = System.nanoTime();
System.err.println("Number of repetitions: " + iter);
double elapsed = (t1 - t0) / ((double)1000 * 1000 * 1000);
System.err.println("Number of seconds required: " + String.format("%.6f", elapsed));
System.out.println(BigDecimal.valueOf(4).divide(sum, digits, RoundingMode.HALF_EVEN));
}
static public void main(String arg[]) {
int digits = Integer.parseInt(arg[0]);
new RamaPi().calc(digits);
}
}
Specify the number of digits in the argument. There may be an error in the last few digits (the error is proportional to the "number of repetitions").
java RamaPi 50
3.14158504007123775601068566340160284951024042742654
3.14159265359762201426827517920229156338712376594028
3.14159265358979322988767708881201158250039662639857
3.14159265358979323846265313702144396820307036308333
3.14159265358979323846264338326814215966511742779251
3.14159265358979323846264338327950289764449222972749
3.14159265358979323846264338327950288419715329619277
3.14159265358979323846264338327950288419716939939455
3.14159265358979323846264338327950288419716939937511
Number of repetitions: 9
Number of seconds required: 0.003002
3.14159265358979323846264338327950288419716939937511
I measured a little in the environment at hand. It takes a long time, but do you think Java's multiple length calculation is fast or slow?
Number of digits | Number of repetitions | Number of seconds required | Μs per time |
---|---|---|---|
10 | 2 | 0.001129 | 564.5 |
100 | 18 | 0.003077 | 170.9 |
1000 | 170 | 0.049033 | 288.4 |
10000 | 1698 | 4.237159 | 2495.3 |
20000 | 3395 | 19.83075 | 5841.1 |
40000 | 6790 | 118.651706 | 17474.4 |
When I tried to run "Super π" in the same environment, "8.38 million digits" was "152 seconds" (200 times the number of digits, which is different).
If you want more speed, use something for multiple precision operations such as GMP (The GNU Multiple Precision Arithmetic Library) instead of Java's BigInteger / BigDecimal. It's definitely better to do it.
Initially, we had a naive factorial implementation, but it took only 20,000 digits and 113 seconds. (At first, I wrote it recursively, but it immediately stack overflowed :hushed :)
.java
/**
*Calculate the factorial.
*/
static final BigInteger facto(final int a) {
BigInteger work = BigInteger.ONE;
for (int i=2 ; i<=a ; i++) {
work = work.multiply(BigInteger.valueOf(i));
}
return work;
}
[Srinivasa Ramanujan (Wikipedia)](https://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A5%E3%83%AA%E3%83%8B%E3 % 83% B4% E3% 82% A1% E3% 83% BC% E3% 82% B5% E3% 83% BB% E3% 83% A9% E3% 83% 9E% E3% 83% 8C% E3% 82 % B8% E3% 83% A3% E3% 83% B3)
[Arbitrary precision calculation](https://ja.wikipedia.org/wiki/%E4%BB%BB%E6%84%8F%E7%B2%BE%E5%BA%A6%E6%BC%94% Introductory article about E7% AE% 97) (Wikipedia)
If you want to perform more proper multiple length calculation, GMP
Recommended Posts