Calculate-Pi/src/formulas/ramanujan.cpp at main · zbob256/Calculate-Pi · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include "formulas/ramanujan.hpp"
#include <mpfr.h>

void ramanujan(mpfr_t pi, long long int n) {
    mpfr_prec_t prec = mpfr_get_prec(pi);

    mpfr_t sum, term, numer, denom, factor, temp;
    mpfr_t fk, f4k;

    mpfr_inits2(prec,
        sum, term, numer, denom, factor, temp,
        fk, f4k,
        (mpfr_ptr)0
    );

    mpfr_set_zero(sum, 1);

    for (long long int k = 0; k < n; k++) {
        mpfr_fac_ui(f4k, 4 * k, MPFR_RNDN);

        mpfr_fac_ui(fk, k, MPFR_RNDN);

        mpfr_set(numer, f4k, MPFR_RNDN);
        mpfr_mul_ui(numer, numer, 1103 + 26390 * k, MPFR_RNDN);

        mpfr_pow_ui(denom, fk, 4, MPFR_RNDN);

        mpfr_set_ui(temp, 396, MPFR_RNDN);
        mpfr_pow_ui(temp, temp, 4 * k, MPFR_RNDN);

        mpfr_mul(denom, denom, temp, MPFR_RNDN);

        mpfr_div(term, numer, denom, MPFR_RNDN);

        mpfr_add(sum, sum, term, MPFR_RNDN);
    }

    mpfr_sqrt_ui(factor, 2, MPFR_RNDN);
    mpfr_mul_ui(factor, factor, 2, MPFR_RNDN);
    mpfr_div_ui(factor, factor, 9801, MPFR_RNDN);

    mpfr_mul(temp, factor, sum, MPFR_RNDN);
    mpfr_ui_div(pi, 1, temp, MPFR_RNDN);

    mpfr_clears(
        sum, term, numer, denom, factor, temp,
        fk, f4k,
        (mpfr_ptr)0
    );
}