// Jason E. Holt <web (at) lunkwill [dot] org>
// Works out how many unit squares can fit into a circle with a certain
// radius, and outputs SVG graphics of circles that contain more squares
//
// This code released into the public domain
// Compile with gcc -lm 
// Make sure you have an svg/ directory if you want svg images

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

#define PI 3.14159

double epsilon = 0.000000001;
double rsquared;
double rsquaredfudge;

int inside(double x, double y) {
//	printf("%f,%f = %f < %f? %d\n", x,y, x*x+y*y, rsquared, 
//		x*x+y*y <= (rsquared+epsilon) );
	return x*x+y*y < rsquaredfudge;
}

// Pack unit squares into circle of radius r, starting at theta
// radians to the left of the y axis.  If makesvg, output an .svg
// graphic of the packing; filename will include, r, theta and
// maxsquares (which you presumably computed previously before deciding
// to output a graphic of the winning combination)
int pack(double r, double theta, int makesvg, int maxsquares) {
			int squares=0;
			double x = -sin(theta)*r;
			double y = sqrt(rsquared-x*x);

			FILE *f;
			if(makesvg) {
				char fn[255];
				sprintf(fn, "svg/%.9f-%d-%.9f.svg", r, maxsquares, theta);

  			f = fopen(fn, "w");
  
    		fprintf(f, "<svg width=\"30cm\" height=\"30cm\" version=\"1.1\" ");
    		fprintf(f, "  xmlns=\"http://www.w3.org/2000/svg\">\n\n");
    		fprintf(f, "<g transform=\"translate(500,500)\">\n");
    		fprintf(f, "<g transform=\"scale(1,-1)\">\n");
    		fprintf(f, "<ellipse stroke=\"0\" ");
  			fprintf(f, " cx=\"0\" cy=\"0\" ry=\"%fcm\" rx=\"%fcm\"/>\n", r, r);
			}
  
			while(inside(x,y)) {
				//printf("theta=%f, (%f,%f) (r=%f)\n", theta, x, y, sqrt(x*x+y*y));
				while(
					((y > 0.5) && inside(x+1.0, y)) ||
					((y <= 0.5) && inside(x+1.0, y-1.0)) ) {

					squares++;

					if(makesvg) {
  					fprintf(f, "<rect x=\"%fcm\" y=\"%fcm\" ", x, y-1);
  					fprintf(f, " style=\"fill:rgb(0,0,255);stroke-width:1\" ");
  					fprintf(f, " width=\"1cm\" height=\"1cm\"/>\n");
					}
  
					x+=1.0;
				}
				y-=1.0;
				x = -sqrt(rsquared - y*y);
				double bottomx = -sqrt(rsquared - (y-1.0)*(y-1.0));
				if(bottomx > x) x = bottomx;
			}

			if(makesvg) {
    		fprintf(f, "</g></g>\n");
    		fprintf(f, "</svg>\n");
    		fclose(f);
			}
  
//			fprintf(stderr,"Theta=%f: %d squares, r=%f\n",
//					theta, squares, r);

	return squares;
}

main() {

	double minr=1.0;
	double r;
	double maxr=20.00;
	double theta=0;
	int lastmaxsquares = 0;
	for(r=minr; r<maxr; r+=1.0) {
		int maxsquares=0;
		double besttheta=0;
		rsquared = r*r;
		rsquaredfudge = rsquared+epsilon;

		for(theta=0; theta < PI/2; theta += 0.0000001) {
			if(cos(theta)*r < r-2.0) {
				fprintf(stderr, "Giving up at theta=%f.\n", theta);
				break;
			}
			int squares = pack(r, theta, 0, 0);
			if(squares > maxsquares) {
				fprintf(stderr,"New record: %d squares, theta=%f, r=%f\n",
					squares, theta, r);
				maxsquares = squares;
				besttheta = theta;
			}
		}
		if(maxsquares > lastmaxsquares) {
			pack(r,besttheta,1, maxsquares);
			lastmaxsquares = maxsquares;
		}
	}
}