#pragma once

#include <array>
#include <cmath>
#include <list>
#include <vector>
#include <fstream>
#include <iomanip>

template <class T>
class Quantity {
private:
  uint64_t N;
  uint64_t n;
  unsigned skip;
  std::list<T> hist;
  double total;
  double total2;
  std::vector<double> C;

public:
  Quantity(unsigned lag, unsigned s = 1) : C(lag) {
    skip = s;
    N = 0;
    n = 0;
    total = 0;
    total2 = 0;
  }

  void read(std::string filename) {
    std::ifstream inFile(filename);
    inFile >> n;
    inFile >> total;
    inFile >> total2;

    for (double& Ci : C) {
      inFile >> Ci;
    }
  }

  void write(std::string filename) const {
    std::ofstream outFile(filename);
    outFile << std::setprecision(15) << n << " " << total << " " << total2 << std::endl;
    for (double Ci : C) {
      outFile << Ci << " ";
    }
    outFile << std::endl;
    outFile.close();
  }

  void reset() {
    total = 0;
    total2 = 0;
    std::fill(C.begin(), C.end(), 0);
    n = 0;
    hist = {};
  }

  void add(const T& x) {
    if (N % skip == 0) {
      hist.push_front(x);
      if (hist.size() > C.size()) {
        hist.pop_back();
        unsigned t = 0;
        for (T a : hist) {
          C[t] += a * x;
          t++;
        }
        double norm = x * x;
        total += sqrt(norm);
        total2 += norm;
        n++;
      }
    }
    N++;
  }

  double avg() const { return total / n; }

  double avg2() const { return total2 / n; }

  std::vector<double> ρ() const {
    double C0 = C.front() / n;
    double avg2 = pow(total / n, 2);

    std::vector<double> ρtmp;

    for (double Ct : C) {
      ρtmp.push_back((Ct / n - avg2) / (C0 - avg2));
    }

    return ρtmp;
  }

  std::array<double, 2> τ() const {
    double τtmp = 0.5;
    unsigned M = 1;
    double c = 8.0;

    std::vector<double> ρ_tmp = this->ρ();

    while (c * τtmp > M && M < C.size()) {
      τtmp += ρ_tmp[M];
      M++;
    }

    return {skip * τtmp, skip * 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / n};
  }

  double σ() const {
    return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
  }

  double serr() const { return sqrt(this->σ()); }

  unsigned num_added() const { return n; }
};