```#PARI/GP, 4127 digits

(10<sup>4127</sup>-1)/9 + 2*10<sup>515</sup>

This is a fairly straightforward search: check only prime digit lengths, then compute the possible primes to use, then iterate through all possibilities. I special-cased the common cases where there are 0 or 1 suitable prime digits to use.

supreme(lim,startAt=3)={
forprime(d=startAt,lim,
my(N=10^d\9, P=select(p->isprime(d+p),[1,2,4,6]), D, n=1);
if(#P==0, next);
if(#P==1,
for(i=0,d-1,
if (ispseudoprime(D=N+n*P), print(D));
n*=10
);
next
);
D=vector(#P-1,i,P[i+1]-P[i]);
for(i=0,d-1,
forstep(k=N+n*P,N+n*P[#P],n*D,
if (ispseudoprime(k), print(k))
);
n*=10
)
)
};
supreme(4200, 4100)

This took 36 minutes to compute on one core of a fairly old machine. It would have no trouble finding such a prime over 5000 digits in an hour, I'm sure, but I'm also impatient.

A better solution would be to use any reasonable language to do everything but the innermost loop, then construct an abc file for <a href="http://openpfgw.sourceforge.net/">primeform</a> which is optimized for that particular sort of calculation. This should be able to push the calculation up to at least 10,000 digits.```