#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <vector>
#include<set>

using namespace std;


// Function returns the index of a stabilizer with coordinates (x,y)
int GetStab(int x, int y, int StabXPos[], int StabYPos[], int NoV)
{
    int s = 999999;

    for (int i =0; i < NoV; i++)
    {
        if (StabXPos[i] == x)
        {
            if (StabYPos[i] == y)
            {
                    s = i;
            }
        }
    }

    return s;
}


// Function to measure star operator about vertex 'Vertex'
// Function to measure star operator about vertex 'Vertex'
bool StabilizerMeasurement(bool Errors[], int Site, int L, int xPos[], int yPos[])
{
    bool State = false;
    //cout << Site << " measures: ";

    /*
    // From each 'Site' 
    // we then use xPos[Site] and yPos[Site] to find their four/six adjacent qubits.

    //       elevenOclock ----- oneOclock
    //       /                                             \
    //     /                                                 \    
    //    /                                                    \
    // nineOclock                                      threeOclock
    //    \                                                     /
    //      \                                                  /  
    //        \                                               /
    //         sevenOclock ----------- fiveOclock

    */

    // ***  9- and 3-O-CLOCK ERROR ***
    // These two errors are adjacent, so we deal with them together



    // Find the correct site.
    // The formula for the correct site depends on the row number modulo 3.

    int ErrorSite = 2*Site;

    if(yPos[Site]%3 == 1)
    {    
        ErrorSite = ErrorSite+2;
    }


    // ***  9-O-CLOCK ERROR ***


    // Do not look at 9-o-clock if it is a four-body stabilizer on the left boundary of the main lattice.

    bool LeftOfLattice = ( xPos[Site] == 0 );
    bool WrongRow = ( yPos[Site]%3 == 2 );

    bool OnLeftBoundary = (LeftOfLattice & WrongRow);

    if( !OnLeftBoundary )
    {
        bool nineOclockError = Errors[ ErrorSite ];
        //cout << ErrorSite << ", ";
        State = State ^ nineOclockError;
    }


    // *** 3-O-CLOCK ERROR ***

    // 3-o-clock error is to the immediate right of the 9-o-clock error
     ErrorSite  = ErrorSite+1;
     
    // Do not check the 3-o-clock error if we are on a boundary tile to the right of the main lattice.

    bool RightOfLattice = 3*xPos[Site] == 3*L-2 - yPos[Site];
    WrongRow = yPos[Site]%3 == 1;

    bool RightOfBoundary = RightOfLattice & WrongRow;

    if( !RightOfBoundary )
    {
        bool threeOclockError = Errors[ ErrorSite  ];
        //cout << ErrorSite << ", ";
        State = State ^ threeOclockError;
    }


    // *** 5- and 7-O-CLOCK ERROR  ***
    // These two errors are easily dealt with together.

    // Boundary tiles on the bottom row do not have 5- and 7-o-clock errors

    bool BottomRow = ( yPos[Site] == 0 );

    if( !BottomRow )
    {

        ErrorSite = 0;

        // The correct formula for the error site depends on the row

        if(yPos[Site] %3 == 0)
            ErrorSite = (int)(2*Site - 2*(L - yPos[Site]/3 ) - 1);


        else if(yPos[Site] %3 == 1)
            ErrorSite = (int)(2*(Site + 1) - 2*(L - (yPos[Site]-1) / 3) - 1);

        else
            ErrorSite = (int)(2*Site - 2*(L - (yPos[Site] +1) / 3) - 1);



        bool fiveOclockError = Errors[ErrorSite+1];
                //cout << ErrorSite+1 << ", ";

        bool sevenOclockError = Errors[ErrorSite];
                //cout << ErrorSite << ", ";

        State = State ^ fiveOclockError ^ sevenOclockError;

    }


    // *** 11- and 1-O-CLOCK ERROR ***

    ErrorSite = 0;


    // 11-O-CLOCK ERROR SITE
    if(yPos[Site]%3 == 0)
    {
        ErrorSite=(int)(2*Site + 2*(L - yPos[Site]/3 )+1);
    }
    else if(yPos[Site]%3 == 1)
    {
        ErrorSite=(int) (2*Site + 2*(L - (yPos[Site]-1)/3 ) + 1);
    }
    else
    {
        ErrorSite=(int) (2*Site + 2*(L - (yPos[Site]-2)/3 ) - 1);
    }




    LeftOfLattice = ( xPos[Site] == 0 );
    WrongRow = ( yPos[Site]%3 == 2 );

    OnLeftBoundary = (LeftOfLattice & WrongRow);

    if( !OnLeftBoundary )
    {
        bool elevenOclockError= Errors[ErrorSite];
                //cout << ErrorSite << ", ";
        State = State^elevenOclockError;
    }


    // 1-O-CLOCK ERROR SITE
    // The site is one on from the 11-o-clock site
    ErrorSite = ErrorSite+1;


    if( !(3*xPos[Site] == 3*L-2 - yPos[Site]  &  yPos[Site] % 3 == 1 )) 
    {
        bool oneOclockError = Errors[ErrorSite];
                //cout << ErrorSite << endl;

        State = State^oneOclockError;

    }

    return State;
}

