|
@@ -2,89 +2,83 @@
|
|
|
#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;
|
|
|
-}
|
|
|
+typedef long double Number;
|
|
|
|
|
|
-double circle(double x) {
|
|
|
- return root(1 - x * x, 25);
|
|
|
+static Number circle(Number x) {
|
|
|
+ return sqrtl(1.0 - x * x);
|
|
|
}
|
|
|
|
|
|
-double monteCarloIntegration() {
|
|
|
+static Number monteCarloIntegration(void) {
|
|
|
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;
|
|
|
+ Number x = rand() / (Number)RAND_MAX;
|
|
|
+ Number y = rand() / (Number)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));
|
|
|
+static Number simpson(Number a, Number b) {
|
|
|
+ return ((b - a) / 6.0) *
|
|
|
+ (circle(a) + 4.0 * circle((a + b) * 0.5) + circle(b));
|
|
|
}
|
|
|
|
|
|
-double simpsonIntegration() {
|
|
|
- double sum = 0.0;
|
|
|
+static Number simpsonIntegration(void) {
|
|
|
+ Number 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;
|
|
|
+static Number lineSum(void) {
|
|
|
+ Number 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);
|
|
|
+ Number x1 = (Number)i / max;
|
|
|
+ Number x2 = (i + 1.0) / max;
|
|
|
+ Number diffX = x2 - x1;
|
|
|
+ Number diffY = circle(x2) - circle(x1);
|
|
|
+ sum += sqrtl(diffX * diffX + diffY * diffY);
|
|
|
}
|
|
|
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);
|
|
|
+static Number lineRecursion(void) {
|
|
|
+ Number x1 = 0.0;
|
|
|
+ Number y1 = 1.0;
|
|
|
+ Number x2 = 1.0;
|
|
|
+ Number y2 = 0.0;
|
|
|
+ Number factor = 2.0;
|
|
|
+ for(int i = 0; i < 26; i++) {
|
|
|
+ Number midX = (x1 + x2) * 0.5;
|
|
|
+ Number midY = (y1 + y2) * 0.5;
|
|
|
+ Number length = sqrtl(midX * midX + midY * midY);
|
|
|
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;
|
|
|
+ Number diffX = x2 - x1;
|
|
|
+ Number diffY = y2 - y1;
|
|
|
+ return sqrtl(diffX * diffX + diffY * diffY) * factor;
|
|
|
}
|
|
|
|
|
|
-void print(double d) {
|
|
|
- double diff = d - M_PI;
|
|
|
- printf("%.20lf %.20lf\n", d, diff < 0.0 ? -diff : diff);
|
|
|
+static void print(Number d) {
|
|
|
+ Number diff = d - M_PI;
|
|
|
+ printf("%.20Lf %.20Lf\n", d, diff < 0.0 ? -diff : diff);
|
|
|
}
|
|
|
|
|
|
-int main() {
|
|
|
+#define STR(x) STR2(x)
|
|
|
+#define STR2(x) #x
|
|
|
+
|
|
|
+int main(void) {
|
|
|
print(monteCarloIntegration());
|
|
|
print(simpsonIntegration());
|
|
|
print(lineSum());
|
|
|
- print(lineRecursion());
|
|
|
- puts("3.14159265358979323846\n");
|
|
|
+ print((double)lineRecursion());
|
|
|
+ puts(STR(M_PI));
|
|
|
return 0;
|
|
|
-}
|
|
|
+}
|