\\ z=a+b*I ; sign of a+b*sqrt(5)
xsign(z) = {
	if (z==0, return (0),
	    real(z) <= 0 && imag(z) <= 0, return (-1),
		real(z) >= 0 && imag(z) >= 0, return (+1),
		real(z) >= 0, sign(real(z)^2 - 5*imag(z)^2),
		              sign(5*imag(z)^2 - real(z)^2)
	);
}

\\ phi^n as complex
xphi(n) = {
	if (n>=0,
		sum (k=0, n, binomial(n,k) * 5^(k\2) * if (k%2==0, 1, I))/2^n,
		sum (k=0, -n, (-1)^(-n-k) * binomial(-n,k) * 5^(k\2) * if (k%2==0, 1, I))/2^-n
	);
}

tophi(n) = {
	my (v=0, p);
	for (i=0, oo,
		p=xphi(i);
		if (xsign(n-p)<=0,
			forstep (j=i, -oo, -1,
				if (n==0,
						return (v),
					xsign(n-p=xphi(j))>=0,
						v+=2^j;
						n-=p;
				);
			);
		);
	);
}

s = 0
unseen = 0
seen(v) = bit test(s, v)
see(v) = s = bit or(s, 2^v); while (seen(unseen), unseen++)

cache = apply(tophi, [0..15 000])

other(p) = {
	my (pp=cache[1+p]);
	for (v=unseen, oo,
		if (!seen(v),
			my (vv=cache[1+v], m=max(denominator(pp), denominator(vv)));
			if (bitand(pp*m, vv*m)==0,
				return (v);
			);
		);
	);
}

{
	for (n=0, 10 000,
		print (n " " v=other(n));
		see(v);
	);
}

quit