allocate mem(2^30) big = 1 000 000 s = 0 S = Set([]) unseen = 1 seen(v) = if (v < big, bit test(s, v), set search(S, v)) see(v) = if (v < big, s = bit or(s, 2^v), S = set union(S, Set([v]))); while (seen(unseen), unseen++) \\ describe a fraction: \\ - integer part \\ - transient part \\ - periodic part frac digits(v, base=10) = { my (int=digits(floor(v), base)); \\ period length my (den=denominator(v), pr=factor(base)[,1]~, p=znorder(Mod(base, den / prod(i=1, #pr, pr[i]^valuation(den,pr[i]))))); \\ transient part my (tra=[], fra=frac(v)); while (fra != frac(fra*base^p), tra = concat(tra, floor(fra*base)); fra = frac(fra*base); ); \\ periodic part my (per = vector(p, x, my (d=floor(fra*base)); fra = frac(fra*base); d)); return ([int,tra,per]); } \\ build fraction back from integer, transient and periodic parts from frac digits(f, base=10) = { fromdigits(f[1],base) + ( fromdigits(f[2],base) + fromdigits(f[3],base)/base^#f[3] / (1-1/base^#f[3]) ) / base^#f[2]; } other(p) = { see(p); my (fp = fracdigits(1/p, 2)); my (mp = floor(2^100/p)); for (v=unseen, oo, if (!seen(v) && !bitand(mp, floor(2^100/v)), my (ok=1); my (fv = fracdigits(1/v, 2)); my (t = max(#fp[2], #fv[2])); if (bitand(floor(2^t/p), floor(2^t/v)), ok = 0, my (w=lcm(#fp[3], #fv[3])); for (i=1, w, if (fp[3][1+((i+#fv[2])%#fp[3])] && fv[3][1+((i+#fp[2])%#fv[3])], ok = 0; break; ); ); ); if (ok, return (v); ); ); ); } v = 1; for (n=1, 2000, print (n " " v); v=other(v)) quit