#include #include #include double root(double x, int steps) { if(x == 0.0) { return 0.0; } double a = x; x /= 2; for(int i = 0; i < steps; i++) { x = (x + a / x) * 0.5; } return x; } double circle(double x) { return root(1 - x * x, 25); } double monteCarloIntegration() { int amount = 512 * 512; int in = 0; for(int i = 0; i < amount; i++) { double x = rand() / (double)RAND_MAX; double y = rand() / (double)RAND_MAX; in += (x * x + y * y < 1); } return 4.0 * in / amount; } double simpson(double a, double b) { return ((b - a) / 6.0) * (circle(a) + 4 * circle((a + b) * 0.5) + circle(b)); } double simpsonIntegration() { double sum = 0.0; for(int i = 0; i < 100; i++) { sum += simpson(i / 100.0, (i + 1) / 100.0); } return 4.0 * sum; } double lineSum() { double sum = 0.0; const int max = 1000000; for(int i = 0; i < max; i++) { double x1 = (double)i / max; double x2 = (i + 1.0) / max; double diffX = x2 - x1; double diffY = circle(x2) - circle(x1); sum += root(diffX * diffX + diffY * diffY, 30); } return sum * 2.0; } double lineRecursion() { double x1 = 0.0; double y1 = 1.0; double x2 = 1.0; double y2 = 0.0; double factor = 2.0; for(int i = 0; i < 24; i++) { double midX = (x1 + x2) * 0.5; double midY = (y1 + y2) * 0.5; double length = root(midX * midX + midY * midY, 10); midX /= length; midY /= length; x2 = midX; y2 = midY; factor *= 2; } double diffX = x2 - x1; double diffY = y2 - y1; return root(diffX * diffX + diffY * diffY, 200) * factor; } void print(double d) { double diff = d - M_PI; printf("%.20lf %.20lf\n", d, diff < 0.0 ? -diff : diff); } int main() { print(monteCarloIntegration()); print(simpsonIntegration()); print(lineSum()); print(lineRecursion()); puts("3.14159265358979323846\n"); return 0; }