2D Ising model

ising2d.c

// http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7
// 上のリンクを参考にして作成
// c99以上が必要
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
#include <math.h>

/*
   dSFMT-src-2.2.3 package is downloaded from
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index-jp.html
*/
#define DSFMT_MEXP 521
#include "dSFMT-src-2.2.3/dSFMT.h"
#include "dSFMT-src-2.2.3/dSFMT.c"

//https://ja.wikipedia.org/wiki/Xorshift
uint32_t xor128(void){ 
    static uint32_t x = 123456789;
    static uint32_t y = 362436069;
    static uint32_t z = 521288629;
    static uint32_t w = 88675123; 
    uint32_t t;

    t = x ^ (x << 11);
    x = y; y = z; z = w;
    return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)); 
}

typedef union{
    uint64_t i;
    double f;
} if64_t;

double xor_1_2(void){
    uint64_t a = xor128();
    uint64_t b = xor128();
    uint64_t c = (a<<32) | b;
    c >>= 12;
    if64_t u;
    u.f = 1.0;
    u.i |= c;
    return u.f;
}

double xor_0_1(void){
    return xor_1_2() - 1.0;
}

/******************************************************************************/

typedef unsigned long ulong;
typedef signed int T;

T* rand_ising2d(ulong n){
    T* v = (T*)malloc(sizeof(T)*n*n);
    for(ulong i=0; i<n*n; i++){
        v[i] = rand()%2 == 0 ? -1 : 1;
    }
    return v;
}

T ising2d_sum_of_adjacent_spins(const T* s, ulong m, ulong n, ulong i, ulong j){
    return s[(i<m-1 ? i+1 : 0)*n+j] + s[(i>0 ? i-1 : m-1)*n+j] +
        s[i*n+(j<n-1 ? j+1 : 0)] + s[i*n+(j>0 ? j-1 : n-1)];
}

void ising2d_sweep(T* s, ulong m, ulong n, double beta, ulong niters){
    double prob[4*2+1];
    dsfmt_t dsfmt;
    dsfmt_init_gen_rand(&dsfmt, rand());
    for(long k=-4; k<=4; k++){
        prob[k+4] = exp(-2*beta*k);
    }
    for(ulong count=0; count<niters/(m*n); count++){
        if(count % 1000 == 0){ 
            fprintf(stderr, "%lu\r", count);
        }
        for(ulong i=0; i<m; i++){
            for(ulong j=0; j<n; j++){
                T s1 = s[i*n+j];
                T k = s1 * ising2d_sum_of_adjacent_spins(s, m, n, i, j);
                s[i*n+j] = dsfmt_genrand_close_open(&dsfmt)<prob[k+4] ? -s1 : s1; //dsfmt
                //s[i*n+j] = xor_0_1()<prob[k+4] ? -s1 : s1; //xorshift128
            }
        }
    }
}

