|
@@ -0,0 +1,90 @@
|
|
|
+#include <math.h>
|
|
|
+#include <stdio.h>
|
|
|
+#include <stdlib.h>
|
|
|
+
|
|
|
+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;
|
|
|
+}
|