uptoqdigits(n) = { my(res = List()); for(i = 12, n, w = withqdigits(i); for(j = 1, #w, listput(res, w[j]); ); print(i" "res); ); res } withqdigits(n) = { my(res = List()); qdfromleft = ceil(n/2) + 1; for(i = 1, qdfromleft, l = uptowithquadratic(i, n, 2*n); for(j = 1, #l, listput(res, l[j]); ); l = uptowithquadratic(i, n, 2*n-1); for(j = 1, #l, listput(res, l[j]); ); ); qdfromright = n - qdfromleft; for(i = 1, qdfromright, l = uptowithmod(i, n); for(j = 1, #l, listput(res, l[j]); ) ); Set(res) } uptowithquadratic(u, m, {qd = 2*m}) = { my(res = List(), llim = 10^(u-1), ulim = 10*llim - 1, pow10shift = 10^(qd - m - u), pow10mod = 10^(qd - u)); for(i = llim, ulim, c = i * 10^(qd - u); k = ceil((pow10shift + sqrt(pow10shift^2 + 4*c))/2); if(k%10 != 0 && k * (k - pow10shift) % pow10mod < pow10shift, listput(res, k) ) ); res } uptowithmod(t, m) = { my(llim = 10^(m - 1), gen = [1..9]); for(i = 1, m, gen = nextgen(gen, t, i); gen = select(x->x<10*llim, gen); ); Set(select(x->llim < x && (x * (x - 10^t))%10^(t + m) < 10^t, gen)); } nextgen(gen, t, qd) = { my(res = List(), pow10 = 10^t, addpow10 = 10^qd); for(i = 1, #gen, for(j = 0, 9, c = gen[i] + j * addpow10; if((c * (c - pow10))%addpow10 < pow10, listput(res, c) ) ) ); res }