domingo, 16 de junio de 2013

Programa Criba

b=6;q=37;c=6;m18=0.82;k0=0;m100=10^12;k1=121*10^7;i=1;d=matrix(primepi(q),c);forprime(p=2,q,a=b%p;d[i,1]=a;for(j=2,c,d[i,j]=a^2%p;a=a^2%p);i++);for(m1=1,primepi(q),for(m2=1,c-1,for(m3=m2+1,c,if(d[m1,m3]==d[m1,m2],d[m1,m3]=0))));e=matrix(primepi(q),c);for(m4=1,primepi(q),for(m5=1,c,e[m4,m5]=prime(m4)-d[m4,m5];if(m5==1&d[m4,m5]==0,e[m4,m5]=0);if(m5>1&d[m4,m5]==0,e[m4,m5]=0)));f=matrix(primepi(q),q);for(m6=1,primepi(q),m7=0;for(m8=1,q,f[m6, m8]=m7;m7++);for(m9=1,c,for(m10=2,prime(m6),if(f[m6,m10]==e[m6,m9],f[m6,m10]=q)));m11=0;for(m12=1,prime(m6),if(f[m6,m12]==q,m11++)));g=matrix(primepi(q),q);for(m13=1,primepi(q),g[m13, ]=vecsort(f[m13, ]));h=matrix(primepi(q),q);for(m15=1,primepi(q),m14=0;for(m16=1,prime(m15),if(g[m15, m16]<prime(m15),m14++));m17=m14/prime(m15);if(m17<m18,h[m15, ]=g[m15, ] ));for(m40=1,3,for(m41=1,prime(m40),if(b%prime(m40)==0,h[m40,m41]=m41)));m19=0;m20=vector(primepi(q));for(m21=1,primepi(q),if(h[m21,2]!=0,m19++;m20[m19]=m21));m22=vector(m19);for(m23=1,m19,m22[m23]=m20[m23]);m24=vector(m19-1);for(m25=2,m19,m26=0;while(h[m22[m25],m26+1]<prime(m22[m25]),m26++);m24[m25-1]=m26);u3=vector(m24[1]*m24[2]);i3=1;for(i1=1,m24[1],for(i2=1,m24[2],u3[i3]=lift(chinese(Mod(h[m22[2],i1],prime(m22[2])),Mod(h[m22[3],i2],prime(m22[3]))));i3++));m27=prime(m22[2])*prime(m22[3]);u5=vector(#(u3)*m24[3]);i4=1;for(i5=1,#(u3),for(i6=1,m24[3],u5[i4]=lift(chinese(Mod(u3[i5],m27),Mod(h[m22[4],i6],prime(m22[4]))));i4++));m28=m27*prime(m22[4]);u7=vector(#(u5)*m24[4]);i7=1;for(i8=1,#(u5),for(i9=1,m24[4],u7[i7]=lift(chinese(Mod(u5[i8],m28),Mod(h[m22[5],i9],prime(m22[5]))));i7++));m29=m28*prime(m22[5]);u9=vector(#(u7));i10=1;for(i11=1,#(u7),u9[i10]=lift(chinese(Mod(u7[i11],m29),Mod(h[1,1],2)));i10++);u10=vecsort(u9);m32=m29*2;u11=vector(#(u9)*m24[5]);i13=1;for(i14=1,#(u9),for(i15=1,m24[5],u11[i13]=lift(chinese(Mod(u9[i14],m32),Mod(h[m22[6],i15],prime(m22[6]))));i13++));u12=vecsort(u11);m34=m32*prime(m22 [6]);u13=vector(#(u11)*m24[6]);i17=1;for(i18=1,#(u11),for(i19=1,m24[6],u13[i17]=lift(chinese(Mod(u11[i18],m34),Mod(h[m22[7],i19],prime(m22[7]))));i17++));u14=vecsort(u13);m35=m34*prime(m22[7]);w1=vector(#(u13));for(i16=1,#(u14)-1,w1[i16]=u14[i16+1]-u14[i16]);w1[#(w1)]=m35-u14[#(u14)]+u14[1];m33=#(u14)/m35;if(m19<8,forstep(k=u14[1]+k0-k0%m35,m100,w1,m30=0;for(n=0,c-1,m=b^2^n+k;if(ispseudoprime(m),m30++;if(m30>c-1,print1(k,", ");write(p6,k," ");if(k>k0+k1,break(2)))))),if(m19=8,u15=vector(#(u13)*m24[7]);i21=1;for(i22=1,#(u13),for(i23=1,m24[7],u15[i21]=lift(chinese(Mod(u13[i22],m35),Mod(h[m22[8],i23],prime(m22[8]))));i21++));u16=vecsort(u15);m36=m35*prime(m22[8]);w2=vector(#(u15));for(i20=1,#(u16)-1,w2[i20]=u16[i20+1]-u16[i20]);w2[#(w2)]=m36-u16[#(u16)]+u16[1];m37=#(u16)/m36);forstep(k=u16[1]+k0-k0%m36,m100,w2,m50=0;for(n=0,c-1,m=b^2^n+k;if(ispseudoprime(m),m50++;if(m50>c-1,print1(k,", ");write(p6,k," ");if(k>k0+k1,break(2))))));if(m19>8,u17=vector(#(u15)*m24[8]);i25=1;for(i26=1,#(u15),for(i27=1,m24[8],u17[i25]=lift(chinese(Mod(u15[i26],m36),Mod(h[m22[9],i27],prime(m22[9]))));i25++));u18=vecsort(u17);m38=m36*prime(m22[9]);w3=vector(#(u17));for(i24=1,#(u18)-1,w3[i24]=u18[i24+1]-u18[i24];w3[#(w3)]=m38-u18[#(u18)]+u18[1];m39=#(u18)/m38);forstep(k=u18[1]+k0-k0%m38,m100,w3,m60=0;for(n=0,c-1,m=b^2^n+k;if(ispseudoprime(m),m60++;if(m60>c-1,print1(k,", ");write(p6,k,"  ");if(k>k0+k1,break(2))))))));print();print();print([readvec(p6)[1],#(readvec(p6)),p6]);c1=7;write(p7, );for(i101=1,#(readvec(p6)),for(n1=c1-1,c1-1,m1=b^2^n1+readvec(p6)[i101];if(ispseudoprime(m1),write(p7,readvec(p6)[i101]," "))));if(#(readvec(p7))>0,print();print([readvec(p7)[1],#(readvec(p7)),p7]));system("del p6");c2=8;write(p8, );for(i102=1,#(readvec(p7)),for(n2=c2-1,c2-1,m2=b^2^n2+readvec(p7)[i102];if(ispseudoprime(m2),write(p8,readvec(p7)[i102]," "))));if(#(readvec(p8))>0,print();print([readvec(p8)[1],#(readvec(p8)),p8]));system("del p7");c3=9;write(p9, );for(i103=1,#(readvec(p8)),for(n3=c3-1,c3-1,m3=b^2^n3+readvec(p8)[i103];if(ispseudoprime(m3),write(p9,readvec(p8)[i103]," "))));if(#(readvec(p9))>0,print();print([readvec(p9)[1],#(readvec(p9)),p9]));system("del p8");