#include <fstream>
#include <cmath>
#include <vector>
#include <deque>
#include <random>

using std::vector;
using std::deque;
using std::ofstream;
using std::random_device;
using std::uniform_real_distribution;
using std::uniform_int_distribution;
using std::mt19937_64;
using std::exp;
using std::abs;
using std::pow;

random_device rd;
mt19937_64 rng(rd());
uniform_real_distribution<double> dds(0,1);

void Wolff(int N, int k, vector<int> &S, vector<vector<int>> A, double p, double &m);

int main()
{
    int N = 128;//lattice length
    double T = 4.0;//given temparature
    
    double p = 1.0-exp(-2.0/T);//flip probability
    double m = 1.0; //magnetization
    int mcs_max = (int)(pow(2,20)+0.5); //max montecarlo step
    int eq_time = 100000; //equilibrium time
   
    vector<int> S(N*N, 1);//spin information array
    vector<vector<int>> A(N*N, vector<int>(4,-1));//the closest neighbor information

    for(int i = 0; i < N; ++i)
    {
        for(int j = 0; j < N; ++j)
        {
            A[i*N + j][0] = i*N + (j+1)%N; //the right particle
            A[i*N + j][1] = i*N + (j-1+N)%N; //the left particle
            A[i*N + j][2] = ((i+1)%N)*N + j; //the above particle
            A[i*N + j][3] = ((i-1+N)%N)*N + j; //the below particle
        }
    }

    ofstream fout;
    fout.open("mag.dat"); //output file

    //equlibrium step
    for(int i = 0; i < eq_time; ++i)
    {
        int loc = ids(rng);
        Wolff(N, loc, S, A, p, m);
    }

    uniform_int_distribution<int> ids(0,(N*N)-1);
    
    //measure the magnetization
    for(int i = 0; i < mcs_max; ++i)
    {
        int loc = ids(rng);
        Wolff(N, loc, S, A, p, m);
        fout << i << "," << m << "\n";
    }
    fout.close();
}

void Wolff(int N, int k, vector<int> &S, vector<vector<int>> A, double p, double &m)
{
    int nc = 1; //# of cluster
    int s0 = S[k]; //the spin information of selected particle
    int ind = -1; //the particle index
    S[k] = -s0; // flip spin

    deque<int> dq; //queue for saving same spin neighbor index
    dq.emplace_back(k); //save selected particle's index
    dq.shrink_to_fit(); //free memories which are not used
    
    while(!dq.empty())
    {
        ind = dq.front(); //extract information of the leftmost queue
        dq.pop_front(); //remove the leftmost queue

        for(int j = 0; j < A[ind].capacity(); ++j)
        {
            if(S[A[ind][j]] == s0) //compare the spins
            {
                if(dds(rng) < p) //proceed in the probability p
                {
                    dq.emplace_back(A[ind][j]);//save the neigbor's index
                    dq.shrink_to_fit();//free memories
                    S[A[ind][j]] = -s0;//flip spin
                    ++nc; //add 1 to # of cluster
                }
            }
        }
    }
    m = 0.0;
    //calculating the magnetization
    for(int i = 0; i < S.capacity(); ++i)
    {
        m += (double)S[i]/(N*N);
    }
    m = abs(m);
}