pi-calculator/Main.c

#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;
}