ising.cpp

#include <sys/time.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>

static unsigned int g_seed;

inline int fastrand() {
  g_seed = (214013*g_seed+2531011);
  return (g_seed>>16)&0x7FFF;
}

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

    int L          = atoi(argv[1]);
    __uint128_t iterations  = atoi(argv[2]);
    float T                 = atof(argv[3]);

    g_seed = static_cast<unsigned int>(time(0));

    __int128_t E = 0;
    __int128_t M = 0;

    unsigned short n[L][L][4][2];
    char grid[L][L];

    for(short x = 0; x < L; x ++) {
        for(short y = 0; y < L; y ++) {

            short xplus = x + 1;
            short yplus = y + 1;

            short xmin = x - 1;
            short ymin = y - 1;

            if (xplus == L) xplus = 0;
            if (yplus == L) yplus = 0;

            if (xmin == -1) xmin = L - 1;
            if (ymin == -1) ymin = L - 1;

            n[x][y][0][0] = x;
            n[x][y][0][1] = ymin;

            n[x][y][1][0] = x;
            n[x][y][1][1] = yplus;

            n[x][y][2][0] = xmin;
            n[x][y][2][1] = y;

            n[x][y][3][0] = xplus;
            n[x][y][3][1] = y;

            grid[x][y] = 1;

            E -= grid[x][y] * (grid[x][yplus] + grid[xplus][y]);
            M += grid[x][y];
        }
    }


    int bf[20];
    rbf[0]  = (int)(exp(8.0 / T)    * 32767);
    bf[4]  = (int)(exp(4.0 / T)    * 32767);
    bf[8]  = (int)(1.0             * 32767);
    bf[12] = (int)(exp( - 4.0 / T) * 32767);
    bf[16] = (int)(exp( - 8.0 / T) * 32767);

    int x = 0;
    int y = 0;

    __int128_t dE = 0;

    __int128_t E_sum  = 0;
    __int128_t E2_sum = 0;

    __int128_t M_sum  = 0;
    __int128_t M2_sum = 0;


    struct timeval start_time;
    struct timeval end_time;

    gettimeofday(&start_time, NULL);

    for (__int128_t N = 0; N < iterations*2; N ++) {

        x = fastrand() % L;
        y = fastrand() % L;


        dE = 2 * grid[x][y] * (grid[n[x][y][0][0]][n[x][y][0][1]] +
                               grid[n[x][y][1][0]][n[x][y][1][1]] +
                               grid[n[x][y][2][0]][n[x][y][2][1]] +
                               grid[n[x][y][3][0]][n[x][y][3][1]]);


        if (dE < 0) {

            grid[x][y] *= -1;
            E += dE;
            M += 2 * grid[x][y];

        } else if (bf[8 + dE] > fastrand()) {

            grid[x][y] *= -1;
            E += dE;
            M += 2 * grid[x][y];
        }


        if (N > iterations) {
            E_sum += E;
            E2_sum += E * E;

            M_sum += M;
            M2_sum +=  M*M;
        }


    }

    gettimeofday(&end_time, NULL);


    double start = start_time.tv_sec*1000000 + (start_time.tv_usec);
    double end   =   end_time.tv_sec*1000000 + (end_time.tv_usec);

    double elapsed = end - start;


    float E_mean  = (1.0 * E_sum)  / (iterations * L*L);
    float E2_mean = (1.0 * E2_sum) / (iterations * L*L*L*L);

    float M_mean  = (1.0 * M_sum)  / (iterations * L*L);
    float M2_mean = (1.0 * M2_sum) / (iterations * L*L*L*L);


    float Cv =  (1.0 /(T * T)) * (E2_mean - E_mean * E_mean);
    float Xi =  1.0 / T * (M2_mean - M_mean * M_mean);


    std::cout << "  time = " << elapsed/1000000.0 << std::endl;
    std::cout << "   T   = " << T       << std::endl;
    std::cout << "  <E>  = " << E_mean  << std::endl;
    std::cout << "  <E2> = " << E2_mean << std::endl;
    std::cout << "  <M>  = " << M_mean  << std::endl;
    std::cout << "  <M2> = " << M2_mean << std::endl;
    std::cout << "   Cv  = " << Cv << std::endl;
    std::cout << "   Xi  = " << Xi << std::endl;


    return 0;

rand(time(NULL));