"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/qfsolve.c" (14 Jan 2021, 32913 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "qfsolve.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

1 /* Copyright (C) 2000-2004  The PARI group.
2
3 This file is part of the PARI/GP package.
4
5 PARI/GP is free software; you can redistribute it and/or modify it under the
7 Foundation; either version 2 of the License, or (at your option) any later
8 version. It is distributed in the hope that it will be useful, but WITHOUT
9 ANY WARRANTY WHATSOEVER.
10
11 Check the License for details. You should have received a copy of it, along
12 with the package; see the file 'COPYING'. If not, write to the Free Software
13 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14
15 /* Copyright (C) 2014 Denis Simon
16  * Adapted from qfsolve.gp v. 09/01/2014
17  *   http://www.math.unicaen.fr/~simon/qfsolve.gp
18  *
19  * Author: Denis SIMON <simon@math.unicaen.fr> */
20
21 #include "pari.h"
22 #include "paripriv.h"
23
24 /* LINEAR ALGEBRA */
25 /* complete by 0s, assume l-1 <= n */
26 static GEN
27 vecextend(GEN v, long n)
28 {
29   long i, l = lg(v);
30   GEN w = cgetg(n+1, t_COL);
31   for (i = 1; i < l; i++) gel(w,i) = gel(v,i);
32   for (     ; i <=n; i++) gel(w,i) = gen_0;
33   return w;
34 }
35
36 /* Gives a unimodular matrix with the last column(s) equal to Mv.
37  * Mv can be a column vector or a rectangular matrix.
38  * redflag = 0 or 1. If redflag = 1, LLL-reduce the n-#v first columns. */
39 static GEN
40 completebasis(GEN Mv, long redflag)
41 {
42   GEN U;
43   long m, n;
44
45   if (typ(Mv) == t_COL) Mv = mkmat(Mv);
46   n = lg(Mv)-1;
47   m = nbrows(Mv); /* m x n */
48   if (m == n) return Mv;
49   (void)ZM_hnfall_i(shallowtrans(Mv), &U, 0);
50   U = ZM_inv(shallowtrans(U),NULL);
51   if (m==1 || !redflag) return U;
52   /* LLL-reduce the m-n first columns */
53   return shallowconcat(ZM_lll(vecslice(U,1,m-n), 0.99, LLL_INPLACE),
54                        vecslice(U, m-n+1,m));
55 }
56
57 /* Compute the kernel of M mod p.
58  * returns [d,U], where
59  * d = dim (ker M mod p)
60  * U in GLn(Z), and its first d columns span the kernel. */
61 static GEN
62 kermodp(GEN M, GEN p, long *d)
63 {
64   long j, l;
65   GEN K, B, U;
66
67   K = FpM_center(FpM_ker(M, p), p, shifti(p,-1));
68   B = completebasis(K,0);
69   l = lg(M); U = cgetg(l, t_MAT);
70   for (j =  1; j < l; j++) gel(U,j) = gel(B,l-j);
71   *d = lg(K)-1; return U;
72 }
73
74 /* INVARIANTS COMPUTATIONS */
75
76 static GEN
77 principal_minor(GEN G, long  i) { return matslice(G,1,i,1,i); }
78 static GEN
79 det_minors(GEN G)
80 {
81   long i, l = lg(G);
82   GEN v = cgetg(l+1, t_VEC);
83   gel(v,1) = gen_1;
84   for (i = 2; i <= l; i++) gel(v,i) = ZM_det(principal_minor(G,i-1));
85   return v;
86 }
87
88 /* Given a symmetric matrix G over Z, compute the Witt invariant
89  *  of G at the prime p (at real place if p = NULL)
90  * Assume that none of the determinant G[1..i,1..i] is 0. */
91 static long
92 qflocalinvariant(GEN G, GEN p)
93 {
94   long i, j, c, l = lg(G);
95   GEN diag, v = det_minors(G);
96   /* Diagonalize G first. */
97   diag = cgetg(l, t_VEC);
98   for (i = 1; i < l; i++) gel(diag,i) = mulii(gel(v,i+1), gel(v,i));
99
100   /* Then compute the product of the Hilbert symbols */
101   /* (diag[i],diag[j])_p for i < j */
102   c = 1;
103   for (i = 1; i < l-1; i++)
104     for (j = i+1; j < l; j++)
105       if (hilbertii(gel(diag,i), gel(diag,j), p) < 0) c = -c;
106   return c;
107 }
108
109 static GEN
110 hilberts(GEN a, GEN b, GEN P, long lP)
111 {
112   GEN v = cgetg(lP, t_VECSMALL);
113   long i;
114   for (i = 1; i < lP; i++) v[i] = hilbertii(a, b, gel(P,i)) < 0;
115   return v;
116 }
117
118 /* G symmetrix matrix or qfb or list of quadratic forms with same discriminant.
119  * P must be equal to factor(-abs(2*matdet(G)))[,1]. */
120 static GEN
121 qflocalinvariants(GEN G, GEN P)
122 {
123   GEN sol;
124   long i, j, l, lP = lg(P);
125
126   /* convert G into a vector of symmetric matrices */
127   G = (typ(G) == t_VEC)? shallowcopy(G): mkvec(G);
128   l = lg(G);
129   for (j = 1; j < l; j++)
130   {
131     GEN g = gel(G,j);
132     if (is_qfb_t(typ(g))) gel(G,j) = gtomat(g);
133   }
134   sol = cgetg(l, t_MAT);
135   if (lg(gel(G,1)) == 3)
136   { /* in dimension 2, each invariant is a single Hilbert symbol. */
137     GEN d = negi(ZM_det(gel(G,1)));
138     for (j = 1; j < l; j++)
139     {
140       GEN a = gcoeff(gel(G,j),1,1);
141       gel(sol,j) = hilberts(a, d, P, lP);
142     }
143   }
144   else /* in dimension n > 2, we compute a product of n Hilbert symbols. */
145     for (j = 1; j <l; j++)
146     {
147       GEN g = gel(G,j), v = det_minors(g), w = cgetg(lP, t_VECSMALL);
148       long n = lg(v);
149       gel(sol,j) = w;
150       for (i = 1; i < lP; i++)
151       {
152         GEN p = gel(P,i);
153         long k = n-2, h = hilbertii(gel(v,k), gel(v,k+1),p);
154         for (k--; k >= 1; k--)
155           if (hilbertii(negi(gel(v,k)), gel(v,k+1),p) < 0) h = -h;
156         w[i] = h < 0;
157       }
158     }
159   return sol;
160 }
161
162 /* QUADRATIC FORM REDUCTION */
163 static GEN
164 qfb(GEN D, GEN a, GEN b, GEN c)
165 {
166   if (signe(D) < 0) return mkqfi(a,b,c);
167   retmkqfr(a,b,c,real_0(DEFAULTPREC));
168 }
169
170 /* Gauss reduction of the binary quadratic form
171  * Q = a*X^2+2*b*X*Y+c*Y^2 of discriminant D (divisible by 4)
172  * returns the reduced form */
173 static GEN
174 qfbreduce(GEN D, GEN Q)
175 {
176   GEN a = gel(Q,1), b = shifti(gel(Q,2),-1), c = gel(Q,3);
177   while (signe(a))
178   {
179     GEN r, q, nexta, nextc;
180     q = dvmdii(b,a, &r); /* FIXME: export as dvmdiiround ? */
181     if (signe(r) > 0 && abscmpii(shifti(r,1), a) > 0)
182     {
183       r = subii(r, absi_shallow(a));
185     }
186     nextc = a; nexta = subii(c, mulii(q, addii(r,b)));
187     if (abscmpii(nexta, a) >= 0) break;
188     c = nextc; b = negi(r); a = nexta;
189   }
190   return qfb(D,a,shifti(b,1),c);
191 }
192
193 /* private version of qfgaussred:
194  * - early abort if k-th principal minor is singular, return stoi(k)
195  * - else return a matrix whose upper triangular part is qfgaussred(a) */
196 static GEN
197 partialgaussred(GEN a)
198 {
199   long n = lg(a)-1, k;
200   a = RgM_shallowcopy(a);
201   for(k = 1; k < n; k++)
202   {
203     GEN ak, p = gcoeff(a,k,k);
204     long i, j;
205     if (isintzero(p)) return stoi(k);
206     ak = row(a, k);
207     for (i=k+1; i<=n; i++) gcoeff(a,k,i) = gdiv(gcoeff(a,k,i), p);
208     for (i=k+1; i<=n; i++)
209     {
210       GEN c = gel(ak,i);
211       if (gequal0(c)) continue;
212       for (j=i; j<=n; j++)
213         gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
214     }
215   }
216   if (isintzero(gcoeff(a,n,n))) return stoi(n);
217   return a;
218 }
219
220 /* LLL-reduce a positive definite qf QD bounding the indefinite G, dim G > 1.
221  * Then finishes by looking for trivial solution */
222 static GEN qftriv(GEN G, GEN z, long base);
223 static GEN
224 qflllgram_indef(GEN G, long base, int *fail)
225 {
226   GEN M, R, g, DM, S, dR;
227   long i, j, n = lg(G)-1;
228
229   *fail = 0;
230   R = partialgaussred(G);
231   if (typ(R) == t_INT) return qftriv(G, R, base);
232   R = Q_remove_denom(R, &dR); /* avoid rational arithmetic */
233   M = zeromatcopy(n,n);
234   DM = zeromatcopy(n,n);
235   for (i = 1; i <= n; i++)
236   {
237     GEN d = absi_shallow(gcoeff(R,i,i));
238     if (dR) {
239       gcoeff(M,i,i) = dR;
240       gcoeff(DM,i,i) = mulii(d,dR);
241     } else {
242       gcoeff(M,i,i) = gen_1;
243       gcoeff(DM,i,i) = d;
244     }
245     for (j = i+1; j <= n; j++)
246     {
247       gcoeff(M,i,j) = gcoeff(R,i,j);
248       gcoeff(DM,i,j) = mulii(d, gcoeff(R,i,j));
249     }
250   }
251   /* G = M~*D*M, D diagonal, DM=|D|*M, g =  M~*|D|*M */
252   g = ZM_transmultosym(M,DM);
253   S = lllgramint(Q_primpart(g));
254   R = qftriv(qf_apply_ZM(G,S), NULL, base);
255   switch(typ(R))
256   {
257     case t_COL: return ZM_ZC_mul(S,R);
258     case t_MAT: *fail = 1; return mkvec2(R, S);
259     default:
260       gel(R,2) = ZM_mul(S, gel(R,2));
261       return R;
262   }
263 }
264
265 /* G symmetric, i < j, let E = E_{i,j}(a), G <- E~*G*E,  U <- U*E.
266  * Everybody integral */
267 static void
268 qf_apply_transvect_Z(GEN G, GEN U, long i, long j, GEN a)
269 {
270   long k, n = lg(G)-1;
271   gel(G, j) =  ZC_lincomb(gen_1, a, gel(G,j), gel(G,i));
272   for (k = 1; k < n; k++) gcoeff(G, j, k) = gcoeff(G, k, j);
275   gel(U, j) =  ZC_lincomb(gen_1, a, gel(U,j), gel(U,i));
276 }
277
278 /* LLL reduction of the quadratic form G (Gram matrix)
279  * where we go on, even if an isotropic vector is found. */
280 static GEN
281 qflllgram_indefgoon(GEN G)
282 {
283   GEN red, U, A, U1,U2,U3,U5,U6, V, B, G2,G3,G4,G5, G6, a, g;
284   long i, j, n = lg(G)-1;
285   int fail;
286
287   red = qflllgram_indef(G,1, &fail);
288   if (fail) return red; /*no isotropic vector found: nothing to do*/
289   /* otherwise a solution is found: */
290   U1 = gel(red,2);
291   G2 = gel(red,1); /* G2[1,1] = 0 */
292   U2 = gel(ZV_extgcd(row(G2,1)), 2);
293   G3 = qf_apply_ZM(G2,U2);
294   U = ZM_mul(U1,U2); /* qf_apply(G,U) = G3 */
295   /* G3[1,] = [0,...,0,g], g^2 | det G */
296   g = gcoeff(G3,1,n);
297   a = diviiround(negi(gcoeff(G3,n,n)), shifti(g,1));
298   if (signe(a)) qf_apply_transvect_Z(G3,U,1,n,a);
299   /* G3[n,n] reduced mod 2g */
300   if (n == 2) return mkvec2(G3,U);
301   V = rowpermute(vecslice(G3, 2,n-1), mkvecsmall2(1,n));
302   A = mkmat2(mkcol2(gcoeff(G3,1,1),gcoeff(G3,1,n)),
303              mkcol2(gcoeff(G3,1,n),gcoeff(G3,2,2)));
304   B = ground(RgM_neg(RgM_mul(RgM_inv(A), V)));
305   U3 = matid(n);
306   for (j = 2; j < n; j++)
307   {
308     gcoeff(U3,1,j) = gcoeff(B,1,j-1);
309     gcoeff(U3,n,j) = gcoeff(B,2,j-1);
310   }
311   G4 = qf_apply_ZM(G3,U3); /* the last column of G4 is reduced */
312   U = ZM_mul(U,U3);
313   if (n == 3) return mkvec2(G4,U);
314
315   red = qflllgram_indefgoon(matslice(G4,2,n-1,2,n-1));
316   if (typ(red) == t_MAT) return mkvec2(G4,U);
317   /* Let U5:=matconcat(diagonal[1,red[2],1])
318    * return [qf_apply_ZM(G5, U5), U*U5] */
319   G5 = gel(red,1);
320   U5 = gel(red,2);
321   G6 = cgetg(n+1,t_MAT);
322   gel(G6,1) = gel(G4,1);
323   gel(G6,n) = gel(G4,n);
324   for (j=2; j<n; j++)
325   {
326     gel(G6,j) = cgetg(n+1,t_COL);
327     gcoeff(G6,1,j) = gcoeff(G4,j,1);
328     gcoeff(G6,n,j) = gcoeff(G4,j,n);
329     for (i=2; i<n; i++) gcoeff(G6,i,j) = gcoeff(G5,i-1,j-1);
330   }
331   U6 = mkvec3(mkmat(gel(U,1)), ZM_mul(vecslice(U,2,n-1),U5), mkmat(gel(U,n)));
332   return mkvec2(G6, shallowconcat1(U6));
333 }
334
335 /* qf_apply_ZM(G,H),  where H = matrix of \tau_{i,j}, i != j */
336 static GEN
337 qf_apply_tau(GEN G, long i, long j)
338 {
339   long l = lg(G), k;
340   G = RgM_shallowcopy(G);
341   swap(gel(G,i), gel(G,j));
342   for (k = 1; k < l; k++) swap(gcoeff(G,i,k), gcoeff(G,j,k));
343   return G;
344 }
345
346 /* LLL reduction of the quadratic form G (Gram matrix)
347  * in dim 3 only, with detG = -1 and sign(G) = [2,1]; */
348 static GEN
349 qflllgram_indefgoon2(GEN G)
350 {
351   GEN red, G2, a, b, c, d, e, f, u, v, r, r3, U2, G3;
352   int fail;
353
354   red = qflllgram_indef(G,1,&fail); /* always find an isotropic vector. */
355   G2 = qf_apply_tau(gel(red,1),1,3); /* G2[3,3] = 0 */
356   r = row(gel(red,2), 3);
357   swap(gel(r,1), gel(r,3)); /* apply tau_{1,3} */
358   a = gcoeff(G2,3,1);
359   b = gcoeff(G2,3,2);
360   d = bezout(a,b, &u,&v);
361   if (!equali1(d))
362   {
363     a = diviiexact(a,d);
364     b = diviiexact(b,d);
365   }
366   /* for U2 = [-u,-b,0;-v,a,0;0,0,1]
367    * G3 = qf_apply_ZM(G2,U2) has known last row (-d, 0, 0),
368    * so apply to principal_minor(G3,2), instead */
369   U2 = mkmat2(mkcol2(negi(u),negi(v)), mkcol2(negi(b),a));
370   G3 = qf_apply_ZM(principal_minor(G2,2),U2);
371   r3 = gel(r,3);
372   r = ZV_ZM_mul(mkvec2(gel(r,1),gel(r,2)),U2);
373
374   a = gcoeff(G3,1,1);
375   b = gcoeff(G3,1,2);
376   c = negi(d); /* G3[1,3] */
377   d = gcoeff(G3,2,2);
378   if (mpodd(a))
379   {
382     e = diviiround(negi(e),c);
383     f = diviiround(negi(a), shifti(c,1));
385   }
386   else
387   {
388     e = diviiround(negi(b),c);
389     f = diviiround(negi(shifti(a,-1)), c);
390     a = addmulii(gel(r,1), f, r3);
391   }
392   b = addmulii(gel(r,2), e, r3);
393   return mkvec3(a,b, r3);
394 }
395
396 /* QUADRATIC FORM MINIMIZATION */
397 /* G symmetric, return ZM_Z_divexact(G,d) */
398 static GEN
399 ZsymM_Z_divexact(GEN G, GEN d)
400 {
401   long i,j,l = lg(G);
402   GEN H = cgetg(l, t_MAT);
403   for(j=1; j<l; j++)
404   {
405     GEN c = cgetg(l, t_COL), b = gel(G,j);
406     for(i=1; i<j; i++) gcoeff(H,j,i) = gel(c,i) = diviiexact(gel(b,i),d);
407     gel(c,j) = diviiexact(gel(b,j),d);
408     gel(H,j) = c;
409   }
410   return H;
411 }
412
413 /* write symmetric G as [A,B;B~,C], A dxd, C (n-d)x(n-d) */
414 static void
415 blocks4(GEN G, long d, long n, GEN *A, GEN *B, GEN *C)
416 {
417   GEN G2 = vecslice(G,d+1,n);
418   *A = principal_minor(G, d);
419   *B = rowslice(G2, 1, d);
420   *C = rowslice(G2, d+1, n);
421 }
422 /* Minimization of the quadratic form G, deg G != 0, dim n >= 2
423  * G symmetric integral
424  * Returns [G',U,factd] with U in GLn(Q) such that G'=U~*G*U*constant
425  * is integral and has minimal determinant.
426  * In dimension 3 or 4, may return a prime p if the reduction at p is
427  * impossible because of local nonsolvability.
428  * P,E = factor(+/- det(G)), "prime" -1 is ignored. Destroy E. */
429 static GEN qfsolvemodp(GEN G, GEN p);
430 static GEN
431 qfminimize(GEN G, GEN P, GEN E)
432 {
433   pari_sp av;
434   GEN d, U, Ker, sol, aux, faE, faP;
435   long n = lg(G)-1, lP = lg(P), i, dimKer, m;
436   faP = vectrunc_init(lP);
437   faE = vecsmalltrunc_init(lP);
438   av = avma;
439   U = NULL;
440   for (i = 1; i < lP; i++)
441   {
442     GEN p = gel(P,i);
443     long vp = E[i];
444     if (!vp || !p) continue;
445
446     if (DEBUGLEVEL >= 4) err_printf("    p^v = %Ps^%ld\n", p,vp);
447     /* The case vp = 1 can be minimized only if n is odd. */
448     if (vp == 1 && n%2 == 0) {
449       vectrunc_append(faP, p);
450       vecsmalltrunc_append(faE, 1);
451       av = avma; continue;
452     }
453     Ker = kermodp(G,p, &dimKer); /* dimKer <= vp */
454     if (DEBUGLEVEL >= 4) err_printf("    dimKer = %ld\n",dimKer);
455     if (dimKer == n)
456     { /* trivial case: dimKer = n */
457       if (DEBUGLEVEL >= 4) err_printf("     case 0: dimKer = n\n");
458       G = ZsymM_Z_divexact(G, p);
459       E[i] -= n;
460       i--; continue; /* same p */
461     }
462     G = qf_apply_ZM(G, Ker);
463     U = U? RgM_mul(U,Ker): Ker;
464
465     /* 1st case: dimKer < vp */
466     /* then the kernel mod p contains a kernel mod p^2 */
467     if (dimKer < vp)
468     {
469       if (DEBUGLEVEL >= 4) err_printf("    case 1: dimker < vp\n");
470       if (dimKer == 1)
471       {
472         long j;
473         gel(G,1) = ZC_Z_divexact(gel(G,1), p);
474         for (j = 1; j<=n; j++) gcoeff(G,1,j) = diviiexact(gcoeff(G,1,j), p);
475         gel(U,1) = RgC_Rg_div(gel(U,1), p);
476         E[i] -= 2;
477       }
478       else
479       {
480         GEN A,B,C, K2 = ZsymM_Z_divexact(principal_minor(G,dimKer),p);
481         long j, dimKer2;
482         K2 = kermodp(K2, p, &dimKer2);
483         for (j = dimKer2+1; j <= dimKer; j++) gel(K2,j) = ZC_Z_mul(gel(K2,j),p);
484         /* Write G = [A,B;B~,C] and apply [K2,0;0,p*Id]/p by blocks */
485         blocks4(G, dimKer,n, &A,&B,&C);
486         A = ZsymM_Z_divexact(qf_apply_ZM(A,K2), sqri(p));
487         B = ZM_Z_divexact(ZM_transmul(B,K2), p);
488         G = shallowmatconcat(mkmat2(mkcol2(A,B),
489                                     mkcol2(shallowtrans(B), C)));
490         /* U *= [K2,0;0,Id] */
491         U = shallowconcat(RgM_Rg_div(RgM_mul(vecslice(U,1,dimKer),K2), p),
492                           vecslice(U,dimKer+1,n));
493         E[i] -= 2*dimKer2;
494       }
495       gerepileall(av, 2, &G, &U);
496       i--; continue; /* same p */
497     }
498
499    /* vp = dimKer
500     * 2nd case: kernel has dim >= 2 and contains an element of norm 0 mod p^2
501     * search for an element of norm p^2... in the kernel */
502     sol = NULL;
503     if (dimKer > 2) {
504       if (DEBUGLEVEL >= 4) err_printf("    case 2.1\n");
505       dimKer = 3;
506       sol = qfsolvemodp(ZsymM_Z_divexact(principal_minor(G,3),p),  p);
507       sol = FpC_red(sol, p);
508     }
509     else if (dimKer == 2)
510     {
511       GEN a = modii(diviiexact(gcoeff(G,1,1),p), p);
512       GEN b = modii(diviiexact(gcoeff(G,1,2),p), p);
513       GEN c = diviiexact(gcoeff(G,2,2),p);
514       GEN di= modii(subii(sqri(b), mulii(a,c)), p);
515       if (kronecker(di,p) >= 0)
516       {
517         if (DEBUGLEVEL >= 4) err_printf("    case 2.2\n");
518         sol = signe(a)? mkcol2(Fp_sub(Fp_sqrt(di,p), b, p), a): vec_ei(2,1);
519       }
520     }
521     if (sol)
522     {
523       long j;
524       sol = FpC_center(sol, p, shifti(p,-1));
525       sol = Q_primpart(sol);
526       if (DEBUGLEVEL >= 4) err_printf("    sol = %Ps\n", sol);
527       Ker = completebasis(vecextend(sol,n), 1);
528       for(j=1; j<n; j++) gel(Ker,j) = ZC_Z_mul(gel(Ker,j), p);
529       G = ZsymM_Z_divexact(qf_apply_ZM(G, Ker), sqri(p));
530       U = RgM_Rg_div(RgM_mul(U,Ker), p);
531       E[i] -= 2;
532       i--; continue; /* same p */
533     }
534     /* Now 1 <= vp = dimKer <= 2 and kernel contains no vector with norm p^2 */
535     /* exchanging kernel and image makes minimization easier ? */
536     m = (n-3)/2;
537     d = ZM_det(G); if (odd(m)) d = negi(d);
538     if ((vp==1 && kronecker(gmod(gdiv(negi(d), gcoeff(G,1,1)),p), p) >= 0)
539      || (vp==2 && odd(n) && n >= 5)
540      || (vp==2 && !odd(n) && kronecker(modii(diviiexact(d,sqri(p)), p),p) < 0))
541     {
542       long j;
543       if (DEBUGLEVEL >= 4) err_printf("    case 3\n");
544       Ker = matid(n);
545       for (j = dimKer+1; j <= n; j++) gcoeff(Ker,j,j) = p;
546       G = ZsymM_Z_divexact(qf_apply_ZM(G, Ker), p);
547       U = RgM_mul(U,Ker);
548       E[i] -= 2*dimKer-n;
549       i--; continue; /* same p */
550     }
551
552     /* Minimization was not possible so far. */
553     /* If n == 3 or 4, this proves the local nonsolubility at p. */
554     if (n == 3 || n == 4)
555     {
556       if (DEBUGLEVEL >= 1) err_printf(" no local solution at %Ps\n",p);
557       return(p);
558     }
559     vectrunc_append(faP, p);
560     vecsmalltrunc_append(faE, vp);
561     av = avma;
562   }
563   if (!U) U = matid(n);
564   else
565   { /* apply LLL to avoid coefficient explosion */
566     aux = lllint(Q_primpart(U));
567     G = qf_apply_ZM(G,aux);
568     U = RgM_mul(U,aux);
569   }
570   return mkvec4(G, U, faP, faE);
571 }
572
573 /* CLASS GROUP COMPUTATIONS */
574
575 /* Compute the square root of the quadratic form q of discriminant D. Not
576  * fully implemented; it only works for detqfb squarefree except at 2, where
577  * the valuation is 2 or 3.
578  * mkmat2(P,zv_to_ZV(E)) = factor(2*abs(det q)) */
579 static GEN
580 qfbsqrt(GEN D, GEN q, GEN P, GEN E)
581 {
582   GEN a = gel(q,1), b = shifti(gel(q,2),-1), c = gel(q,3), mb = negi(b);
583   GEN m,n, aux,Q1,M, A,B,C;
584   GEN d = subii(mulii(a,c), sqri(b));
585   long i;
586
587   /* 1st step: solve m^2 = a (d), m*n = -b (d), n^2 = c (d) */
588   m = n = mkintmod(gen_0,gen_1);
589   E[1] -= 3;
590   for (i = 1; i < lg(P); i++)
591   {
592     GEN p = gel(P,i), N, M;
593     if (!E[i]) continue;
594     if (dvdii(a,p)) {
595       aux = Fp_sqrt(c, p);
596       N = aux;
597       M = Fp_div(mb, aux, p);
598     } else {
599       aux = Fp_sqrt(a, p);
600       M = aux;
601       N = Fp_div(mb, aux, p);
602     }
603     n = chinese(n, mkintmod(N,p));
604     m = chinese(m, mkintmod(M,p));
605   }
606   m = centerlift(m);
607   n = centerlift(n);
608   if (DEBUGLEVEL >=4) err_printf("    [m,n] = [%Ps, %Ps]\n",m,n);
609
610   /* 2nd step: build Q1, with det=-1 such that Q1(x,y,0) = G(x,y) */
611   A = diviiexact(subii(sqri(n),c), d);
612   B = diviiexact(addii(b, mulii(m,n)), d);
613   C = diviiexact(subii(sqri(m), a), d);
614   Q1 = mkmat3(mkcol3(A,B,n), mkcol3(B,C,m), mkcol3(n,m,d));
616
617   /* 3rd step: reduce Q1 to [0,0,-1;0,1,0;-1,0,0] */
618   M = qflllgram_indefgoon2(Q1);
619   if (signe(gel(M,1)) < 0) M = ZC_neg(M);
620   a = gel(M,1);
621   b = gel(M,2);
622   c = gel(M,3);
623   if (mpodd(a))
624     return qfb(D, a, shifti(b,1), shifti(c,1));
625   else
626     return qfb(D, c, shifti(negi(b),1), shifti(a,1));
627 }
628
629 /* \prod gen[i]^e[i] as a Qfb, e in {0,1}^n nonzero */
630 static GEN
631 qfb_factorback(GEN D, GEN gen, GEN e)
632 {
633   GEN q = NULL;
634   long j, l = lg(gen), n = 0;
635   for (j = 1; j < l; j++)
636     if (e[j]) { n++; q = q? qfbcompraw(q, gel(gen,j)): gel(gen,j); }
637   return (n <= 1)? q: qfbreduce(D, q);
638 }
639
640 /* unit form, assuming 4 | D */
641 static GEN
642 id(GEN D)
643 { return mkmat2(mkcol2(gen_1,gen_0),mkcol2(gen_0,shifti(negi(D),-2))); }
644
645 /* Shanks/Bosma-Stevenhagen algorithm to compute the 2-Sylow of the class
646  * group of discriminant D. Only works for D = fundamental discriminant.
647  * When D = 1(4), work with 4D.
648  * P2D,E2D = factor(abs(2*D))
649  * Pm2D = factor(-abs(2*D))[,1].
650  * Return a form having Witt invariants W at Pm2D */
651 static GEN
652 quadclass2(GEN D, GEN P2D, GEN E2D, GEN Pm2D, GEN W, int n_is_4)
653 {
654   GEN gen, Wgen, U2;
655   long i, n, r, m, vD;
656
657   if (mpodd(D))
658   {
659     D = shifti(D,2);
660     E2D = shallowcopy(E2D);
661     E2D[1] = 3;
662   }
663   if (zv_equal0(W)) return id(D);
664
665   n = lg(Pm2D)-1; /* >= 3 since W != 0 */
666   r = n-3;
667   m = (signe(D)>0)? r+1: r;
668   /* n=4: look among forms of type q or 2*q, since Q can be imprimitive */
669   U2 = n_is_4? mkmat(hilberts(gen_2, D, Pm2D, lg(Pm2D))): NULL;
670   if (U2 && zv_equal(gel(U2,1),W)) return gmul2n(id(D),1);
671
672   gen = cgetg(m+1, t_VEC);
673   for (i = 1; i <= m; i++) /* no need to look at Pm2D[1]=-1, nor Pm2D[2]=2 */
674   {
675     GEN p = gel(Pm2D,i+2), d;
676     long vp = Z_pvalrem(D,p, &d);
677     gel(gen,i) = qfb(D, powiu(p,vp), gen_0, negi(shifti(d,-2)));
678   }
679   vD = Z_lval(D,2);  /* = 2 or 3 */
680   if (vD == 2 && smodis(D,16) != 4)
681   {
682     GEN q2 = qfb(D, gen_2,gen_2, shifti(subsi(4,D),-3));
683     m++; r++; gen = vec_append(gen, q2);
684   }
685   if (vD == 3)
686   {
687     GEN q2 = qfb(D, gen_2,gen_0, negi(shifti(D,-3)));
688     m++; r++; gen = vec_append(gen, q2);
689   }
690   if (!r) return id(D);
691   Wgen = qflocalinvariants(gen,Pm2D);
692   for(;;)
693   {
694     GEN Wgen2, gen2, Ker, indexim = gel(Flm_indexrank(Wgen,2), 2);
695     long dKer;
696     if (lg(indexim)-1 >= r)
697     {
698       GEN W2 = Wgen, V;
699       if (lg(indexim) < lg(Wgen)) W2 = vecpermute(Wgen,indexim);
700       if (U2) W2 = shallowconcat(W2,U2);
701       V = Flm_Flc_invimage(W2, W,2);
702       if (V) {
703         GEN Q = qfb_factorback(D, vecpermute(gen,indexim), V);
704         Q = gtomat(Q);
705         if (U2 && V[lg(V)-1]) Q = gmul2n(Q,1);
706         return Q;
707       }
708     }
709     Ker = Flm_ker(Wgen,2); dKer = lg(Ker)-1;
710     gen2 = cgetg(m+1, t_VEC);
711     Wgen2 = cgetg(m+1, t_MAT);
712     for (i = 1; i <= dKer; i++)
713     {
714       GEN q = qfb_factorback(D, gen, gel(Ker,i));
715       q = qfbsqrt(D,q,P2D,E2D);
716       gel(gen2,i) = q;
717       gel(Wgen2,i) = gel(qflocalinvariants(q,Pm2D), 1);
718     }
719     for (; i <=m; i++)
720     {
721       long j = indexim[i-dKer];
722       gel(gen2,i) = gel(gen,j);
723       gel(Wgen2,i) = gel(Wgen,j);
724     }
725     gen = gen2; Wgen = Wgen2;
726   }
727 }
728
730 /* is x*y = -1 ? */
731 static int
732 both_pm1(GEN x, GEN y)
733 { return is_pm1(x) && is_pm1(y) && signe(x) == -signe(y); }
734
735 /* Try to solve G = 0 with small coefficients. This is proved to work if
736  * -  det(G) = 1, dim <= 6 and G is LLL reduced
737  * Returns G if no solution is found.
738  * Exit with a norm 0 vector if one such is found.
739  * If base == 1 and norm 0 is obtained, returns [H~*G*H,H] where
740  * the 1st column of H is a norm 0 vector */
741 static GEN
742 qftriv(GEN G, GEN R, long base)
743 {
744   long n = lg(G)-1, i;
745   GEN s, H;
746
747   /* case 1: A basis vector is isotropic */
748   for (i = 1; i <= n; i++)
749     if (!signe(gcoeff(G,i,i)))
750     {
751       if (!base) return col_ei(n,i);
752       H = matid(n); swap(gel(H,1), gel(H,i));
753       return mkvec2(qf_apply_tau(G,1,i),H);
754     }
755   /* case 2: G has a block +- [1,0;0,-1] on the diagonal */
756   for (i = 2; i <= n; i++)
757     if (!signe(gcoeff(G,i-1,i)) && both_pm1(gcoeff(G,i-1,i-1),gcoeff(G,i,i)))
758     {
759       s = col_ei(n,i); gel(s,i-1) = gen_m1;
760       if (!base) return s;
761       H = matid(n); gel(H,i) = gel(H,1); gel(H,1) = s;
762       return mkvec2(qf_apply_ZM(G,H),H);
763     }
764   if (!R) return G; /* fail */
765   /* case 3: a principal minor is 0 */
766   s = ZM_ker(principal_minor(G, itos(R)));
767   s = vecextend(Q_primpart(gel(s,1)), n);
768   if (!base) return s;
769   H = completebasis(s, 0);
770   gel(H,n) = ZC_neg(gel(H,1)); gel(H,1) = s;
771   return mkvec2(qf_apply_ZM(G,H),H);
772 }
773
774 /* p a prime number, G 3x3 symmetric. Finds X!=0 such that X^t G X = 0 mod p.
775  * Allow returning a shorter X: to be completed with 0s. */
776 static GEN
777 qfsolvemodp(GEN G, GEN p)
778 {
779   GEN a,b,c,d,e,f, v1,v2,v3,v4,v5, x1,x2,x3,N1,N2,N3,s,r;
780
781   /* principal_minor(G,3) = [a,b,d; b,c,e; d,e,f] */
782   a = modii(gcoeff(G,1,1), p);
783   if (!signe(a)) return mkcol(gen_1);
784   v1 = a;
785   b = modii(gcoeff(G,1,2), p);
786   c = modii(gcoeff(G,2,2), p);
787   v2 = modii(subii(mulii(a,c), sqri(b)), p);
788   if (!signe(v2)) return mkcol2(Fp_neg(b,p), a);
789   d = modii(gcoeff(G,1,3), p);
790   e = modii(gcoeff(G,2,3), p);
791   f = modii(gcoeff(G,3,3), p);
792   v4 = modii(subii(mulii(c,d), mulii(e,b)), p);
793   v5 = modii(subii(mulii(a,e), mulii(d,b)), p);
794   v3 = subii(mulii(v2,f), addii(mulii(v4,d), mulii(v5,e))); /* det(G) */
795   v3 = modii(v3, p);
796   N1 =  Fp_neg(v2,  p);
797   x3 = mkcol3(v4, v5, N1);
798   if (!signe(v3)) return x3;
799
800   /* now, solve in dimension 3... reduction to the diagonal case: */
801   x1 = mkcol3(gen_1, gen_0, gen_0);
802   x2 = mkcol3(negi(b), a, gen_0);
803   if (kronecker(N1,p) == 1) return ZC_lincomb(Fp_sqrt(N1,p),gen_1,x1,x2);
804   N2 = Fp_div(Fp_neg(v3,p), v1, p);
805   if (kronecker(N2,p) == 1) return ZC_lincomb(Fp_sqrt(N2,p),gen_1,x2,x3);
806   N3 = Fp_mul(v2, N2, p);
807   if (kronecker(N3,p) == 1) return ZC_lincomb(Fp_sqrt(N3,p),gen_1,x1,x3);
808   r = gen_1;
809   for(;;)
810   {
811     s = Fp_sub(gen_1, Fp_mul(N1,Fp_sqr(r,p),p), p);
812     if (kronecker(s, p) <= 0) break;
813     r = randomi(p);
814   }
815   s = Fp_sqrt(Fp_div(s,N3,p), p);
817 }
818
819 /* assume G square integral */
820 static void
821 check_symmetric(GEN G)
822 {
823   long i,j, l = lg(G);
824   for (i = 1; i < l; i++)
825     for(j = 1; j < i; j++)
826       if (!equalii(gcoeff(G,i,j), gcoeff(G,j,i)))
827         pari_err_TYPE("qfsolve [not symmetric]",G);
828 }
829
830 /* Given a square matrix G of dimension n >= 1, */
831 /* solves over Z the quadratic equation X^tGX = 0. */
832 /* G is assumed to have integral coprime coefficients. */
833 /* The solution might be a vectorv or a matrix. */
834 /* If no solution exists, returns an integer, that can */
835 /* be a prime p such that there is no local solution at p, */
836 /* or -1 if there is no real solution, */
837 /* or 0 in some rare cases. */
838 static  GEN
839 qfsolve_i(GEN G)
840 {
841   GEN M, signG, Min, U, G1, M1, G2, M2, solG2, P, E;
842   GEN solG1, sol, Q, d, dQ, detG2, fam2detG;
843   long n, np, codim, dim;
844   int fail;
845
846   if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
847   n = lg(G)-1;
848   if (n == 0) pari_err_DOMAIN("qfsolve", "dimension" , "=", gen_0, G);
849   if (n != nbrows(G)) pari_err_DIM("qfsolve");
850   G = Q_primpart(G); RgM_check_ZM(G, "qfsolve");
851   check_symmetric(G);
852
853   /* Trivial case: det = 0 */
854   d = ZM_det(G);
855   if (!signe(d))
856   {
857     if (n == 1) return mkcol(gen_1);
858     sol = ZM_ker(G);
859     if (lg(sol) == 2) sol = gel(sol,1);
860     return sol;
861   }
862
863   /* Small dimension: n <= 2 */
864   if (n == 1) return gen_m1;
865   if (n == 2)
866   {
867     GEN t, a =  gcoeff(G,1,1);
868     if (!signe(a)) return mkcol2(gen_1, gen_0);
869     if (signe(d) > 0) return gen_m1; /* no real solution */
870     if (!Z_issquareall(negi(d), &t)) return gen_m2;
871     return mkcol2(subii(t,gcoeff(G,1,2)), a);
872   }
873
874   /* 1st reduction of the coefficients of G */
875   M = qflllgram_indef(G,0,&fail);
876   if (typ(M) == t_COL) return M;
877   G = gel(M,1);
878   M = gel(M,2);
879
880   /* real solubility */
881   signG = ZV_to_zv(qfsign(G));
882   {
883     long r =  signG[1], s = signG[2];
884     if (!r || !s) return gen_m1;
885     if (r < s) { G = ZM_neg(G); signG = mkvecsmall2(s,r);  }
886   }
887
888   /* factorization of the determinant */
889   fam2detG = absZ_factor(d);
890   P = gel(fam2detG,1);
891   E = ZV_to_zv(gel(fam2detG,2));
892   /* P,E = factor(|det(G)|) */
893
894   /* Minimization and local solubility */
895   Min = qfminimize(G, P, E);
896   if (typ(Min) == t_INT) return Min;
897
898   M = RgM_mul(M, gel(Min,2));
899   G = gel(Min,1);
900   P = gel(Min,3);
901   E = gel(Min,4);
902   /* P,E = factor(|det(G))| */
903
904   /* Now, we know that local solutions exist (except maybe at 2 if n==4)
905    * if n==3, det(G) = +-1
906    * if n==4, or n is odd, det(G) is squarefree.
907    * if n>=6, det(G) has all its valuations <=2. */
908
909   /* Reduction of G and search for trivial solutions. */
910   /* When |det G|=1, such trivial solutions always exist. */
911   U = qflllgram_indef(G,0,&fail);
912   if(typ(U) == t_COL) return Q_primpart(RgM_RgC_mul(M,U));
913   G = gel(U,1);
914   M = RgM_mul(M, gel(U,2));
915   /* P,E = factor(|det(G))| */
916
917   /* If n >= 6 is even, need to increment the dimension by 1 to suppress all
918    * squares from det(G) */
919   np = lg(P)-1;
920   if (n < 6 || odd(n) || !np)
921   {
922     codim = 0;
923     G1 = G;
924     M1 = NULL;
925   }
926   else
927   {
928     GEN aux;
929     long i;
930     codim = 1; n++;
931     /* largest square divisor of d */
932     aux = gen_1;
933     for (i = 1; i <= np; i++)
934       if (E[i] == 2) { aux = mulii(aux, gel(P,i)); E[i] = 3; }
935     /* Choose sign(aux) so as to balance the signature of G1 */
936     if (signG[1] > signG[2])
937     {
938       signG[2]++;
939       aux = negi(aux);
940     }
941     else
942       signG[1]++;
943     G1 = shallowmatconcat(diagonal_shallow(mkvec2(G,aux)));
944     /* P,E = factor(|det G1|) */
945     Min = qfminimize(G1, P, E);
946     G1 = gel(Min,1);
947     M1 = gel(Min,2);
948     P = gel(Min,3);
949     E = gel(Min,4);
950     np = lg(P)-1;
951   }
952
953   /* now, d is squarefree */
954   if (!np)
955   { /* |d| = 1 */
956      G2 = G1;
957      M2 = NULL;
958   }
959   else
960   { /* |d| > 1: increment dimension by 2 */
961     GEN factdP, factdE, W;
962     long i, lfactdP;
963     codim += 2;
964     d = ZV_prod(P); /* d = abs(matdet(G1)); */
965     if (odd(signG[2])) togglesign_safe(&d); /* d = matdet(G1); */
966     /* solubility at 2 (this is the only remaining bad prime). */
967     if (n == 4 && smodis(d,8) == 1 && qflocalinvariant(G,gen_2) == 1)
968       return gen_2;
969
970     P = shallowconcat(mpodd(d)? mkvec2(NULL,gen_2): mkvec(NULL), P);
971     /* build a binary quadratic form with given Witt invariants */
972     W = const_vecsmall(lg(P)-1, 0);
973     /* choose signature of Q (real invariant and sign of the discriminant) */
974     dQ = absi(d);
975     if (signG[1] > signG[2]) togglesign_safe(&dQ); /* signQ = [2,0]; */
976     if (n == 4 && smodis(dQ,4) != 1) dQ = shifti(dQ,2);
977     if (n >= 5) dQ = shifti(dQ,3);
978
980     if (n == 4)
981     {
982       GEN t = qflocalinvariants(ZM_neg(G1),P);
983       for (i = 3; i < lg(P); i++) W[i] = ucoeff(t,i,1);
984     }
985     else
986     {
987       long s = signe(dQ) == signe(d)? 1: -1;
988       GEN t;
989       if (odd((n-3)/2)) s = -s;
990       t = s > 0? utoipos(8): utoineg(8);
991       for (i = 3; i < lg(P); i++)
992         W[i] = hilbertii(t, gel(P,i), gel(P,i)) > 0;
993     }
994     /* for p = 2, the choice is fixed from the product formula */
995     W[2] = Flv_sum(W, 2);
996
997     /* Construction of the 2-class group of discriminant dQ until some product
998      * of the generators gives the desired invariants. */
999     factdP = vecsplice(P, 1); lfactdP =  lg(factdP);
1000     factdE = cgetg(lfactdP, t_VECSMALL);
1001     for (i = 1; i < lfactdP; i++) factdE[i] = Z_pval(dQ, gel(factdP,i));
1002     factdE[1]++;
1003     /* factdP,factdE = factor(2|dQ|), P = factor(-2|dQ|)[,1] */
1004     Q = quadclass2(dQ, factdP,factdE, P, W, n == 4);
1005     /* Build a form of dim=n+2 potentially unimodular */
1006     G2 = shallowmatconcat(diagonal_shallow(mkvec2(G1,ZM_neg(Q))));
1007     /* Minimization of G2 */
1008     detG2 = mulii(d, ZM_det(Q));
1009     for (i = 1; i < lfactdP; i++) factdE[i] = Z_pval(detG2, gel(factdP,i));
1010     /* factdP,factdE = factor(|det G2|) */
1011     Min = qfminimize(G2, factdP,factdE);
1012     M2 = gel(Min,2);
1013     G2 = gel(Min,1);
1014   }
1015   /* |det(G2)| = 1, find a totally isotropic subspace for G2 */
1016   solG2 = qflllgram_indefgoon(G2);
1017   /* G2 must have a subspace of solutions of dimension > codim */
1018   dim = codim+2;
1019   while(gequal0(principal_minor(gel(solG2,1), dim))) dim ++;
1020   solG2 = vecslice(gel(solG2,2), 1, dim-1);
1021
1022   if (!M2)
1023     solG1 = solG2;
1024   else
1025   { /* solution of G1 is simultaneously in solG2 and x[n+1] = x[n+2] = 0*/
1026     GEN K;
1027     solG1 = RgM_mul(M2,solG2);
1028     K = ker(rowslice(solG1,n+1,n+2));
1029     solG1 = RgM_mul(rowslice(solG1,1,n), K);
1030   }
1031   if (!M1)
1032     sol = solG1;
1033   else
1034   { /* solution of G1 is simultaneously in solG2 and x[n] = 0 */
1035     GEN K;
1036     sol = RgM_mul(M1,solG1);
1037     K = ker(rowslice(sol,n,n));
1038     sol = RgM_mul(rowslice(sol,1,n-1), K);
1039   }
1040   sol = Q_primpart(RgM_mul(M, sol));
1041   if (lg(sol) == 2) sol = gel(sol,1);
1042   return sol;
1043 }
1044 GEN
1045 qfsolve(GEN G)
1046 {
1047   pari_sp av = avma;
1048   return gerepilecopy(av, qfsolve_i(G));
1049 }
1050
1051 /* G is a symmetric 3x3 matrix, and sol a solution of sol~*G*sol=0.
1052  * Returns a parametrization of the solutions with the good invariants,
1053  * as a matrix 3x3, where each line contains
1054  * the coefficients of each of the 3 quadratic forms.
1055  * If fl!=0, the fl-th form is reduced. */
1056 GEN
1057 qfparam(GEN G, GEN sol, long fl)
1058 {
1059   pari_sp av = avma;
1060   GEN U, G1, G2, a, b, c, d, e;
1061   long n, tx = typ(sol);
1062
1063   if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
1064   if (!is_vec_t(tx)) pari_err_TYPE("qfsolve", G);
1065   if (tx == t_VEC) sol = shallowtrans(sol);
1066   n = lg(G)-1;
1067   if (n == 0) pari_err_DOMAIN("qfsolve", "dimension" , "=", gen_0, G);
1068   if (n != nbrows(G) || n != 3 || lg(sol) != 4) pari_err_DIM("qfsolve");
1069   G = Q_primpart(G); RgM_check_ZM(G,"qfsolve");
1070   check_symmetric(G);
1071   sol = Q_primpart(sol); RgV_check_ZV(sol,"qfsolve");
1072   /* build U such that U[,3] = sol, and |det(U)| = 1 */
1073   U = completebasis(sol,1);
1074   G1 = qf_apply_ZM(G,U); /* G1 has a 0 at the bottom right corner */
1075   a = shifti(gcoeff(G1,1,2),1);
1076   b = shifti(negi(gcoeff(G1,1,3)),1);
1077   c = shifti(negi(gcoeff(G1,2,3)),1);
1078   d = gcoeff(G1,1,1);
1079   e = gcoeff(G1,2,2);
1080   G2 = mkmat3(mkcol3(b,gen_0,d), mkcol3(c,b,a), mkcol3(gen_0,c,e));
1081   sol = ZM_mul(U,G2);
1082   if (fl)
1083   {
1084     GEN v = row(sol,fl);
1085     int fail;
1086     a = gel(v,1);
1087     b = gmul2n(gel(v,2),-1);
1088     c = gel(v,3);
1089     U = qflllgram_indef(mkmat2(mkcol2(a,b),mkcol2(b,c)), 1, &fail);
1090     U = gel(U,2);
1091     a = gcoeff(U,1,1); b = gcoeff(U,1,2);
1092     c = gcoeff(U,2,1); d = gcoeff(U,2,2);
1093     U = mkmat3(mkcol3(sqri(a),mulii(a,c),sqri(c)),