N:= 30: # for terms of up to 2*N digits
R:= {1}: T:= {1, 9}:
for d from 2 to N do
T:= select(t -> convert(convert(t^2 mod 10^d, base, 10), set) subset {0, 1, 4, 9}, map(t -> (t, t + 10^(d-1), t + 4*10^(d-1), t + 9*10^(d-1)), T));
R:= R union select(t -> convert(convert(t^2, base, 10), set) subset {0, 1, 4, 9}, T);
od:
R2:= map(t -> t^2, R):
Res:= map(t -> seq(t*10^(2*i), i=0..(2*N-ilog10(t)-1)/2), R2) union {0}:
