\\ Experimental prime cluster search, save as, e.g., paricode.gp
\\ Language PARI/GP, freeware from  http://pari.math.u-bordeaux.fr/
\\ Call PARI/GP by gp, execute the file by \r paricode.gp  
\\ or just by the call main(...), once the procedures are loaded

{loop(np,ind,class,b0,nbl) = if(ind<=np, for(k=1,len[ind],
   class[ind+1]=(class[ind]+term[ind,k])%P;
   loop(np,ind+1,class,b0,nbl)), test(class[ind]%P,b0,nbl))}

{test(class,b0,nbl) = sieve=vector(nbl,i1,1);
forprime(p=pn+1,q,
   for(i1=1,nc,st=component(Mod(-b0-(class+c[i1])/P,p),2);
   forstep(i2=st+1,nbl,p,sieve[i2]=0)
));
for(i2=1,nbl,
   if(sieve[i2]==1,x=class+(b0+i2-1)*P; v=1; i1=1;
      while(v&i1<=nc, v=ispseudoprime(abs(x+c[i1])); i1=i1+1);
      if(v, count=count+1; data[count]=x)
   )
)}

{main(d,np,q,b0,nbl) =
\\ d: difference pattern (vector);  np: number of primes in Chinese RMT
\\ q: limit for primes in sieve; b0, nbl determine the search interval
\\ [b0*P, (b0+nbl)*P] where P=2*3*5*...*p_np
c=concat(0,d); nc=length(c); for(j=2,nc, c[j]=c[j-1]+c[j]); data=vector(300);
\\ Find vacant residue classes rem and length vector len
pp=primes(np); len=0*pp; pn=pp[np]; col=pn-nc+3; rem=matrix(np,col,i,j,0);
i=0; forprime(p=2,pn, i=i+1; l=0; for(j=0,p-1,\
     if(prod(k=1,nc, c[k]%p-j)!=0, l=l+1;rem[i,l]=p-j)); len[i]=l);
\\ Chinese remainders
P=prod(j=1,np,pp[j]); p1=vector(np,j,P/pp[j]);
s=vector(np,j,component(Mod(1/p1[j],pp[j]),2)); term=0*rem;
for(j=1,np,\
   for(k=1,len[j], f=(s[j]*rem[j,k])%pp[j]; term[j,k]=f*p1[j]);\
); print("Pattern ",d,",   Interval  ",b0*P,"   ",(b0+nbl)*P);
\\ Nested loops
count=0; loop(np,1,vector(np+1),b0,nbl);
res=vecsort(vector(count,k, data[k])); res}

\\ Test examples, calling statement (timings: 3 GHz processor)
main([2,4,2],           7,641,    0,    1)   \\  .17 sec, 101 sets
main([4,2,4,2,4],       8, 53,    0,   31)   \\  .44 sec, 159 sets
main([2,4,2,4,6,2],     8, 53, -103,  206)   \\  .71 sec, 102 sets
main([2,4,6,2,6,4,2],   8, 53,    0, 1031)   \\  4.1 sec, 101 sets
main([2,4,2,4,6,2,6],   8, 61,-1031, 2062)   \\  3.1 sec,  80 sets
main([2,4,2,4,6,2,6,4], 9,101,-1345, 2690)   \\ 42.6 sec, 101 sets
main([4,2,4,6,2,6,4,2], 9,101,-1345, 2690)   \\ 85.5 sec, 189 sets
main([2,4,2,10,2,10,2,4,2], 9,199,0, 2690)   \\ 19.8 sec,   6 sets

\\ Next 3 calls yield 71542018620258822095351
main([2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2], 8,997, 7375701555437217,1) \\0.01 sec
main([2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2],10,313,   11058023321495,1) \\0.26 sec
main([2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2],11,101,     356710429725,1) \\2.87 sec

\\ Yield 49380784584721484163799421  and  52739555754897822763982917
main(  [2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2],11,101,246213920561889,1) \\ 2.9 sec
main([4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2],12, 71,  7107049837468,1) \\35.8 sec

\\ Long examples, (-73,...,-13) and (13,...,83), > 1 hour per block
\\main([2,4,6,2,6,6,4,2,4,6,2,6,4,2,4],     14, 541,      0, 30575)
\\main([4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4], 16, 541, -10000, 20000)