/*****************************  * 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; }