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
| // g++ -std=c++0x -o pi pi.cpp -DDIGITS=5 && ./pi
#if DIGITS < 1
# error DIGITS < 1
#endif
#include <iostream>
#include <iomanip>
#include <array>
constexpr double make_target_prec(int i) {
return (i > 1) ? 10.0 * make_target_prec(i-1) : 10.0;
}
const int PART_SIZE = 60;
const int TOTAL_RECURSION = PART_SIZE * PART_SIZE * PART_SIZE;
const double PI = 3.141592654;
const double TARGET_PREC = make_target_prec(DIGITS);
const double TARGET = (((int)(PI*TARGET_PREC))/TARGET_PREC);
#define PI_STEP(x) ((((x) % 2) == 0 ? 4.0: -4.0) / (2.0*(x)+1.0))
constexpr double pi_part2(int s, int c) {
return
(s <= 0 || c == PART_SIZE) ?
0.0 :
PI_STEP(s) + pi_part2(s-1, c+1);
}
constexpr double pi_part1(int s, int c) {
return
(s <= 0 || c == PART_SIZE) ?
0.0 :
pi_part2(s, 0) + pi_part1(s-PART_SIZE, c+1);
}
constexpr double pi(int s) {
return
(s <= 0) ?
4.0 :
pi_part1(s, 0) + pi(s-PART_SIZE*PART_SIZE);
}
constexpr int find_part2(int s, int c, double p) {
return
(c == PART_SIZE) ?
-1 :
(((int)((p+PI_STEP(s))*TARGET_PREC))/TARGET_PREC) == TARGET ?
s :
find_part2(s+1, c+1, p+PI_STEP(s));
}
constexpr int find_part1(int s, int c, double p) {
return
(c == PART_SIZE) ?
-1 :
find_part2(s, 0, p) != -1 ?
find_part2(s, 0, p) :
find_part1(s+PART_SIZE,
c+1,
((s == 0) ? pi(0): 0.0) + p+pi_part2(s+PART_SIZE-1, 1));
}
constexpr int find_k(int s=0, int c=0, double p=0.0) {
return
(c == PART_SIZE) ?
-1 :
find_part1(s, 0, p) != -1 ?
find_part1(s, 0, p) :
find_k(s+PART_SIZE*PART_SIZE,
c+1,
((s == 0) ? pi(0): 0.0) + p+pi_part1(s+PART_SIZE*PART_SIZE-1, 0));
}
int main() {
constexpr int k = find_k();
// être sur d'avoir une constante de compilation
std::array<int, ((unsigned int)(k>>30)) + 1> a;
std::cout << "k=" << k << " target=" << TARGET
<< " pi(k)=" << std::setprecision(15) << pi(k) << std::endl;
return 0;
} |
Partager