pi-calculator/Main.c

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef long double Number;

static Number circle(Number x) {
    return sqrtl(1.0 - x * x);
}

static Number monteCarloIntegration(void) {
    int amount = 512 * 512;
    int in = 0;
    for(int i = 0; i < amount; i++) {
        Number x = rand() / (Number)RAND_MAX;
        Number y = rand() / (Number)RAND_MAX;
        in += (x * x + y * y < 1);
    }
    return 4.0 * in / amount;
}

static Number simpson(Number a, Number b) {
    return ((b - a) / 6.0) *
           (circle(a) + 4.0 * circle((a + b) * 0.5) + circle(b));
}

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

static Number lineSum(void) {
    Number sum = 0.0;
    const int max = 1000000;
    for(int i = 0; i < max; i++) {
        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;
}

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;
    }
    Number diffX = x2 - x1;
    Number diffY = y2 - y1;
    return sqrtl(diffX * diffX + diffY * diffY) * factor;
}

static void print(Number d) {
    Number diff = d - M_PI;
    printf("%.20Lf %.20Lf\n", d, diff < 0.0 ? -diff : diff);
}

#define STR(x) STR2(x)
#define STR2(x) #x

int main(void) {
    print(monteCarloIntegration());
    print(simpsonIntegration());
    print(lineSum());
    print((double)lineRecursion());
    puts(STR(M_PI));
    return 0;
}