\\ find best suffix index findsuffix(m) = { forstep (k=#M, 1, -1, if (m%M[k]==P[k], return (k); ); ); print ("bug"); quit } \\ initialize for "n" init(n) = { my (w=#binary(n)); M = vector(1+w, k, 2^(k-1)); P = vector(1+w, k, n\2^(w-k+1)); P0 = apply(v -> findsuffix(2*v+0), P); P0[#P]=#P; P1 = apply(v -> findsuffix(2*v+1), P); P1[#P]=#P; } \\ add "0" or "1" any(v) = { my (w=vector(#P)); for (k=1, #P, w[P0[k]] += v[k]; w[P1[k]] += v[k] ); w } \\ find best suffix (or #P if contains "n") findbest(b) = { my (m=0); for (k=1, #b, m = 2*m+b[k]; if (findsuffix(m)==#P, return (#P); ); ); findsuffix(m); } \\ count up to "lim" count(lim) = { my (nb=0, b=binary(lim)); for (k=1, #b, if (b[k], b[k]=0; my (v=vector(#P)); v[findbest(b[1..k])]++; for (x=k+1, #b, v=any(v); ); nb += v[#P]; b[k]=1 ); ); if (findbest(b)==#P, nb++; ); nb } a(n) = { init(n); my (k=n, i=1); while (2*i < k, i = count(k = 2*k-2*i); ); return (k); } for (n=1, 256, print (n " " a(n))) quit