/*
 * Copyright ©️ 2024 Mario Forzanini <mf@marioforzanini.com>
 *
 * This file is part of my exercises for LSN.
 *
 * My exercises for LSN are free software: you can redistribute them
 * and/or modify them under the terms of the GNU General Public License
 * as published by the Free Software Foundation, either version 3 of
 * the License, or (at your option) any later version.
 *
 * My exercises for LSN are distributed in the hope that they will be
 * useful, but WITHOUT ANY WARRANTY; without even the implied
 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with my exercises for LSN. If not, see
 * <https://www.gnu.org/licenses/>.
 *
 */

#include <errno.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>

#include "random.h"
#include "std.h"

/* Total number of "throws" of the MonteCarlo algorithm. */
#define TOTAL 10000
/* Argument to select uniform sampling */
#define UNIFORM_NAME "UNIFORM"
/* Argument to select importance sampling */
#define IMPORTANCE_NAME "IMPORTANCE"

enum method {
	METHOD_UNKNOWN,
	METHOD_UNIFORM,
	METHOD_IMPORTANCE,
};

const double M_PI_SQUARE = M_PI * M_PI;

#define M_PI_2 (M_PI / 2.0)

double
integrand(double x)
{
	return M_PI_2 * cos(M_PI_2 * x);
}

double
d(double x)
{
	return 3. / 2. * (1 - x * x);
}

double
sample(Random rnd[static 1])
{
	double x, r;
	do {
	        x = random_rannyu(rnd);
	        r = random_rannyu(rnd);
	} while (r >= 1 - x * x);
	return r;
}

void
usage(char *argv0)
{
	fprintf(stderr, "USAGE: %s [UNIFORM|IMPORTANCE]\n", argv0);
	exit(1);
}

/* Numbers of blocks to be tested. */
static int Ns[] = { 5, 10, 20, 50, 100, 200 };

int
main(int argc, char *argv[])
{
	FILE *primes, *input;
	Random rnd = { 0 };
	if (!(primes = fopen("Primes", "r"))) {
		STD_ERROR("Error opening Primes: %s\n", strerror(errno));
		return 1;
	} else if (!(input = fopen("seed.in", "r"))) {
		STD_ERROR("Error opening seed.in: %s\n", strerror(errno));
		return 1;
	}
	if (!random_init(&rnd, primes, input)) {
		return 1;
	}
	enum method method = METHOD_UNKNOWN;
	if (argc < 2) {
		usage(argv[0]);
	} else {
		if (strcmp(argv[1], UNIFORM_NAME) == 0) {
			method = METHOD_UNIFORM;
		} else if (strcmp(argv[1], IMPORTANCE_NAME) == 0) {
			method = METHOD_IMPORTANCE;
		}
	}
	for (size_t i = 0; i < LEN(Ns); i++) {
		int nblock = Ns[i];
		int block_steps = TOTAL / nblock;
		double avg = 0.0, avg2 = 0.0, sum = 0.0;
		if (method == METHOD_IMPORTANCE) {
			for (int block = 0; block < nblock; block++) {
				for (int step = 0; step < block_steps;
				        step++) {
					double x = sample(&rnd);
					sum += integrand(x) / d(x);
				}
				sum /= (double)(block_steps);
				avg += sum;
				avg2 += sum * sum;
			}
		} else if (method == METHOD_UNIFORM) {
			for (int block = 0; block < nblock; block++) {
				for (int step = 0; step < block_steps;
				        step++) {
					sum += integrand(random_rannyu(&rnd));
				}
				sum /= (double)(block_steps);
				avg += sum;
				avg2 += sum * sum;
			}
		}
		avg /= (double)(nblock);
		printf("%d %lf %lf\n", nblock, avg,
		        sqrt(avg2 / (double)(nblock)-avg * avg));
	}
	random_save_seed(&rnd, "seed.out");
	return 0;
}