x-or/unkgen.cpp

## unkgen.cpp
///// Kiss My Glowing Pixel! (demo)
///// Copyright (c) 2013, Sergey Shishmintzev
///// License: Apache/GPLv3
///// Example: http://youtu.be/xVpsqkSMgOI
///// Compiler: gcc-4.7.3 (cygwin)
///// Compile command: g++ -O3  unkgen.cpp -lpng
///// P.S. Sorry for my English. It isn't my mother tongue.

#include <assert.h>
#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <png.h>
#include <iostream>
#include <string>
#include <vector>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/random.hpp>
#include <boost/generator_iterator.hpp>

using namespace boost::numeric::ublas;

inline int mod(int a, int b) { // http://stackoverflow.com/questions/4003232/how-to-code-a-modulo-operator-in-c-c-obj-c-that-handles-negative-numbers
	if(b < 0) //you can check for b == 0 separately and do what you want
		return mod(-a, -b);
	int ret = a % b;
	if (ret < 0)
		ret += b;
	return ret;
}


inline int make_rgb(int r, int g, int b) {
	return (b&0xFF) | ((g&0xFF)<<8) | ((r&0xFF)<<16);
}


inline int get_r(int rgb) {
	return (rgb>>16) & 0xFF;
}


inline int get_g(int rgb) {
	return (rgb>>8) & 0xFF;
}


inline int get_b(int rgb) {
	return (rgb) & 0xFF;
}


inline void setRGB(png_byte *ptr, int val) {
	ptr[0] = get_r(val);
	ptr[1] = get_g(val);
	ptr[2] = get_b(val);
}


