#include <cmath>
#include <map>
#include <limits>
#include <algorithm>
#include <vector>
#include <random>

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

std::vector<double> DATA;

void setup_data(double avg, double sigma, size_t count)
{
	const int seed = 0; // 再現性のため固定
	std::mt19937 rand_src(seed);
	std::normal_distribution<> nd(avg, sigma);

	DATA.clear();

	double cur = 1.0;
	for (size_t i=0;i<count;i++) {
		double nxt = cur * nd(rand_src);
		DATA.push_back(nxt);
		cur = nxt;
	}
}

double find_cdf_pos(const std::vector<double> &data, double cdf_target) {
	size_t max = data.size();
	size_t idx = static_cast<size_t>( std::floor((max-1) * cdf_target ) );
	if (idx == (max-1)) {
		return data[idx];
	}
	double p = data[idx];
	double n = data[idx+1];

	double ip_range = 1.0 / (max-1);
	double frac = cdf_target - (static_cast<double>(idx)/(max-1));
	return (p*(ip_range-frac)+n*frac)/ip_range; // 線形補間
}

void print_stat(int interval) {

	std::vector<double> work_data;

	size_t max = DATA.size() - interval;
	double sum = 0.0;
	double sum_sq = 0.0;
	for (size_t i=0; i<max; ++i) {
		double ratio = DATA[i+interval] / DATA[i];
		work_data.push_back(ratio);
		sum += ratio;
		sum_sq += ratio*ratio;
	}
	double avg = sum / max;
	double sq_avg = sum_sq / max;
	double sigma = std::sqrt(sq_avg - (avg*avg));

	std::sort(work_data.begin(), work_data.end());
	double low_rank_16 = find_cdf_pos(work_data, 0.16);
	double high_rank_16 = find_cdf_pos(work_data, 1-0.16);
	
	printf("%d\t%f\t%f\t%f\t%f\n", interval, sigma, high_rank_16, avg, low_rank_16);
}

int main(int argc, char **argv) {

	setup_data(1.006512, 0.052017, 12*1000);

	printf("interval\tsigma\thigh_rank_16\taverage\tlow_rank_16\n");
	for (int interval=1; interval<12; ++interval ) {
		print_stat(interval);
	}

	for (int interval=12; interval<=240; interval+=12 ) {
		print_stat(interval);
	}

	return EXIT_SUCCESS;
}