\\ \\ c12.gp \\ \\ calculate c1 of LLL-arigorithm \\ \\ Refference: \\ [1]Nigel P. Smart,"The Algorithmic Resolution of Diophantine Equations", \\ Camdridge Univsersity Press, 1998. \\ \\ 2021.03.08 modify functions: eeb1(),eeb2(),eeb3() \\ \\ change ee*abs(e.omega[1])*sqrt(A0*imtau)/(u0^2*sqrt(3*Pi)) \\ to ee*abs(e.omega[1])*sqrt(A0*imtau)/(abs(u0)*sqrt(3*Pi)); \\ \\ exchange e.omgega[1] and e.omega[2] \\ such that |tau| >=1, Imag(tau)>0, -1/2< Real(tau) <=1/2 \\ with Real(tau) >= 0 if |tau|=1; \\ \\ 2021.03.17 calculate omega1,omega2 using AGM function agm() \\ \\ \\The Gram-Schmidt Process: [1] Theorem V.7(p.67) \\ gramschmidt(b)= { local(bs,s,r,m); bs=b; m=b; s=matsize(b)[1]; r=matsize(b)[2]; for(i=2,r, for(j=1,r, bs[j,i]=1.0*b[j,i]-sum(k=1,(i-1),((sum(l=1,s,b[l,i]*bs[l,j])/sum(l=1,s,bs[l,j]*bs[l,j]))*bs[k,j])) ) ); return(bs) } \\ \\ c0=log|u| if v <= 0 \\ =log|u|+(1/2)*v/(X0-v) if v > 0 \\ cc0(X0,u,v)=if(v<=0,log(abs(u)),log(abs(u))+(1/2)*v/(X0-v)); \\ \\ c1 \\ cc1(b)= { local(bs,s,r,nb1,nb,nbmin); bs=gramschmidt(b); s=matsize(bs)[1]; r=matsize(bs)[2]; nb1=sum(k=1,s,bs[k,1]^2); nbmin=nb1; for(i=2,r, nb=sum(k=1,s,bs[k,i]^2); if(nbmin gamma_2 > gamma_3 \\ gamma_i=4*e_i \\ \\ omega_1=2*\pi/agm(sqrt(gamma_1-gamma_3),sqrt(gamma_1-gamma_2)); \\ omega_2=2*\pi*I/agm(sqrt(gamma_1-gamma_3),sqrt(gamma_2-gamma_3)); \\ omega_pe(a,b)= { local(d,om1,om2,ww,g1,g2,g3,gg); \\ local(AA,BB,gg1,gg2,gg3); d=poldisc(x^3+a*x+b); \\ print("d=",d); if(d==0, return([])); ww=polroots(x^3+a*x+b); \\ print(ww); if(d>0, g1=max3(real(ww[1]),real(ww[2]),real(ww[3])), if(imag(ww[1])==0, g1=real(ww[1]), g1=real(ww[2]) ) ); if(d>0, g3=min3(real(ww[1]),real(ww[2]),real(ww[3])); g2=real(ww[1])+real(ww[2])+real(ww[3])-g1-g3; om1=2*Pi/agm(sqrt(g1-g3),sqrt(g1-g2)); om2=2*Pi*I/agm(sqrt(g1-g3),sqrt(g2-g3)); [om1,om2], gg=sqrt(3*g1^2+a); om1=Pi/agm(sqrt(gg),sqrt(3*g1+2*gg)/2); om2=Pi*I/agm(sqrt(gg),sqrt(-3*g1+2*gg)/2); [om1+om2,om1-om2] ) } \\ \\ calclate om1,om2 by omega_pe() \\ \\ \\ A_i=max3(ellheight(e,p_i),he,2*Pi*abs(u_i)^2/(abs(e.omega[1])^2*Im(tau)) \\ AA(e,pi,he,om)= { local(ui,om1,om2,imtau,X0); X0=10^20; om1=om[1]; om2=om[2]; ui=om1*phi(e,pi,X0); imtau=imag(om2/om1); if(imtau<0, imtau=abs(imtau) ); return(max3(ellheight(e,pi), he, 2*Pi*abs(ui)^2/(abs(om1)^2*imtau))) } \\ \\ Epsiron=min_{i=0,1,2,3}{exp(1)*|\omega_1|*sqrt(Ai*Imag(tau))/(|u_i|*sqrt(3*Pi))} \\ \\ where \\ |tau| >=1, Imag(tau)>0, -1/2< Real(tau) <=1/2 \\ with Real(tau) >= 0 if |tau|=1 \\ eeb3(e,p1,p2,p3,A0,A1,A2,A3,om)= { local(u0,u1,u2,u3,tau,om1,om2,X0); X0=10^20; om1=om[1]; om2=om[2]; u0=om1; u1=om1*phi(e,p1,X0); u2=om1*phi(e,p2,X0); u3=om1*phi(e,p3,X0); imtau=imag(om2/om1); if(imtau<0, imtau=abs(imtau) ); ee=exp(1); return(min4(ee*abs(om1)*sqrt(A0*imtau)/(abs(u0)*sqrt(3*Pi)), ee*abs(om1)*sqrt(A1*imtau)/(abs(u1)*sqrt(3*Pi)), ee*abs(om1)*sqrt(A2*imtau)/(abs(u2)*sqrt(3*Pi)), ee*abs(om1)*sqrt(A3*imtau)/(abs(u3)*sqrt(3*Pi)))) } eeb2(e,p1,p2,A0,A1,A2,om)= { local(u0,u1,tau,om1,om2,X0); X0=10^20; om1=om[1]; om2=om[2]; u0=om1; u1=om1*phi(e,p1,X0); u2=om1*phi(e,p2,X0); imtau=imag(om2/om1); if(imtau<0, imtau=abs(imtau) ); ee=exp(1); return(min3(ee*abs(om1)*sqrt(A0*imtau)/(abs(u0)*sqrt(3*Pi)), ee*abs(om1)*sqrt(A1*imtau)/(abs(u1)*sqrt(3*Pi)), ee*abs(om1)*sqrt(A2*imtau)/(abs(u2)*sqrt(3*Pi)))) } eeb1(e,p1,A0,A1,om)= { local(u0,u1,tau,om1,om2,X0); X0=10^20; om1=om[1]; om2=om[2]; u0=abs(om1); u1=abs(om1)*phi(e,p1,X0); imtau=imag(om2/om1); if(imtau<0, imtau=abs(imtau) ); ee=exp(1); return(min(ee*abs(om1)*sqrt(A0*imtau)/(abs(u0)*sqrt(3*Pi)), ee*abs(om1)*sqrt(A1*imtau)/(abs(u1)*sqrt(3*Pi)))) } \\ \\ c4=2*10^(7*r+5)*(2/e)^(4*r^2+18*r+14)*log(eps)^(-2*r-3)*\prod_{i=0}^r{A_i} \\ cc4(r,eps,ProdAi)=2*10^(7*r+5)*(2/exp(1))^(4*r^2+18*r+14)*log(eps)^(-2*r-3)*ProdAi; \\ \\ c5=log(eps) \\ cc5(eps)=log(eps); \\ \\ c6=c5+he; \\ cc6(c5,he)=c5+he; \\ \\ c7=c6+log(t0*r)+(2*t0-1)/(16*t0*r); \\ cc7(c6,r,t0)=c6+log(t0*r)+(2*t0-1)/(16*t0*r); \\ \\ c8=c6+(log(t0*r)+(2*t0-1)/(16*t0*r))/log(16) \\ cc8(c6,r,t0)=c6+(log(t0*r)+(2*t0-1)/(16*t0*r))/log(16); \\ \\ ||b_1|| of bs \\ nr(bs)= { local(l,s); l=matsize(bs)[1]; sqrt(sum(i=1,l,bs[i,1]^2)) } \\ \\ case Im(a_i)=0 for all i \\ \\ Lemma VI.1([1] p.83) \\ hmax(C,c2,c3,X,b)= { local(S,T,l,c4,h); l=matsize(b)[1]; S=(l-1)*X; T=(1+l*X)/2; c4=nr(b)/sqrt(cc1(b)); \\ print("c4=",c4); if(c4^2>=(T^2+S), \\ print("OK"); h=(log(C*c2)-log(sqrt(c4^2-S)-T))/c3, \\ print("NG"); h=-1 ); return(h) } \\ \\ general case: Re(a_n)Im(a_{n-1})-Re(a_{n-1})R(a_{n}) != 0 \\ \\ Lemmna VI.2([1] p.84) \\ hmax2(C,c2,c3,X,b)= { local(S,T,l,m,c4,h); l=matsize(b)[1]; m=matsize(b)[2]; S=(l-2)*X^2; T=(1+l*X)/sqrt(2); c4=nr(b)/sqrt(cc1(b)); \\ print("c4=",c4); if(c4^2>=(T^2+S), \\ print("OK"); h=(log(C*c2)-log(sqrt(c4^2-S)-T))/c3, \\ print("NG"); h=-1 ); return(h) } \\ \\ check |x*log(2)+y*log(3)+z*log(5)-log(7)| <= exp(-max(|x|,|y|,|z|)) \\ check(X)= { local(p,q,m); for(x=-X,X, for(y=-X,X, m=max(abs(x),abs(y)); for(z=-X,X, p=abs(x*log(2)+y*log(3)+z*log(5)-log(7)); q=max(m,abs(z)); if(p <= exp(-q),print([x,y,z])) ) ) ) } \\ \\ Zagier's Algorithm : elliptic logarithm \\ \\ \phi(P)=(1/{\omega}) \int_O^P{dt/sqrt{f(t)}} \\ =\sum_{i=0}^{\infinity}{c_i/2^(i+1)} \\ where c_i=0 if y([2^i]P) >= 0 or [2^i]P=O \\ c_i=1 if y([2^i]P) < 0 \\ phi(e,p,x0)= { local(m,n,c,s,q,x,y); if(p==[0],return(1)); m=1; s=0; n=10*x0; q=[p[1]*1.0,p[2]*1.0]; while(n>0, q=ellpow(e,q,2); c=if((q==[0] || q[2]>=0),0,1); s=s*2+c; m=2*m; n=n\2 ); return((s/m)*1.0) } \\ \\ H=(1/2)[]; U=[\delta_{i,j}]; ff=|H-r*U|; \\ HH2(e,p1,p2)= { local(h); local(h11,h12,h22); h11=ellbil(e,p1,p1)/2; h12=ellbil(e,p1,p2)/2; h22=ellbil(e,p2,p2)/2; h=[h11,h12;h12,h22]; return(h) } HH3(e,p1,p2,p3)= { local(h); local(h11,h12,h13,h22,h23,h33); h11=ellbil(e,p1,p1)/2; h12=ellbil(e,p1,p2)/2; h13=ellbil(e,p1,p3)/2; h22=ellbil(e,p2,p2)/2; h23=ellbil(e,p2,p3)/2; h33=ellbil(e,p3,p3)/2; h=[h11,h12,h13; h12,h22,h23; h13,h23,h33]; return(h) } \\ \\ K0 >= (2^(r/2)*t*K3*sqrt(r^2+r))^(r+1) \\ KK0(r,t,K3)=(2^(r/2)*t*K3*sqrt(r^2+r))^(r+1); \\ \\ check whether i*P1+j*P2+k*P3 is an integral point or not \\ \\ max(|i|,|j|,|k|)<= X \\ check3(X,e,p1,p2,p3)= { local(p,q,r); local(count); count=0; for(i=-X,X, p=ellpow(e,p1,i); for(j=-X,X, q=elladd(e,p,ellpow(e,p2,j)); for(k=-X,X, r=elladd(e,q,ellpow(e,p3,k)); if(isintegral(r), count++; print(r) ) ) ) ); print("There are ",count," integral points.") } \\ \\ check whether i*P1+j*P2 is an integral point or not \\ \\ max(|i|,|j|)<= X check2(X,e,p1,p2)= { local(p,q); local(count); count=0; for(i=-X,X, p=ellpow(e,p1,i); for(j=-X,X, q=elladd(e,p,ellpow(e,p2,j)); if(isintegral(q), count++; print(q) ) ) ); print("There are ",count," integral points.") } \\ \\ check whether i*P1 is an integral point or not \\ \\ max(|i|)<= X \\ check1(X,e,p1)= { local(p); local(count); count=0; for(i=-X,X, p=ellpow(e,p1,i); if(isintegral(p), count++; print(p) ) ); print("There are ",count," integral points.") } \\ \\ check whether i*P1+T is an integral point or not \\ \\ max(|i|)<= X \\ check1t(X,e,p1,t)= { local(p); local(count); count=0; for(i=-X,X, p=elladd(e,ellpow(e,p1,i),t); if(isintegral(p), count++; print(p) ) ); print("There are ",count," integral points.") } \\ \\ check whether i*P1+j*P2+T is an integral point or not \\ \\ max(|i|,|j|)<= X \\ check2t(X,e,p1,p2,t)= { local(p,q); local(count); count=0; for(i=-X,X, p=elladd(e,ellpow(e,p1,i),t); for(j=-X,X, q=elladd(e,p,ellpow(e,p2,j)); if(isintegral(q), count++; print(q) ) ) ); print("There are ",count," integral points.") } \\ \\ check whether i*P1+j*P2+k*P3+T is an integral point or not \\ \\ max(|i|,|j|,|k|)<= X \\ check3t(X,e,p1,p2,p3,t)= { local(p,q,r); local(count); count=0; for(i=-X,X, p=elladd(e,ellpow(e,p1,i),t); for(j=-X,X, q=elladd(e,p,ellpow(e,p2,j)); for(k=-X,X, r=elladd(e,q,ellpow(e,p3,k)); if(isintegral(r), count++; print(r) ) ) ) ); print("There are ",count," integral points.") } \\ \\ if p is an integral point \\ then 1 else 0 \\ isintegral(p)= { local(x,y); if(p==[0],return(0)); x=p[1]; y=p[2]; return(denominator(x)==1 && denominator(y)==1) }