/*
 * estimate pi using numerical integration
 * input is a convergence threshold:
 * program repeats the calculation for powers-of-2 numbers of steps
 * until it "converges" (successive iterations are "close", as defined
 * by the input threshold), or until we reach largest "int" power of 2.
 */
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
/* copied from not-strictly-standard part of math.h */
#define M_PI		3.14159265358979323846

double estimate_pi(int steps);

int main(void) {
	double threshold;
	printf("threshold?\n");
	if ((scanf("%lf", &threshold) != 1) || (threshold <= 0)) {
		printf("invalid input\n");
		return EXIT_FAILURE;
	}
	int steps = 1;
	double pi = 0.0; /* initial estimate */
	double pi_prev;
	do {
		pi_prev = pi;
		steps *= 2;
		pi = estimate_pi(steps);

	} while ((fabs(pi_prev-pi) > threshold) && (steps < (INT_MAX/2)));

	printf("results for threshold %g:\n", threshold);
	printf("number of steps = %d\n",steps);
	printf("computed pi = %g  (%17.15f)\n",pi, pi);
	printf("difference between computed pi and math.h M_PI = %g\n", 
			fabs(pi-M_PI));
	if (fabs(pi_prev-pi) > threshold) {
		printf("computation failed to converge\n");
	}
	return EXIT_SUCCESS;
}

/* returns result of computation with "steps" steps */
double estimate_pi(int steps) {
	double sum = 0;
	double step = 1.0/(double) steps;
	for (int i = 0 ; i < steps;  ++i) {
		double x = (i+0.5)*step;
		sum += 4.0/(1.0+x*x);
	}
	return step * sum;
}