\\ 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)