// Function to locate qubits about vertex 'Vertex'
void GetSurrQubitLocs(bool Errors[], int Site, int L, int xPos[], int yPos[], int QubXPos[], int QubYPos)
{
    bool State = false;
    //cout << Site << " measures: ";

    /*
    // From each 'Site' 
    // we then use xPos[Site] and yPos[Site] to find their four/six adjacent qubits.

    //       elevenOclock ----- oneOclock
    //       /                                             \
    //     /                                                 \    
    //    /                                                    \
    // nineOclock                                      threeOclock
    //    \                                                     /
    //      \                                                  /  
    //        \                                               /
    //         sevenOclock ----------- fiveOclock

    */

    // ***  9- and 3-O-CLOCK ERROR ***
    // These two errors are adjacent, so we deal with them together



    // Find the correct site.
    // The formula for the correct site depends on the row number modulo 3.

    int ErrorSite = 2*Site;

    if(yPos[Site]%3 == 1)
    {    
        ErrorSite = ErrorSite+2;
    }


    // ***  9-O-CLOCK ERROR ***


    // Do not look at 9-o-clock if it is a four-body stabilizer on the left boundary of the main lattice.

    bool LeftOfLattice = ( xPos[Site] == 0 );
    bool WrongRow = ( yPos[Site]%3 == 2 );

    bool OnLeftBoundary = (LeftOfLattice & WrongRow);

    if( !OnLeftBoundary )
    {
        bool nineOclockError = Errors[ ErrorSite ];
        //cout << ErrorSite << ", ";
        State = State ^ nineOclockError;
    }


    // *** 3-O-CLOCK ERROR ***

    // 3-o-clock error is to the immediate right of the 9-o-clock error
     ErrorSite  = ErrorSite+1;
     
    // Do not check the 3-o-clock error if we are on a boundary tile to the right of the main lattice.

    bool RightOfLattice = 3*xPos[Site] == 3*L-2 - yPos[Site];
    WrongRow = yPos[Site]%3 == 1;

    bool RightOfBoundary = RightOfLattice & WrongRow;

    if( !RightOfBoundary )
    {
        bool threeOclockError = Errors[ ErrorSite  ];
        //cout << ErrorSite << ", ";
        State = State ^ threeOclockError;
    }


    // *** 5- and 7-O-CLOCK ERROR  ***
    // These two errors are easily dealt with together.

    // Boundary tiles on the bottom row do not have 5- and 7-o-clock errors

    bool BottomRow = ( yPos[Site] == 0 );

    if( !BottomRow )
    {

        ErrorSite = 0;

        // The correct formula for the error site depends on the row

        if(yPos[Site] %3 == 0)
            ErrorSite = (int)(2*Site - 2*(L - yPos[Site]/3 ) - 1);


        else if(yPos[Site] %3 == 1)
            ErrorSite = (int)(2*(Site + 1) - 2*(L - (yPos[Site]-1) / 3) - 1);

        else
            ErrorSite = (int)(2*Site - 2*(L - (yPos[Site] +1) / 3) - 1);



        bool fiveOclockError = Errors[ErrorSite+1];
                //cout << ErrorSite+1 << ", ";

        bool sevenOclockError = Errors[ErrorSite];
                //cout << ErrorSite << ", ";

        State = State ^ fiveOclockError ^ sevenOclockError;

    }


    // *** 11- and 1-O-CLOCK ERROR ***

    ErrorSite = 0;


    // 11-O-CLOCK ERROR SITE
    if(yPos[Site]%3 == 0)
    {
        ErrorSite=(int)(2*Site + 2*(L - yPos[Site]/3 )+1);
    }
    else if(yPos[Site]%3 == 1)
    {
        ErrorSite=(int) (2*Site + 2*(L - (yPos[Site]-1)/3 ) + 1);
    }
    else
    {
        ErrorSite=(int) (2*Site + 2*(L - (yPos[Site]-2)/3 ) - 1);
    }




    LeftOfLattice = ( xPos[Site] == 0 );
    WrongRow = ( yPos[Site]%3 == 2 );

    OnLeftBoundary = (LeftOfLattice & WrongRow);

    if( !OnLeftBoundary )
    {
        bool elevenOclockError= Errors[ErrorSite];
                //cout << ErrorSite << ", ";
        State = State^elevenOclockError;
    }


    // 1-O-CLOCK ERROR SITE
    // The site is one on from the 11-o-clock site
    ErrorSite = ErrorSite+1;


    if( !(3*xPos[Site] == 3*L-2 - yPos[Site]  &  yPos[Site] % 3 == 1 )) 
    {
        bool oneOclockError = Errors[ErrorSite];
                //cout << ErrorSite << endl;

        State = State^oneOclockError;

    }

}


