Topic: Calculate pi

The following code estimates the value of pi by estimating the area of the upper right quadrant of a circle of radius R whose center is at (0,0) and using the known area of pi*R*R/4.  It estimates the area by dividing the quadrant into horizontal slices, each of which has a top edge that is shorter than the bottom edge.  The program first takes y to be from 0 to R by 1.  This is essentially calculating the area of each slice to be a rectangle whose width is the bottom edge.  Since this is the longest width of the true area, this overestimates the true area of each slice.  The program then takes y to be from 1 to R by 1.  This is essentially calculating the area of each slice to be a rectangle whose width is the top edge.  Since this is the shortest width of the true area, this underestimates the true area of each slice.  The program then takes y to be from 1 to R by 2.  This is essentially calculating the area of each slice to be a rectangle whose width is the width of a line in the middle of the slice.  Hence, although this calulates the area of half as many slices, the improved estimate of each slice causes the estimate of the entire area to be better than either of the prior two.  Finally, the program then takes y to be from 1 to 2*R by 2.  This is again calculating the area of each slice to be a rectangle whose width is the width of a line in the middle of the slice but is calulating the area of as many slices as before, further improving the estimate.

#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
using namespace std;

double calcpi1(long R, long start) {
    long R2 = R*R;
    double A = 0.0;
    for (long y = start; y < R; y++) {
        double x = sqrt(R2 - y*y);
        A += x;
    }
    return 4 * A / R2;
}
double calcpi2(long R, long start) {
    long R2 = R*R;
    double A = 0.0;
    for (long y = start; y < R; y += 2) {
        double x = sqrt(R2 - y*y);
        A += x;
    }
    return 8 * A / R2;
}
int main(int argc, char *argv[]) {
    long r = 1000;
    if (argc > 1) r = atol(argv[1]);
    printf("R = %ld\n", r);
    double pi0 = calcpi1(r, 0);
    double pi1 = calcpi1(r, 1);
    double pi2 = calcpi2(r, 1);
    double pi3 = calcpi2(r*2, 1);
    double pi = M_PI;
    printf("R = 0 to %4ld by 1: pi = %-12.10g  %g %% diff\n", r-1, pi0, 100 * (pi0 - pi) / pi);
    printf("R = 1 to %4ld by 1: pi = %-12.10g  %g %% diff\n", r-1, pi1, 100 * (pi1 - pi) / pi);
    printf("R = 1 to %4ld by 2: pi = %-12.10g  %g %% diff\n", r-1, pi2, 100 * (pi2 - pi) / pi);
    printf("R = 1 to %4ld by 2: pi = %-12.10g  %g %% diff\n", r*2-1, pi3, 100 * (pi3 - pi) / pi);
    printf("      Actual value: pi = %-12.10g\n", pi);
    return 0;
}

Following is the output from the program:

R = 1000
R = 0 to  999 by 1: pi = 3.143555467   0.0624783 % diff
R = 1 to  999 by 1: pi = 3.139555467   -0.0648457 % diff
R = 1 to  999 by 2: pi = 3.141623457   0.000980497 % diff
R = 1 to 1999 by 2: pi = 3.141603545   0.000346682 % diff
      Actual value: pi = 3.141592654

As shown above, the program uses a method that overestimates pi by 0.062 percent, then one that underestimates it by 0.064 percent, then one that overestimates it by only 0.00098 percent, and then a final one that overestimates by just 0.00035 percent.  All methods calculate the area of R slices except the third which only calculates the area of R/2 slices.