/**************************************************************** Rutinas hall para compilar con gp2c-run 22-7-02 ****************************************************************/ \\ residuo x%m centrado res(x,m, r)= r=x%m;if(r>m\2,r=r-m,);r; /*** Raices cúbicas modulo p (cuberoots_p) y mod m (cuberoots). 14-6-02 ***/ /**** Metodo Niven ... "A Introduction to NT" p.115 ejer. 8 **/ { cuberoots_p(x,p)= local(h,z,zm,z2,ex,m,k,kk,kkk,r,n,nn,c,cc1,cc2,nc1,nc2,nc11,nc22,i); x%=p; if(x==0,return([0])); if(p<=3,return([x])); /* Caso p=3k+2 */ if(p%3==2,if(x==1,return([1]),return([lift(Mod(x,p)^((2*p-1)/3))]))); /* Caso p=3k+1 */ h=Mod(1,p);ex=(p-1)/3;until(zm!=Mod(1,p),h++;zm=h^ex); z=lift(zm);z2=z^2%p; \\print("z="z" zm="zm); /* Caso p=9k+7 */ if(p%9==7,r=lift(Mod(x,p)^((p+2)/9)); if(r^3%p!=x,return([]),return(vecsort([r,r*z%p,r*z2%p])))); /* Algoritmo en Niven ... p.115 ejer. 8*/ m=p-1;k=0;until(m%3!=0,m/=3;k++); if(m%3==2,r=Mod(x,p)^((m+1)/3);n=Mod(x,p)^m, r=Mod(x,p)^((2*m+1)/3);n=Mod(x,p)^(2*m)); c=h^m; /* c tiene orden 3^k */ kk=0;nn=n;while(nn!=Mod(1,p),kk++;nn=nn^3); /* n tiene orden 3^kk */ if(k==kk,return([])); while(n!=Mod(1,p),cc1=c;for(i=1,k-kk-1,cc1=cc1^3);nc1=nc11=n*cc1^3; cc2=c^2;for(i=1,k-kk-1,cc2=cc2^3);nc2=nc22=n*cc2^3; if(k==kk,return([])); kkk=0; while(1, if(nc11==Mod(1,p),n=nc1;c=cc1^3;r=r*cc1;break(),nc11=nc11^3); if(nc22==Mod(1,p),n=nc2;c=cc2^3;r=r*cc2;break(),nc22=nc22^3); kkk++; ); k=kk;kk=kkk; ); r=lift(r); return([r,r*z%p,r*z2%p]); } { cuberoots(x,m,modo=0)= local(f,j,k,np,np1,np3,ip1,ip3,r1,r3,p,e,r,raiz3,a,d,fa,rr,c,i); /*print([x,factor(m)]);*/ if(modo==0, /* modo default */ np1=np3=0;ip1=ip3=1; /* Contadores de r1(raices simples) y r3(raices triples) */ /* Primero se trata el caso p=3, el único no-singular. Si p=3,e=1 tenemos una raiz tanto si m=1,2. Si e>1 tenemos raices triples*/ e=0;while(m%3==0,m/=3;e++); if(e!=0, if(e==1,raiz3=[Mod(x%3,3)];np1++ , if(x%9==1|x%9==8,/* Caso de que exista raiz se eleva a mod 3^e */ a=x%3;c=1;d=3; for(i=2,e-1,c*=3;d*=3; raiz3=[a,(a+c)%d,(a+2*c)%d]; j=1;until(fa==0,a=raiz3[j];fa=(a^3-x)%(3*d);j++) ); c*=3;d*=3;raiz3=[Mod(a,d),Mod(a+c,d),Mod(a+2*c,d)];np3++; ,return([]) ) ) ,raiz3=0 ); f=factor(m);np=matsize(f)[1];/* Cuantas raices simples y triples? */ for(i=1,np,p=f[i,1];if(p%3==2,np1++,np3++)); r1=vector(np1);r3=vector(np3*3); /* Cargar raiz3 si la hubiere */ if(raiz3!=0, if(length(raiz3)==1,r1[ip1]=raiz3[1];ip1++ ,for(i=1,3,r3[i]=raiz3[i]);ip3+=3 )); for(i=1,np,p=f[i,1];e=f[i,2];r=cuberoots_p(x,p); if(r==[],return(r)); if(p%3==2,a=r[1];d=3*a^2%p;d=bezout(d,p)[1];c=p^e; /* raices simples en r1*/ for(k=2,e,a-=((a^3-x)*d)%c);r1[ip1]=Mod(a,c);ip1++ ); if(p%3==1, for(j=1,3,a=r[j];d=3*a^2%p;d=bezout(d,p)[1];c=p^e; /* raices triples en r3 */ for(k=2,e,a-=(((a^3-x)*d)%c) );r3[ip3+j-1]=Mod(a,c); );ip3+=3; ); ); rr=vector(3^np3); /* En rr todas las raices mod m */ /* Resto chino de r1 */ if(np1!=0,c=r1[1];for(i=2,np1,c=chinese(c,r1[i])),c=Mod(0,1)); /* Resto chino de r1+r3 */ rr[1]=c; for(i=0,np3-1,forstep(j=3,1,-1,forstep(k=3^i,1,-1,rr[3^i*(j-1)+k]=chinese(rr[k],r3[i*3+j]) )));rr=vecsort(lift(rr)); ,rr=[];for(i=0,m-1,k=i^3%m;if(k==x,rr=concat(rr,i)))); /* Modo 1 (naive) */ rr; } hall_k(x, y,k)=y=round(sqrt(x^3));k=x^3-y^2; hall_bb(a,b, bb)=bb=2*b;if(a%2==0,bb/=2);if(b%3==0,bb/=3);bb; hall_alfa(a,b, alfa)=alfa=res(a^2,b^2);alfa; { hall_a(b,C,i)= local(n,A,a,alfa,k,aa); if(gcd(2*C,b)!=1,return([])); /* a_0 (y también a) con b es reducible */ aa=[]; /* Casos */ if(b%2==0, if(b%3!=0, A=cuberoots(2*C,b^2);if(A==[],return([])); for(n=1,length(A),a=A[n];alfa=res(a^2,b^2); k=bezout(-3*alfa,2*b)[1]*(3*i*a+(2*a^3-3*alfa*a+2*C)/b^2)%(2*b); a+=k*b^2;a%=(2*b^3); aa=concat(aa,[[a,2]]) ); return(aa) ,A=cuberoots(2*C,3*b^2);if(A==[],return([])); for(n=1,length(A)/3,a=A[n];alfa=res(a^2,b^2); k=bezout(-alfa,2*b/3)[1]*(i*a+(2*a^3-3*alfa*a+2*C)/(3*b^2))%(2*b/3); a+=k*b^2;a%=(2*b^3/3); aa=concat(aa,[[a,2/3]]) ); return(aa) ) ,if(type(C)=="t_INT", if(b%3!=0, A=cuberoots(2*C,2*b^2);if(A==[],return([])); for(n=1,length(A),a=A[n];alfa=res(a^2,b^2); if((i+alfa)%2==1, k=bezout(-3*alfa,b)[1]*(3*i*a/2+(a^3-3*alfa*a/2+C)/b^2)%b; a+=2*k*b^2;a%=(2*b^3); aa=concat(aa,[[a,2]]); , k=bezout(-3*alfa/2,b)[1]*(3*i*a+(2*a^3-3*alfa*a+2*C)/b^2)/2%b; a+=k*b^2;a%=(b^3); aa=concat(aa,[[a,1]]); ) ); return(aa) , A=cuberoots(2*C,6*b^2);if(A==[],return([])); for(n=1,length(A)/3,a=A[n];alfa=res(a^2,b^2); if((i+alfa)%2==1, k=bezout(-alfa,b/3)[1]*(i*a/2+(2*a^3-3*alfa*a+2*C)/(6*b^2))%(b/3); a+=2*k*b^2;a%=(2*b^3/3); aa=concat(aa,[[a,2/3]]); , a%=(3*b^2); k=bezout(-alfa/2,b/3)[1]*(i*a/2+(2*a^3-3*alfa*a+2*C)/(6*b^2))%(b/3); a+=k*b^2;a%=(b^3/3); aa=concat(aa,[[a,1/3]]); ) ); return(aa) ) ,if(b%3!=0, A=cuberoots(2*C,2*b^2);if(A==[],return([])); for(n=1,length(A),a=A[n];alfa=res(a^2,b^2); if((alfa+i)%2==1, k=bezout(-3*alfa,b)[1]*(3*i*a+(2*a^3-3*alfa*a+2*C)/b^2)/2%(b); a+=k*2*b^2;a%=(2*b^3); aa=concat(aa,[[a,2]]) ) ); return(aa) , A=cuberoots(2*C,6*b^2);if(A==[],return([])); for(n=1,length(A)/3,a=A[n];alfa=res(a^2,b^2); if((alfa+i)%2==1, k=bezout(-alfa,b/3)[1]*(i*a/2+(2*a^3-3*alfa*a+2*C)/(6*b^2))%(b/3); a+=2*k*b^2;a%=(2*b^3/3); aa=concat(aa,[[a,2/3]]) ) ); return(aa) ) ) ) } /**************************************************************/ /******* Busqueda a partir de valores de b,C, i=w=0 !!! ******/ /**************************************************************/ { hall_search_bC00(bi,bf,guarda=1)= local(tiempo,b,b2,b3,C,C2,iC2,as,j,a0,s,alfa,n,t,a,x0,x,x2,y,kx,rk,dat,dats); print("x r t C tiempo(seg.)"); tiempo=0;gettime(); for(b=bi,bf,if(b%100000==0,print("b=",b));b2=b^2;b3=b*b2; if(b%2==0,iC2=2,iC2=1); forstep(C2=1,sqrtint(sqrtint(b)),iC2, /* Ojo!, poner limite */ if(gcd(C2,b)==1,C=C2/2; as=hall_a(b,C,0); for(j=1,length(as),a0=as[j][1];s=as[j][2]; alfa=hall_alfa(a0,b); n=round((3/8*alfa^2/C-a0)/s/b3); a=a0+n*s*b3;t=a/b;x0=round(t^2);x=x0; x2=sqrt(x);y=round(x2*x);kx=x^3-y^2; if(kx!=0&abs(kx)> Hall00.dat"); if(guarda,system(dats) ) ) ) ) ) ) } /******* Busqueda a partir de valores de b,C,i ******/ /*****************************************************/ { hall_search_bCi(bi,bf,guarda=1)= local(wi,wf,tiempo,b,c,b2,b3,C,C2,iC2,as,i,j,a0,s,alfa,Cs,A,B,u,CC,r,bb,w,rw,a, n,t,x0,x,x2,y,kx,rk,dat,dats); wi=-4;wf=4; /* Rango de w */ print("x r t C i w time(sec.)"); tiempo=0;gettime(); { forstep(b=bi,bf,1,if(b%3==0,c=3,c=1);if(b%1000==0,print("b=",b));b2=b^2;b3=b*b2; if(b%2==0,iC2=2,iC2=1); forstep(C2=1,sqrtint(sqrtint(b)),iC2, /* Rango de C hasta b^(1/4) */ if(gcd(C2,b)==1,C=C2/2; for(i=0,2*b/c-1, /* Rango de i hasta 2*b' */ as=hall_a(b,C,i); for(j=1,length(as),a0=as[j][1];s=as[j][2]; alfa=hall_alfa(a0,b);Cs=3./(8.*C*s); A=b3*Cs;u=b2*i-alfa;B=2*u*Cs;CC=(3*u^2/8/C-a0)/(s*b3); if(a0%2==1,r=2/c,r=1/c);bb=b*r;CC*=1.; for(w=wi,wf,rw=r*w*1.; n=round(A*rw^2+B*rw+CC); a=a0+n*s*b3;t=a/b;x0=round(t^2);x=x0+i+w*bb; x2=sqrt(x);y=round(x2*x);kx=x^3-y^2; if(kx!=0&abs(kx)> Hall.dat"); if(guarda,system(dats) ) ) ) ) ) ) ) ) } /* Please, comment from here to the end of the file before use 'gp2c-run' */ Xs=[2, 5234, 8158,93844, 367806, 421351,720114, 939787, 28187351,\ 110781386, 154319269, 384242766, 390620082, 3790689201, 65589428378,952764389446,\ 12438517260105, 35495694227489, 53197086958290, 5853886516781223, 12813608766102806, 23415546067124892, 38115991067861271, 322001299796379844, 471477085999389882, 810574762403977064,9870884617163518770,42532374580189966073,51698891432429706382,44648329463517920535,231411667627225650649,601724682280310364065,4996798823245299750533,14038790674256691230847, 372193377967238474960883]; default(realprecision,100);default(format,"g0.5"); print("\nGood Examples of Hall's conjecture\n(you can find them in the vector 'Xs')"); print("\n# x r"); print("=================================="); for(i=1,length(Xs),x=Xs[i];k=hall_k(x);print(i," ",x," ",sqrt(x)/k)); { print("\nTo search for 'examples', you can use the routine \n \thall_search_bC00(bi,bf,save=1)\n\n Use, for example, \n \thall_search_bC00(2,10^3,0)\n Put an 1 in the third entry to save the findings in file 'hall00.dat'."); }