\\ \\ Auteur: \\ Denis SIMON -> simon@math.unicaen.fr \\ adresse du fichier: \\ www.math.unicaen.fr/~simon/qfsolve.gp \\ \\ ********************************************* \\ * VERSION TRES PROVISOIRE A USAGE PERSONNEL * \\ ********************************************* \\ * VERSION 09/03/2003 * \\ ********************************************* \\ * SANS AUCUNE GARANTIE * \\ ********************************************* \\ \\ Programme de resolution des equations quadratiques G(x,y,z)=0 \\ langage: GP \\ pour l'utiliser, lancer gp, puis taper \\ \r qfsolve.gp \\ \\ Ce fichier contient 4 principales fonctions: \\ \\ - diagonalsolve(a,b,c): pour resoudre l'equation de Legendre \\ ax^2+by^2+cz^2=0 ou a,b et c sont entiers, premiers entre eux \\ et sans carres. On suppose qu'une solution existe. \\ \\ - Qfsolve(G): pour resoudre l'equation quadratique X^t*G*X = 0 \\ G doit etre une matrice symetrique 3*3, a coefficients dans Z, \\ et det(G) <> 0. S'il n'existe pas de solution, la response est 0. \\ \\ - IndefiniteLLL(G,c): pour resoudre ou reduire la forme quadratique \\ G a coefficients entiers. Il s'agit d'un algorithme type LLL, avec la \\ constante 1/40, la 'fl'eme forme quadratique est reduite. \\ \\ global(DEBUGLEVEL); DEBUGLEVEL = 0; {diagonalsolve(a,b,c)= local(fa,xa,fb,xb,fc,xc,bb,u,v,Q,QQ,QQ2,sol); if (DEBUGLEVEL >=5, print("entree dans diagonalsolve")); \\ cette fonction resout l'equation aX^2+bY^2+cZ^2=0 \\ la solution est [X,Y,Z]~. \\ on suppose qu'une solution existe (conditions locales). \\ on suppose que a,b et c sont premiers entre eux (et sans facteurs carres) fa=factor(abs(a))[,1];fb=factor(abs(b))[,1];fc=factor(abs(c))[,1]; xa=Mod(1,1); for(i=1,length(fa), xa=chinese(xa,sqrt(Mod(-c/b,fa[i])))); xa=centerlift(xa); xb=Mod(1,1); for(i=1,length(fb), xb=chinese(xb,sqrt(Mod(-c/a,fb[i])))); xb=centerlift(xb); xc=Mod(1,1); for(i=1,length(fc), xc=chinese(xc,sqrt(Mod(-b/a,fc[i])))); xc=centerlift(xc); bb=bezout(b,c);u=bb[1];v=bb[2]; Q=[b*c,a*b*u*xc,xc*b*u*xa+xb*c*v;0,a,xa;0,0,1]; QQ=(Q~*[a,0,0;0,b,0;0,0,c]*Q)/(a*b*c); QQ2=IndefiniteLLL(QQ); sol=Q*QQ2; if (DEBUGLEVEL >=5, print("test = ",sol[1]^2*a+sol[2]^2*b+sol[3]^2*c)); return(sol); } {QfbReduce(M)= \\ Reduction de la forme quadratique binaire \\ qf=(a,b,c)=a*X^2+2*b*X*Y+c*Y^2 \\ Renvoit la matrice de reduction. \\ On suppose que la forme ne s'annule pas \\ ie -det(M) n'est pas un carre. \\ M=[a,b;b;c]; \\ Coefficients entiers local(a,b,c,H,test,di,q,r,nexta,nextb,nextc,aux,compteur); compteur=0; \\print(log(vecmax(abs(M))),";",.5*log(abs(matdet(M)))); a=M[1,1];b=M[1,2];c=M[2,2]; H=matid(2);test=1; while(test, compteur++; di=divrem(b,a); q=di[1];r=di[2]; if(2*r>abs(a),r-=abs(a);q+=sign(a)); H[,2]-=q*H[,1]; nextc=a;nextb=-r;nexta=(nextb-b)*q+c; if( test = abs(nexta)=2, print("Minimization")); for( i = 1, length(factdetG[,1]), p=factdetG[i,1]; dU = kermodp(G,p); if(dU[1]==3, G/=p;U*=p; if(factdetG[i,2]-=3,i--); next ); if(dU[1]==2, UU=dU[2]*[1,0,0;0,1,0;0,0,p]; G=UU~*G*UU/p; U=U*UU; if(factdetG[i,2]--,i--); next ); if(dU[1]==1 && factdetG[i,2]>1, UU=dU[2]*[1,0,0;0,p,0;0,0,p]; G=UU~*G*UU/p^2; U=U*UU; if(factdetG[i,2]-=2,i--); next ); \\ maintenant : dU[1]=1 && factdetG[i,2]=1 G=dU[2]~*G*dU[2]; U=U*dU[2]; if(!(G[2,2]%p), UU=[1,0,0;0,1,0;0,0,p] , d = Mod(G[3,2]^2-G[2,2]*G[3,3],p); if(!issquare(d), if (DEBUGLEVEL >=2, print(" Pas de solution locale en ",p)); return(0)); d=sqrt(d); UU=[1,0,0;0,p,centerlift((-G[2,3]+d)/G[2,2]);0,0,1]; ); G=UU~*G*UU/p; U=U*UU; ); \\ 2eme etape: on reduit la nouvelle forme G. H = IndefiniteLLL(G,1,1); if(length(H)!=3, if (DEBUGLEVEL >=2, print("Pas de solution reelle")); return(0)); \\ Enfin, on reconstruit une solution de G initial return(U*H); } {Qfparam(G,sol,fl=3)= \\ G est une matrice symetrique 3*3, et sol une solution de sol~*G*sol=0. \\ Renvoit une parametrisation des solutions avec de bons invariants, \\ sous la forme d'une matrice 3*3, dont chaque ligne contient \\ les coefficients de chacune des 3 formes quadratiques. \\ si fl!=0, la 'fl'eme forme quadratique est reduite. local(G1,G2,U); \\print("entree dans Qfparam"); sol/=content(sol); \\ construction de U telle que U[,3] est sol, et det(U)=+-1 U=(mathnf(Mat(sol~),1)[2]~)^-1; G1=U~*G*U; \\ G1 a un 0 en bas a droite. G2=[-2*G1[1,3],-2*G1[2,3],0; 0,-2*G1[1,3],-2*G1[2,3]; G1[1,1],2*G1[1,2],G1[2,2]]; sol=U*G2; if(fl && !issquare(matdet(U=[sol[fl,1],sol[fl,2]/2; sol[fl,2]/2,sol[fl,3]])), U=QfbReduce(U); U=[U[1,1]^2,2*U[1,2]*U[1,1],U[1,2]^2; U[2,1]*U[1,1],U[2,2]*U[1,1]+U[2,1]*U[1,2],U[1,2]*U[2,2]; U[2,1]^2,2*U[2,1]*U[2,2],U[2,2]^2]; sol=sol*U; ); \\print("sortie de Qfparam"); return(sol); } \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ \\ DEBUT DE LA PARTIE DEVELOPPEMENT \\ \\ \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ {Qfparamold(G,sol,fl=3)= \\ G est une matrice symetrique 3*3, et sol une solution de sol~*G*sol=0. \\ Renvoit une parametrisation des solutions avec de bons invariants, \\ sous la forme d'une matrice 3*3, dont chaque ligne contient \\ les coefficients de chacune des 3 formes quadratiques. \\ si fl!=0, la 'fl'eme forme quadratique est reduite. local(U,X0,q0,q1,M1,M2,solparam); \\print("entree dans Qfparamold"); \\ construction de U telle que U[,3] est sol, et det(U)=+-1 sol/=content(sol); U=(mathnf(Mat(sol~),1)[2]~)^-1; X0=vecextract(U,3); q0=Mat([X0~*G*X0]); q1=X0~*G*sol; M1=matrix(2,2);M1[,1]=q1;M1=Mat([M1+M1~]); M2=matrix(2,2);M2[,2]=q1;M2=Mat([M2+M2~]); \\ Rem: Ces changements de types permettent de multiplier terme a terme les \\ coeff d'un vecteur par une matrice. solparam=Mat(sol)*q0-Mat(X0[,1])*M1-Mat(X0[,2])*M2; solparam=solparam[,1]; if(fl && !issquare(-matdet(solparam[fl])), U=QfbReduce(solparam[fl]); solparam=vector(3,i,U~*solparam[i]*U); ); sol = matrix(3,3); for(i=1,3, sol[i,1]=solparam[i][1,1]; sol[i,2]=solparam[i][1,2]*2; sol[i,3]=solparam[i][2,2] ); \\print("sortie de Qfparamold"); return(sol); } {LLLgoon(G,c=1)= local(red,U1,G2,b,U2,G3,U3,cc); \\ reduction LLL de la forme quadratique G (matrice de Gram) \\ en dim 3 seulement. \\ ou l'on continue meme si on trouve un vecteur isotrope red = IndefiniteLLL(G,c); \\ si on ne trouve pas de vecteur isotrope, il n'y a rien a faire. if(length(red)==2,return(red)); \\ si on trouve un vecteur isotrope red /= content(red); U1 = (mathnf(Mat(red~),1)[2]~)^-1; \\ la 3eme colonne de U1 est red, et det(U1)=+-1 G2 = U1~*G*U1; \\ la matrice G2 a un 0 dans le coin en bas a droite. b = bezout(G2[3,1],G2[3,2]); U2 = [b[1],G2[3,2]/b[3],0;b[2],-G2[3,1]/b[3],0;0,0,-1]; G3 = U2~*G2*U2; \\ la matrice G3 a des 0 sous l'anti-diagonale. cc = G3[1,1]%2; U3 = [1,0,0; cc,1,0; round(-(G3[1,1]+cc*(2*G3[1,2]+G3[2,2]*cc))/2/G3[1,3]), round(-(G3[1,2]+cc*G3[2,2])/G3[1,3]),1]; return([U3~*G3*U3,U1*U2*U3]); } {LLLgoon2(G,c=1)= local(red,U1,G2,l,U2,G3,U3,G4,U,V,B,U4,G5,U5,G6); \\ reduction LLL de la forme quadratique G (matrice de Gram) \\ ou l'on continue meme si on trouve un vecteur isotrope red = IndefiniteLLL(G,c); \\ si on ne trouve pas de vecteur isotrope, il n'y a rien a faire. if(length(red)==2,return(red)); \\ si on trouve un vecteur isotrope red /= content(red); U1 = (mathnf(Mat(red~),1)[2]~)^-1; \\ la derniere colonne de U1 est red, et det(U1)=+-1 G2 = U1~*G*U1; \\ la matrice G2 a un 0 dans le coin en bas a droite. U2 = mathnf(Mat(G2[l=length(G),]),1)[2]; G3 = U2~*G2*U2; \\ la matrice de G3 a des 0 sur toute la 1ere ligne, sauf un 'g' au bout a droite, avec g^2| det G. U3 = matid(l); U3[1,l] = round(-G3[l,l]/G3[1,l]/2); G4 = U3~*G3*U3; \\ le coeff G4[l,l] est reduit modulo 2g U = vecextract(G4,[1,l],[1,l]); V = vecextract(G4,[1,l],1<<(l-1)-2); B = round(-U^-1*V); U4 = matid(l); for(j=2,l-1,U4[1,j]=B[1,j-1];U4[l,j]=B[2,j-1]); G5 = U4~*G4*U4; \\ la dernire colonne de G5 est reduite if(l>=4, red = LLLgoon2(matrix(l-2,l-2,i,j,G5[i+1,j+1]),c); U5 = matdiagonalblock([Mat([1]),red[2],Mat([1])]); G6 = U5~*G5*U5; return([G6,U1*U2*U3*U4*U5]) , return([G5,U1*U2*U3*U4])) } {setim(G,c=1)= local(red,U1,G2,U2,G3,U3,G4,U4,G5,v3,s1,s2,s3,sol); \\ recherche d'un setim de dim 3 de la forme quadratique entiere G \\ (matrice de Gram) en dim 6 seulement. \\ on suppose que det(G)=-1 if( DEBUGLEVEL >= 4, print("entree dans setim avec ")); if( DEBUGLEVEL >= 5, print("G=",G)); red = IndefiniteLLL(G,c); red /= content(red); U1 = (mathnf(Mat(red~),1)[2]~)^-1; G2 = U1~*G*U1; \\ la matrice G2 a un 0 en bas a droite. if( DEBUGLEVEL >= 5, print("G2=",G2)); U2 = mathnf(vecextract(G2,32,31),1)[2]; G3 = vecextract(U2~*vecextract(G2,31,31)*U2,15,15); if( DEBUGLEVEL >= 5, print("G3=",G3)); if( DEBUGLEVEL >= 5, print("det3=",matdet(G3))); red = IndefiniteLLL(G3,c); red /= content(red); U3 = (mathnf(Mat(red~),1)[2]~)^-1; G4 = U3~*G3*U3; \\ la matrice G4 a un 0 en bas a droite. if( DEBUGLEVEL >= 5, print("G4=",G4)); U4 = mathnf(vecextract(G4,8,7),1)[2]; G5 = vecextract(U4~*vecextract(G4,7,7)*U4,3,3); if( DEBUGLEVEL >= 5, print("G5=",G5)); if( DEBUGLEVEL >= 5, print("det5=",matdet(G5))); if (!G5[1,1], v3 = [1,0,0]~ , v3 = [-G5[1,2]+1,G5[1,1],0]~); s1 = U1[,6]; if( DEBUGLEVEL >= 5, print("s1=",s1)); s2 = U1*concat(U2*concat(U3[,4],[0]~),[0]~); if( DEBUGLEVEL >= 5, print("s2=",s2)); s3 = U1*concat(U2*concat(U3*concat(U4*v3,[0]~),[0]~),[0]~); if( DEBUGLEVEL >= 5, print("s3=",s3)); if( DEBUGLEVEL >= 4, print("fin de setim")); sol = matrix(6,3);sol[,1]=s1;sol[,2]=s2;sol[,3]=s3; if( DEBUGLEVEL >= 3, print("test setim = ",if(sol~*G*sol,"NOT ",""),"ok")); return(sol); } {QfIndiv(G)= \\ pseudo-diagonalisation de G : \\ avec des 0 partout sauf sur et autour de la diagonale n=length(G); for(i=1,n-1, for(j=i+2,n, if(G[i,j], be=bezout(G[i,i+1],G[i,j]); M=matid(n); M[i+1,i+1]=be[1];M[j,i+1]=be[2]; M[i+1,j]=-G[i,j];M[j,j]=G[i,i+1]; G=M~*G*M ) ) ); return(G); } {QfbSolve(M)= \\ Suppose que la forme binaire quadratique M (donnee par sa matrice \\ de Gram) a une solution (equivaut a -det(M) = carre). \\ trouve une solution. local(a,b,c,di); \\ Coefficients entiers \\ Existe-t-il une solution ? a=M[1,1];b=M[1,2];c=M[2,2]; if(!a, return([1,0]~)); if(!c, return([0,1]~)); di=b^2-a*c; print("solution ici = ",[-b+sqrtint(di),a]*M*[-b+sqrtint(di),a]~); return([-b+sqrtint(di),a]~); } {legendresolve(a, b) = local(aa,bb,na,nb,mat,result,t,fact,l,f,xx,alpha,oldnb,newnb,verif); if (DEBUGLEVEL >=5, print("entree dans legendresolve")); \\ cette fonction resout l'equation X^2-aY^2=bZ^2 \\ la solution est [X,Y,Z]~. \\ on suppose qu'une solution existe (conditions locales). if(DEBUGLEVEL >=3, print("(a,b)=(",a,",",b,")")); aa = a; bb = b; \\ cas degeneres aa ou bb = 0 if ( aa == 0, if (DEBUGLEVEL >= 5, print("fin de legendresolve")); return ([0,1,0]~)); if ( bb == 0, if (DEBUGLEVEL >= 5, print("fin de legendresolve")); return ([0,0,1]~)); na = abs(aa); nb = abs(bb); mat = matid(3); while (1, if(DEBUGLEVEL >=4, print("[a,b]=",[aa,bb])); if(DEBUGLEVEL >=5, print("***"nb"*** ",nb)); \\ cas faciles if ( aa == 1 ,result = [1,1,0]~; break); if ( bb == 1 ,result = [1,0,1]~; break); if ( aa+bb == 1 ,result = [1,1,1]~; break); if ( aa+bb == 0 ,result = [0,1,1]~; break); if ( aa == bb , t = aa*mat[,1]; mat[,1] = mat[,3]; mat[,3] = t; aa = -1; na = 1); if ( issquare(aa) , result=[sqrtint(aa),1,0]~; break); if ( issquare(bb) , result=[sqrtint(bb),0,1]~; break); if ( na > nb , t = aa; aa = bb; bb = t; t = na; na = nb; nb = t; t = mat[,3]; mat[,3] = mat[,2]; mat[,2] = t); print("factoring",nb); fact = factor(nb); l = length(fact[,1]); for (i = 1, l, if ( (f=(fact[i,2]>>1)) , fact[i,2]-=(f<<1); t = fact[i,1]^f; mat[,3] /= t; t=t^2; nb /= t; bb/=t)); if(DEBUGLEVEL >=4, print("bb=",bb)); if ( na >= nb , next ); l=length(fact[,1]); xx=Mod(1,1); for (i = 1, l, if (fact[i,2] , xx = chinese(xx,sqrt(Mod(aa,fact[i,1]))))); xx=centerlift(xx); alpha = xx^2-aa; if ( alpha == 0, result=[xx,1,0]~; break); t = alpha/bb; oldnb = nb; newnb = abs(t); while ( nb > newnb , mat[,3] *= t; bb = t; nb = newnb; t = xx*mat[,1]+mat[,2]; mat[,2] = aa*mat[,1] + xx*mat[,2]; mat[,1] = t; xx = centerlift(Mod(-xx,bb)); alpha = xx^2-aa; t = alpha/bb; newnb = abs(t) ) ); result = mat*result; if (DEBUGLEVEL >= 5, print("result with content=",result)); if (DEBUGLEVEL >= 5, print("content=",content(result))); result /= content(result); if (DEBUGLEVEL >= 5, print("result=",result)); verif = (result[1]^2-a*result[2]^2-b*result[3]^2==0); if ( !verif && DEBUGLEVEL >= 0, print("***** erreur dans legendresolve : mauvaise solution *******",[a,b]);1/0); if (DEBUGLEVEL >= 3, print("fin de legendresolve")); return (result); } { legendresolve2(a,b,fa,fb)= \\ for the moment, only valid for \\ a,b squarefree and coprime. \\ \\ assume the local conditions (as usual...) \\ uses the factorization fa of |a| and fb of |b|. local(xa,xb,qinit,be,M,newq); xa=Mod(1,1); for (i = 1, length(fb[,1]), xa = chinese(xa,sqrt(Mod(a,fb[i,1])))); xa=centerlift(xa); xb=Mod(1,1); for (i = 1, length(fa[,1]), xb = chinese(xb,sqrt(Mod(b,fa[i,1])))); xb=centerlift(xb); qinit=[1,0,0;0,-a,0;0,0,-b]; be=bezout(a,b); M=[a*b,xa*a*be[1],xb*b*be[2];0,1,0;0,0,1]; newq=M~*qinit*M; newq/=a*b; \\ The symetrix matrix newq should have det=+-1 \\ Then use lll to find a solution. sol=M*IndefiniteLLL(newq); return(sol~*qinit*sol); } {Qfinvariantsdet1(G,p)= \\ pour une forme quadratique unimodulaire (det=+-1) local(n,r,s); n=length(G); s=qfsign(G);r=s[1];s=s[2]; return([n,r-s,vector(n,i,G[i,i]%2)]); } {Qfinvariantc(G,p)= \\ calcule l'invariant c d'une forme quadratique, p-adique \\ on calcule d'abord les coefficients diagonaux (modulo les carres) local(n); n=length(G); diag = vector(n,i,0); for(i=1,n, diag[i] = matdet(vecextract(G,1<=3, print("p=",p,"^",factdetG[i,2])); \\ pour l'instant, on ne regarde que la valuation >=2 if((vp=factdetG[i,2])<2,next); ker=kermodp(G,p); dimker=ker[1];ker=ker[2]; if (DEBUGLEVEL >=3, print("dimker = ",dimker)); G=ker~*G*ker; U=U*ker; \\ 1er cas: la dimension du noyau est plus petite que la valuation if(dimker=3, print("cas 1")); ker2=kermodp(vecextract(G,1<=3, print("fin cas 1")); i--;next; ); \\ maintenant, vp = dimker <= n \\ 2eme cas: la dimension du noyau est tres grande if(2*dimker>n, if (DEBUGLEVEL >=3, print("cas 2")); ker=matid(n); for(j=dimker+1,n,ker[j,j]=p); G=ker~*G*ker/p; U=U*ker; factdetG[i,2]-=2*dimker-n; if (DEBUGLEVEL >=3, print("fin cas 2")); i--;next; ); \\ maintenant, vp = dimker <= n/2 \\ 3eme cas: la dimension du noyau est >=2 et contient un elt de norme p^2... if( dimker >2 || (dimker==2 && issquare(di=Mod((G[1,2]^2-G[1,1]*G[2,2])/p^2,p))), \\ on cherche dans le noyau un elt de norme p^2... if(dimker>2, if (DEBUGLEVEL >=3, print("cas 3.1")); dimker=3; sol=Qfsolvemodp(vecextract(G,7,7)/p,p) , if (DEBUGLEVEL >=3, print("cas 3.2")); if( !(G[1,1]%p^2), if (DEBUGLEVEL >=3, print("cas 3.2.1")); sol=[1,0]~ , if (DEBUGLEVEL >=3, print("cas 3.2.2")); sol=[-G[1,2]/p+sqrt(di),Mod(G[1,1]/p,p)] ) ); sol=centerlift(sol); sol/=content(sol); if (DEBUGLEVEL >=3, print("sol=",sol)); ker=matrix(1,n,j,k,0);for(j=1,dimker,ker[1,j]=sol[j]); \\ on complete avec des 0 ker=(mathnf(ker,1)[2]~)^-1; \\ on trouve une base dont le dernier vecteur soit ker for(j=1,n-1,ker[,j]*=p); G=ker~*G*ker/p^2; U=U*ker; factdetG[i,2]-=2; if (DEBUGLEVEL >=3, print("fin cas 3")); i--;next; ); \\ maintenant, vp = dimker <= 2 \\ et le noyau ne contient pas de vecteur de norme p^2... \\ RESTE A FAIRE ); return([G,U]) } {mymat(qfb)=qfb=Vec(qfb);[qfb[1],qfb[2]/2;qfb[2]/2,qfb[3]] } {Qfbsqrtgauss(G,factdetG)= \\ pour l'instant, ca ne marche que pour detG pair local(a,b,c,d,m,n,p,aux,Q1,M); if (DEBUGLEVEL >=3, print("entree dans Qfbsqrtgauss avec",G)); G = Vec(G); a = G[1]; b = G[2]/2; c = G[3]; d = a*c-b^2; \\ 1ere etape: on resout m^2 = a (d), m*n = -b (d), n^2 = c (d) m = n = Mod(1,1); factdetG[1,2]-=3; for( i = 1, length(factdetG[,1]), if(!factdetG[i,2],next); p = factdetG[i,1]; if(gcd(a,p)==1, aux = lift(polrootspadic(x^2-a,p,factdetG[i,2])[1]); aux = Mod(aux,p^factdetG[i,2]); m = chinese(m,aux); n = chinese(n,Mod(-b,p^factdetG[i,2])/aux); \\ n = chinese(n,-b/aux); , aux = lift(polrootspadic(x^2-c,p,factdetG[i,2])[1]); aux = Mod(aux,p^factdetG[i,2]); n = chinese(n,aux); m = chinese(m,Mod(-b,p^factdetG[i,2])/aux); \\ m = chinese(m,-b/aux); ); ); m = centerlift(m); n = centerlift(n); if (DEBUGLEVEL >=4, print("m=",m);print("n=",n)); \\ 2eme etape: on construit Q1 de det -1 tq Q1(x,y,0)=G(x,y) Q1 = [(n^2-c)/d, (m*n+b)/d, n ; (m*n+b)/d, (m^2-a)/d, m ; n, m, d ]; Q1 = -matadjoint(Q1); \\ 3eme etape: on reduit Q1 jusqu'a [0,0,-1;0,1,0;-1,0,0] M = LLLgoon(Q1)[2][3,]; if(M[1]<0,M=-M); if (DEBUGLEVEL >=3, print("fin de Qfbsqrtgauss ")); if( M[1]%2, return(Qfb(M[1],2*M[2],2*M[3])) , return(Qfb(M[3],-2*M[2],2*M[1]))); } {class2(D,factdetG)= \\ On suit l'algorithme de Shanks/Bosma-Stevenhagen \\ pour construire tout le 2-groupe de classes \\ seulement pour D discriminant fondamental \\ lorsque D=1(4), on travaille avec 4D. \\ Si on connait la factorisation de abs(2*D), \\ on peut la donner dans factdetG, et dans ce cas le reste \\ de l'algorithme est polynomial. local(factD,n,rang,m,listgen,vD,p,vp,aux,invgen,im,ker,kerim,listgen2,G2,struct); if( DEBUGLEVEL >= 1, print("Construction du 2-groupe de classe de discriminant ",D)); if( D%4 == 2 || D%4 == 3, print("Discriminant not congruent to 0,1 mod 4");return(0)); if(D==-4, return([[1],[Qfb(1,0,1)]])); if(!factdetG, factdetG = factor(2*abs(D))); factD = concat([-1],factdetG[,1]~); if( D%4 == 1, D*=4; factdetG[1,2]+=2); n = length(factD); rang = n-3; if(D>0, m = rang+1, m = rang); if( DEBUGLEVEL >= 3, print("factD = ",factD)); listgen = vector(m); if( vD = valuation(D,2), E = Qfb(1,0,-D/4); , E = Qfb(1,1,(1-D)/4) ); if( DEBUGLEVEL >= 3, print("E = ",E)); for( i = 1, m, \\ on ne regarde pas factD[1]=-1, ni factD[2]=2 p = factD[i+2]; vp = valuation(D,p); aux = p^vp; if (vD, listgen[i] = Qfb(aux,0,-D/4/aux) , listgen[i] = Qfb(aux,aux,(aux-D/aux)/4)); ): if( vD == 2 && D%16 != 4, m++;rang++; listgen = concat(listgen,[Qfb(2,2,(4-D)/8)])); if( vD == 3, m++;rang++; listgen = concat(listgen,[Qfb(2^(vD-2),0,-D/2^vD)])); if( DEBUGLEVEL >= 3, print("listgen = ",listgen)); if( DEBUGLEVEL >= 2, print("rang = ",rang)); if( !rang, return([[1],[E]])); invgen = matrix(n,m); for( i = 1, m, invgen[,i] = Qflisteinvariants(mymat(listgen[i]),factD)[,2]*Mod(1,2) ); if( DEBUGLEVEL >= 3, printp("invgen = ",lift(invgen))); struct = vector(m,i,2); compteur = 0; im = lift(matinverseimage(invgen,matimage(invgen))); while( (length(im) < rang), if( DEBUGLEVEL >= 4, print("passage ",compteur++)); ker = lift(matker(invgen)); kerim = concat(ker,im); listgen2 = vector(m); for( i = 1, m, listgen2[i] = E; for( j = 1, m, if( kerim[j,i], listgen2[i] = qfbcompraw(listgen2[i],listgen[j]))); if( norml2(kerim[,i]) > 1, red = QfbReduce(aux=mymat(listgen2[i])); aux = red~*aux*red; listgen2[i] = Qfb(aux[1,1],2*aux[1,2],aux[2,2])) ); listgen = listgen2; invgen = invgen*kerim; if( DEBUGLEVEL >= 4, print("listgen = ",listgen)); if( DEBUGLEVEL >= 4, printp("invgen = ",lift(invgen))); for( i = 1, length(ker), G2 = Qfbsqrtgauss(listgen[i],factdetG); struct[i]<<=1: listgen[i] = G2; invgen[,i] = Qflisteinvariants(mymat(G2),factD)[,2]*Mod(1,2) ); if( DEBUGLEVEL >= 3, print("listgen = ",listgen)); if( DEBUGLEVEL >= 3, printp("invgen = ",lift(invgen))); if( DEBUGLEVEL >= 3, print("struct = ",struct)); im = lift(matinverseimage(invgen,matimage(invgen))); ); listgen2 = vector(rang); for( i = 1, rang, listgen2[i] = E; for( j = 1, m, if( im[j,i], listgen2[i] = qfbcompraw(listgen2[i],listgen[j]))); if( norml2(im[,i]) > 1, red = QfbReduce(aux=mymat(listgen2[i])); aux = red~*aux*red; listgen2[i] = Qfb(aux[1,1],2*aux[1,2],aux[2,2])) ); listgen = listgen2; \\ listgen = vector(rang,i,listgen[m-rang+i]); struct = vector(rang,i,struct[m-rang+i]); if( DEBUGLEVEL >= 2, print("listgen = ",listgen)); if( DEBUGLEVEL >= 2, print("struct = ",struct)); return([struct,listgen]); } { {Qfsolve4(G,factD)= \\ Resolution de la forme quadratique quaternaire X^tGX=0 \\ On suppose que G est entiere et primitive. \\ Si on connait la factorisation de -abs(2*matdet(G)), \\ on peut lui donner pour gagner du temps. local(M,signG,vp,p,kerp,dimkerp,Q1,Q3,d,U,mask,Winvariants,signQ,clgp2,V,Q,G6,detG6,M6,setimG6,solG); \\ reduction des coefficients de G \\ et recherche d'une solution triviale M=IndefiniteLLL(G); if(length(M)==4, if( DEBUGLEVEL >= 1, print("solution triviale",M)); return(M)); G=M[1];M=M[2]; \\ factorisation du determinant if( DEBUGLEVEL >= 1, print1("factorisation du determinant")); if(!factD, factD=factor(-abs(2*matdet(G)))); if( DEBUGLEVEL >= 1, print(factD)); \\ cas particulier ou det(G)=0 if(!factD[1,1], if( DEBUGLEVEL >= 2, print("det(G) = 0")); return(matker(G)[,1])); factD[1,2]=0; factD[2,2]--; if( DEBUGLEVEL >= 3, print("factD=",factD)); \\ solubilite reelle signG=qfsign(G); if(signG[1]==4 || signG[1]==0, if( DEBUGLEVEL >= 1, print("pas de solution reelle")); return(0)); if( DEBUGLEVEL >= 1, print("reduction au cas sans facteur carre")); \\ on se ramene d'abord au cas ou det(G) est sans facteur carre \\ et on detecte la solubilite locale. for(i=2,matsize(factD)[1], vp=factD[i,2]; if(vp<2,next); p=factD[i,1]; if( DEBUGLEVEL >= 4, print("p = ",p)); kerp=kermodp(G,p);dimkerp=kerp[1];kerp=kerp[2]; \\ 1er cas: dimkerp = 3 if(dimkerp==3, if( DEBUGLEVEL >= 4, print("cas 1")); kerp[,4]*=p; M=M*kerp; G=(kerp~*G*kerp)/p; factD[i,2]-=2; i--;next; ); \\ 2eme cas: dimkerp = 1 if(dimkerp==1, if( DEBUGLEVEL >= 4, print("cas 2")); kerp[,1]/=p; M=M*kerp; G=(kerp~*G*kerp); factD[i,2]-=2; i--;next; ); \\ maintenant dimkerp = 2 \\ on se ramene au cas ou les 2 premieres lignes et colonnes \\ sont divisibles par p. M=M*kerp; G=kerp~*G*kerp; Q1=Mod(vecextract(G,3,3)/p,p); Q3=Mod(vecextract(G,12,12),p); \\ 3eme cas: dim kerp = 2 et -det(Q1) est un carre mod p d=-matdet(Q1); if(issquare(d), if( DEBUGLEVEL >= 4, print("cas 3")); U=matid(4); if(!Q1[1,1], U[1,1]=1/p; , U[1,2]=centerlift((-Q1[1,2]+sqrt(d))/Q1[1,1])/p; U[2,2]=1/p; ); M=M*U; G=U~*G*U; factD[i,2]-=2; i--;next; ); \\ 4eme cas: dim kerp = 2 et -det(Q3) est un carre mod p d=-matdet(Q3); if(issquare(d), if( DEBUGLEVEL >= 4, print("cas 4")); U=matid(4); if(!Q3[1,1], U[4,4]=p; , U[3,4]=centerlift((-Q3[1,2]+sqrt(d))/Q3[1,1]); U[3,3]=p; ); M=M*U; G=(U~*G*U)/p; factD[i,2]-=2; i--;next; ); \\ 5eme cas: dim kerp = 2 et -det(Q1) et -det(Q3) ne sont pas des carres \\ alors il n'y a pas de solution locale en p if( DEBUGLEVEL >= 1, print("pas de solution locale en ",p)); return(0) ); if( DEBUGLEVEL >= 4, print("G reduite = ",G)); \\ maintenant, on sait qu'il y a des solutions locales \\ (sauf peut-etre en 2), \\ et det(G) est sans facteur carre. \\ reduction de G et recherche de solution triviales \\ par exemple quand det G=+-1, il y en a toujours. U=IndefiniteLLL(G); if(length(U)==4, if( DEBUGLEVEL >= 1, print("solution triviale",M*U)); return(M*U)); G=U[1];M=M*U[2]; mask=[1,2]; for(i=3,matsize(factD)[1], if(factD[i,2],mask=concat(mask,[i]))); factD=vecextract(factD[,1],mask); if( DEBUGLEVEL >= 4, print("factD=",factD)); Winvariants=Qflisteinvariants(G,factD); if( DEBUGLEVEL >= 3, print("Witt invariants de G =",Winvariants)); Winvariants[,2]*=Mod(1,2); d=matdet(G); if( DEBUGLEVEL >= 3, print("-d=",-d)); \\ solubilite locale en 2 if(issquare(Mod(d,8)) && !Winvariants[2,2], if( DEBUGLEVEL >= 1, print("pas de solution locale en 2")); return(0) , if( DEBUGLEVEL >= 1, print("equation localement soluble partout")) ); \\ construction d'une forme quadratique binaire ayant les bons invariants. signQ=signG-[1,1]; if( DEBUGLEVEL >= 4, print("signQ=",signQ)); for(i=1,length(factD), if(hilbert(-1,-d,factD[i])==-1,Winvariants[i,2]++)): if( DEBUGLEVEL >= 1, print("Recherche d'un forme binaire d'invariants de Witt = ",lift(Winvariants))); if(d%4!=1, d*=4); factd = matrix(length(factD)-1,2); for(i=1,length(factD)-1, factd[i,1]=factD[i+1]; factd[i,2]=valuation(d,factd[i,1])); factd[1,2]++; clgp2=class2(d,factd); if( DEBUGLEVEL >= 4, print("clgp=",clgp2)); clgp2=clgp2[2]; U=matrix(length(factD),2+length(clgp2)); if( DEBUGLEVEL >= 4, print("U=", matsize(U))); for(j=1,length(clgp2), U[,j]=Qflisteinvariants(mymat(clgp2[j]),factD)[,2]); for(i=1,length(factD),U[i,length(clgp2)+1]=(1-hilbert(-1,d,factD[i]))/2); for(i=1,length(factD),U[i,length(clgp2)+2]=(1-hilbert( 2,d,factD[i]))/2); if( DEBUGLEVEL >= 4, printp("U=",U)); U=U*Mod(1,2); V=lift(matinverseimage(U,Winvariants[,2])); if( DEBUGLEVEL >= 4, print("V=",V)); if(!length(V), print("****** ERREUR: pas de forme quadratique ayant les bons invariants !!!!"); return(0)); \\ bug de pari a corriger !! \\ Q=prod(i=1,length(V)-2,clgp2[i]^V[i]);print("Q=",Q); Q=clgp2[1]^V[1];for(i=2,length(V)-2,Q=qfbcompraw(Q,clgp2[i]^V[i]));print("Q=",Q); Q=mymat(Q); if(V[length(V)-1],Q*=-1); if(V[length(V) ],Q*= 2); if( DEBUGLEVEL >= 2, print("Q=",Q)); if( DEBUGLEVEL >= 3, print("Witt invariants de Q=",Qflisteinvariants(Q,factD))); \\ construction d'une forme de dim=6 potentiellement unimodulaire G6=matrix(6,6); for(i=1,4,for(j=1,4,G6[i,j]=G[i,j])); G6[5,5]=-Q[1,1];G6[5,6]=G6[6,5]=-Q[1,2];G6[6,6]=-Q[2,2]; if( DEBUGLEVEL >= 4, print("G6=",G6)); if( DEBUGLEVEL >= 1, print("reduction de la forme quadratique de dimension 6")); \\ minimisation de G6 detG6=matdet(G6); M6=matid(6); if( DEBUGLEVEL >= 3, print("det(G6)=",factor(detG6))); for(i=2,length(factD), p=factD[i]; if(detG6%p,next); kerp=kermodp(G6,p); if(kerp[1]!=2, print("****** ERREUR: NOYAU DE MAUVAISE DIMENSION !! "); return(0)); M6=M6*kerp[2]; G6=kerp[2]~*G6*kerp[2]; Q1=Mod(vecextract(G6,3,3)/p,p); d=-matdet(Q1); if(!issquare(d), print("****** ERREUR: NON RESIDU QUADRATIQUE "); return(0)); U=matid(6); if(!Q1[1,1], U[1,1]=1/p; , U[1,2]=centerlift((-Q1[1,2]+sqrt(d))/Q1[1,1])/p; U[2,2]=1/p; ); M6=M6*U; G6=U~*G6*U; detG6/=p^2; i--;next; ); \\ maintenant det(G6)=-1 if( DEBUGLEVEL >= 4, print("G6=",G6)); if( DEBUGLEVEL >= 4, print("det=",matdet(G6))); \\ on construit un setim de G6 de dimension 3 if( DEBUGLEVEL >= 1, print("recherche d'un Setim de dimension 3")); setimG6 = setim(G6,1); if( DEBUGLEVEL >= 3, print("setimG6=",setimG6)); \\ ce setim rencontre le sous-espace matdiagonal([1,1,1,1,0,0]); if( DEBUGLEVEL >= 1, print("Reconstruction de la solution")); solG = matintersect(matrix(6,4,i,j,i==j),M6*setimG6)[,1]; if( DEBUGLEVEL >= 3, print("solG = ",solG)); solG = vecextract(solG,15); solG /= content(solG); if( DEBUGLEVEL >= 3, print("solG = ",solG)); if( DEBUGLEVEL >= 3, print("test solG = ", if(solG~*G*solG,"NOT ",""),"ok")); if( DEBUGLEVEL >= 1, print("fin de Qfsolve4")); return(M*solG); } { matdiagonalblock(v)= local(lv,lt); lv=length(v); lt=sum(i=1,lv,length(v[i])); M=matrix(lt,lt); lt=0; for(i=1,lv, for(j=1,length(v[i]), for(k=1,length(v[i]), M[lt+j,lt+k]=v[i][j,k])); lt+=length(v[i]); ); return(M) } \\ \\\\\\\\\\\\\\\\\ \\ \\ ESSAIS DIRECTS \\ \\ \\ \\\\\\\\\\\\\\\\\ \\diagonal4(301,210,-362,1113) \\ time = 55,170 ms (citron) \\ time = 28,760 ms (citron) \\diagonal4(30,43,159,-2534) \\ time = 8,300 ms (citron) \\ time = 4,250 ms (citron) \\M=matrix(5,5,i,j,random);M=(M+M~)\4;factor(matdet(M)) \\factor(matdet(Qfminim(M)[1])) \\ for(i=1,10,M=matrix(4,4,i,j,random(50));M=(M+M~)\4;print("M=",M);print("fa=",fa=factor(matdet(M)));print("sign=",qfsign(M));if(vecmax(fa[,2])==1 && vecmax(qfsign(M))<4,inva=Qflisteinvariants(M);inva[,2]+=vectorv(length(inva[,2]),j,hilbert(-1,-matdet(M),inva[j,1]));print(Qfbwithgiveninvariants(inva,-matdet(M))))) \\for(i=1,10,M=matrix(2,2,i,j,random(50));M=(M+M~)\4;print("M=",M);print("fa=",fa=factor(matdet(M)));if(vecmax(fa[,2])==1 , inva=Qflisteinvariants(M);print(Qfbwithgiveninvariants(inva,matdet(M))))); \\for(i=1,100,N=1<<5;a=rr(N);b=rr(N);c=rr(N);q=[a,b;b,c];if(issquarefree(matdet(q)),print(qq=Qfbwithgiveninvariants(li=Qflisteinvariants(q),matdet(q)));print("fini");if(li!=Qflisteinvariants(qq),print("ERREUR *****************",q);1/0))) \\rr(N)=random(2*N)-N; \\G=matrix(2,2,i,j,rr(10))/2;print(M=G+G~); \\M=[-12, -13; -13, 3]; \\inva=Qflisteinvariants(M); \\sol=Qfbwithgiveninvariants(inva,-4*matdet(M),20) \\addprimes(vp=[10^100+949,10^100+1293,10^100+2809,10^100+6637,10^100+22261]); \\D=prod(i=1,5,vp[i]); \\listgroups=[1,173,61,137,1129,433,-43,-8,-311,-359,-2711,-1663]; \\gettime();for(i=1,length(listgroups),print(listgroups[i],class2(D*listgroups[i])[1],gettime())); \\D2=541*557; \\listgroups2=[1,5,13,17,-3,-7,-11,-15,-19]; \\gettime();for(i=1,length(listgroups2),print(listgroups2[i],class2(D2*listgroups2[i])[1],gettime())); \\ infinite loop ???? \\diagonal4(23,prime(1000),prime(1009),-prime(194)) \\diagonal4(23,7919,8011,-1181) \\Qfsolve4(M) \\rr(N)=random(2*N)-N; \\n=21;gettime();timefact=0;timesolve=0;for(k=1,100,G=matrix(4,4,i,j,rr(2^n));print(G=(G+G~)\2);G/=content(G);fa=factor(-abs(2*matdet(G)));timefact+=gettime();sol=Qfsolve4(G,fa);timesolve+=gettime();if(sol && sol~*G*sol !=0,1/0,print(sol)));print("n=",n,"; solve = ",timesolve,"; fact = ",timefact); \\n=4; solve = 2830: fact = 70 \\n=5; solve = 4240; fact = 70 \\n=6; solve = 6230; fact = 20 \\n=7; solve = 7450; fact = 70 \\n=8; solve = 10220; fact = 60 \\n=9; solve = 11190; fact = 120 \\n=10; solve = 12550; fact = 120 \\n=11; solve = 11640; fact = 110 \\n=12; solve = 13880; fact = 110 \\n=13; solve = 15510; fact = 150 \\n=14; solve = 16130; fact = 120 \\n=15; solve = 16870; fact = 120 \\n=16; solve = 18410; fact = 320 \\n=17; solve = 18120; fact = 390 \\n=18; solve = 18540; fact = 430 \\n=19; solve = 21160; fact = 750 \\n=20; solve = 21680; fact = 890 \\n=21; solve = 23600; fact = 970 \\n=22; solve = 23700; fact = 1640 \\n=23; solve = 23730; fact = 1500 \\n=24; solve = 25410; fact = 1930 \\n=25; solve = 28360; fact = 1970 \\n=26; solve = 27420; fact = 2880 \\n=27; solve = 29150; fact = 3550 \\n=28; solve = 29610; fact = 4690 \\n=29; solve = 31700; fact = 6350 \\n=30; solve = 34770; fact = 5580 \\n=31; solve = 33980; fact = 7530 \\n=32; solve = 35670; fact = 10160 \\n=33; solve = 37130; fact = 14700 \\n=34; solve = 39960; fact = 19330 \\n=35; solve = 40000; fact = 29980 \\n=36; solve = 40340; fact = 25930 \\n=37; solve = 42470; fact = 44350 \\n=38; solve = 39720; fact = 62030 \\n=39; solve = 44500; fact = 34110 \\n=40; solve = 44970; fact = 99110 \\n=41; solve = 45150; fact = 154480 \\n=42; solve = 46530; fact = 154390 \\n=43; solve = 50330; fact = 220970 \\n=44; solve = 51600; fact = 259870 \\n=45; solve = 49800; fact = 363540 \\n=46; solve = 52290; fact = 356630 \\n=47; solve = 56650; fact = 757780 \\n=48; solve = 53240; fact = 666500 \\n=49; solve = 59600; fact = 1170430 \\n=50; solve = 56320; fact = 1572480 \\n=51; solve = 58540; fact = 1800840 \\n=52; solve = 61240; fact = 3005540 \\n=53; solve = 59490; fact = 3201520 \\n=54; solve = 60370; fact = 5133560 \\n=55; solve = 62260; fact = 5014180 \\n=56; solve = 66720; fact = 8054360 \\n=57; solve = 65450; fact = 11194240 \\n=58; solve = 70650; fact = 7931900 \\n=59; solve = 71410; fact = 19075400 \\ \\ \\ \\{v=[2830,4240,6230,7450,10220,11190 \\,12550,11640,13880,15510,16130,16870,18410,18120,18540,21160 \\,21680,23600,23700,23730,25410,28360,27420,29150,29610,31700 \\,34770,33980,35670,37130,39960,40000,40340,42470,39720,44500 \\,44970,45150,46530,50330,51600,49800,52290,56650,53240,59600 \\,56320,58540,61240,59490,60370,62260,66720,65450,70650,71410 \\];plothraw(vector(length(v),i,3+i),v,1)} \\ \\