Bellard's formula always generates 3.14506

I've been working on a simple C++ program to approximate pi with Bellard's formula. No matter what number of digits it tries calculating to, it results in the number 3.145063.

I tried changing it from storing as a float to a double, because I thought it might just be rounding to fit in the smaller data type, but it still isn't working.

Here's my code:

#include <iostream>
#include <cmath>

double approximate(int digits) {
    double summed = 0;
    for (int n = 0; n < digits; n++) {
        summed += (pow(-1, n) / pow(2, n * 10)) * ((-(pow(2, 5)) / (4 * n + 1)) - (1 / (4 * n + 3)) + (pow(2, 8) / (10 * n + 1)) - (pow(2, 6) / (10 * n + 3)) - (pow(2, 2) / (10 * n + 5)) - (pow(2, 2) / (10 * n + 7)) + (1 / (10 * n + 9)));
    }
    return (1 / pow(2, 6)) * summed;
}

int main() {
    std::cout << "How many digits?: ";
    std::string its;
    std::getline(std::cin, its);
    std::string approx = std::to_string(approximate(std::stoi(its)));
    std::cout << approx << std::endl;
}

Solution

As mentioned, if you had broken the line up into several statements, you would have seen that you are performing integer division instead of floating point division.

(1 / (4 * n + 3)) and (1 / (10 * n + 9))

are not correct. The simple correction is to use 1.0 to force a floating point division to occur:

#include <iostream>
#include <cmath>

double approximate(int digits) {
    double summed = 0;
    for (int n = 0; n < digits; n++) {
        summed += (pow(-1, n) / pow(2, n * 10)) * 
                  ((-(pow(2, 5)) / (4 * n + 1)) - (1.0 / (4 * n + 3)) + // <-- 1.0
                  (pow(2, 8) / (10 * n + 1)) - 
                  (pow(2, 6) / (10 * n + 3)) - 
                  (pow(2, 2) / (10 * n + 5)) - 
                  (pow(2, 2) / (10 * n + 7)) + (1.0 / (10 * n + 9))); // <-- 1.0
    }
    return (1 / pow(2, 6)) * summed;
}

int main() {
    std::string approx = std::to_string(approximate(10));
    std::cout << approx << std::endl;
}

Output:

3.141593

The other things in your code that can also be improved, more specifically, are the calls to pow. You know that pow(2,2), pow(2,5), pow(2,8), pow(2,6), (1.0 / pow(2,6)), etc. are constants. They do not be need to be computed each and every time the forloop iterates, or at thereturn` point. Just make these floating point constants and use them.

Edit:

Here is a compile-time version of your code, which can calculate to a maximum of 6 digits of precision.

Note that the assembly output doesn't have any trace of the approximate function call anywhere, as the entire value was computed at compile-time.