
/*
	ulam.c - Generates a PBM of an Ulam spiral, up to 46340 x 46340 pixels
	Copyright (C) 2007  Guillem Cantallops Ramis <guillem-at-cantallops-dot-net>

	This program is free software; you can redistribute it and/or
	modify it under the terms of the GNU General Public License
	as published by the Free Software Foundation; either version 2
	of the License, or (at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

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

unsigned int dim, x, y, s, i, m, step, cont, max;
unsigned char *table;
#define T(X,Y) table[(X-1)+dim*(Y-1)]

inline unsigned char isprime (unsigned int n)
{
	if (n <= 1) return 0;
	if (n == 2) return 1;
	if (n % 2 == 0) return 0;
	m = (unsigned int)round (sqrt (n));
	for (i = 3; i <= m; i += 2)
	{
		if (n % i == 0)
		{
			return 0;
		}
	}
	return 1;
}

int main (int argc, char **argv)
{

	if (argc != 2) return 1;
	dim = atoi (argv[1]);
	if ((dim < 0) || (dim > 46340)) return 1;
	dim |= 1;
	max = dim * dim;
	table = calloc (max, sizeof(unsigned char));

	x = y = dim / 2 + 1;
	step = cont = 1;
	while (step <= dim)
	{
		for (s = 0; cont <= max && s < step; s++)
		{
			T(x,y) = isprime (cont);
			x++;
			cont++;
		}
		for (s = 0; cont <= max && s < step; s++)
		{
			T(x,y) = isprime (cont);
			y--;
			cont++;
		}
		step++;
		for (s = 0; cont <= max && s < step; s++)
		{
			T(x,y) = isprime (cont);
			x--;
			cont++;
		}
		for (s = 0; cont <= max && s < step; s++)
		{
			T(x,y) = isprime (cont);
			y++;
			cont++;
		}
		step++;
	}

	printf ("P1\n%d %d\n", dim, dim);
	for (s = 0; s < max; s++)
	{
		printf ("%1d", table[s]);
		if (! (s % 69))
		{
			printf ("\n");
		}
	}


	free (table);
	return 0;
}

