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]))
		);
	);
}

\\ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

neighbours = Set()

{
	nmx = 50;
	for (x=-nmx, +nmx,
		my (ym=sqrt int(nmx^2 - x^2));
		neighbours = set union(
			neighbours,
			Set(vector(2*ym+1, y, z = x+(y-ym)*I; [norm(z), z]))
		);
	);
}

ppos = 0
pos = 0
nor = 0

other(n) = {
	my (old = pos - ppos);
	ppos = pos;

	if (n==1,	pos = 0,
	    n==2,	pos = 1,
		      	my (best = 0, best n = 0, best i = 0);
			  	for (ne = 2, #neighbours,
					my (new = neighbours[ne]);
					if (best n && best n < new[1],
						break;
					);
					my (ang = new[2] / old);
					if (real(ang)>=0 && imag(ang)>0,
						if (!seen(pos + new[2]),
							if (best n==0 || best i > imag(ang),
								best   = new[2];
								best n = new[1];
								best i = imag(ang);
							);
						);
					);
			  	);
			  	if (best,
					nor = best n;
					pos += best,
					print ("** oops"); quit
			  	);
	);

	see(pos);
}

{
	for (n=1, 12 000,
		other(n);
		print (n " " real(pos));
	)
}

quit