uptoqdigits(n) = { my(res = List([2,3]), i, j, gen); gen = [1,3,7,9]; for(i = 1, n, for(j = 1, #gen, if(gen[j] % 3 != 0, dc = digits(gen[j]^2); if(dc == vecsort(dc, , 4), if(isprime(gen[j]), listput(res, gen[j]); ); ) ) ); gen = nxtgen(gen, i); ); Set(res) } nxtgen(v, qd) = { my(res = List(), pow10 = 10^qd, c, dc, i, j); for(i = 0, 9, for(j = 1, #v, c = (i*pow10 + v[j]); dc = digits(c^2 % (10*pow10)); if((c^2 % (10*pow10) > pow10) && dc == vecsort(dc,,4), listput(res, c) ) ) ); res } is(n) = my(d = digits(n^2)); d == vecsort(d, , 4) && isprime(n)