/*
 * Program to repeatedly estimate pi using numerical integration,
 * two ways.  (One way does the sum going from left to right and
 * the other from right to left.)
 *
 * Input is a file containing numbers of "slices".
 *
 * Output is a file containing one line per input value, with the
 *   input value and two computed values representing how close the
 *   two methods come to the best-available value of pi (so, each
 *   value will the absolute value of the difference between the 
 *   estimate and M_PI).
 *
 * Names for the two files are given as command-line arguments.
 */
#include <stdio.h>
#include <math.h>

/* best(?) available value of pi, to check computed approximation */
/* surprisingly, no library-defined constant for this in standard C! */
#define M_PI acos(-1.0)

double curve(double x);
double estimate_lr(long num_slices);
double estimate_rl(long num_slices);

int main(int argc, char *argv[]) {

    if (argc < 3) {
        printf("usage %s infile outfile\n", argv[0]);
        return 1;
    }
    FILE * infile = fopen(argv[1], "r");
    if (infile == NULL) {
        printf("cannot open input file\n");
        return 1;
    }
    FILE * outfile = fopen(argv[2], "w");
    if (outfile == NULL) {
        printf("cannot open output file\n");
        return 1;
    }

    long num_slices;
    int rval = 0;
    while ((fscanf(infile, "%ld\n", &num_slices) == 1) && (rval == 0)) {
        if (num_slices > 0) {
            fprintf(outfile, "%ld %22.20g %22.20g\n", 
                    num_slices,
                    fabs(estimate_lr(num_slices) - M_PI),
                    fabs(estimate_rl(num_slices) - M_PI));
        }
        else {
            printf("input values must be positive\n");
            rval = 1;
        }
    }
    if ((rval == 0) && !feof(infile)) {
        printf("input values must be integers\n");
        rval = 1;
    }
    fclose(infile);
    fclose(outfile);
    return rval;
}

double curve(double x) {
	return 4 / (1 + x*x) ;
}

double estimate_lr(long num_slices) {

	double work = 0.0; /* sum of rectangle heights */
	double width = 1.0/num_slices;

	for (long i = 0; i < num_slices; ++i) {
		work += curve(width*(i+0.5));
	}
	return work * width;
}

double estimate_rl(long num_slices) {

	double work = 0.0; /* sum of rectangle heights */
	double width = 1.0/num_slices;

	for (long i = num_slices; i > 0; --i) {
		work += curve(width*(i-0.5));
	}
	return work * width;
}