//Function indicates whether site1 and site 2 are in two different restricted lattices (and which ones) or not (returns 0)
int CreaseCross(int site1, int site2, int NoOfVertices, int prev)
{
    //The function returns a value = blueflip + 2*redflip + 4*greenflip. This value is able to store information about all three RLs
    // RL0 = RB
    // RL1 = GB
    // RL2 = RG

    bool blueflip = prev % 2;
    bool redflip = (prev / 2 ) % 2;
    bool greenflip =  (prev / 4 ) % 2;

    int lat1 = site1/NoOfVertices; //returns 0, 1 or 2
    int lat2 = site2/NoOfVertices;

        if (lat1 != lat2)
        {
            if (lat1 + lat2 == 1) // if crosses blue edge
            {
                blueflip = blueflip^1;
            }

            if (lat1 + lat2 == 2) // if crosses red edge
            {
                redflip = redflip^1;
            }

            if (lat1 + lat2 == 3) // if crosses green edge
            {
                greenflip = greenflip^1;
            }
        }

    return blueflip + 2*int(redflip) + 4*int(greenflip);
}


// Function to complete logical measurement
int LogicalMeasurement(bool QubitsX[], int L)
{
    
    int NoOfQubits = 1 + 3*L*(L+1);

    //there are 3 edges on the mobius strip
    //RB-RG edge (red edge)
    //RB-GB edge (blue edge)
    //GB-RG edge (green edge)


    //The function returns a value = blueflip + 2*redflip + 4*greenflip. This value is able to store information about all three edges
    //the qubits on the edges of the triangle (there are 2L+1 qubits on each edge) can be identified using quadratic equations dependant on i and L

    bool blueflip = 0;
    bool redflip = 0;
    bool greenflip = 0;

    //measuring greenflip (qubits on green edge are along the bottom of the triangle)
    for (int i = 0; i < 2*L+1; i++)
    {
        greenflip = QubitsX[i]^greenflip;
    }

    //measuring redflip (along the left of the triangle)
    for (int i = 0; i < 2*L+1; i++)
    {
        if (i%2 == 0)
        {
            int site = -0.75*i*i + 3*(2*L+1)*i/2;
            redflip = QubitsX[site]^redflip;
        }
        else
        {
            int site = -0.75*i*i - L - 0.25 + (3*L+2)*i;
            redflip = QubitsX[site]^redflip;           
        }
    }


    //measuring blueflip (along the right of the triangle)
    for (int i = 0; i < 2*L+1; i++)
    {
        if (i%2 == 0)
        {
            int site = -0.75*i*i + (6*L+1)*i/2 + 2*L;
            //cout << site << "  ";
            blueflip = QubitsX[site]^blueflip;
        }
        else
        {
            int site = -0.75*i*i + 3*L*i + 3*L - 0.25;
            //cout << site << "  ";
            blueflip = QubitsX[site]^blueflip;           
        }

    }
    //cout << endl;


    return blueflip + 2*redflip + 4*greenflip;

}



