numberphile-codes/0735-new-pi.cpp at main · vinaydes/numberphile-codes · 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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#include <cstdio>
#include <cstdint>
#include <gmp.h>

constexpr int NUM_BITS = 2048;

const char ref_pi_str[] = "3.141592653589793238462643383279502884197169399375"
                          "10582097494459230781640628620899862803482534211706"
                          "79821480865132823066470938446095505822317253594081"
                          "28481117450284102701938521105559644622948954930381"
                          "96442881097566593344612847564823378678316527120190"
                          "91456485669234603486104543266482133936072602491412"
                          "73724587006606315588174881520920962829254091715364"
                          "36789259036001133053054882046652138414695194151160"
                          "94330572703657595919530921861173819326117931051185"
                          "480744623799627495673518857527248912279381830119491";


/*
 * pi = 4 + sum( 1/n! * ( 1/n+lambda - 4/2n+1) * ((2n+1)^2/4(n+lambda) - n)_n-1
 *
 */

void calculate_nth_term(uint32_t n, mpf_t term, mpf_t lambda) {
  mpf_t t1, t2;
  mpf_init(t1);
  mpf_init(t2);
  mpf_set_ui(t1, 1);

  for (uint32_t i = 1; i <= n; i++) {
    mpf_mul_ui(t1, t1, i);
  }
  mpf_ui_div(t1, 1, t1);
  mpf_set(term, t1);

  mpf_set(t1, lambda);
  mpf_add_ui(t1, t1, n);
  mpf_ui_div(t1, 1, t1);

  mpf_set_ui(t2, (2*n+1));
  mpf_ui_div(t2, 4, t2);

  mpf_sub(t1, t1, t2);

  mpf_mul(term, term, t1);

  mpf_set_ui(t1, (2*n+1)*(2*n+1));
  mpf_set(t2, lambda);
  mpf_add_ui(t2, t2, n);
  mpf_mul_ui(t2, t2, 4);
  mpf_div(t2, t1, t2);
  mpf_sub_ui(t2, t2, n);
  mpf_set_ui(t1, 1);

  for (uint32_t i = 0; i < n-1; i++) {
    mpf_mul(t1, t1, t2);
    mpf_add_ui(t2, t2, 1);
  }

  mpf_mul(term, term, t1);
  mpf_clear(t1);
  mpf_clear(t2);
}


void estimate_pi(mpf_t pi_est, mpf_t lambda, uint32_t num_terms) {
  mpf_t term;
  mpf_init(term);

  mpf_set_ui(pi_est, 4);
  for (uint32_t i = 1; i < num_terms; i++) {
    calculate_nth_term(i, term, lambda);
    mpf_add(pi_est, pi_est, term);
  }

  mpf_clear(term);
}

int main() {
  mpf_set_default_prec(NUM_BITS);
  mpf_t REF_PI;
  mpf_init(REF_PI);
  mpf_set_str(REF_PI, ref_pi_str, 10);

  mpf_t pi_est1, pi_est2;
  mpf_inits(pi_est1, pi_est2, NULL);

  mpf_t lambda, lower, upper, diff1, diff2;
  mpf_inits(lambda, lower, upper, diff1, diff2, NULL);

  mpf_set_d(lower,  2.0);
  mpf_set_d(upper,  3.0);

  estimate_pi(pi_est1, lower, 4);
  estimate_pi(pi_est2, upper, 4);

  mpf_sub(diff1, pi_est1, REF_PI);
  mpf_sub(diff2, pi_est2, REF_PI);
  gmp_printf("diff1 = %.80Fg\n", diff1);
  gmp_printf("diff2 = %.80Fg\n", diff2);

  mpf_mul(diff1, diff1, diff2);

  if (mpf_cmp_d(diff1, 0.0) > 0) {
    printf("Incorrect inital values of lower and upper\n");
    return 0;
  }

  for (int itr = 0; itr < 300; itr++) {
    mpf_add(lambda, lower, upper);
    mpf_div_ui(lambda, lambda, 2);

    estimate_pi(pi_est1, lambda, 4);
    mpf_sub(diff1, pi_est1, REF_PI);
    mpf_sub(diff2, upper, lower);
    gmp_printf("itr = %d, lambda = %.20Fg, window = %.20Fg, ", itr, lambda, diff2);
    //gmp_printf("Estimated value of Pi = %.100Ff\n", pi_est1);
    gmp_printf("diff1 = %.20Fg\n", diff1);

    if ( mpf_sgn(diff1) == mpf_sgn(upper)) {
      mpf_set(upper, lambda);
    }
    else {
      mpf_set(lower, lambda);
    }
  }
  //mpf_out_str(stdout, 10, 500, REF_PI);
  mpf_clear(REF_PI);
  mpf_clears(pi_est1, pi_est2, NULL);
  mpf_clears(lambda, lower, upper, diff1, diff2, NULL);
  return 0;
}