upto(n) = {my(res = List(), dk = 10); seed = 5; dk = 10; vals = List(); for(i = 1, n, c = is(i, seed); if(c > 0, listput(res, i); listput(vals, c); if(#vals > 3, dk = ceil((vals[#vals] - vals[#vals - 1]) / (res[#res] - res[#res - 1])) ); seed = c + dk; ); seed += dk ; ); res } is(n, {k = (n*log(n)/log(2)+5)\1}) = my(b=1. + 1 / n); while(floor(b^k) > k && k>1, k--); return((floor(b^k)==k)*k)