/*
 * 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"
#define TOTAL 10000
// Number of blocks
#define NBLOCKS 100
#define BLOCK_STEPS (TOTAL / NBLOCKS)

/*! Array in which results are stored for each block of the blocking
 *  method. */
static double blocks[NBLOCKS] = { 0 };
/*! Average to calculate the standard deviation of the mean. */
static double prog_avg[NBLOCKS] = { 0 };
/*! Squared average to calculate the standard deviation of the mean */
static double prog_avg2[NBLOCKS] = { 0 };
/*! Successive values of the variance as the number of blocks considered
 * enlarges. */
static double variances[NBLOCKS] = { 0 };

int
main(void)
{
	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;
	}
	for (int block = 0; block < NBLOCKS; block++) {
		double sum = 0.0;
		for (int i = 0; i < BLOCK_STEPS; i++) {
			sum += pow(random_rannyu(&rnd) - 0.5, 2);
		}
		blocks[block] = sum / (double)(BLOCK_STEPS);
	}
	for (int block = 0; block < NBLOCKS; block++) {
		for (int i = 0; i < block + 1; i++) {
			double Ai = blocks[i];
			prog_avg[block] += Ai;
			prog_avg2[block] += Ai * Ai;
		}
		prog_avg[block] /= block + 1;
		prog_avg2[block] /= block + 1;
	}
	for (int block = 1; block < NBLOCKS; block++) {
		variances[block] = sqrt(
		        (prog_avg2[block] - prog_avg[block] * prog_avg[block])
		        / (double)block);
	}
	for (int block = 0; block < NBLOCKS; block++) {
		printf("%lf %lf\n", prog_avg[block] - 1. / 12.,
		        variances[block]);
	}
	random_save_seed(&rnd, "seed.out");
	return 0;
}