void print_pbm(const char* s, const T* v, ulong m, ulong n){
    // PNM形式で出力
    // https://ja.wikipedia.org/wiki/PNM_(%E7%94%BB%E5%83%8F%E3%83%95%E3%82%A9%E3%83%BC%E3%83%9E%E3%83%83%E3%83%88)
    FILE* fp = fopen(s, "w");
    fprintf(fp, "P1\n");
    fprintf(fp, "%lu %lu\n", n , m);
    for(ulong i=0; i<m; i++){
        for(ulong j=0; j<n; j++){
            fprintf(fp, "%d ", (int)(v[i*n+j]==1 ? 1 : 0));
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
}

int main(void){
    srand((unsigned)time(NULL));
    //雑に初期化
    int k=rand();
    for(int i=0; i<k; i++) xor128();
    double beta_crit=log(1+sqrt(2))/2;
    const int n=100;
    T* s = rand_ising2d(n);
    print_pbm("start.pbm", s, n, n);
    time_t st = clock();
    ising2d_sweep(s, n, n, beta_crit, 1000000000);
    time_t et = clock();
    print_pbm("result.pbm", s, n, n);
    printf("%.3f\n", (double)(et-st)/CLOCKS_PER_SEC);
}

ising2d.cpp

// http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7
// 上のリンクを参考にして作成
// c++11以上が必要 (g++, clang++なら -std=c++11)

#include <iostream>
#include <fstream>
#include <random>
#include <vector>
#include <chrono>
#include <math.h>

typedef unsigned long ulong;

template <typename T> void rand_ising2d(std::vector<T>& v, ulong n, std::mt19937& engine){
    std::uniform_int_distribution<> dist(0, 1);
    v.resize(n*n);
    for(ulong i=0; i<n*n; i++){
        v[i] = dist(engine) == 0 ? -1 : 1;
    }
}

template <typename T> T ising2d_sum_of_adjacent_spins(const std::vector<T>& s, ulong m, ulong n, ulong i, ulong j){
    return s[(i<m-1 ? i+1 : 0)*n+j] + s[(i>0 ? i-1 : m-1)*n+j] +
        s[i*n+(j<n-1 ? j+1 : 0)] + s[i*n+(j>0 ? j-1 : n-1)];
}

template <typename T> void ising2d_sweep(std::vector<T>& s, ulong m, ulong n, double beta, ulong niters, std::mt19937& engine){
    constexpr std::size_t bits = std::numeric_limits<double>::digits;
    std::vector<double> prob(4*2+1);
    for(long k=-4; k<=4; k++){
        prob[k+4] = exp(-2*beta*k);
    }
    for(ulong count=0; count<niters/(m*n); count++){
        if(count % 1000 == 0){ 
            std::cerr << count << "\r" << std::flush;
        }
        for(ulong i=0; i<m; i++){
            for(ulong j=0; j<n; j++){
                T s1 = s[i*n+j];
                T k = s1 * ising2d_sum_of_adjacent_spins(s, m, n, i, j);
                s[i*n+j] = std::generate_canonical<double, bits>(engine)<prob[k+4] ? -s1 : s1;
            }
        }
    }
}

template <typename T> void print_pbm(const char* s, std::vector<T>& v, ulong m, ulong n){
    // PNM形式で出力
    // https://ja.wikipedia.org/wiki/PNM_(%E7%94%BB%E5%83%8F%E3%83%95%E3%82%A9%E3%83%BC%E3%83%9E%E3%83%83%E3%83%88)
    std::ofstream file;
    file.open(s, std::ofstream::out);
    file << "P1" << std::endl;
    file << n << " " << m << std::endl;
    for(ulong i=0; i<m; i++){
        for(ulong j=0; j<n; j++){
            file << (v[i*n+j]==1 ? 1 : 0) << " ";
        }
        file << "\n";
    }
    file.close();
}

int main(void){
    std::random_device seed_gen;
    std::mt19937 engine(seed_gen());
    constexpr double beta_crit=log(1+sqrt(2))/2;
    std::vector<signed char> s;
    rand_ising2d(s, 100, engine);
    print_pbm("start.pbm", s, 100, 100);
    auto begin = std::chrono::high_resolution_clock::now();
    ising2d_sweep(s, 100, 100, beta_crit, 1000000000, engine);
    auto end = std::chrono::high_resolution_clock::now();
    print_pbm("result.pbm", s, 100, 100);
    auto elapsed_time = std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin);
    std::cerr << elapsed_time.count()*1e-9 << "s" << std::endl;
}

ising2d_dSFMT.cpp

// http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7
// 上のリンクを参考にして作成
// c++11以上が必要 (g++, clang++なら -std=c++11)

#include <iostream>
#include <fstream>
#include <random>
#include <vector>
#include <chrono>
#include <math.h>

/*
   dSFMT-src-2.2.3 package is downloaded from
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index-jp.html
*/
#define DSFMT_MEXP 521
#include "dSFMT-src-2.2.3/dSFMT.h"
#include "dSFMT-src-2.2.3/dSFMT.c"

typedef unsigned long ulong;

template <typename T> void rand_ising2d(std::vector<T>& v, ulong n, std::mt19937& engine){
    std::uniform_int_distribution<> dist(0, 1);
    v.resize(n*n);
    for(ulong i=0; i<n*n; i++){
        v[i] = dist(engine) == 0 ? -1 : 1;
    }
}

template <typename T> T ising2d_sum_of_adjacent_spins(const std::vector<T>& s, ulong m, ulong n, ulong i, ulong j){
    return s[(i<m-1 ? i+1 : 0)*n+j] + s[(i>0 ? i-1 : m-1)*n+j] +
        s[i*n+(j<n-1 ? j+1 : 0)] + s[i*n+(j>0 ? j-1 : n-1)];
}

template <typename T> void ising2d_sweep(std::vector<T>& s, ulong m, ulong n, double beta, ulong niters, std::mt19937& engine){
    dsfmt_t dsfmt;
    dsfmt_init_gen_rand(&dsfmt, engine());
    std::vector<double> prob(4*2+1);
    for(long k=-4; k<=4; k++){
        prob[k+4] = exp(-2*beta*k);
    }
    for(ulong count=0; count<niters/(m*n); count++){
        if(count % 1000 == 0){ 
            std::cerr << count << "\r" << std::flush;
        }
        for(ulong i=0; i<m; i++){
            for(ulong j=0; j<n; j++){
                T s1 = s[i*n+j];
                T k = s1 * ising2d_sum_of_adjacent_spins(s, m, n, i, j);
                s[i*n+j] = dsfmt_genrand_close_open(&dsfmt)<prob[k+4] ? -s1 : s1;
            }
        }
    }
}

template <typename T> void print_pbm(const char* s, std::vector<T>& v, ulong m, ulong n){
    // PNM形式で出力
    // https://ja.wikipedia.org/wiki/PNM_(%E7%94%BB%E5%83%8F%E3%83%95%E3%82%A9%E3%83%BC%E3%83%9E%E3%83%83%E3%83%88)
    std::ofstream file;
    file.open(s, std::ofstream::out);
    file << "P1" << std::endl;
    file << n << " " << m << std::endl;
    for(ulong i=0; i<m; i++){
        for(ulong j=0; j<n; j++){
            file << (v[i*n+j]==1 ? 1 : 0) << " ";
        }
        file << "\n";
    }
    file.close();
}

int main(void){
    std::random_device seed_gen;
    std::mt19937 engine(seed_gen());
    double beta_crit=log(1+sqrt(2))/2;
    std::vector<signed char> s;
    rand_ising2d(s, 100, engine);
    print_pbm("start.pbm", s, 100, 100);
    auto begin = std::chrono::high_resolution_clock::now();
    ising2d_sweep(s, 100, 100, beta_crit, 1000000000, engine);
    auto end = std::chrono::high_resolution_clock::now();
    print_pbm("result.pbm", s, 100, 100);
    auto elapsed_time = std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin);
    std::cerr << elapsed_time.count()*1e-9 << "s" << std::endl;
}