int save_matrix(const matrix<int> &m, const std::string &filename, const std::vector<int> &palette) {
	// http://www.labbookpages.co.uk/software/imgProc/libPNG.html
	int code = 0;
	FILE *fp = NULL;
	png_structp png_ptr = NULL;
	png_infop info_ptr = NULL;
	png_bytep row = NULL;

	// Open file for writing (binary mode)
	fp = fopen(filename.c_str(), "wb");
	if (fp == NULL) {
		fprintf(stderr, "Could not open file %s for writing\n", filename.c_str());
		code = 1;
		goto finalise;
	}

	// Initialize write structure
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
	if (png_ptr == NULL) {
		fprintf(stderr, "Could not allocate write struct\n");
		code = 1;
		goto finalise;
	}

	// Initialize info structure
	info_ptr = png_create_info_struct(png_ptr);
	if (info_ptr == NULL) {
		fprintf(stderr, "Could not allocate info struct\n");
		code = 1;
		goto finalise;
	}

	// Setup Exception handling
	if (setjmp(png_jmpbuf(png_ptr))) {
		fprintf(stderr, "Error during png creation\n");
		code = 1;
		goto finalise;
	}

	png_init_io(png_ptr, fp);

	// Write header (8 bit colour depth)
	png_set_IHDR(png_ptr, info_ptr, m.size2(), m.size1(),
			8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
			PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

	png_write_info(png_ptr, info_ptr);

	// Allocate memory for one row (3 bytes per pixel - RGB)
	row = (png_bytep) malloc(3 * m.size2() * sizeof(png_byte));

	// Write image data
	unsigned int i, j;
	for (i = 0; i < m.size1(); i++) {
		for (j = 0; j < m.size2(); j++) {
			setRGB(&(row[j*3]), palette[m(i, j)]);
		}
		png_write_row(png_ptr, row);
	}

	// End write
	png_write_end(png_ptr, NULL);

finalise:
	if (fp != NULL) fclose(fp);
	if (info_ptr != NULL) png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
	if (png_ptr != NULL) png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
	if (row != NULL) free(row);

	return code;
}

// Gradient palette routines
// http://jtauber.com/blog/2008/05/18/creating_gradients_programmatically_in_python/
int linear_gradient(int start_value, int stop_value, float start_offset, float stop_offset, float offset) {
	return floor(start_value + ((offset - start_offset) / (stop_offset - start_offset) * (stop_value - start_value)));
}

int gradient(int start, int end, float initial_offset, float stop_offset, float offset) {
	// make gradient color
	int r = linear_gradient(get_r(start), get_r(end), initial_offset, stop_offset, offset);
	int g = linear_gradient(get_g(start), get_g(end), initial_offset, stop_offset, offset);
	int b = linear_gradient(get_b(start), get_b(end), initial_offset, stop_offset, offset);
	return make_rgb(r, g, b);
}


typedef struct point_t {
	int value; // point weight
	int i, j; // coordinates
	int di, dj; // coordinate increment
	int edge; // circle edge
} point_t;


int main() {
	using std::cout;
	using std::endl;
	using std::min;
	using std::max;
	using std::vector;

	// random generators seed
	const int seed = time(0);
	const int modulo = 1024; // arbitrary even number
	// output image size
	const int width = 1920;
	const int height = 1080;
	const int points_count = 24;

	cout << "seed = " << seed << endl;
	cout << "modulo = " << modulo << endl;
	cout << "size = " << width << "x" << height << endl;
	cout << "points count = " << points_count << endl;

	// make double gradient palette
	vector<int> palette(modulo);

	//const int start_color1 = make_rgb(171, 214, 121); // green
	const int start_color1 = make_rgb(204, 255, 145); // green
	const int stop_color1 = make_rgb(65, 66, 225); // blue
	//const int start_color2 = make_rgb(171, 214, 121); // green
	const int start_color2 = start_color1; // green
	const int stop_color2 = make_rgb(209, 35, 42); // red

	for(int i = 0; i < modulo/2; i++) {
		palette[i] = gradient(start_color1, stop_color1, 0.0, 1.0, min(i, modulo/2-1)*1.0/(modulo/2-1));
		palette[modulo-i-1] = gradient(start_color2, stop_color2, 0.0, 1.0, min(i, modulo/2-1)*1.0/(modulo/2-1));
	}
#if 1
	const int intermediate_color = make_rgb(145, 35, 145); // red
	for(int i = 0; i < modulo/4; i++) {
		palette[modulo/4+i] = gradient(palette[modulo/4], intermediate_color, 0.0, 1.0, i*1.0/(modulo/4-1));
		palette[-modulo/4+modulo-i-1] = gradient(palette[-modulo/4+modulo-1], intermediate_color, 0.0, 1.0, i*1.0/(modulo/4-1));
	}
#endif
	vector<point_t> point(points_count);

	// particles random stuff
	typedef boost::mt19937 RNGType;
	RNGType rng1(seed);
	RNGType rng2(seed^3);
	RNGType rng3(seed^7);
	RNGType rng4(seed^15);
	RNGType rng5(seed^31);
	RNGType rng01(seed^1023);

	boost::uniform_int<> i_range( 0, height-1 );
	boost::uniform_int<> j_range( 0, width-1 );
	boost::uniform_int<> direction_range( 0, 4 );
	boost::uniform_int<> value_range( 1, 25 );
	boost::uniform_int<> edge_range( 100, 200 );
	boost::variate_generator< RNGType, boost::uniform_int<> > random_i(rng1, i_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_j(rng2, j_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_direction_index(rng3, direction_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_value(rng4, value_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_edge(rng5, edge_range);
	boost::uniform_01<RNGType> zeroone(rng01);

	const int direction[] = {-3, 3, 0, -4, 4}; // slow and fast moves

	// generate random positions for points
	for (int p = 0; p < points_count; p++) {
		point[p].edge = random_edge();
		int i = random_i();
		int j = random_j();
		assert(i >= 0);
		assert(j >= 0);
		assert(i < height);
		assert(j < width);
		point[p].i = i;
		point[p].j = j;
		point[p].di = direction[random_direction_index()];
		point[p].dj = direction[random_direction_index()];
		point[p].value = 15*point[p].edge/100;
		if (zeroone() < 0.5)
			point[p].value = mod(-point[p].value, modulo);
	}

	const int frames = 10000; // count of the images to produce

        matrix<int> state(height, width);
        matrix<int> future_state(height, width);

	const int bacground_cell_size = trunc(sqrt(height)); // background pixelization

	for (int frame = 0; frame < frames; frame++) {
		cout << "frame #" << frame << endl;
		future_state.clear();

		// circle shaped background (final circle formula is x*(x+1)/2 + y*(y+1)/2 == v (mod modulo))
		for (int i = 0; i < height; i+=bacground_cell_size)
				future_state(i, 0) = mod((2*i)/bacground_cell_size, modulo);
		for (int j = 0; j < width; j+=bacground_cell_size)
				future_state(0, j) = mod((2*j)/bacground_cell_size, modulo);

		for (int p = 0; p < points_count; p++) {
			int i, j;
			int di = point[p].di;
			int dj = point[p].dj;
			while (1) { // failback loop
				i = point[p].i + di;
				j = point[p].j + dj;

				// canvas edge detection
				if ((i >= height) || (i < 0)) {
					di = direction[random_direction_index()];
					continue;
				}
				if ((j >= width) || (j < 0)) {
					dj = direction[random_direction_index()];
					continue;
				}

				// chaotic moving
				if (zeroone() < 0.01)
					di = direction[random_direction_index()];
				if (zeroone() < 0.01)
					dj = direction[random_direction_index()];
				if (di == 0 && dj == 0)
					continue;
				break;
			} // retry ...
			assert(i >= 0);
			assert(j >= 0);
			assert(i < height);
			assert(j < width);

			// new point position
			point[p].i = i;
			point[p].j = j;
			// new point position increment
			point[p].di = di;
			point[p].dj = dj;

			int ci, cj;
			const int edge = point[p].edge;
			int v;
			// draw a diamond. avoid glitches
			for (int c = 0; c < edge; c++) {
				v = c/13; // smooth colors

				// left top edge
				ci = i + c;
				cj = j + edge - c - 1;
				if (ci >= 0 && ci < height && cj >= 0 && cj < width)
					future_state(ci, cj) = mod(future_state(ci, cj) + (v), modulo);

				// left bottom edge
				ci = i + 2*edge - c - 1;
				cj = j + edge - c - 1;
				if (ci >= 0 && ci < height && cj >= 0 && cj < width)
					future_state(ci, cj) = mod(future_state(ci, cj) + (-v), modulo);

				// right top edge
				ci = i + c;
				cj = j + edge + c;
				if (ci >= 0 && ci < height && cj >= 0 && cj < width)
					future_state(ci, cj) = mod(future_state(ci, cj) + (-v), modulo);

				// right bottom edge
				ci = i + 2*edge - c - 1;
				cj = j + edge + c;
				if (ci >= 0 && ci < height && cj >= 0 && cj < width)
					future_state(ci, cj) = mod(future_state(ci, cj) + (v), modulo);
			}
		} // p = ...

		state.clear();

		// Backwards an automata rule: future_state(i0, j0) == state(i0, j0) - state(i0, j1) - state(i1, j0) + state(i1, j1) (mod modulo).
		// It is the rule 15 [Khan1997] when modulo == 2.
		// It is the unambiguous when we know initial values of the state (state(0, k) == 0 and state(k, 0)  == 0).
		// Reference:
		//   [Khan1997] A.R. Khan, et.al., "VLSI architecture of a cellular automata machine", 1997,
		//              http://dx.doi.org/10.1016/S0898-1221(97)00021-7
		unsigned int i0, i1, j0, j1, solved_state;
		for (i0 = 0; i0 < state.size1()-1; i0++) {
			i1 = i0 + 1;
			for (j0 = 0; j0 < state.size2()-1; j0++) {
				j1 = j0 + 1;
				solved_state = mod(-state(i0, j0) + state(i0, j1) + state(i1, j0) + future_state(i0, j0), modulo);
				assert(solved_state >= 0);
				assert(solved_state < modulo);
				state(i1, j1) = solved_state;
			}
		}
		// produce image
		char fn[100];
		sprintf(fn, "rand%d-%dx%d-frame%05d.png", modulo, width, height, frame);
		save_matrix(state, fn, palette);
	}
	return 0;
}
	///// Kiss My Glowing Pixel! (demo)
	///// Copyright (c) 2013, Sergey Shishmintzev
	///// License: Apache/GPLv3
	///// Example: http://youtu.be/xVpsqkSMgOI
	///// Compiler: gcc-4.7.3 (cygwin)
	///// Compile command: g++ -O3 unkgen.cpp -lpng
	///// P.S. Sorry for my English. It isn't my mother tongue.

	#include <assert.h>
	#include <stdio.h>
	#include <malloc.h>
	#include <math.h>
	#include <png.h>
	#include <iostream>
	#include <string>
	#include <vector>
	#include <boost/numeric/ublas/matrix.hpp>
	#include <boost/random.hpp>
	#include <boost/generator_iterator.hpp>

	using namespace boost::numeric::ublas;

	inline int mod(int a, int b) { // http://stackoverflow.com/questions/4003232/how-to-code-a-modulo-operator-in-c-c-obj-c-that-handles-negative-numbers
	if(b < 0) //you can check for b == 0 separately and do what you want
	return mod(-a, -b);
	int ret = a % b;
	if (ret < 0)
	ret += b;
	return ret;
	}


	inline int make_rgb(int r, int g, int b) {
	return (b&0xFF) \| ((g&0xFF)<<8) \| ((r&0xFF)<<16);
	}


	inline int get_r(int rgb) {
	return (rgb>>16) & 0xFF;
	}


	inline int get_g(int rgb) {
	return (rgb>>8) & 0xFF;
	}


	inline int get_b(int rgb) {
	return (rgb) & 0xFF;
	}


	inline void setRGB(png_byte *ptr, int val) {
	ptr[0] = get_r(val);
	ptr[1] = get_g(val);
	ptr[2] = get_b(val);
	}


	int save_matrix(const matrix<int> &m, const std::string &filename, const std::vector<int> &palette) {
	// http://www.labbookpages.co.uk/software/imgProc/libPNG.html
	int code = 0;
	FILE *fp = NULL;
	png_structp png_ptr = NULL;
	png_infop info_ptr = NULL;
	png_bytep row = NULL;

	// Open file for writing (binary mode)
	fp = fopen(filename.c_str(), "wb");
	if (fp == NULL) {
	fprintf(stderr, "Could not open file %s for writing\n", filename.c_str());
	code = 1;
	goto finalise;
	}

	// Initialize write structure
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
	if (png_ptr == NULL) {
	fprintf(stderr, "Could not allocate write struct\n");
	code = 1;
	goto finalise;
	}

	// Initialize info structure
	info_ptr = png_create_info_struct(png_ptr);
	if (info_ptr == NULL) {
	fprintf(stderr, "Could not allocate info struct\n");
	code = 1;
	goto finalise;
	}

	// Setup Exception handling
	if (setjmp(png_jmpbuf(png_ptr))) {
	fprintf(stderr, "Error during png creation\n");
	code = 1;
	goto finalise;
	}

	png_init_io(png_ptr, fp);

	// Write header (8 bit colour depth)
	png_set_IHDR(png_ptr, info_ptr, m.size2(), m.size1(),
	8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
	PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

	png_write_info(png_ptr, info_ptr);

	// Allocate memory for one row (3 bytes per pixel - RGB)
	row = (png_bytep) malloc(3 * m.size2() * sizeof(png_byte));

	// Write image data
	unsigned int i, j;
	for (i = 0; i < m.size1(); i++) {
	for (j = 0; j < m.size2(); j++) {
	setRGB(&(row[j*3]), palette[m(i, j)]);
	}
	png_write_row(png_ptr, row);
	}

	// End write
	png_write_end(png_ptr, NULL);

	finalise:
	if (fp != NULL) fclose(fp);
	if (info_ptr != NULL) png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
	if (png_ptr != NULL) png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
	if (row != NULL) free(row);

	return code;
	}

	// Gradient palette routines
	// http://jtauber.com/blog/2008/05/18/creating_gradients_programmatically_in_python/
	int linear_gradient(int start_value, int stop_value, float start_offset, float stop_offset, float offset) {
	return floor(start_value + ((offset - start_offset) / (stop_offset - start_offset) * (stop_value - start_value)));
	}

	int gradient(int start, int end, float initial_offset, float stop_offset, float offset) {
	// make gradient color
	int r = linear_gradient(get_r(start), get_r(end), initial_offset, stop_offset, offset);
	int g = linear_gradient(get_g(start), get_g(end), initial_offset, stop_offset, offset);
	int b = linear_gradient(get_b(start), get_b(end), initial_offset, stop_offset, offset);
	return make_rgb(r, g, b);
	}


	typedef struct point_t {
	int value; // point weight
	int i, j; // coordinates
	int di, dj; // coordinate increment
	int edge; // circle edge
	} point_t;


	int main() {
	using std::cout;
	using std::endl;
	using std::min;
	using std::max;
	using std::vector;

	// random generators seed
	const int seed = time(0);
	const int modulo = 1024; // arbitrary even number
	// output image size
	const int width = 1920;
	const int height = 1080;
	const int points_count = 24;

	cout << "seed = " << seed << endl;
	cout << "modulo = " << modulo << endl;
	cout << "size = " << width << "x" << height << endl;
	cout << "points count = " << points_count << endl;

	// make double gradient palette
	vector<int> palette(modulo);

	//const int start_color1 = make_rgb(171, 214, 121); // green
	const int start_color1 = make_rgb(204, 255, 145); // green
	const int stop_color1 = make_rgb(65, 66, 225); // blue
	//const int start_color2 = make_rgb(171, 214, 121); // green
	const int start_color2 = start_color1; // green
	const int stop_color2 = make_rgb(209, 35, 42); // red

	for(int i = 0; i < modulo/2; i++) {
	palette[i] = gradient(start_color1, stop_color1, 0.0, 1.0, min(i, modulo/2-1)*1.0/(modulo/2-1));
	palette[modulo-i-1] = gradient(start_color2, stop_color2, 0.0, 1.0, min(i, modulo/2-1)*1.0/(modulo/2-1));
	}
	#if 1
	const int intermediate_color = make_rgb(145, 35, 145); // red
	for(int i = 0; i < modulo/4; i++) {
	palette[modulo/4+i] = gradient(palette[modulo/4], intermediate_color, 0.0, 1.0, i*1.0/(modulo/4-1));
	palette[-modulo/4+modulo-i-1] = gradient(palette[-modulo/4+modulo-1], intermediate_color, 0.0, 1.0, i*1.0/(modulo/4-1));
	}
	#endif
	vector<point_t> point(points_count);

	// particles random stuff
	typedef boost::mt19937 RNGType;
	RNGType rng1(seed);
	RNGType rng2(seed^3);
	RNGType rng3(seed^7);
	RNGType rng4(seed^15);
	RNGType rng5(seed^31);
	RNGType rng01(seed^1023);

	boost::uniform_int<> i_range( 0, height-1 );
	boost::uniform_int<> j_range( 0, width-1 );
	boost::uniform_int<> direction_range( 0, 4 );
	boost::uniform_int<> value_range( 1, 25 );
	boost::uniform_int<> edge_range( 100, 200 );
	boost::variate_generator< RNGType, boost::uniform_int<> > random_i(rng1, i_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_j(rng2, j_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_direction_index(rng3, direction_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_value(rng4, value_range);
	boost::variate_generator< RNGType, boost::uniform_int<> > random_edge(rng5, edge_range);
	boost::uniform_01<RNGType> zeroone(rng01);

	const int direction[] = {-3, 3, 0, -4, 4}; // slow and fast moves

	// generate random positions for points
	for (int p = 0; p < points_count; p++) {
	point[p].edge = random_edge();
	int i = random_i();
	int j = random_j();
	assert(i >= 0);
	assert(j >= 0);
	assert(i < height);
	assert(j < width);
	point[p].i = i;
	point[p].j = j;
	point[p].di = direction[random_direction_index()];
	point[p].dj = direction[random_direction_index()];
	point[p].value = 15*point[p].edge/100;
	if (zeroone() < 0.5)
	point[p].value = mod(-point[p].value, modulo);
	}

	const int frames = 10000; // count of the images to produce

	matrix<int> state(height, width);
	matrix<int> future_state(height, width);

	const int bacground_cell_size = trunc(sqrt(height)); // background pixelization

	for (int frame = 0; frame < frames; frame++) {
	cout << "frame #" << frame << endl;
	future_state.clear();

	// circle shaped background (final circle formula is x(x+1)/2 + y(y+1)/2 == v (mod modulo))
	for (int i = 0; i < height; i+=bacground_cell_size)
	future_state(i, 0) = mod((2*i)/bacground_cell_size, modulo);
	for (int j = 0; j < width; j+=bacground_cell_size)
	future_state(0, j) = mod((2*j)/bacground_cell_size, modulo);

	for (int p = 0; p < points_count; p++) {
	int i, j;
	int di = point[p].di;
	int dj = point[p].dj;
	while (1) { // failback loop
	i = point[p].i + di;
	j = point[p].j + dj;

	// canvas edge detection
	if ((i >= height) \|\| (i < 0)) {
	di = direction[random_direction_index()];
	continue;
	}
	if ((j >= width) \|\| (j < 0)) {
	dj = direction[random_direction_index()];
	continue;
	}

	// chaotic moving
	if (zeroone() < 0.01)
	di = direction[random_direction_index()];
	if (zeroone() < 0.01)
	dj = direction[random_direction_index()];
	if (di == 0 && dj == 0)
	continue;
	break;
	} // retry ...
	assert(i >= 0);
	assert(j >= 0);
	assert(i < height);
	assert(j < width);

	// new point position
	point[p].i = i;
	point[p].j = j;
	// new point position increment
	point[p].di = di;
	point[p].dj = dj;

	int ci, cj;
	const int edge = point[p].edge;
	int v;
	// draw a diamond. avoid glitches
	for (int c = 0; c < edge; c++) {
	v = c/13; // smooth colors

	// left top edge
	ci = i + c;
	cj = j + edge - c - 1;
	if (ci >= 0 && ci < height && cj >= 0 && cj < width)
	future_state(ci, cj) = mod(future_state(ci, cj) + (v), modulo);

	// left bottom edge
	ci = i + 2*edge - c - 1;
	cj = j + edge - c - 1;
	if (ci >= 0 && ci < height && cj >= 0 && cj < width)
	future_state(ci, cj) = mod(future_state(ci, cj) + (-v), modulo);

	// right top edge
	ci = i + c;
	cj = j + edge + c;
	if (ci >= 0 && ci < height && cj >= 0 && cj < width)
	future_state(ci, cj) = mod(future_state(ci, cj) + (-v), modulo);

	// right bottom edge
	ci = i + 2*edge - c - 1;
	cj = j + edge + c;
	if (ci >= 0 && ci < height && cj >= 0 && cj < width)
	future_state(ci, cj) = mod(future_state(ci, cj) + (v), modulo);
	}
	} // p = ...

	state.clear();

	// Backwards an automata rule: future_state(i0, j0) == state(i0, j0) - state(i0, j1) - state(i1, j0) + state(i1, j1) (mod modulo).
	// It is the rule 15 [Khan1997] when modulo == 2.
	// It is the unambiguous when we know initial values of the state (state(0, k) == 0 and state(k, 0) == 0).
	// Reference:
	// [Khan1997] A.R. Khan, et.al., "VLSI architecture of a cellular automata machine", 1997,
	// http://dx.doi.org/10.1016/S0898-1221(97)00021-7
	unsigned int i0, i1, j0, j1, solved_state;
	for (i0 = 0; i0 < state.size1()-1; i0++) {
	i1 = i0 + 1;
	for (j0 = 0; j0 < state.size2()-1; j0++) {
	j1 = j0 + 1;
	solved_state = mod(-state(i0, j0) + state(i0, j1) + state(i1, j0) + future_state(i0, j0), modulo);
	assert(solved_state >= 0);
	assert(solved_state < modulo);
	state(i1, j1) = solved_state;
	}
	}
	// produce image
	char fn[100];
	sprintf(fn, "rand%d-%dx%d-frame%05d.png", modulo, width, height, frame);
	save_matrix(state, fn, palette);
	}
	return 0;
	}