#define GRID_WIDTH  1024
#define GRID_HEIGHT 1024
#define MAX_STEPS 1024

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

double grid[GRID_HEIGHT][GRID_WIDTH];
double gridNewArray[GRID_HEIGHT][GRID_WIDTH];
double gridOldArray[GRID_HEIGHT][GRID_WIDTH];

void printAddress(void *p)
{
    printf("address %p\n", p);
}

void test(double (*grid)[GRID_WIDTH])
{
    double (*gridOld)[GRID_WIDTH] = gridOldArray;
    double (*gridNew)[GRID_WIDTH] = gridNewArray;
    printAddress(&gridNew[0][0]);
    printAddress(&gridOld[0][0]);

    // copy initial state
    for (int y = 0; y < GRID_HEIGHT; ++y) {
        memcpy(&gridOld[y][0], &grid[y][0], GRID_WIDTH * sizeof(double));
        memset(&gridNew[y][0], 0, GRID_WIDTH * sizeof(double));
    }

    // update matrices
    for (int step = 0; step < MAX_STEPS; ++step) {
        for (int y = 1; y < GRID_HEIGHT-1; ++y) 
            for (int x = 1; x < GRID_WIDTH-1; ++x)
                gridNew[y][x] = 
                    (gridOld[y-1][x  ] + 
                     gridOld[y  ][x-1] + 
                     gridOld[y  ][x  ] + 
                     gridOld[y  ][x+1] + 
                     gridOld[y+1][x  ]) * 0.2;
        double (*tmp)[GRID_WIDTH] = gridOld;
        gridOld = gridNew;
        gridNew = tmp;
    }

    // copy result back
    for (int y = 0; y < GRID_HEIGHT; ++y)
        memcpy(&grid[y][0], &gridOld[y][0], GRID_WIDTH * sizeof(double));
}

void setupGrid()
{
    for (int y = 0; y < GRID_HEIGHT; ++y)
        for (int x = 0; x < GRID_WIDTH; ++x)
            grid[y][x] = 0;

    for (int y = 10; y < 20; ++y)
        for (int x = 10; x < 20; ++x)
            grid[y][x] = 1;
}

int main(int argc, char** argv)
{
    setupGrid();
    test(grid);
    printf("res: %f\n", grid[10][10]); // prevent dead code elimination
    return 0;
}