\\ 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..12 000]) other(p) = { see(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=if (n==0, 0, other(v))); ); } quit