//Function to help show which cases cause the decoder to fail
/*
 string PrintFail(bool QubitsX[], int NoOfQubits)
{
    int counter = 0;
    string outstr = "Qubits: ";
    for (int i = 0; i < NoOfQubits; i++)
    {
        if (QubitsX[i] == 1)
        {
            outstr = outstr + to_string(i) + "  ";
            counter++;
        }
    }
    outstr = outstr + "\n";
    return outstr;
    //cout << "No of Qubits Failed: " << counter << endl;
}
*/

int CountFail(bool QubitsX[], int NoOfQubits)
{
    int counter = 0;
    for (int i = 0; i < NoOfQubits; i++)
    {
        if (QubitsX[i] == 1)
        {
            counter++;
        }
    }
    return counter;    //cout << "No of Qubits Failed: " << counter << endl;
}


//Function to print a representation of the Mobius Strip
void MobiusPrinter(vector<int> WrappedMobius, int V)
{

    //V is the number of stabilizers in the MS
    //WrappedMobius is a vector of length V*V


    //Converting the vector into a matrix
    int MobiusMatrix[V][V];

    for (int i = 0; i < V; i++)
    {
        for (int j = 0; j < V; j++)
        {
            MobiusMatrix[i][j] = WrappedMobius[i*V+j];
        }
    }

    cout << endl;

    int NoOfVertices = V/3;


    //The representation format used is that first RB stabilizers are stored, then GB, then RG

    for (int i = 0; i < V; i++)
    {
        if ( i == NoOfVertices/2)
            cout << "RB" << "";
        if ( i == NoOfVertices/2 +NoOfVertices)
            cout << "GB" << "";
        if ( i == NoOfVertices/2 +2*NoOfVertices)
            cout << "RG" << "";
        else
            cout << "   ";

        if ((i+1) % NoOfVertices == 0)
            cout << " ";
    }
    cout << endl;


    for (int i = 0; i < V; i++)
    {
        if ( i < 10)
            cout << i << "  ";
        else
            cout << i << " ";

        if ((i+1) % NoOfVertices == 0)
            cout << "   ";
    }
    cout << endl;

    for (int j = 0; j < V+2; j++)
            {
                cout << "__ ";
            }

     cout << endl;


    for (int i = 0; i < V; i++)
    {
        for (int j = 0; j< V; j++)
        {
            if (MobiusMatrix[i][j] < 99999)
                cout << MobiusMatrix[i][j] << "  ";
            else 
                cout << "   ";

            if ((j+1) % NoOfVertices == 0)
                cout << "|  ";
        }
        cout << endl;

        if ((i+1) % NoOfVertices == 0)
        {
            for (int j = 0; j < V+2; j++)
            {
                cout << "__ ";
            }

            cout << endl;
        }

    }
    cout << endl;


}

