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

#define RADIUS  3000
#define HIT_ZONE  ((double) RADIUS * RADIUS)
#define HIT_ZONE_FIXPOINT  ((RADIUS * RADIUS) << 2)

typedef double (*calculate_pi)(void);


double calculate_pi_a(void)
{
    uint16_t i, j;
    double x, y;
    uint32_t hit_count = 0;

    for (i = 0; i < RADIUS; ++i) {
        x = i + 0.5;
        x *= x;
        for (j = 0; j < RADIUS; ++j) {
            y = j + 0.5;
            y *= y;
            if (x + y <= HIT_ZONE) {
                ++hit_count;
            }
        }
    }
    return 4 * hit_count / HIT_ZONE;
}


double calculate_pi_b(void)
{
    uint16_t i, j;
    double x, *values;
    uint32_t hit_count = 0;

    values = malloc(RADIUS * sizeof(*values));
    for (i = 0; i < RADIUS; ++i) {
        x = i + 0.5;
        values[i] = x * x;
    }

    for (i = 0; i < RADIUS; ++i) {
        for (j = 0; j < RADIUS; ++j) {
            if (values[i] + values[j] <= HIT_ZONE) {
                ++hit_count;
            }
        }
    }

    free(values);
    return 4 * hit_count / HIT_ZONE;
}


double calculate_pi_c(void)
{
    uint16_t i, j;
    double x, *values;
    uint32_t hit_count = 0;

    values = malloc(RADIUS * sizeof(*values));
    for (i = 0; i < RADIUS; ++i) {
        x = i + 0.5;
        values[i] = x * x;
    }

    for (i = 0; i < RADIUS; ++i) {
        for (j = i + 1; j < RADIUS; ++j) {
            if (values[i] + values[j] <= HIT_ZONE) {
                ++hit_count;
            }
        }
    }
    hit_count *= 2;

    for (i = 0; i < RADIUS; ++i) {
        if (values[i] * 2 <= HIT_ZONE) {
            ++hit_count;
        }
    }

    free(values);
    return 4 * hit_count / HIT_ZONE;
}


double calculate_pi_d(void)
{
    uint16_t i, j;
    uint32_t *values;   /* Fixed point values with two bits fractional part. */
    uint32_t hit_count = 0;

    values = malloc(RADIUS * sizeof(*values));
    /* 
     * `j` is a fixed point value with one bit fractional part here.
     * So it starts at 0.5 and is incremented by 1.0 each loop iteration.
     */
    for (i = 0, j = 1; i < RADIUS; ++i, j += 2) {
        values[i] = j * j;
    }

    for (i = 0; i < RADIUS; ++i) {
        for (j = i + 1; j < RADIUS; ++j) {
            if (values[i] + values[j] <= HIT_ZONE_FIXPOINT) {
                ++hit_count;
            }
        }
    }
    hit_count *= 2;

    for (i = 0; i < RADIUS; ++i) {
        if (values[i] * 2 <= HIT_ZONE_FIXPOINT) {
            ++hit_count;
        }
    }

    free(values);
    return 4.0 * hit_count / (HIT_ZONE_FIXPOINT >> 2);
}


int main(int argc, char *argv[])
{
    uint i;
    calculate_pi functions[] = {
        calculate_pi_a, calculate_pi_b, calculate_pi_c, calculate_pi_d
    };
    uint function_count = sizeof(functions) / sizeof(calculate_pi);

    if (argc > 1 && sscanf(argv[1], "%u", &i) == 1 && i < function_count) {
        printf("method %c: %f\n", 'a' + i, functions[i]());
    } else {
        printf("usage: %s [0-%u]\n", argv[0], function_count - 1);
    }
    return 0;
}