#include <tr1/random>
#include <boost/random.hpp>

#include <boost/multi_array.hpp>
#include <iostream>
#include <fstream>

using namespace std;
using namespace boost;

typedef mt19937 boost_rng_t;
typedef tr1::mt19937 tr1_rng_t;

typedef boost::multi_array<int, 2> array_t;

void wall(int& x, int& y, size_t n)
{
    x = max(min(x, int(n-1)), 0);
    y = max(min(y, int(n-1)), 0);
}


int main()
{
    size_t n = 100;

    array_t bcount(extents[n][n]);
    array_t tcount(extents[n][n]);


    boost_rng_t brng(static_cast<boost_rng_t::result_type>(42));
    uniform_int<int> brandom(-1, 1);

    tr1_rng_t trng(static_cast<tr1_rng_t::result_type>(42));
    tr1::uniform_int<int> trandom(-1, 1);

    // start 10 random walks starting from the middle of the lattice
    for (size_t j = 0; j < 10; j++)
    {
        int xt, yt, xb, yb;
        xt = yt = xb = yb = n/2;

        for (size_t i = 0; i < 10000000; ++i)
        {
            xb += brandom(brng);
            yb += brandom(brng);
            wall(xb, yb, n);
            bcount[xb][yb]++;

            xt += trandom(trng);
            yt += trandom(trng);
            wall(xt, yt, n);
            tcount[xt][yt]++;
        }
    }

    ofstream f;

    f.open("bcount");
    for (size_t i = 0; i < n; ++i)
        for (size_t j = 0; j < n; ++j)
            f << i << " " << j << " " << bcount[i][j] << endl;
    f.close();

    f.open("tcount");
    for (size_t i = 0; i < n; ++i)
        for (size_t j = 0; j < n; ++j)
            f << i << " " << j << " " << tcount[i][j] << endl;
    f.close();

    return 1;
}