/*****************************
 * Dyadic transformation Plot
 * (C)2010 Claudio Rocchini
 * [email protected]
 *****************************/

#include <stdio.h>
#include <assert.h>
#include <vector>
#include <set>

typedef unsigned char byte;

class image {			// image definition
public:
	byte * mbuf;		// Image memory
	int    dimx,dimy;	// Image size
	 image( int ndimx, int ndimy ) : dimx(ndimx),dimy(ndimy) {
		 mbuf = new byte[dimx*dimy*3];
		 memset(mbuf,255,dimx*dimy*3);
	 }
	~image() { delete[] mbuf; }
	bool save_ppm( const char * filename ) const {
		FILE * f = fopen(filename,"wb"); if(f==0) return false;
		fprintf(f,"P6\n%d %d\n255\n",dimx,dimy);
		fwrite(mbuf,3,dimx*dimy,f);
		fclose(f);
		return true;
	}
	void set_black_pixel( int x, int y, double trasp )
	{
		if(x<0 || x>=dimx) return; if(y<0 || y>=dimy) return;
		byte * p = mbuf 3*(x y*dimx);
		p[0] = byte(p[0]*(1-trasp) /*  (trasp)*0*/);
		p[1] = byte(p[1]*(1-trasp) /*  (trasp)*0*/);
		p[2] = byte(p[2]*(1-trasp) /*  (trasp)*0*/);
	}
	void dot( double x, double y, double trasp )	// float precision put pixel
	{
		int ix = int(x); double rx = x-ix;
		int iy = int(y); double ry = y-iy;
		set_black_pixel(ix 0,iy 0,trasp*(1-rx)*(1-ry));
		set_black_pixel(ix 1,iy 0,trasp*(  rx)*(1-ry));
		set_black_pixel(ix 0,iy 1,trasp*(1-rx)*(  ry));
		set_black_pixel(ix 1,iy 1,trasp*(  rx)*(  ry));
	}
};

	// rational definition
typedef __int64 bigint;

bigint gcd(bigint m, bigint n) {
	while (n != 0){
		bigint t = m % n;
		m = n; n = t;
	}
	return m;
}

class rational {
public:
	bigint p,q;

	rational(){}
	rational( int np, int nq ) { p=np; q=nq; normalize(); }
	bool operator< ( const rational & x ) const { return p*x.q < q*x.p; }
	void normalize()  { bigint g = gcd(p,q); p /= g; q /= g; }
	void doubled()    {
		if( q%2==0 ) { assert(q!=1);        q /= 2; }		// assert for overflow checking
		else         { assert(p< (1<<30) ); p *= 2; }
	}
	void modulo()     { while( p >= q ) p -= q; normalize(); }
	operator double() { return double(p)/q; }
};

int main() {
	const int SX = 600;
	const int SY = 600;
	image im(SX,SY);
	std::set< rational > ss;					// fraction done
	for(int i=2;i<150;  i) {					// enumerate the fractions
		for(int j=1;j<i;  j) {
			rational w(j,i);				// current start fractions
			if( ss.find(w)!=ss.end() ) continue;		// just do it!
			ss.insert(w);
			double xx = SX*w;
			std::set< rational > tt;			// row set
			for(;;){
				if( tt.find(w)!=tt.end() ) break;	// just do it: is a loop
				tt.insert(w);
				double yy = SX*w;
				im.dot(xx,yy,0.5);			// dot this value
				w.doubled();				// ww = (w*2) mod 1
				w.modulo();
			}
		}
	}
	im.save_ppm("c:\\temp\\dyadic.ppm");
	return 0;
}