big limit = 100 000
big values = Set([])
small values = vector(4, i, [])

see index(z) = {
	my (s=1, x=real(z), y=imag(z));
	if (x<1,
		s += 1;
		x = 1-x;
	);
	if (y<0,
		s += 2;
		y = -y;
	);
	return ([s,x,y]);
}

see(z) = {
	my (index = see index(z));
	if (max(index[2],index[3]) > big limit,
		big values = set union(big values, Set([z])),
		if (#small values[index[1]] < index[2],
			small values[index[1]] = concat(small values[index[1]], vector(index[2] - #small values[index[1]]))
		);
		small values[index[1]][index[2]] = bit or(small values[index[1]][index[2]], 2^index[3])
	);
}

seen(z) = {
	my (index = see index(z));
	if (max(index[2],index[3]) > big limit,
		return (set search(big values, z)),
		if (#small values[index[1]] < index[2],
			return (0),
			return (bit test(small values[index[1]][index[2]], index[3]))
		);
	);
}

cell(z) = (real(z) - 2*floor(imag(z)/4))%14 + I*(imag(z)%4);

forbidden = [ [ 1,      1 ],		\
              [ 2,      I ],		\
              [ 2+I,    -1+I ],		\
              [ 1+2*I,  -1 ],		\
              [ 2*I,    -I ],		\
              [ I,      1-I ],		\
			  [ 6,      1 ],		\
			  [ 9+2*I,	1 ],		\
			  [ 10+2*I, I ],		\
			  [ 10+3*I, -1+I],		\
			  [ 9+2*I, -1+I ],		\
			  [ 8+3*I, +I ] ];
			  
forbidden = Set(concat(forbidden, apply(zd -> [cell(zd[1]+zd[2]), -zd[2]], forbidden)));

is hex center(z) = real(z)%2==1 && imag(z)%2==1

accept(z, d) = {
	if (is hex center(z) || is hex center(z+d), return (0),
		set search(forbidden, [cell(z), d]),  return (0),
		return (1))
}

dz = [ +1, +I, -1+I, -1, -I, +1-I ];
neighbours(z) = apply(d -> z+d, select(d -> accept(z,d), dz));

\\ coordinates for rendering
flat(z) = 2*real(z) + (1+2*I) * imag(z);

\\  0 -> A320495
\\  I -> A320496
\\ -I -> A320497
\\ -1 -> A320498

current = Set(I)

{
	for (r=0, 1 000,
		print (r " " #current);
		apply(see, current);

		current = Set(select(z -> !seen(z), Set(concat(apply(neighbours, current)))));
	);
}

quit