Ticket #6955: qfsolve.gp

File qfsolve.gp, 31.1 KB (added by John Cremona, 13 years ago)
Line 
1\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
2\\       Copyright (C) 2007 Denis Simon
3\\
4\\ Distributed under the terms of the GNU General Public License (GPL)
5\\
6\\    This code is distributed in the hope that it will be useful,
7\\    but WITHOUT ANY WARRANTY; without even the implied warranty of
8\\    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
9\\    General Public License for more details.
10\\
11\\ The full text of the GPL is available at:
12\\
13\\                 http://www.gnu.org/licenses/
14\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
15
16\\
17\\ Auteur:
18\\ Denis SIMON -> simon@math.unicaen.fr
19\\ adresse du fichier:
20\\ www.math.unicaen.fr/~simon/qfsolve.gp
21\\
22\\  *********************************************         
23\\  *          VERSION 02/10/2009               *
24\\  *********************************************         
25\\
26\\ Programme de resolution des equations quadratiques
27\\ langage: GP
28\\ pour l'utiliser, lancer gp, puis taper
29\\ \r qfsolve.gp
30\\
31\\ Ce fichier contient 4 principales fonctions:
32\\
33\\ - Qfsolve(G,factD): pour resoudre l'equation quadratique X^t*G*X = 0
34\\ G doit etre une matrice symetrique n*n, a coefficients dans Z.
35\\ S'il n'existe pas de solution, la reponse est un entier
36\\ indiquant un corps local dans lequel aucune solution n'existe
37\\ (-1 pour les reels, p pour Q_p).
38\\ Si on connait la factorisation de -abs(2*matdet(G)),
39\\ on peut la passer par le parametre factD pour gagner du temps.
40\\
41\\ - Qfparam(G,sol,fl): pour parametrer les solutions de la forme
42\\ quadratique ternaire G, en utilisant la solution particuliere sol.
43\\ si fl>0, la 'fl'eme forme quadratique est reduite.
44\\
45\\ - IndefiniteLLL(G,c): pour resoudre ou reduire la forme quadratique
46\\ G a coefficients entiers. Il s'agit d'un algorithme type LLL, avec la
47\\ constante 1/4<c<1.
48\\
49\\ - class2(d,factd): determine le 2-Sylow du (narrow) groupe de classes de
50\\ discriminant d, ou d est un discriminant fondamental.
51\\ Si on connait la factorisation de abs(2*d),
52\\ on peut la donner dans factd, et dans ce cas le reste
53\\ de l'algorithme est polynomial.
54\\
55
56\\
57DEBUGLEVEL_qfsolve = 0;
58
59\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
60\\
61\\  DEBUT DES SOURCES                \\
62\\
63\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
64{QfbReduce(M) =
65\\ M = [a,b;b;c] est a coefficients entiers.
66\\ Reduction de la forme quadratique binaire
67\\   qf = (a,b,c)=a*X^2+2*b*X*Y+c*Y^2
68\\ Renvoit la matrice de reduction de det = +1.
69
70local(a,b,c,H,test,di,q,r,nexta,nextb,nextc,aux);
71
72if( DEBUGLEVEL_qfsolve >= 5, print("entree dans QfbReduce avec ",M));
73
74  a = M[1,1]; b = M[1,2]; c = M[2,2];
75 
76  H = matid(2); test = 1;
77  while( test && a,
78    di = divrem(b,a); q = di[1]; r = di[2];
79    if( 2*r > abs(a), r -= abs(a); q += sign(a));
80    H[,2] -= q*H[,1];
81    nextc = a; nextb = -r; nexta= (nextb-b)*q+c;
82 
83    if( test = abs(nexta) < abs(a),
84      c = nextc; b = nextb; a = nexta;
85      aux = H[,1]; H[,1] = -H[,2]; H[,2] = aux
86    )
87  );
88
89if( DEBUGLEVEL_qfsolve >= 5, print("sortie de QfbReduce avec ",H));
90return(H);
91}
92{IndefiniteLLL(G,c=1,base=0) =
93\\ Performs first a LLL reduction on a positive definite
94\\ quadratic form QD bounding the indefinite G.
95\\ Then finishes the reduction with IndefiniteLLL2.
96
97local(n,M,QD,M1,S,red);
98
99  n = length(G);
100  M = matid(n);
101  QD = G;
102  for( i = 1, n-1,
103    if( !QD[i,i],
104return(IndefiniteLLL2(G,c,base))
105    );
106    M1 = matid(n);
107    M1[i,] = -QD[i,]/QD[i,i];
108    M1[i,i] = 1;
109    M = M*M1;
110    QD = M1~*QD*M1
111  );
112  M = M^(-1);
113  QD = M~*abs(QD)*M;
114  S = qflllgram(QD/content(QD));
115  red = IndefiniteLLL2(S~*G*S,c,base);
116  if( type(red) == "t_COL",
117return(S*red));
118  if( length(red) == 3,
119return([red[1],S*red[2],S*red[3]]));
120return([red[1],S*red[2]]);
121}
122{IndefiniteLLL2(G,c=1,base=0) =
123\\ following Cohen's book p. 86
124\\ but without b and bstar: works on G
125\\ returns [H~*G*H,H] where det(H) = 1 and H~*G*H is reduced.
126\\ Exit with a norm 0 vector if one such is found.
127\\ If base == 1 and norm 0 is obtained, returns [H~*G*H,H,sol] where
128\\   sol is a norm 0 vector and is the 1st column of H.
129
130local(n,H,M,A,aux,sol,k,nextk,swap,q,di,HM,aux1,aux2,Mkk1,bk1new,Mkk1new,newG);
131
132  n = length(G);
133if( DEBUGLEVEL_qfsolve >= 3, print("LLL dim ",n," avec G=",log(vecmax(abs(G)))/log(10)));
134if( DEBUGLEVEL_qfsolve >= 4, print("LLL avec ");printp(G));
135
136 if( n <= 1, return([G,matid(n)]));
137
138  H = M = matid(n); A = matrix(n,n);
139
140\\ compute Gram-Schmidt 
141
142  for( i = 1, n,
143    if( !(A[i,i] = G[i,i]),
144      if( base,
145        aux = H[,1]; H[,1] = H[,i]; H[,i] = -aux;
146        return([H~*G*H,H,H[,1]])
147      , return(M[,i])));
148    for( j = 1, i-1,
149      A[i,j] = G[i,j] - sum( k = 1, j-1, M[j,k]*A[i,k]);
150      M[i,j] = A[i,j]/A[j,j];
151      A[i,i] -= M[i,j]*A[i,j];
152      if( !A[i,i],
153        sol = (M^(-1))~[,i]; sol /= content(sol);
154        if( base,
155          H = completebasis(sol);
156          aux = H[,1]; H[,1] = H[,n]; H[,n]= -aux;
157          return([H~*G*H,H,H[,1]])
158        , return(sol)))
159    )
160  );
161
162\\ LLL loop
163
164  k = 2; nextk = 1;
165  while( k <= n,
166
167    swap = 1;
168    while( swap,
169      swap = 0;
170
171\\ red(k,k-1);
172      if( q = round(M[k,k-1]),
173        for( i = 1, k-2, M[k,i] -= q*M[k-1,i]);
174        M[k,k-1] -= q;
175        for( i = 1, n,
176          A[k,i] -= q*A[k-1,i];
177          H[i,k] -= q*H[i,k-1]
178        )
179      );
180
181\\ preparation du swap(k,k-1)
182
183      if( issquare( di = -A[k-1,k-1]*A[k,k]),
184\\ di est le determinant de matr
185\\ On trouve une solution
186        HM = (M^(-1))~;
187        aux1 = sqrtint(numerator(di));
188        aux2 = sqrtint(denominator(di));
189        sol = aux1*HM[,k-1]+aux2*A[k-1,k-1]*HM[,k];
190        sol /= content(sol);
191        if( base,
192          H = H*completebasis(sol,1);
193          aux = H[,1]; H[,1] = H[,n]; H[,n] = -aux;
194          return([H~*G*H,H,H[,1]])
195        , return(H*sol)
196        )
197      );
198
199\\ On reduit [k,k-1].
200      Mkk1 = M[k,k-1];
201      bk1new = Mkk1^2*A[k-1,k-1] + A[k,k];
202      if( swap = abs(bk1new) < c*abs(A[k-1,k-1]),
203        Mkk1new = -Mkk1*A[k-1,k-1]/bk1new
204      );
205
206\\ Calcul des nouvelles matrices apres le swap.
207      if( swap,
208        for( j = 1, n,
209          aux = H[j,k-1]; H[j,k-1] = H[j,k]; H[j,k] = -aux);
210        for( j = 1, k-2,
211          aux = M[k-1,j]; M[k-1,j] = M[k,j]; M[k,j] = -aux);
212        for( j = k+1, n,
213          aux = M[j,k]; M[j,k] = -M[j,k-1]+Mkk1*aux; M[j,k-1] = aux+Mkk1new*M[j,k]);
214        for( j = 1, n, if( j != k && j != k-1,
215          aux = A[k-1,j]; A[k-1,j] = A[k,j]; A[k,j] =- aux;
216          aux = A[j,k-1]; A[j,k-1] = Mkk1*aux+A[j,k]; A[j,k] = -aux-Mkk1new*A[j,k-1]));
217
218        aux1 = A[k-1,k-1];
219        aux2 = A[k,k-1];
220        A[k,k-1]  =-A[k-1,k] - Mkk1*aux1;
221        A[k-1,k-1]= A[k,k]   + Mkk1*aux2;
222        A[k,k]    = aux1     - Mkk1new*A[k,k-1];
223        A[k-1,k]  =-aux2     - Mkk1new*A[k-1,k-1];
224
225        M[k,k-1] = Mkk1new;
226
227if( DEBUGLEVEL_qfsolve >=4, newG=H~*G*H;print(vector(n,i,matdet(vecextract(newG,1<<i-1,1<<i-1)))));
228
229        if( k != 2, k--)
230      )
231    );
232
233    forstep( l = k-2, 1, -1,
234\\ red(k,l)
235      if( q = round(M[k,l]),
236        for( i = 1, l-1, M[k,i] -= q*M[l,i]);
237        M[k,l] -= q;
238        for( i = 1, n,
239          A[k,i] -= q*A[l,i];
240          H[i,k] -= q*H[i,l]
241        )
242      )
243    );
244    k++
245  );
246return([H~*G*H,H]);
247}
248{kermodp(M,p) =
249\\ Compute the kernel of M mod p.
250\\ returns [d,U], where
251\\ d = dim (ker M mod p)
252\\ U in SLn(Z), and its first d columns span the kernel.
253
254local(n,U,d);
255
256  n = length(M);
257  U = centerlift(matker(M*Mod(1,p)));
258  d = length(U);
259  U = completebasis(U);
260  U = matrix(n,n,i,j,U[i,n+1-j]);
261return([d,U]);
262}
263{Qfparam(G,sol,fl=3) =
264\\ G est une matrice symetrique 3*3, et sol une solution de sol~*G*sol=0.
265\\ Renvoit une parametrisation des solutions avec de bons invariants,
266\\ sous la forme d'une matrice 3*3, dont chaque ligne contient
267\\ les coefficients de chacune des 3 formes quadratiques.
268
269\\ si fl!=0, la 'fl'eme forme quadratique est reduite.
270
271local(U,G1,G2);
272
273if( DEBUGLEVEL_qfsolve >= 5, print("entree dans Qfparam"));
274  sol /= content(sol);
275\\ construction de U telle que U[,3] est sol, et det(U) = +-1
276  U = completebasis(sol,1);
277  G1 = U~*G*U; \\ G1 a un 0 en bas a droite.
278  G2 = [-2*G1[1,3],-2*G1[2,3],0;
279      0,-2*G1[1,3],-2*G1[2,3];
280      G1[1,1],2*G1[1,2],G1[2,2]];
281  sol = U*G2;
282  if(fl && !issquare(matdet( U = [sol[fl,1],sol[fl,2]/2;
283                                  sol[fl,2]/2,sol[fl,3]])),
284    U = QfbReduce(U);
285    U = [U[1,1]^2,2*U[1,2]*U[1,1],U[1,2]^2;
286         U[2,1]*U[1,1],U[2,2]*U[1,1]+U[2,1]*U[1,2],U[1,2]*U[2,2];
287         U[2,1]^2,2*U[2,1]*U[2,2],U[2,2]^2];
288    sol = sol*U
289  );
290if( DEBUGLEVEL_qfsolve >= 5, print("sortie de Qfparam"));
291return(sol);
292}
293{LLLgoon3(G,c=1) =
294\\ reduction LLL de la forme quadratique G (matrice de Gram)
295\\ en dim 3 seulement avec detG = -1 et sign(G) = [1,2];
296
297local(red,U1,G2,bez,U2,G3,cc,U3);
298
299  red = IndefiniteLLL(G,c,1);
300\\ On trouve toujours un vecteur isotrope.
301  U1 = [0,0,1;0,1,0;1,0,0];
302  G2 = U1~*red[1]*U1;
303\\ la matrice G2 a un 0 dans le coin en bas a droite.
304  bez = bezout(G2[3,1],G2[3,2]);
305  U2 = [bez[1],G2[3,2]/bez[3],0;bez[2],-G2[3,1]/bez[3],0;0,0,-1];
306  G3 = U2~*G2*U2;
307\\ la matrice G3 a des 0 sous l'anti-diagonale.
308  cc = G3[1,1]%2;
309  U3 = [1,0,0;  cc,1,0;
310        round(-(G3[1,1]+cc*(2*G3[1,2]+G3[2,2]*cc))/2/G3[1,3]),
311        round(-(G3[1,2]+cc*G3[2,2])/G3[1,3]),1];
312return([U3~*G3*U3,red[2]*U1*U2*U3]);
313}
314{completebasis(v,redflag=0) =
315\\ Donne une matrice unimodulaire dont la derniere colonne est v.
316\\ Si redflag <> 0, alors en plus on reduit le reste.
317
318local(U,n,re);
319
320  U = (mathnf(Mat(v~),1)[2]~)^-1;
321  n = length(v~);
322  if( n==1 || !redflag, return(U));
323  re = qflll(vecextract(U,1<<n-1,1<<(n-1)-1));
324return( U*matdiagonalblock([re,Mat(1)]));
325}
326{LLLgoon(G,c=1) =
327\\ reduction LLL de la forme quadratique G (matrice de Gram)
328\\ ou l'on continue meme si on trouve un vecteur isotrope
329
330local(red,U1,G2,U2,G3,n,U3,G4,U,V,B,U4,G5,U5,G6);
331
332  red = IndefiniteLLL(G,c,1);
333\\ si on ne trouve pas de vecteur isotrope, il n'y a rien a faire.
334  if( length(red) == 2, return(red));
335\\ sinon :
336  U1 = red[2];
337  G2 = red[1]; \\ On a G2[1,1] = 0
338  U2 = mathnf(Mat(G2[1,]),4)[2];
339  G3 = U2~*G2*U2;
340\\ la matrice de G3 a des 0 sur toute la 1ere ligne,
341\\ sauf un 'g' au bout a droite, avec g^2| det G.
342  n = length(G);
343  U3 = matid(n); U3[1,n] = round(-G3[n,n]/G3[1,n]/2);
344  G4 = U3~*G3*U3;
345\\ le coeff G4[n,n] est reduit modulo 2g
346  U = vecextract(G4,[1,n],[1,n]);
347  if( n == 2,
348    V = matrix(2,0)
349  , V = vecextract(G4,[1,n],1<<(n-1)-2));
350  B = round(-U^-1*V);
351  U4 = matid(n);
352  for( j = 2, n-1,
353    U4[1,j] = B[1,j-1];
354    U4[n,j] = B[2,j-1]
355  );
356  G5 = U4~*G4*U4;
357\\ la derniere colonne de G5 est reduite
358  if( n < 4, return([G5,U1*U2*U3*U4]));
359
360  red = LLLgoon(matrix(n-2,n-2,i,j,G5[i+1,j+1]),c);
361  U5 = matdiagonalblock([Mat(1),red[2],Mat(1)]);
362  G6 = U5~*G5*U5;
363return([G6,U1*U2*U3*U4*U5]);
364}
365{QfWittinvariant(G,p) =
366\\ calcule l'invariant c (=invariant de Witt) d'une forme quadratique,
367\\ p-adique (reel si p = -1)
368
369local(n,det,diag,c);
370
371  n = length(G);
372\\ On diagonalise d'abord G
373  det  = vector( n+1, i, matdet(matrix(i-1,i-1,j,k,G[j,k])));
374  diag = vector( n, i, det[i+1]/det[i]);
375
376\\ puis on calcule les symboles de Hilbert
377  c = prod( i = 1, n,
378        prod( j = i+1, n,
379          hilbert( diag[i], diag[j], p)));
380return(c);
381}
382{Qflisteinvariants(G,fa=[]) =
383\\ G est une forme quadratique, ou une matrice symetrique,
384\\ ou un vecteur contenant des formes quadratiques de meme discriminant.
385\\ fa = factor(-abs(2*matdet(G)))[,1] est facultatif.
386
387local(l,sol,n,det);
388
389if( DEBUGLEVEL_qfsolve >= 4, print("entree dans Qflisteinvariants",G));
390  if( type(G) != "t_VEC", G = [G]);
391  l = length(G);
392  for( j = 1, l,   
393    if( type(G[j]) == "t_QFI" || type(G[j]) == "t_QFR",
394      G[j] = mymat(G[j])));
395
396  if( !length(fa),
397    fa = factor(-abs(2*matdet(G[1])))[,1]);
398
399  if( length(G[1]) == 2,
400\\ En dimension 2, chaque invariant est un unique symbole de Hilbert.
401    det = -matdet(G[1]);
402    sol = matrix(length(fa),l,i,j,hilbert(G[j][1,1],det,fa[i])<0);
403if( DEBUGLEVEL_qfsolve >= 4, print("sortie de Qflisteinvariants"));
404    return([fa,sol])
405  );
406
407  sol = matrix(length(fa),l);
408  for( j = 1, l,
409    n = length(G[j]);
410\\ En dimension n, on calcule un produit de n symboles de Hilbert.
411    det = vector(n+1, i, matdet(matrix(i-1,i-1,k,m,G[j][k,m])));
412    for( i = 1, length(fa),
413      sol[i,j] = prod( k = 1, n-1, hilbert(-det[k],det[k+1],fa[i]))*hilbert(det[n],det[n+1],fa[i]) < 0;
414    )
415  );
416if( DEBUGLEVEL_qfsolve >= 4, print("sortie de Qflisteinvariants"));
417return([fa,sol]);
418}
419{Qfsolvemodp(G,p) =
420\\ p a prime number.
421\\ finds a solution mod p for the quadatic form G
422\\ such that det(G) !=0 mod p and dim G = n>=3;
423local(n,vdet,G2,sol,x1,x2,x3,N1,N2,N3,s,r);
424
425  n = length(G);
426  vdet = [0,0,0];
427  for( i = 1, 3,
428    G2 = vecextract(G,1<<i-1,1<<i-1)*Mod(1,p);
429    if( !(vdet[i] = matdet(G2)),
430      sol = kermodp(lift(G2),p)[2][,1];
431      sol = vectorv(n, j, if( j <= i, sol[j]));
432      return(sol)
433    )
434  );
435\\ maintenant, on resout en dimension 3...
436\\ d'abord, on se ramene au cas diagonal:
437  x1 = [1,0,0]~;
438  x2 = [-G2[1,2],G2[1,1],0]~;
439  x3 = [G2[2,2]*G2[1,3]-G2[2,3]*G2[1,2],G2[1,1]*G2[2,3]-G2[1,3]*G2[1,2],G2[1,2]^2-G2[1,1]*G2[2,2]]~;
440  while(1,
441    if( issquare( N1 = -vdet[2]),                 s = sqrt(N1); sol = s*x1+x2; break);
442    if( issquare( N2 = -vdet[3]/vdet[1]),         s = sqrt(N2); sol = s*x2+x3; break);
443    if( issquare( N3 = -vdet[2]*vdet[3]/vdet[1]), s = sqrt(N3); sol = s*x1+x3; break);
444    r = 1;
445    while( !issquare( s = (1-N1*r^2)/N3), r = random(p));
446    s = sqrt(s); sol = x1+r*x2+s*x3; break
447  );
448  sol = vectorv(n, j, if( j <= 3, sol[j]));
449return(sol);
450}
451{Qfminim(G,factdetG=0) =
452\\ Minimisation de la forme quadratique entiere non degeneree G,
453\\   de dimension n >=2. On suppose que G est symetrique a coefficients entiers.
454\\ Renvoit [G',U,factd] avec U \in GLn(Q) telle que G'=U~*G*U*constante est entiere
455\\   de determinant minimal.
456\\ Renvoit p si la reduction est impossible a cause de la non-solubilite locale
457\\   en p (seulement en dimension 3-4).
458\\ factd est la factorisation de abs(det(G'))
459\\ Si on connait la factorisation de matdet(G), on peut la donner en parametre
460\\   pour gagner du temps.
461local(n,U,factd,detG,i,vp,Ker,dimKer,Ker2,dimKer2,sol,aux,p,di,m);
462
463  n = length(G);
464  U = matid(n);
465  factd = matrix(0,2);
466  detG = matdet(G);
467  if( !factdetG, factdetG = factor(detG));
468  i = 1;
469  while(i <= length(factdetG[,1]),
470    p = factdetG[i,1];
471    if( p == -1, i++; next);
472if( DEBUGLEVEL_qfsolve >= 4, print("p=",p,"^",factdetG[i,2]));
473    vp = factdetG[i,2];
474    if( vp == 0, i++; next);
475\\ Le cas vp=1 n'est minimisable que si n est impair.
476    if( vp == 1 && n%2 == 0,
477      factd = concat(factd~, Mat([p,1])~)~;
478      i++; next
479    );
480    Ker = kermodp(G,p); dimKer = Ker[1]; Ker = Ker[2];
481\\ Rem: on a toujours dimKer <= vp
482if( DEBUGLEVEL_qfsolve >= 4, print("dimKer = ",dimKer));
483\\ cas trivial: G est divisible par p.
484    if( dimKer == n,
485      G /= p;
486      factdetG[i,2] -= n;
487      next
488    );
489    G = Ker~*G*Ker;
490    U = U*Ker;
491\\ 1er cas: la dimension du noyau est plus petite que la valuation
492\\ alors le noyau mod p contient un noyau mod p^2
493    if( dimKer < vp,
494if( DEBUGLEVEL_qfsolve >= 4, print("cas 1"));
495      Ker2 = kermodp(matrix(dimKer,dimKer,j,k,G[j,k]/p),p);
496      dimKer2 = Ker2[1]; Ker2 = Ker2[2];
497      for( j = 1, dimKer2, Ker2[,j] /= p);
498      Ker2 = matdiagonalblock([Ker2,matid(n-dimKer)]);
499      G = Ker2~*G*Ker2;
500      U = U*Ker2;
501      factdetG[i,2] -= 2*dimKer2;
502if( DEBUGLEVEL_qfsolve >= 4, print("fin cas 1"));
503      next
504    );
505
506\\ Maintenant, vp = dimKer
507\\ 2eme cas: la dimension du noyau est >=2 et contient un element de norme 0 mod p^2
508    if( dimKer > 2 ||
509       (dimKer == 2 && issquare(di=Mod((G[1,2]^2-G[1,1]*G[2,2])/p^2,p))),
510\\ on cherche dans le noyau un elt de norme p^2...
511      if( dimKer > 2,
512if( DEBUGLEVEL_qfsolve >= 4, print("cas 2.1"));
513        dimKer = 3;
514        sol = Qfsolvemodp(matrix(3,3,j,k,G[j,k]/p),p)
515      , 
516if( DEBUGLEVEL_qfsolve >= 4, print("cas 2.2"));
517        if( G[1,1]%p^2 == 0,
518          sol = [1,0]~
519        , sol = [-G[1,2]/p+sqrt(di),Mod(G[1,1]/p,p)]~
520        )
521      );
522      sol = centerlift(sol); sol /= content(sol);
523if( DEBUGLEVEL_qfsolve >= 4, print("sol=",sol));
524      Ker = vectorv(n, j, if( j<= dimKer, sol[j], 0)); \\ on complete avec des 0
525      Ker = completebasis(Ker,1);
526      Ker[,n] /= p;
527      G = Ker~*G*Ker;
528      U = U*Ker;
529      factdetG[i,2] -= 2;
530if( DEBUGLEVEL_qfsolve >= 4, print("fin cas 2"));
531      next
532    );
533
534\\ Maintenant, vp = dimKer <= 2
535\\   et le noyau ne contient pas de vecteur de norme p^2...
536
537\\ Dans certains cas, on peut echanger le noyau et l'image
538\\   pour continuer a reduire.
539
540    m = (n-1)\2-1;
541    if( ( vp == 1 && issquare(Mod(-(-1)^m*matdet(G)/G[1,1],p)))
542     || ( vp == 2 && n%2 == 1 && n >= 5)
543     || ( vp == 2 && n%2 == 0 && !issquare(Mod((-1)^m*matdet(G)/p^2,p)))
544    ,
545if( DEBUGLEVEL_qfsolve >= 4, print("cas 3"));
546      Ker = matid(n);
547      for( j = dimKer+1, n, Ker[j,j] = p);
548      G = Ker~*G*Ker/p;
549      U = U*Ker;
550      factdetG[i,2] -= 2*dimKer-n;
551if( DEBUGLEVEL_qfsolve >= 4, print("fin cas 3"));
552      next
553    );
554
555\\ On n'a pas pu minimiser.
556\\ Si n == 3 ou n == 4 cela demontre la non-solubilite locale en p.
557    if( n == 3 || n == 4,
558if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution locale en ",p));
559      return(p));
560
561if( DEBUGLEVEL_qfsolve >= 4, print("plus de minimisation possible"));
562    factd = concat(factd~,Mat([p,vp])~)~;
563    i++
564  );
565\\print("Un=",log(vecmax(abs(U))));
566  aux = qflll(U);
567\\print("Ur=",log(vecmax(abs(U*aux))));
568return([aux~*G*aux,U*aux,factd]);
569}
570{mymat(qfb) = qfb = Vec(qfb);[qfb[1],qfb[2]/2;qfb[2]/2,qfb[3]];
571}
572{Qfbsqrtgauss(G,factdetG) =
573\\ pour l'instant, ca ne marche que pour detG sans facteur carre
574\\ sauf en 2, ou la valuation est 2 ou 3.
575\\ factdetG contient la factorisation de 2*abs(disc G)
576local(a,b,c,d,m,n,p,aux,Q1,M);
577if( DEBUGLEVEL_qfsolve >=3, print("entree dans Qfbsqrtgauss avec",G,factdetG));
578  G = Vec(G);
579  a = G[1]; b = G[2]/2; c = G[3]; d = a*c-b^2;
580
581\\ 1ere etape: on resout m^2 = a (d), m*n = -b (d), n^2 = c (d)
582  m = n = Mod(1,1);
583  factdetG[1,2] -= 3;
584  for( i = 1, length(factdetG[,1]),
585    if( !factdetG[i,2], next);
586    p = factdetG[i,1];
587    if( gcd(a,p) == 1,
588      aux = sqrt(Mod(a,p));
589      m = chinese(m,aux);
590      n = chinese(n,-b/aux)
591    ,
592      aux = sqrt(Mod(c,p));
593      n = chinese(n,aux);
594      m = chinese(m,-b/aux)
595    )
596  );
597  m = centerlift(m);  n = centerlift(n);
598if( DEBUGLEVEL_qfsolve >=4, print("m=",m); print("n=",n));
599 
600\\ 2eme etape: on construit Q1 de det -1 tq Q1(x,y,0) = G(x,y)
601  Q1 = [(n^2-c)/d, (m*n+b)/d, n ;
602        (m*n+b)/d, (m^2-a)/d, m ;
603        n,         m,         d ];
604  Q1 = -matadjoint(Q1);
605
606\\ 3eme etape: on reduit Q1 jusqu'a [0,0,-1;0,1,0;-1,0,0]
607  M = LLLgoon3(Q1)[2][3,];
608  if( M[1] < 0, M = -M);
609if( DEBUGLEVEL_qfsolve >=3, print("fin de Qfbsqrtgauss "));
610  if( M[1]%2,
611    return(Qfb(M[1],2*M[2],2*M[3]))
612  , return(Qfb(M[3],-2*M[2],2*M[1])));
613}
614{class2(D,factdetG,Winvariants,U2) =
615\\ On suit l'algorithme de Shanks/Bosma-Stevenhagen
616\\ pour construire tout le 2-groupe de classes
617\\ seulement pour D discriminant fondamental.
618\\ lorsque D = 1(4), on travaille avec 4D.
619\\ Si on connait la factorisation de abs(2*D),
620\\ on peut la donner dans factdetG, et dans ce cas le reste
621\\ de l'algorithme est polynomial.
622\\ Si Winvariants est donne, l'algorithme s'arrete des qu'un element d'invariants=Winvariants
623\\ est trouve.
624
625local(factD,n,rang,m,listgen,vD,p,vp,aux,invgen,im,Ker,Kerim,listgen2,G2,struct,E,compteur,red);
626
627if( DEBUGLEVEL_qfsolve >= 1, print("Construction du 2-groupe de classe de discriminant ",D));
628  if( D%4 == 2 || D%4 == 3, print("Discriminant not congruent to 0,1 mod 4");return(0));
629
630  if( D==-4, return([[1],[Qfb(1,0,1)]]));
631
632  if( !factdetG, factdetG = factor(2*abs(D)));
633  factD = concat([-1],factdetG[,1]~);
634  if( D%4 == 1, D *= 4; factdetG[1,2] += 2);
635 
636  n = length(factD); rang = n-3;
637  if(D>0, m = rang+1, m = rang);
638if( DEBUGLEVEL_qfsolve >= 3, print("factD = ",factD));
639  listgen = vector(max(0,m));
640 
641  if( vD = valuation(D,2),
642    E = Qfb(1,0,-D/4)
643  , E = Qfb(1,1,(1-D)/4)
644  );
645if( DEBUGLEVEL_qfsolve >= 3, print("E = ",E));
646
647  if( type(Winvariants) == "t_COL" && (Winvariants == 0 || length(matinverseimage(U2*Mod(1,2),Winvariants))>0), return([[1],[E]]));
648 
649  for( i = 1, m, \\ on  ne regarde pas factD[1]=-1, ni factD[2]=2
650    p = factD[i+2];
651    vp = valuation(D,p);
652    aux = p^vp;
653    if( vD,
654      listgen[i] = Qfb(aux,0,-D/4/aux)
655    , listgen[i] = Qfb(aux,aux,(aux-D/aux)/4))
656  );
657  if( vD == 2 && D%16 != 4,
658    m++; rang++; listgen = concat(listgen,[Qfb(2,2,(4-D)/8)]));
659  if( vD == 3,
660    m++; rang++; listgen = concat(listgen,[Qfb(2^(vD-2),0,-D/2^vD)]));
661 
662if( DEBUGLEVEL_qfsolve >= 3, print("listgen = ",listgen));
663if( DEBUGLEVEL_qfsolve >= 2, print("rang = ",rang));
664
665  if( !rang, return([[1],[E]]));
666
667 invgen = Qflisteinvariants(listgen,factD)[2]*Mod(1,2);
668if( DEBUGLEVEL_qfsolve >= 3, printp("invgen = ",lift(invgen)));
669
670  struct = vector(m,i,2);
671  im = lift(matinverseimage(invgen,matimage(invgen)));
672  while( (length(im) < rang)
673  || (type(Winvariants) == "t_COL" && length(matinverseimage(concat(invgen,U2),Winvariants) == 0)),
674    Ker = lift(matker(invgen));
675    Kerim = concat(Ker,im);
676    listgen2 = vector(m);
677    for( i = 1, m,
678      listgen2[i] = E;
679      for( j = 1, m,
680        if( Kerim[j,i],
681          listgen2[i] = qfbcompraw(listgen2[i],listgen[j])));
682      if( norml2(Kerim[,i]) > 1,
683        red = QfbReduce(aux=mymat(listgen2[i]));
684        aux = red~*aux*red;
685        listgen2[i] = Qfb(aux[1,1],2*aux[1,2],aux[2,2]))
686    );
687    listgen = listgen2;
688    invgen = invgen*Kerim;
689
690if( DEBUGLEVEL_qfsolve >= 4, print("listgen = ",listgen));
691if( DEBUGLEVEL_qfsolve >= 4, printp("invgen = ",lift(invgen)));
692
693    for( i = 1, length(Ker),
694      G2 = Qfbsqrtgauss(listgen[i],factdetG);
695      struct[i] <<= 1;
696      listgen[i] = G2;
697      invgen[,i] = Qflisteinvariants(G2,factD)[2][,1]*Mod(1,2)
698    );
699
700if( DEBUGLEVEL_qfsolve >= 3, print("listgen = ",listgen));
701if( DEBUGLEVEL_qfsolve >= 3, printp("invgen = ",lift(invgen)));
702if( DEBUGLEVEL_qfsolve >= 3, print("struct = ",struct));
703
704    im = lift(matinverseimage(invgen,matimage(invgen)))
705  );
706 
707  listgen2 = vector(rang);
708  for( i = 1, rang,
709    listgen2[i] = E;
710    for( j = 1, m,
711      if( im[j,i],
712        listgen2[i] = qfbcompraw(listgen2[i],listgen[j])));
713    if( norml2(im[,i]) > 1,
714      red = QfbReduce(aux=mymat(listgen2[i]));
715      aux = red~*aux*red;
716      listgen2[i] = Qfb(aux[1,1],2*aux[1,2],aux[2,2]))
717  );
718  listgen = listgen2;
719\\ listgen = vector(rang,i,listgen[m-rang+i]);
720  struct = vector(rang,i,struct[m-rang+i]);
721
722if( DEBUGLEVEL_qfsolve >= 2, print("listgen = ",listgen));
723if( DEBUGLEVEL_qfsolve >= 2, print("struct = ",struct));
724
725return([struct,listgen]);
726}
727{Qfsolve(G,factD) =
728\\ Resolution de la forme quadratique X^tGX=0 de dimension n >= 3.
729\\ On suppose que G est entiere et primitive.
730\\ La solution peut etre un vectorv ou une matrice.
731\\ S'il n'y a pas de solution, alors renvoit un p tel
732\\ qu'il n'existe pas de solution locale en p.
733\\
734\\ Si on connait la factorisation de -abs(2*matdet(G)),
735\\ on peut la passer par le parametre factD pour gagner du temps.
736
737local(n,M,signG,d,Min,U,codim,aux,G1,detG1,M1,subspace1,G2,subspace2,M2,solG2,Winvariants,dQ,factd,U2,clgp2,V,detG2,dimseti,solG1,sol,Q);
738
739if( DEBUGLEVEL_qfsolve >= 1, print("entree dans Qfsolve"));
740\\
741\\ 1ere reduction des coefficients de G
742\\
743
744  n = length(G);
745  M = IndefiniteLLL(G);
746  if( type(M) == "t_COL",
747if( DEBUGLEVEL_qfsolve >= 1, print("solution ",M));
748    return(M));
749  G = M[1]; M = M[2];
750
751\\ Solubilite reelle
752  signG = qfsign(G);
753  if( signG[1] == 0 || signG[2] == 0,
754if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution reelle"));
755    return(-1));
756  if( signG[1] < signG[2], G = -G; signG = signG*[0,1;1,0]);
757
758\\ Factorisation du determinant
759  d = matdet(G);
760  if( !factD,
761if( DEBUGLEVEL_qfsolve >= 1, print("factorisation du determinant"));
762    factD = factor(-abs(2*d));
763if( DEBUGLEVEL_qfsolve >= 1, print(factD))
764  );
765  factD[1,2] = 0;
766  factD[2,2] --;
767
768\\
769\\ Minimisation et solubilite locale.
770\\
771
772if( DEBUGLEVEL_qfsolve >= 1, print("minimisation du determinant"));
773  Min = Qfminim(G,factD);
774  if( type(Min) == "t_INT",
775if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution locale en ",Min));
776    return(Min));
777
778  M = M*Min[2];
779  G = Min[1];
780\\  Min[3] contient la factorisation de abs(matdet(G));
781
782if( DEBUGLEVEL_qfsolve >= 4, print("G minime = ",G));
783if( DEBUGLEVEL_qfsolve >= 4, print("d=",d));
784
785\\ Maintenant, on sait qu'il y a des solutions locales
786\\ (sauf peut-etre en 2 si n==4),
787\\ si n==3, det(G) = +-1
788\\ si n==4, ou si n est impair, det(G) est sans facteur carre.
789\\ si n>=6, det(G) toutes les valuations sont <=2.
790
791\\ Reduction de G et recherche de solution triviales
792\\ par exemple quand det G=+-1, il y en a toujours.
793
794if( DEBUGLEVEL_qfsolve >= 1, print("reduction"));
795  U = IndefiniteLLL(G);
796  if(type(U) == "t_COL",
797if( DEBUGLEVEL_qfsolve >= 1, print("solution ",M*U));
798    return(M*U));
799  G = U[1]; M = M*U[2];
800
801\\
802\\ Quand n >= 6 est pair, il faut passer en dimension n+1
803\\ pour supprimer tous les carres de det(G).
804\\
805
806  if( n >= 6 && n%2 == 0 && matsize(Min[3])[1] != 0,
807if( DEBUGLEVEL_qfsolve >= 1, print("On passe en dimension ",n+1));
808    codim = 1; n++;
809\\ On calcule le plus grand diviseur carre de d.   
810    aux = prod( i = 1, matsize(Min[3])[1], if( Min[3][i,1] == 2, Min[3][i,1], 1));
811\\ On choisit le signe de aux pour que la signature de G1
812\\ soit la plus equilibree possible.
813    if( signG[1] > signG[2],
814      signG[2] ++; aux = -aux
815    , signG[1] ++
816    );
817    G1 = matdiagonalblock([G,Mat(aux)]);
818    detG1 = 2*matdet(G1);
819    for( i = 2, length(factD[,1]),
820      factD[i,2] = valuation(detG1,factD[i,1]));
821    factD[2,2]--;
822    Min = Qfminim(G1,factD);
823    G1 = Min[1];
824    M1 = Min[2];
825    subspace1 = matrix(n,n-1,i,j, i == j)
826  , codim = 0;
827    G1 = G;
828    subspace1 = M1 = matid(n)
829  );
830\\ Maintenant d est sans facteur carre.
831
832\\
833\\ Si d n'est pas +-1, il faut passer en dimension n+2
834\\
835
836  if( matsize(Min[3])[1] == 0, \\ if( abs(d) == 1,
837if( DEBUGLEVEL_qfsolve >= 2, print(" detG2 = 1"));
838     G2 = G1;
839     subspace2 = M2 = matid(n);
840     solG2 = LLLgoon(G2,1)
841  ,
842if( DEBUGLEVEL_qfsolve >= 1, print("On passe en dimension ",n+2));
843    codim += 2;
844    subspace2 = matrix( n+2, n, i, j, i == j);
845    d = prod( i = 1, matsize(Min[3])[1],Min[3][i,1]);    \\ d = abs(matdet(G1));
846    if( signG[2]%2 == 1, d = -d);                        \\ d = matdet(G1);
847    if( Min[3][1,1] == 2, factD = [-1], factD = [-1,2]); \\ si d est pair ...
848    factD = concat(factD, Min[3][,1]~);
849if( DEBUGLEVEL_qfsolve >= 4, print("factD=",factD));
850
851\\ Solubilite locale en 2 (c'est le seul cas qui restait a traiter !!)
852    if( n == 4 && d%8 == 1,
853      if( QfWittinvariant(G,2) == 1,
854if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution locale en 2"));
855        return(2)));
856 
857\\
858\\ Construction d'une forme Q de dim 2, a partir de ces invariants.
859\\
860    Winvariants = vectorv(length(factD));
861
862\\ choix de la signature de Q.
863\\ invariant reel et signe du discriminant.
864    dQ = abs(d);
865    if( signG[1] ==signG[2], dQ = dQ; Winvariants[1] = 0); \\ signQ = [1,1];
866    if( signG[1] > signG[2], dQ =-dQ; Winvariants[1] = 0); \\ signQ = [2,0];
867    if( n == 4 && dQ%4 != 1, dQ *= 4);
868    if( n >= 5, dQ *= 8);
869   
870\\ invariants p-adiques
871\\ pour p = 2, on ne choisit pas.
872    if( n == 4,
873if( DEBUGLEVEL_qfsolve >= 1, print("calcul des invariants locaux de G1"));
874      aux = Qflisteinvariants(-G1,factD)[2][,1];
875      for( i = 3, length(factD), Winvariants[i] = aux[i])
876    ,
877      aux = (-1)^((n-3)/2)*dQ/d; \\ ici aux = 8 ou -8
878      for( i = 3, length(factD), Winvariants[i] = hilbert(aux,factD[i],factD[i]) > 0)
879    );
880    Winvariants[2] = sum( i = 1, length(factD), Winvariants[i])%2;
881
882if( DEBUGLEVEL_qfsolve >= 1,
883  print("Recherche d'un forme binaire de discriminant = ",dQ);
884  print("et d'invariants de Witt = ",Winvariants));
885
886\\ On construit le 2-groupe de classes de disc dQ,
887\\ jusqu'a obtenir une combinaison des generateurs qui donne les bons invariants.
888\\ En dim 4, il faut chercher parmi les formes du type q, 2*q
889\\ car Q peut etre imprimitive.
890
891    factd = matrix(length(factD)-1,2);
892    for( i = 1, length(factD)-1,
893      factd[i,1] = factD[i+1];
894      factd[i,2] = valuation(dQ,factd[i,1]));
895    factd[1,2]++;
896    U2 = matrix(length(factD), n == 4, i,j, hilbert(2,dQ,factD[i])<0);
897    clgp2 = class2(dQ,factd,Winvariants,U2);
898if( DEBUGLEVEL_qfsolve >= 4, print("clgp2=",clgp2));
899
900    clgp2 = clgp2[2];
901    U = Qflisteinvariants(clgp2,factD)[2];
902    if( n == 4, U = concat(U,U2));
903if( DEBUGLEVEL_qfsolve >= 4, printp("U=",U));
904
905    V = lift(matinverseimage(U*Mod(1,2),Winvariants*Mod(1,2)));
906    if( !length(V), next);
907if( DEBUGLEVEL_qfsolve >= 4, print("V=",V));
908
909    if( dQ%2 == 1, Q = qfbprimeform(4*dQ,1), Q = qfbprimeform(dQ,1));
910    for( i = 1, length(clgp2),
911      if( V[i], Q = qfbcompraw(Q,clgp2[i])));
912    Q = mymat(Q);
913    if( norml2(V) > 1, aux = QfbReduce(Q); Q = aux~*Q*aux);
914    if( n == 4 && V[length(V)], Q*=  2);
915if( DEBUGLEVEL_qfsolve >= 2, print("Q=",Q));
916if( DEBUGLEVEL_qfsolve >= 3, print("invariants de Witt de Q=",Qflisteinvariants(Q,factD)));
917
918\\
919\\ Construction d'une forme de dim=n+2 potentiellement unimodulaire
920\\
921
922    G2 = matdiagonalblock([G1,-Q]);
923if( DEBUGLEVEL_qfsolve >= 4, print("G2=",G2));
924
925if( DEBUGLEVEL_qfsolve >= 2, print("minimisation de la forme quadratique de dimension ",length(G2)));
926\\ Minimisation de G2
927    detG2 = matdet(G2);
928    factd = matrix(length(factD)-1,2);
929    for( i = 1, length(factD)-1,
930      factd[i,2] = valuation(detG2, factd[i,1] = factD[i+1]));
931if( DEBUGLEVEL_qfsolve >= 3, print("det(G2) = ",factd));
932    Min = Qfminim(G2,factd);
933    M2 = Min[2]; G2 = Min[1];
934if( abs(matdet(G2)) > 2, print("******* ERREUR dans Qfsolve: det(G2) <> +-1 *******",matdet(G2));return(0));
935if( DEBUGLEVEL_qfsolve >= 4, print("G2=",G2));
936
937\\ Maintenant det(G2) = +-1
938
939\\ On construit un seti pour G2 (Sous-Espace Totalement Isotrope)
940if( DEBUGLEVEL_qfsolve >= 2, print("recherche d'un espace de solutions pour G2"));
941    solG2 = LLLgoon(G2,1);
942    if( matrix(codim+1,codim+1,i,j,solG2[1][i,j]) != 0,
943if( DEBUGLEVEL_qfsolve >= 2, print(" pas assez de solutions dans G2"));return(0))
944  );
945
946\\ G2 doit avoir un espace de solutions de dimension > codim
947  dimseti = 0;
948  while( matrix(dimseti+1,dimseti+1,i,j,solG2[1][i,j]) == 0, dimseti ++);
949  if( dimseti <= codim,
950print("dimseti = ",dimseti," <= codim = ",codim);
951print("************* ERREUR : pas assez de solutions pour G2"); return(0));   
952  solG2 = matrix(length(G2),dimseti,i,j,solG2[2][i,j]);
953if( DEBUGLEVEL_qfsolve >= 3, print("solG2=",solG2));
954
955\\ La solution de G1 se trouve a la fois dans solG2 et dans subspace2
956if( DEBUGLEVEL_qfsolve >= 1, print("Reconstruction de la solution de G1"));
957  solG1 = matintersect(subspace2,M2*solG2);
958  solG1 = subspace2~*solG1;
959if( DEBUGLEVEL_qfsolve >= 3, print("solG1 = ",solG1));
960
961\\ La solution de G se trouve a la fois dans solG et dans subspace1
962if( DEBUGLEVEL_qfsolve >= 1, print("Reconstruction de la solution de G"));
963  sol = matintersect(subspace1,M1*solG1);
964  sol = subspace1~*sol;
965  sol = M*sol;
966  sol /= content(sol);
967  if( length(sol) == 1, sol = sol[,1]);
968if( DEBUGLEVEL_qfsolve >= 3, print("sol = ",sol));
969if( DEBUGLEVEL_qfsolve >= 1, print("fin de Qfsolve"));
970  return(sol);
971}
972{matdiagonalblock(v) =
973local(lv,lt,M);
974  lv = length(v);
975  lt = sum( i = 1, lv, length(v[i]));
976  M = matrix(lt,lt);
977  lt = 0;
978  for( i = 1, lv,
979    for( j = 1, length(v[i]),
980      for( k = 1, length(v[i]),
981        M[lt+j,lt+k] = v[i][j,k]));
982    lt += length(v[i])
983  );
984return(M);
985}
986
987
988
989
990
991