//Function to fill Stabilizer and Qubit locations on triangular lattice
void XYPosFiller(int StabXPos[], int StabYPos[], int QubXPos[], int QubYPos[], int L, int NoOfVertices)
{
    int NoOfQubits = 2*NoOfVertices + 1;
    
    for(int j = 0; j < NoOfVertices; j++)
    {

        StabXPos[j] = 0;
        StabYPos[j] = 0;

    }
    for(int j = 0; j < NoOfQubits; j++)
    {

        QubXPos[j] = 0;
        QubYPos[j] = 0;

    }

    int StartPoint = 0;

        for(int Rows = 0; Rows < L; ++Rows)
        {
            int RowLength = L - Rows;

            for(int j = 0; j < RowLength; ++j)
            {
                StabXPos[StartPoint + j] = j;
                StabXPos[StartPoint + RowLength + j ] = j;
                StabXPos[StartPoint + 2*RowLength + j ] = j;


                StabYPos[StartPoint+j] = (3*Rows);
                StabYPos[StartPoint + RowLength + j ] = (3*Rows + 1);
                StabYPos[StartPoint + 2*RowLength + j] = (3*Rows + 2);
            }

            StartPoint = StartPoint + 3*RowLength;

        }


        int QubStartPoint = 0;
        int StabStartPoint = 0;

        for(int Rows = 0; Rows < 3*L +1; Rows++)
        {
            int RowLength = 3*L+1 - Rows;
            int CColor = Rows % 3;

            for(int j = 0; j < RowLength; j++)
            {

                    if ((CColor == 0 && j % 3 != 1) || (CColor == 1 && j % 3 != 2) || (CColor == 2 && j % 3 != 0))
                    {
                        QubYPos[QubStartPoint] = Rows;
                        QubXPos[QubStartPoint] = j;
                        QubStartPoint++;
                    }
                    else
                    {
                        //cout << StabStartPoint << " : (" << TriStabXPos[StabStartPoint] << "," << TriStabYPos[StabStartPoint] << ") sheared : (" << StabXPos60[StabStartPoint] << "," << StabYPos60[StabStartPoint] << ")" << endl;
                        StabStartPoint++;
                    }              

            }

        }

}



//Function to copy stabilizer measurements on triangular lattice over to RLs on Mobius Strip
void ExpandToMobius(bool MobiusStabs[], int StabXPos[], int StabYPos[], bool StabilizersZ[], int NoOfVertices)
{
        for (int i = 0; i < NoOfVertices; i++)
        {
            MobiusStabs[i] = 0;
            MobiusStabs[i+NoOfVertices] = 0;
            MobiusStabs[i+2*NoOfVertices] = 0;

            if (StabYPos[i]%3 == 0) // if blue add to RB and GB
            {
                MobiusStabs[i] = StabilizersZ[i];
                MobiusStabs[i+NoOfVertices] = StabilizersZ[i]; 
            }
            if (StabYPos[i]%3 == 1) // if red add to RB and RG
            {
                MobiusStabs[i] = StabilizersZ[i];
                MobiusStabs[i+2*NoOfVertices] = StabilizersZ[i]; 
            }
            if (StabYPos[i]%3 == 2) // if green add to RG and GB
            {
                MobiusStabs[i+2*NoOfVertices] = StabilizersZ[i];
                MobiusStabs[i+NoOfVertices] = StabilizersZ[i]; 
            }
        }
}


//find lexographically next combination
bool next_combination(vector<int>& a, int n) {
    int k = (int)a.size();
    for (int i = k - 1; i >= 0; i--) {
        if (a[i] < n - k + i + 1) {
            a[i]++;
            for (int j = i + 1; j < k; j++)
                a[j] = a[j - 1] + 1;
            return true;
        }
    }
    return false;
}