// 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-minimum-r/ 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-minimum-r/%.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() {

	int target=5;

	double minincrement=0.0000001;

	for( ; target < 50; target++) {
		double r;
		double increment;
		for(r=1.0; r<10; r+=0.01) {
			/* Thanks to David W. Cantrell for the approximation */
			double predicted = PI*r*r - 2.9*r + 1.5;
			if(predicted > target) {
				increment = r/10.0;
				break;
			}
		}
		fprintf(stderr, "Target: %d.  Starting at r=%f\n", target, r);

		int maxsquares;
		double besttheta;
		while(increment > minincrement) {
			maxsquares = 0;
			besttheta = 0;

  		rsquared = r*r;
  		rsquaredfudge = rsquared+epsilon;
  
			double theta=0;
  		for(theta=0; theta < PI/2; theta += 0.000001) {
  			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 < target) {
				r += increment;
			} else {
				r -= increment;
			}

			fprintf(stderr, "   target=%d, r=%f, maxsquares=%d\n",
				target, r, maxsquares);

			increment /= 2.0;

			if((increment < minincrement) && (maxsquares < target)) {
				increment *= 2.0;
			}
  	}
		fprintf(stderr, "target=%d, mininum r=%f\n", target, r);
		pack(r,besttheta,1, maxsquares);
	}
}