Line data Source code
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
6 : terms of the GNU General Public License as published by the Free Software
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 : #define DEBUGLEVEL DEBUGLEVEL_qfsolve
25 :
26 : /* LINEAR ALGEBRA */
27 : /* complete by 0s, assume l-1 <= n */
28 : static GEN
29 117346 : vecextend(GEN v, long n)
30 : {
31 117346 : long i, l = lg(v);
32 117346 : GEN w = cgetg(n+1, t_COL);
33 362228 : for (i = 1; i < l; i++) gel(w,i) = gel(v,i);
34 388391 : for ( ; i <=n; i++) gel(w,i) = gen_0;
35 117346 : return w;
36 : }
37 :
38 : /* Gives a unimodular matrix with the last column(s) equal to Mv.
39 : * Mv can be a column vector or a rectangular matrix.
40 : * redflag = 0 or 1. If redflag = 1, LLL-reduce the n-#v first columns. */
41 : static GEN
42 410926 : completebasis(GEN Mv, long redflag)
43 : {
44 : GEN U;
45 : long m, n;
46 :
47 410926 : if (typ(Mv) == t_COL) Mv = mkmat(Mv);
48 410926 : n = lg(Mv)-1;
49 410926 : m = nbrows(Mv); /* m x n */
50 410926 : if (m == n) return Mv;
51 410401 : (void)ZM_hnfall_i(shallowtrans(Mv), &U, 0);
52 410401 : U = ZM_inv(shallowtrans(U),NULL);
53 410401 : if (m==1 || !redflag) return U;
54 : /* LLL-reduce the m-n first columns */
55 124041 : return shallowconcat(ZM_lll(vecslice(U,1,m-n), 0.99, LLL_INPLACE),
56 124041 : vecslice(U, m-n+1,m));
57 : }
58 :
59 : /* Return U in GLn(Z) whose first d columns span Ker (M mod p). */
60 : static GEN
61 286622 : kermodp(GEN M, GEN p, long *d)
62 : {
63 : long j, l;
64 : GEN K, B, U;
65 :
66 286622 : K = FpM_center(FpM_ker(M, p), p, shifti(p,-1));
67 286622 : B = completebasis(K,0);
68 286622 : l = lg(M); U = cgetg(l, t_MAT);
69 1366953 : for (j = 1; j < l; j++) gel(U,j) = gel(B,l-j);
70 286622 : *d = lg(K)-1; return U;
71 : }
72 :
73 : /* INVARIANTS COMPUTATIONS */
74 :
75 : static GEN
76 33253 : principal_minor(GEN G, long i) { return matslice(G,1,i,1,i); }
77 : static GEN
78 707 : det_minors(GEN G)
79 : {
80 707 : long i, l = lg(G);
81 707 : GEN v = cgetg(l+1, t_VEC);
82 707 : gel(v,1) = gen_1;
83 3535 : for (i = 2; i <= l; i++) gel(v,i) = ZM_det(principal_minor(G,i-1));
84 707 : return v;
85 : }
86 :
87 : static GEN
88 7280 : hilberts(GEN a, GEN b, GEN P)
89 : {
90 7280 : long i, lP = lg(P);
91 7280 : GEN v = cgetg(lP, t_VECSMALL);
92 42217 : for (i = 1; i < lP; i++) v[i] = hilbertii(a, b, gel(P,i)) < 0;
93 7280 : return v;
94 : }
95 :
96 : /* 4 | disc(q); special case of gtomat */
97 : static GEN
98 1666 : qfbmat(GEN q)
99 : {
100 1666 : GEN a = gel(q,1), b = shifti(gel(q,2), -1), c = gel(q,3);
101 1666 : retmkmat22(a, b, b, c);
102 : }
103 : /* 2*qfbmat(q) */
104 : static GEN
105 21 : qfbmat2(GEN q)
106 : {
107 21 : GEN a = shifti(gel(q,1), 1), b = gel(q,2), c = shifti(gel(q,3), 1);
108 21 : retmkmat22(a, b, b, c);
109 : }
110 : /* Given a symmetric matrix G over Z, compute the Witt invariant of G at p
111 : * v = det_minors(G) [G diagonalized]; assume that none of the v[i] is 0. */
112 : static long
113 2128 : witt(GEN v, GEN p)
114 : {
115 2128 : long k = lg(v)-2, h = hilbertii(gel(v,k), gel(v,k+1), p);
116 8512 : for (k--; k >= 1; k--)
117 6384 : if (hilbertii(negi(gel(v,k)), gel(v,k+1),p) < 0) h = -h;
118 2128 : return h;
119 : }
120 :
121 : /* QUADRATIC FORM REDUCTION */
122 : /* private version of qfgaussred:
123 : * - early abort if k-th principal minor is singular, return stoi(k)
124 : * - else return a matrix whose upper triangular part is qfgaussred(a) */
125 : static GEN
126 44881 : partialgaussred(GEN a)
127 : {
128 44881 : long n = lg(a)-1, k;
129 44881 : a = RgM_shallowcopy(a);
130 148963 : for(k = 1; k < n; k++)
131 : {
132 113587 : GEN ak, p = gcoeff(a,k,k);
133 : long i, j;
134 113587 : if (isintzero(p)) return stoi(k);
135 104082 : ak = row(a, k);
136 385290 : for (i=k+1; i<=n; i++) gcoeff(a,k,i) = gdiv(gcoeff(a,k,i), p);
137 385290 : for (i=k+1; i<=n; i++)
138 : {
139 281208 : GEN c = gel(ak,i);
140 281208 : if (gequal0(c)) continue;
141 844598 : for (j=i; j<=n; j++)
142 616113 : gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
143 : }
144 : }
145 35376 : if (isintzero(gcoeff(a,n,n))) return stoi(n);
146 35355 : return a;
147 : }
148 :
149 : /* LLL-reduce a positive definite qf QD bounding the indefinite G, dim G > 1.
150 : * Then finishes by looking for trivial solution */
151 : static GEN qftriv(GEN G, GEN z, long base);
152 : static GEN
153 44881 : qflllgram_indef(GEN G, long base, int *fail)
154 : {
155 : GEN M, R, g, DM, S, dR;
156 44881 : long i, j, n = lg(G)-1;
157 :
158 44881 : *fail = 0;
159 44881 : R = partialgaussred(G);
160 44881 : if (typ(R) == t_INT) return qftriv(G, R, base);
161 35355 : R = Q_remove_denom(R, &dR); /* avoid rational arithmetic */
162 35355 : M = zeromatcopy(n,n);
163 35355 : DM = zeromatcopy(n,n);
164 170074 : for (i = 1; i <= n; i++)
165 : {
166 134719 : GEN d = absi_shallow(gcoeff(R,i,i));
167 134719 : if (dR) {
168 110505 : gcoeff(M,i,i) = dR;
169 110505 : gcoeff(DM,i,i) = mulii(d,dR);
170 : } else {
171 24214 : gcoeff(M,i,i) = gen_1;
172 24214 : gcoeff(DM,i,i) = d;
173 : }
174 395797 : for (j = i+1; j <= n; j++)
175 : {
176 261078 : gcoeff(M,i,j) = gcoeff(R,i,j);
177 261078 : gcoeff(DM,i,j) = mulii(d, gcoeff(R,i,j));
178 : }
179 : }
180 : /* G = M~*D*M, D diagonal, DM=|D|*M, g = M~*|D|*M */
181 35355 : g = ZM_transmultosym(M,DM);
182 35355 : S = lllgramint(Q_primpart(g));
183 35355 : R = qftriv(qf_ZM_apply(G,S), NULL, base);
184 35355 : switch(typ(R))
185 : {
186 6216 : case t_COL: return ZM_ZC_mul(S,R);
187 21440 : case t_MAT: *fail = 1; return mkvec2(R, S);
188 7699 : default:
189 7699 : gel(R,2) = ZM_mul(S, gel(R,2));
190 7699 : return R;
191 : }
192 : }
193 :
194 : /* G symmetric, i < j, let E = E_{i,j}(a), G <- E~*G*E, U <- U*E.
195 : * Everybody integral */
196 : static void
197 4205 : qf_apply_transvect_Z(GEN G, GEN U, long i, long j, GEN a)
198 : {
199 4205 : long k, n = lg(G)-1;
200 4205 : gel(G, j) = ZC_lincomb(gen_1, a, gel(G,j), gel(G,i));
201 23089 : for (k = 1; k < n; k++) gcoeff(G, j, k) = gcoeff(G, k, j);
202 4205 : gcoeff(G,j,j) = addmulii(gcoeff(G,j,j), a,
203 4205 : addmulii(gcoeff(G,i,j), a,gcoeff(G,i,i)));
204 4205 : gel(U, j) = ZC_lincomb(gen_1, a, gel(U,j), gel(U,i));
205 4205 : }
206 :
207 : /* LLL reduction of the quadratic form G (Gram matrix)
208 : * where we go on, even if an isotropic vector is found. */
209 : static GEN
210 11095 : qflllgram_indefgoon(GEN G)
211 : {
212 : GEN red, U, A, U1,U2,U3,U5,U6, V, B, G2,G3,G4,G5, G6, a, g;
213 11095 : long i, j, n = lg(G)-1;
214 : int fail;
215 :
216 11095 : red = qflllgram_indef(G,1, &fail);
217 11095 : if (fail) return red; /*no isotropic vector found: nothing to do*/
218 : /* otherwise a solution is found: */
219 11074 : U1 = gel(red,2);
220 11074 : G2 = gel(red,1); /* G2[1,1] = 0 */
221 11074 : U2 = gel(ZV_extgcd(row(G2,1)), 2);
222 11074 : G3 = qf_ZM_apply(G2,U2);
223 11074 : U = ZM_mul(U1,U2); /* qf_apply(G,U) = G3 */
224 : /* G3[1,] = [0,...,0,g], g^2 | det G */
225 11074 : g = gcoeff(G3,1,n);
226 11074 : a = diviiround(negi(gcoeff(G3,n,n)), shifti(g,1));
227 11074 : if (signe(a)) qf_apply_transvect_Z(G3,U,1,n,a);
228 : /* G3[n,n] reduced mod 2g */
229 11074 : if (n == 2) return mkvec2(G3,U);
230 10409 : V = rowpermute(vecslice(G3, 2,n-1), mkvecsmall2(1,n));
231 10409 : A = mkmat22(gcoeff(G3,1,1),gcoeff(G3,1,n),gcoeff(G3,1,n),gcoeff(G3,2,2));
232 10409 : B = ground(RgM_neg(QM_mul(QM_inv(A), V)));
233 10409 : U3 = matid(n);
234 51835 : for (j = 2; j < n; j++)
235 : {
236 41426 : gcoeff(U3,1,j) = gcoeff(B,1,j-1);
237 41426 : gcoeff(U3,n,j) = gcoeff(B,2,j-1);
238 : }
239 10409 : G4 = qf_ZM_apply(G3,U3); /* the last column of G4 is reduced */
240 10409 : U = ZM_mul(U,U3);
241 10409 : if (n == 3) return mkvec2(G4,U);
242 :
243 8141 : red = qflllgram_indefgoon(matslice(G4,2,n-1,2,n-1));
244 8141 : if (typ(red) == t_MAT) return mkvec2(G4,U);
245 : /* Let U5:=matconcat(diagonal[1,red[2],1])
246 : * return [qf_ZM_apply(G5, U5), U*U5] */
247 8141 : G5 = gel(red,1);
248 8141 : U5 = gel(red,2);
249 8141 : G6 = cgetg(n+1,t_MAT);
250 8141 : gel(G6,1) = gel(G4,1);
251 8141 : gel(G6,n) = gel(G4,n);
252 47299 : for (j=2; j<n; j++)
253 : {
254 39158 : gel(G6,j) = cgetg(n+1,t_COL);
255 39158 : gcoeff(G6,1,j) = gcoeff(G4,j,1);
256 39158 : gcoeff(G6,n,j) = gcoeff(G4,j,n);
257 258958 : for (i=2; i<n; i++) gcoeff(G6,i,j) = gcoeff(G5,i-1,j-1);
258 : }
259 8141 : U6 = mkvec3(mkmat(gel(U,1)), ZM_mul(vecslice(U,2,n-1),U5), mkmat(gel(U,n)));
260 8141 : return mkvec2(G6, shallowconcat1(U6));
261 : }
262 :
263 : /* qf_ZM_apply(G,H), where H = matrix of \tau_{i,j}, i != j */
264 : static GEN
265 13324 : qf_apply_tau(GEN G, long i, long j)
266 : {
267 13324 : long l = lg(G), k;
268 13324 : G = RgM_shallowcopy(G);
269 13324 : swap(gel(G,i), gel(G,j));
270 80556 : for (k = 1; k < l; k++) swap(gcoeff(G,i,k), gcoeff(G,j,k));
271 13324 : return G;
272 : }
273 :
274 : /* LLL reduction of the quadratic form G (Gram matrix)
275 : * in dim 3 only, with detG = -1 and sign(G) = [2,1]; */
276 : static GEN
277 2275 : qflllgram_indefgoon2(GEN G)
278 : {
279 : GEN red, G2, a, b, c, d, e, f, u, v, r, r3, U2, G3;
280 : int fail;
281 :
282 2275 : red = qflllgram_indef(G,1,&fail); /* always find an isotropic vector. */
283 2275 : G2 = qf_apply_tau(gel(red,1),1,3); /* G2[3,3] = 0 */
284 2275 : r = row(gel(red,2), 3);
285 2275 : swap(gel(r,1), gel(r,3)); /* apply tau_{1,3} */
286 2275 : a = gcoeff(G2,3,1);
287 2275 : b = gcoeff(G2,3,2);
288 2275 : d = bezout(a,b, &u,&v);
289 2275 : if (!equali1(d))
290 : {
291 0 : a = diviiexact(a,d);
292 0 : b = diviiexact(b,d);
293 : }
294 : /* for U2 = [-u,-b,0;-v,a,0;0,0,1]
295 : * G3 = qf_ZM_apply(G2,U2) has known last row (-d, 0, 0),
296 : * so apply to principal_minor(G3,2), instead */
297 2275 : U2 = mkmat22(negi(u),negi(b),negi(v),a);
298 2275 : G3 = qf_ZM_apply(principal_minor(G2,2),U2);
299 2275 : r3 = gel(r,3);
300 2275 : r = ZV_ZM_mul(mkvec2(gel(r,1),gel(r,2)),U2);
301 :
302 2275 : a = gcoeff(G3,1,1);
303 2275 : b = gcoeff(G3,1,2);
304 2275 : c = negi(d); /* G3[1,3] */
305 2275 : d = gcoeff(G3,2,2);
306 2275 : if (mpodd(a))
307 : {
308 1204 : e = addii(b,d);
309 1204 : a = addii(a, addii(b,e));
310 1204 : e = diviiround(negi(e),c);
311 1204 : f = diviiround(negi(a), shifti(c,1));
312 1204 : a = addmulii(addii(gel(r,1),gel(r,2)), f,r3);
313 : }
314 : else
315 : {
316 1071 : e = diviiround(negi(b),c);
317 1071 : f = diviiround(negi(shifti(a,-1)), c);
318 1071 : a = addmulii(gel(r,1), f, r3);
319 : }
320 2275 : b = addmulii(gel(r,2), e, r3);
321 2275 : return mkvec3(a,b, r3);
322 : }
323 :
324 : /* QUADRATIC FORM MINIMIZATION */
325 : /* G symmetric, return ZM_Z_divexact(G,d) */
326 : static GEN
327 20328 : ZsymM_Z_divexact(GEN G, GEN d)
328 : {
329 20328 : long i,j,l = lg(G);
330 20328 : GEN H = cgetg(l, t_MAT);
331 79149 : for (j = 1; j < l; j++)
332 : {
333 58821 : GEN c = cgetg(l, t_COL), b = gel(G,j);
334 115486 : for (i=1; i < j; i++) gcoeff(H,j,i) = gel(c,i) = diviiexact(gel(b,i),d);
335 58821 : gel(c,j) = diviiexact(gel(b,j),d);
336 58821 : gel(H,j) = c;
337 : }
338 20328 : return H;
339 : }
340 : /* by blocks, in place: G[1,1] /= d, G[2,2] *= d */
341 : static void
342 92235 : ZsymM_Z_divexact_partial(GEN G, long n, GEN d)
343 : {
344 92235 : long i, j, l = lg(G);
345 185058 : for(j = 1; j <= n; j++)
346 : {
347 93411 : for (i = 1; i < j; i++)
348 588 : gcoeff(G,i,j) = gcoeff(G,j,i) = diviiexact(gcoeff(G,i,j), d);
349 92823 : gcoeff(G,j,j) = diviiexact(gcoeff(G,j,j), d);
350 : }
351 309829 : for (; j < l; j++)
352 : {
353 454722 : for (i = n+1; i < j; i++)
354 237128 : gcoeff(G,i,j) = gcoeff(G,j,i) = mulii(gcoeff(G,i,j), d);
355 217594 : gcoeff(G,j,j) = mulii(gcoeff(G,j,j), d);
356 : }
357 92235 : }
358 :
359 : /* write symmetric G as [A,B;B~,C], A dxd, C (n-d)x(n-d) */
360 : static void
361 2289 : blocks4(GEN G, long d, long n, GEN *A, GEN *B, GEN *C)
362 : {
363 2289 : GEN G2 = vecslice(G,d+1,n);
364 2289 : *A = principal_minor(G, d);
365 2289 : *B = rowslice(G2, 1, d);
366 2289 : *C = rowslice(G2, d+1, n);
367 2289 : }
368 : static GEN qfsolvemodp(GEN G, GEN p);
369 : static void
370 72676 : update_fm(GEN f, GEN a, long i)
371 : {
372 72676 : if (!odd(i))
373 51169 : gel(f,1) = a;
374 : else
375 : {
376 21507 : long v = vals(i+1), k;
377 21507 : GEN b = gel(f,1), u = QM_mul(b, a);
378 21507 : gel(f,1) = gen_0;
379 21507 : if (v+2 >= lg(f)) pari_err_BUG("update_fm");
380 26152 : for (k = 1; k < v; k++)
381 : {
382 4645 : u = QM_mul(gel(f, k+1), u);
383 4645 : gel(f,k+1) = gen_0; /* for (void)gc_all */
384 : }
385 21507 : gel(f,v+1) = u;
386 : }
387 72676 : }
388 : static GEN
389 41286 : prod_fm(GEN f, long i)
390 : {
391 41286 : long v = vals(i), k;
392 41286 : GEN u = gel(f, ++v);
393 47819 : for (i >>= v, k = v+1; i; i >>= 1, k++)
394 6533 : if (odd(i)) u = QM_mul(gel(f,k), u);
395 41286 : return u;
396 : }
397 :
398 : /* Minimization of the quadratic form G, deg G != 0, dim n >= 2
399 : * G symmetric integral
400 : * Returns [G',U,factd] with U in GLn(Q) such that G'=U~*G*U*constant
401 : * is integral and has minimal determinant.
402 : * In dimension 3 or 4, may return a prime p if the reduction at p is
403 : * impossible because of local nonsolvability.
404 : * P,E = factor(+/- det(G)), d = det(G) "prime" -1 is ignored,
405 : * Either E or d should be NULL, but not both */
406 : static GEN
407 16183 : qfminimize_fact(GEN G, GEN P, GEN E, GEN d, long loc)
408 : {
409 16183 : GEN U = NULL, Ker = NULL, faE, faP;
410 16183 : long n = lg(G)-1, lP = lg(P), i;
411 16183 : faP = vectrunc_init(lP);
412 16183 : faE = vecsmalltrunc_init(lP);
413 157597 : for (i = 1; i < lP; i++)
414 : {
415 : long Ei, vp, wp;
416 141897 : GEN p = gel(P,i);
417 141897 : if (is_pm1(p)) continue;
418 141295 : Ei = E ? E[i]: Z_pval(d, p); vp = Ei; wp = vp;
419 141295 : if (!vp) continue;
420 140833 : if (DEBUGLEVEL >= 3) err_printf("qfminimize: for %Ps^%ld:", p,vp);
421 347576 : while (vp)
422 : {
423 248818 : long idx = 0, dimKer = 0; /* -Wall */
424 248818 : GEN d, sol = NULL, FU = zerovec(2*expu(vp)+100);
425 248818 : pari_sp av = avma;
426 321494 : while (vp) /* loop until vp <= n */
427 : {
428 291900 : if (DEBUGLEVEL>=3 && vp <= wp)
429 0 : { err_printf(" %ld%%", (Ei-vp)*100/Ei); wp -= Ei/100; }
430 : /* The case vp = 1 can be minimized only if n is odd. */
431 291900 : if (vp == 1 && !odd(n)) break;
432 284333 : Ker = kermodp(G,p, &dimKer); /* dimKer <= vp */
433 284333 : if (DEBUGLEVEL >= 4) err_printf(" dimKer = %ld\n",dimKer);
434 284333 : if (dimKer == n) break;
435 284333 : G = qf_ZM_apply(G, Ker);
436 : /* 1st case: dimKer < vp */
437 : /* then the kernel mod p contains a kernel mod p^2 */
438 284333 : if (dimKer >= vp) break;
439 :
440 72676 : if (DEBUGLEVEL >= 4) err_printf(" case 1: dimker < vp\n");
441 72676 : if (dimKer == 1)
442 : {
443 : long j;
444 70387 : gcoeff(G,1,1) = diviiexact(gcoeff(G,1,1), sqri(p));
445 215641 : for (j = 2; j <= n; j++)
446 145254 : gcoeff(G,1,j) = gcoeff(G,j,1) = diviiexact(gcoeff(G,j,1), p);
447 70387 : gel(Ker,1) = RgC_Rg_div(gel(Ker,1), p);
448 70387 : vp -= 2;
449 : }
450 : else
451 : {
452 : GEN A, B, C, K2;
453 : long j, dimKer2;
454 2289 : blocks4(G, dimKer,n, &A,&B,&C); A = ZsymM_Z_divexact(A, p);
455 2289 : K2 = kermodp(A, p, &dimKer2);
456 : /* Write G = [pA,B;B~,C] and apply [K2/p,0;0,Id] by blocks */
457 2289 : A = qf_ZM_apply(A,K2); ZsymM_Z_divexact_partial(A, dimKer2, p);
458 2289 : B = ZM_transmul(B,K2);
459 5131 : for (j = 1; j <= dimKer2; j++) gel(B,j) = ZC_Z_divexact(gel(B,j), p);
460 2289 : G = shallowmatconcat(mkmat22(A,shallowtrans(B),B,C));
461 : /* Ker *= [K2,0;0,Id] */
462 2289 : B = ZM_mul(vecslice(Ker,1,dimKer),K2);
463 5131 : for (j = 1; j <= dimKer2; j++) gel(B,j) = RgC_Rg_div(gel(B,j), p);
464 2289 : Ker = shallowconcat(B, vecslice(Ker,dimKer+1,n));
465 2289 : vp -= 2*dimKer2;
466 : }
467 72676 : update_fm(FU, Ker, idx++);
468 72676 : if (gc_needed(av, 1))
469 : {
470 0 : if (DEBUGMEM >= 2) pari_warn(warnmem,"qfminimize");
471 0 : (void)gc_all(av, 2, &G, &FU);
472 : }
473 : }
474 248818 : if (idx)
475 : {
476 41286 : GEN PU = prod_fm(FU, idx);
477 41286 : U = U ? QM_mul(U, PU): PU;
478 : }
479 260816 : if (vp == 0) break;
480 219224 : if (vp == 1 && !odd(n))
481 : {
482 7567 : vectrunc_append(faP, p);
483 7567 : vecsmalltrunc_append(faE, 1);
484 7567 : break;
485 : }
486 211657 : if (dimKer == n)
487 : { /* trivial case: dimKer = n */
488 0 : if (DEBUGLEVEL >= 4) err_printf(" case 0: dimKer = n\n");
489 0 : G = ZsymM_Z_divexact(G, p);
490 206743 : vp -= n; continue;
491 : }
492 211657 : U = U ? QM_mul(U, Ker): Ker;
493 : /* vp = dimKer
494 : * 2nd case: kernel has dim >= 2 and contains an elt of norm 0 mod p^2,
495 : * find it */
496 211657 : if (dimKer > 2) {
497 18039 : if (DEBUGLEVEL >= 4) err_printf(" case 2.1\n");
498 18039 : dimKer = 3;
499 18039 : sol = qfsolvemodp(ZsymM_Z_divexact(principal_minor(G,3),p), p);
500 18039 : sol = FpC_red(sol, p);
501 : }
502 193618 : else if (dimKer == 2)
503 : {
504 98807 : GEN a = modii(diviiexact(gcoeff(G,1,1),p), p);
505 98807 : GEN b = modii(diviiexact(gcoeff(G,1,2),p), p);
506 98807 : GEN c = modii(diviiexact(gcoeff(G,2,2),p), p);
507 98807 : GEN D = modii(subii(sqri(b), mulii(a,c)), p);
508 98807 : if (kronecker(D,p) >= 0)
509 : {
510 98758 : if (DEBUGLEVEL >= 4) err_printf(" case 2.2\n");
511 98758 : sol = signe(a)? mkcol2(Fp_sub(Fp_sqrt(D,p), b, p), a): vec_ei(2,1);
512 : }
513 : }
514 211657 : if (sol)
515 116797 : {
516 : long j;
517 116797 : sol = FpC_center(sol, p, shifti(p,-1));
518 116797 : sol = Q_primpart(sol);
519 116797 : if (DEBUGLEVEL >= 4) err_printf(" sol = %Ps\n", sol);
520 116797 : Ker = completebasis(vecextend(sol,n), 1);
521 116797 : G = qf_ZM_apply(G, Ker);
522 513197 : for (j = 1; j < n; j++)
523 396400 : gcoeff(G,n,j) = gcoeff(G,j,n) = diviiexact(gcoeff(G,j,n), p);
524 116797 : gcoeff(G,n,n) = diviiexact(gcoeff(G,n,n), sqri(p));
525 116797 : U = QM_mul(U,Ker); gel(U,n) = RgC_Rg_div(gel(U,n), p);
526 116797 : vp -= 2; continue;
527 : }
528 : /* Now 0 < vp = dimKer < 3 and kernel contains no vector with norm p^2 */
529 : /* exchanging kernel and image makes minimization easier ? */
530 94860 : d = ZM_det(G); if (odd((n-3) / 2)) d = negi(d);
531 94860 : if ((vp==1 && kronecker(gmod(gdiv(negi(d), gcoeff(G,1,1)),p), p) >= 0)
532 4949 : || (vp==2 && odd(n) && n >= 5)
533 4935 : || (vp==2 && !odd(n) && kronecker(modii(diviiexact(d,sqri(p)), p),p) < 0))
534 89946 : {
535 : long j;
536 89946 : if (DEBUGLEVEL >= 4) err_printf(" case 3\n");
537 89946 : ZsymM_Z_divexact_partial(G, dimKer, p);
538 305678 : for (j = dimKer+1; j <= n; j++) gel(U,j) = RgC_Rg_mul(gel(U,j), p);
539 89946 : vp -= 2*dimKer-n; continue;
540 : }
541 :
542 : /* Minimization was not possible so far. */
543 : /* If n == 3 or 4, this proves the local nonsolubility at p. */
544 4914 : if (loc && (n == 3 || n == 4))
545 : {
546 483 : if (DEBUGLEVEL >= 1) err_printf(" no local solution at %Ps\n",p);
547 483 : return p;
548 : }
549 4431 : vectrunc_append(faP, p);
550 4431 : vecsmalltrunc_append(faE, vp);
551 4431 : break;
552 : }
553 140350 : if (DEBUGLEVEL >= 3) err_printf("\n");
554 : }
555 15700 : if (!U) U = matid(n);
556 : else
557 : { /* apply LLL to avoid coefficient explosion */
558 14629 : GEN u = ZM_lll(Q_primpart(U), .99, LLL_IM);
559 14629 : G = qf_ZM_apply(G, u);
560 14629 : U = QM_mul(U, u);
561 : }
562 15700 : return mkvec4(G, U, faP, faE);
563 : }
564 :
565 : /* assume G square integral */
566 : static void
567 19976 : check_symmetric(GEN G)
568 : {
569 19976 : long i,j, l = lg(G);
570 90432 : for (i = 1; i < l; i++)
571 176241 : for(j = 1; j < i; j++)
572 105785 : if (!equalii(gcoeff(G,i,j), gcoeff(G,j,i)))
573 7 : pari_err_TYPE("qfsolve [not symmetric]",G);
574 19969 : }
575 : /* assume G symmetric and integral */
576 : static void
577 14 : symmetric_non0_coeff(GEN G, long *pi, long *pj)
578 : {
579 14 : long i, j, l = lg(G);
580 14 : *pi = *pj = 0;
581 14 : for (i = 1; i < l; i++)
582 14 : for (j = 1; j <= i; j++)
583 14 : if (signe(gcoeff(G,i,j))) { *pi = i; *pj = j; return; }
584 : }
585 :
586 : GEN
587 14 : qfminimize(GEN G)
588 : {
589 14 : pari_sp av = avma;
590 : GEN c, d, F, H, U;
591 14 : long i, j, n = lg(G)-1;
592 14 : if (typ(G) != t_MAT) pari_err_TYPE("qfminimize", G);
593 14 : if (n == 0) pari_err_DOMAIN("qfminimize", "dimension", "=", gen_0, G);
594 14 : if (n != nbrows(G)) pari_err_DIM("qfminimize");
595 14 : G = Q_primpart(G); RgM_check_ZM(G, "qfminimize");
596 14 : check_symmetric(G);
597 14 : d = ZM_det(G);
598 14 : if (!signe(d)) pari_err_DOMAIN("qfminimize", "det", "=", gen_0, gen_0);
599 14 : F = absZ_factor(d);
600 14 : H = qfminimize_fact(G, gel(F,1), ZV_to_zv(gel(F,2)), NULL, 0);
601 14 : symmetric_non0_coeff(G, &i, &j);
602 14 : U = gel(H,2); H = gel(H,1);
603 14 : c = gdiv(gcoeff(H,i,j), RgV_dotproduct(gel(U,i), RgM_RgC_mul(G, gel(U,j))));
604 14 : return gc_GEN(av, mkvec3(H, U, c));
605 : }
606 :
607 : /* CLASS GROUP COMPUTATIONS */
608 :
609 : /* Compute the square root of the quadratic form q of discriminant D = 4 * md
610 : * Not fully implemented; only works for D squarefree except at 2, where the
611 : * valuation is 2 or 3. Finally, [P,E] = factor(2*abs(D)) if valuation is 3 and
612 : * factor(abs(D / 4)) otherwise */
613 : static GEN
614 2275 : qfbsqrt(GEN D, GEN md, GEN q, GEN P)
615 : {
616 2275 : GEN a = gel(q,1), b = shifti(gel(q,2),-1), c = gel(q,3), B = negi(b);
617 2275 : GEN m, n, Q, M, N, d = negi(md); /* ac - b^2 */
618 2275 : long i, lP = lg(P);
619 :
620 : /* 1) solve m^2 = a, m*n = -b, n^2 = c in Z/dZ => q(n,m) = 0 mod d */
621 2275 : M = cgetg(lP, t_VEC);
622 2275 : N = cgetg(lP, t_VEC);
623 9716 : for (i = 1; i < lg(P); i++)
624 : {
625 7441 : GEN p = gel(P,i);
626 7441 : if (dvdii(a,p)) { n = Fp_sqrt(c, p); m = Fp_div(B, n, p); }
627 5628 : else { m = Fp_sqrt(a, p); n = Fp_div(B, m, p); }
628 7441 : gel(M, i) = m;
629 7441 : gel(N, i) = n;
630 : }
631 2275 : m = ZV_chinese_center(M, P, NULL);
632 2275 : n = ZV_chinese_center(N, P, NULL);
633 :
634 : /* 2) build Q, with det=-1 such that Q(x,y,0) = G(x,y) */
635 2275 : N = diviiexact(addii(mulii(a,n), mulii(b,m)), d);
636 2275 : M = diviiexact(addii(mulii(b,n), mulii(c,m)), d);
637 2275 : Q = diviiexact(subiu(addii(mulii(m,M), mulii(n,N)), 1), d); /*(q(n,m)-d)/d^2 */
638 2275 : Q = mkmat3(mkcol3(a,b,N), mkcol3(b,c,M), mkcol3(N,M,Q)); /* det = -1 */
639 :
640 : /* 3) reduce Q to [0,0,-1; 0,1,0; -1,0,0] */
641 2275 : M = qflllgram_indefgoon2(Q);
642 2275 : if (signe(gel(M,1)) < 0) M = ZC_neg(M);
643 2275 : a = gel(M,1);
644 2275 : b = gel(M,2);
645 2275 : c = gel(M,3);
646 2275 : if (!mpodd(a)) { swap(a, c); togglesign_safe(&b); }
647 2275 : return mkqfb(a, shifti(b,1), shifti(c,1), D);
648 : }
649 :
650 : /* \prod gen[i]^e[i] as a Qfb, e in {0,1}^n nonzero */
651 : static GEN
652 3962 : qfb_factorback(GEN gen, GEN e, GEN isqrtD)
653 : {
654 3962 : GEN q = NULL;
655 3962 : long j, l = lg(gen), n = 0;
656 13496 : for (j = 1; j < l; j++)
657 9534 : if (e[j]) { n++; q = q? qfbcompraw(q, gel(gen,j)): gel(gen,j); }
658 3962 : return (n <= 1)? q: qfbred0(q, 0, isqrtD, NULL);
659 : }
660 :
661 : /* unit form discriminant 4d */
662 : static GEN
663 1001 : id(GEN d)
664 1001 : { retmkmat22(gen_1, gen_0, gen_0, negi(d)); }
665 :
666 : /* Shanks/Bosma-Stevenhagen algorithm to compute the 2-Sylow of the class
667 : * group of discriminant D. Only works for D = fundamental discriminant.
668 : * When D = 1(4), work with 4D.
669 : * P2D,E2D = factor(abs(2*D))
670 : * Pm2D = factor(-abs(2*D))[,1].
671 : * Return a form having Witt invariants W at Pm2D */
672 : static GEN
673 2688 : quadclass2(GEN D, GEN P2D, GEN E2D, GEN Pm2D, GEN W, int n_is_4)
674 : {
675 2688 : GEN U2 = NULL, gen, Wgen, isqrtD, d;
676 : long i, r, m;
677 2688 : int splice2 = mpodd(D);
678 :
679 2688 : if (!splice2) d = shifti(D,-2); else { d = D; D = shifti(D,2); }
680 : /* D = 4d */
681 2688 : if (zv_equal0(W)) return id(d);
682 1806 : r = lg(Pm2D) - 4; /* >= 0 since W != 0 */
683 1806 : m = (signe(D) > 0)? r+1: r;
684 1806 : if (n_is_4)
685 : { /* n = 4: look among forms of type q or 2*q, since Q can be imprimitive */
686 539 : U2 = hilberts(gen_2, d, Pm2D);
687 539 : if (zv_equal(U2,W)) return gmul2n(id(d),1);
688 : }
689 :
690 1687 : gen = cgetg(m+2, t_VEC);
691 4571 : for (i = 1; i <= m; i++)
692 : { /* no need to look at P2D[1]=2*/
693 2884 : GEN q = powiu(gel(P2D,i+1), E2D[i+1]);
694 2884 : gel(gen,i) = mkqfb(q, gen_0, negi(diviiexact(d,q)), D);
695 : }
696 1687 : if (!mpodd(d))
697 : {
698 1463 : gel(gen, ++m) = mkqfb(gen_2, gen_0, negi(shifti(d,-1)), D);
699 1463 : r++;
700 : }
701 224 : else if (Mod4(d) != 1)
702 : {
703 119 : gel(gen, ++m) = mkqfb(gen_2, gen_2, shifti(subsi(1,d),-1), D);
704 119 : r++;
705 : }
706 105 : else setlg(gen, m+1);
707 1687 : if (!r) return id(d);
708 : /* remove 2^3; leave alone 2^4 */
709 1687 : if (splice2) P2D = vecsplice(P2D, 1);
710 1687 : Wgen = cgetg(m+1, t_MAT);
711 6153 : for (i = 1; i <= m; i++) gel(Wgen,i) = hilberts(gmael(gen,i,1), d, Pm2D);
712 1687 : isqrtD = signe(D) > 0? sqrti(D) : NULL;
713 : for(;;)
714 2065 : {
715 3752 : GEN Wgen2, gen2, Ker, indexim = gel(Flm_indexrank(Wgen,2), 2);
716 : long dKer;
717 3752 : if (lg(indexim)-1 >= r)
718 : {
719 1687 : GEN W2 = Wgen, V;
720 1687 : if (lg(indexim) < lg(Wgen)) W2 = vecpermute(Wgen,indexim);
721 1687 : if (U2) W2 = vec_append(W2,U2);
722 1687 : V = Flm_Flc_invimage(W2, W, 2);
723 1687 : if (V)
724 : {
725 1687 : GEN Q = qfb_factorback(vecpermute(gen,indexim), V, isqrtD);
726 1687 : return (U2 && V[lg(V)-1])? qfbmat2(Q): qfbmat(Q);
727 : }
728 : }
729 2065 : Ker = Flm_ker(Wgen,2); dKer = lg(Ker)-1;
730 2065 : gen2 = cgetg(m+1, t_VEC);
731 2065 : Wgen2 = cgetg(m+1, t_MAT);
732 4340 : for (i = 1; i <= dKer; i++)
733 : {
734 2275 : GEN q = qfb_factorback(gen, gel(Ker,i), isqrtD);
735 2275 : gel(gen2,i) = q = qfbsqrt(D, d, q, P2D);
736 2275 : gel(Wgen2,i) = hilberts(gel(q,1), d, Pm2D);
737 : }
738 4445 : for (; i <=m; i++)
739 : {
740 2380 : long j = indexim[i-dKer];
741 2380 : gel(gen2,i) = gel(gen,j);
742 2380 : gel(Wgen2,i) = gel(Wgen,j);
743 : }
744 2065 : gen = gen2; Wgen = Wgen2;
745 : }
746 : }
747 :
748 : /* QUADRATIC EQUATIONS */
749 : /* is x*y = -1 ? */
750 : static int
751 19674 : both_pm1(GEN x, GEN y)
752 19674 : { return is_pm1(x) && is_pm1(y) && signe(x) == -signe(y); }
753 :
754 : /* Try to solve G = 0 with small coefficients. This is proved to work if
755 : * - det(G) = 1, dim <= 6 and G is LLL reduced
756 : * Returns G if no solution is found.
757 : * Exit with a norm 0 vector if one such is found.
758 : * If base == 1 and norm 0 is obtained, returns [H~*G*H,H] where
759 : * the 1st column of H is a norm 0 vector */
760 : static GEN
761 44881 : qftriv(GEN G, GEN R, long base)
762 : {
763 44881 : long n = lg(G)-1, i;
764 : GEN s, H;
765 :
766 : /* case 1: A basis vector is isotropic */
767 142258 : for (i = 1; i <= n; i++)
768 116479 : if (!signe(gcoeff(G,i,i)))
769 : {
770 19102 : if (!base) return col_ei(n,i);
771 11049 : H = matid(n); swap(gel(H,1), gel(H,i));
772 11049 : return mkvec2(qf_apply_tau(G,1,i),H);
773 : }
774 : /* case 2: G has a block +- [1,0;0,-1] on the diagonal */
775 85378 : for (i = 2; i <= n; i++)
776 63389 : if (!signe(gcoeff(G,i-1,i)) && both_pm1(gcoeff(G,i-1,i-1),gcoeff(G,i,i)))
777 : {
778 3790 : s = col_ei(n,i); gel(s,i-1) = gen_m1;
779 3790 : if (!base) return s;
780 2995 : H = matid(n); gel(H,i) = gel(H,1); gel(H,1) = s;
781 2995 : return mkvec2(qf_ZM_apply(G,H),H);
782 : }
783 21989 : if (!R) return G; /* fail */
784 : /* case 3: a principal minor is 0 */
785 549 : s = ZM_ker(principal_minor(G, itos(R)));
786 549 : s = vecextend(Q_primpart(gel(s,1)), n);
787 549 : if (!base) return s;
788 263 : H = completebasis(s, 0);
789 263 : gel(H,n) = ZC_neg(gel(H,1)); gel(H,1) = s;
790 263 : return mkvec2(qf_ZM_apply(G,H),H);
791 : }
792 :
793 : /* p a prime number, G 3x3 symmetric. Finds X!=0 such that X^t G X = 0 mod p.
794 : * Allow returning a shorter X: to be completed with 0s. */
795 : static GEN
796 18039 : qfsolvemodp(GEN G, GEN p)
797 : {
798 : GEN a,b,c,d,e,f, v1,v2,v3,v4,v5, x1,x2,x3,N1,N2,N3,s,r;
799 :
800 : /* principal_minor(G,3) = [a,b,d; b,c,e; d,e,f] */
801 18039 : a = modii(gcoeff(G,1,1), p);
802 18039 : if (!signe(a)) return mkcol(gen_1);
803 15234 : v1 = a;
804 15234 : b = modii(gcoeff(G,1,2), p);
805 15234 : c = modii(gcoeff(G,2,2), p);
806 15234 : v2 = modii(subii(mulii(a,c), sqri(b)), p);
807 15234 : if (!signe(v2)) return mkcol2(Fp_neg(b,p), a);
808 12319 : d = modii(gcoeff(G,1,3), p);
809 12319 : e = modii(gcoeff(G,2,3), p);
810 12319 : f = modii(gcoeff(G,3,3), p);
811 12319 : v4 = modii(subii(mulii(c,d), mulii(e,b)), p);
812 12319 : v5 = modii(subii(mulii(a,e), mulii(d,b)), p);
813 12319 : v3 = subii(mulii(v2,f), addii(mulii(v4,d), mulii(v5,e))); /* det(G) */
814 12319 : v3 = modii(v3, p);
815 12319 : N1 = Fp_neg(v2, p);
816 12319 : x3 = mkcol3(v4, v5, N1);
817 12319 : if (!signe(v3)) return x3;
818 :
819 : /* now, solve in dimension 3... reduction to the diagonal case: */
820 10620 : x1 = mkcol3(gen_1, gen_0, gen_0);
821 10620 : x2 = mkcol3(negi(b), a, gen_0);
822 10620 : if (kronecker(N1,p) == 1) return ZC_lincomb(Fp_sqrt(N1,p),gen_1,x1,x2);
823 4256 : N2 = Fp_div(Fp_neg(v3,p), v1, p);
824 4256 : if (kronecker(N2,p) == 1) return ZC_lincomb(Fp_sqrt(N2,p),gen_1,x2,x3);
825 2129 : N3 = Fp_mul(v2, N2, p);
826 2129 : if (kronecker(N3,p) == 1) return ZC_lincomb(Fp_sqrt(N3,p),gen_1,x1,x3);
827 1065 : r = gen_1;
828 : for(;;)
829 : {
830 2046 : s = Fp_sub(gen_1, Fp_mul(N1,Fp_sqr(r,p),p), p);
831 2046 : if (kronecker(s, p) <= 0) break;
832 981 : r = randomi(p);
833 : }
834 1065 : s = Fp_sqrt(Fp_div(s,N3,p), p);
835 1065 : return ZC_add(x1, ZC_lincomb(r,s,x2,x3));
836 : }
837 :
838 : /* Given a square rational matrix G of dimension n >= 1, solves over Z the
839 : * quadratic equation X^tGX = 0. The solution is a t_VEC (a solution) or a
840 : * t_MAT (totally isotropic subspace). If no solution exists, returns an
841 : * integer: a prime p (no local solution at p), or -1 (no real solution), or
842 : * -2 (n = 2 and -deg(G) not a square). */
843 : static GEN
844 12718 : qfsolve_i(GEN G)
845 : {
846 12718 : GEN M, signG, Min, U, G1, M1, G2, M2, solG2, P = NULL, E = NULL;
847 : GEN solG1, sol, Q, d, dQ, detG2;
848 : long n, np, codim, dim;
849 : int fail;
850 :
851 12718 : if (typ(G) == t_VEC && lg(G)==3)
852 : {
853 8420 : P = gel(G,2);
854 8420 : G = gel(G,1);
855 8420 : if (typ(P)==t_MAT)
856 : {
857 7 : if (!is_Z_factornon0(P)) pari_err_TYPE("qfsolve", P);
858 7 : P = gel(P,1);
859 : } else
860 8413 : if (!is_vec_t(typ(P)) || !RgV_is_ZVnon0(P))
861 0 : pari_err_TYPE("qfsolve", P);
862 : }
863 12718 : if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
864 12718 : n = lg(G)-1;
865 12718 : if (n == 0) pari_err_DOMAIN("qfsolve", "dimension", "=", gen_0, G);
866 12718 : if (n != nbrows(G)) pari_err_DIM("qfsolve");
867 12718 : G = Q_primpart(G); RgM_check_ZM(G, "qfsolve");
868 12718 : check_symmetric(G);
869 :
870 : /* Trivial case: det = 0 */
871 12711 : d = ZM_det(G);
872 12711 : if (!signe(d))
873 : {
874 7 : if (n == 1) return mkcol(gen_1);
875 0 : sol = ZM_ker(G);
876 0 : if (lg(sol) == 2) sol = gel(sol,1);
877 0 : return sol;
878 : }
879 :
880 : /* Small dimension: n <= 2 */
881 12704 : if (n == 1) return gen_m1;
882 12697 : if (n == 2)
883 : {
884 21 : GEN t, a = gcoeff(G,1,1);
885 21 : if (!signe(a)) return mkcol2(gen_1, gen_0);
886 14 : if (signe(d) > 0) return gen_m1; /* no real solution */
887 14 : if (!Z_issquareall(negi(d), &t)) return gen_m2;
888 7 : return mkcol2(subii(t,gcoeff(G,1,2)), a);
889 : }
890 :
891 : /* 1st reduction of the coefficients of G */
892 12676 : M = qflllgram_indef(G,0,&fail);
893 12676 : if (typ(M) == t_COL) return M;
894 12165 : G = gel(M,1);
895 12165 : M = gel(M,2);
896 :
897 : /* real solubility */
898 12165 : signG = ZV_to_zv(qfsign(G));
899 : {
900 12165 : long r = signG[1], s = signG[2];
901 12165 : if (!r || !s) return gen_m1;
902 12095 : if (r < s) { G = ZM_neg(G); signG = mkvecsmall2(s,r); }
903 : }
904 :
905 : /* factorization of the determinant */
906 12095 : if (!P)
907 : {
908 4172 : GEN F = absZ_factor(d);
909 4172 : P = gel(F,1); E = ZV_to_zv(gel(F,2));
910 : }
911 :
912 : /* Minimization and local solubility */
913 12095 : Min = qfminimize_fact(G, P, E, d, 1);
914 12095 : if (typ(Min) == t_INT) return Min;
915 :
916 11612 : M = QM_mul(M, gel(Min,2));
917 11612 : G = gel(Min,1);
918 11612 : P = gel(Min,3);
919 11612 : E = gel(Min,4);
920 : /* P,E = factor(|det(G))| */
921 :
922 : /* Now, we know that local solutions exist (except maybe at 2 if n==4)
923 : * if n==3, det(G) = +-1
924 : * if n==4, or n is odd, det(G) is squarefree.
925 : * if n>=6, det(G) has all its valuations <=2. */
926 :
927 : /* Reduction of G and search for trivial solutions. */
928 : /* When |det G|=1, such trivial solutions always exist. */
929 11612 : U = qflllgram_indef(G,0,&fail);
930 11612 : if (typ(U) == t_COL) return Q_primpart(RgM_RgC_mul(M,U));
931 2989 : G = gel(U,1);
932 2989 : M = QM_mul(M, gel(U,2));
933 : /* P,E = factor(|det(G))| */
934 :
935 : /* If n >= 6 is even, need to increment the dimension by 1 to suppress all
936 : * squares from det(G) */
937 2989 : np = lg(P)-1;
938 2989 : if (n < 6 || odd(n) || !np)
939 : {
940 1603 : codim = 0;
941 1603 : G1 = G;
942 1603 : M1 = NULL;
943 : }
944 : else
945 : {
946 : GEN aux;
947 : long i;
948 1386 : codim = 1; n++;
949 1386 : aux = gen_1;
950 6524 : for (i = 1; i <= np; i++)
951 5138 : if (E[i] == 2) { aux = mulii(aux, gel(P,i)); E[i] = 3; }
952 : /* aux = largest square divisor of d; choose sign(aux) so as to balance
953 : * the signature of G1 */
954 1386 : if (signG[1] > signG[2])
955 : {
956 546 : signG[2]++;
957 546 : aux = negi(aux);
958 : }
959 : else
960 840 : signG[1]++;
961 1386 : G1 = shallowmatconcat(diagonal_shallow(mkvec2(G,aux)));
962 : /* P,E = factor(|det G1|) */
963 1386 : Min = qfminimize_fact(G1, P, E, NULL, 1);
964 1386 : G1 = gel(Min,1);
965 1386 : M1 = gel(Min,2);
966 1386 : P = gel(Min,3);
967 1386 : E = gel(Min,4);
968 1386 : np = lg(P)-1;
969 : }
970 :
971 : /* now, d is squarefree */
972 2989 : if (!np) { G2 = G1; M2 = NULL; } /* |d| = 1 */
973 : else
974 : { /* |d| > 1: increment dimension by 2 */
975 2723 : GEN factdP, factdE, W, v = NULL;
976 : long i, lfactdP, lP;
977 2723 : codim += 2;
978 2723 : d = ZV_prod(P); /* d = abs(matdet(G1)); */
979 2723 : if (odd(signG[2])) togglesign_safe(&d); /* d = matdet(G1); */
980 2723 : if (n == 4)
981 : { /* solubility at 2, the only remaining bad prime */
982 707 : v = det_minors(G1);
983 707 : if (Mod8(d) == 1 && witt(v, gen_2) == 1) return gen_2;
984 : }
985 2688 : P = shallowconcat(mpodd(d)? mkvec2(NULL,gen_2): mkvec(NULL), P);
986 : /* build a binary quadratic form with given Witt invariants */
987 2688 : lP = lg(P); W = const_vecsmall(lP-1, 0);
988 : /* choose signature of Q (real invariant and sign of the discriminant) */
989 2688 : dQ = absi(d);
990 2688 : if (signG[1] > signG[2]) togglesign_safe(&dQ); /* signQ = [2,0]; */
991 2688 : if (n == 4 && Mod4(dQ) != 1) dQ = shifti(dQ,2);
992 2688 : if (n >= 5) dQ = shifti(dQ,3);
993 :
994 : /* p-adic invariants */
995 2688 : if (n == 4)
996 : {
997 672 : togglesign(gel(v,2));
998 672 : togglesign(gel(v,4)); /* v = det_minors(-G1) */
999 2681 : for (i = 3; i < lP; i++) W[i] = witt(v, gel(P,i)) < 0;
1000 : }
1001 : else
1002 : {
1003 2016 : long s = signe(dQ) == signe(d)? 1: -1;
1004 : GEN t;
1005 2016 : if (odd((n-3)/2)) s = -s;
1006 2016 : t = s > 0? utoipos(8): utoineg(8);
1007 6055 : for (i = 3; i < lP; i++) W[i] = hilbertii(t, gel(P,i), gel(P,i)) > 0;
1008 : }
1009 2688 : W[2] = Flv_sum(W, 2); /* for p = 2, use product formula */
1010 :
1011 : /* Construction of the 2-class group of discriminant dQ until some product
1012 : * of the generators gives the desired invariants. */
1013 2688 : factdP = vecsplice(P, 1); lfactdP = lg(factdP);
1014 2688 : factdE = cgetg(lfactdP, t_VECSMALL);
1015 11424 : for (i = 1; i < lfactdP; i++) factdE[i] = Z_pval(dQ, gel(factdP,i));
1016 2688 : factdE[1]++;
1017 : /* factdP,factdE = factor(2|dQ|), P = factor(-2|dQ|)[,1] */
1018 2688 : Q = quadclass2(dQ, factdP,factdE, P, W, n == 4);
1019 : /* Build a form of dim=n+2 potentially unimodular */
1020 2688 : G2 = shallowmatconcat(diagonal_shallow(mkvec2(G1,ZM_neg(Q))));
1021 : /* Minimization of G2 */
1022 2688 : detG2 = mulii(d, ZM_det(Q));
1023 : /* factdP,factdE = factor(|det G2|) */
1024 2688 : Min = qfminimize_fact(G2, factdP, NULL, detG2, 1);
1025 2688 : M2 = gel(Min,2);
1026 2688 : G2 = gel(Min,1);
1027 : }
1028 : /* |det(G2)| = 1, find a totally isotropic subspace for G2 */
1029 2954 : solG2 = qflllgram_indefgoon(G2);
1030 : /* G2 must have a subspace of solutions of dimension > codim */
1031 2954 : dim = codim+1;
1032 7273 : while(gequal0(principal_minor(gel(solG2,1), dim))) dim ++;
1033 2954 : if (dim == codim+1)
1034 14 : pari_err_IMPL("qfsolve, dim >= 10");
1035 2940 : solG2 = vecslice(gel(solG2,2), 1, dim-1);
1036 :
1037 2940 : if (!M2) solG1 = solG2;
1038 : else
1039 : { /* solution of G1 is simultaneously in solG2 and x[n+1] = x[n+2] = 0*/
1040 : GEN K;
1041 2681 : solG1 = QM_mul(M2,solG2);
1042 2681 : K = QM_ker(rowslice(solG1,n+1,n+2));
1043 2681 : solG1 = QM_mul(rowslice(solG1,1,n), K);
1044 : }
1045 2940 : if (!M1) sol = solG1;
1046 : else
1047 : { /* solution of G1 is simultaneously in solG2 and x[n] = 0 */
1048 : GEN K;
1049 1386 : sol = QM_mul(M1,solG1);
1050 1386 : K = QM_ker(rowslice(sol,n,n));
1051 1386 : sol = QM_mul(rowslice(sol,1,n-1), K);
1052 : }
1053 2940 : sol = Q_primpart(QM_mul(M, sol));
1054 2940 : if (lg(sol) == 2) sol = gel(sol,1);
1055 2940 : return sol;
1056 : }
1057 : GEN
1058 12718 : qfsolve(GEN G)
1059 12718 : { pari_sp av = avma; return gc_GEN(av, qfsolve_i(G)); }
1060 :
1061 : /* G is a symmetric 3x3 matrix, and sol a solution of sol~*G*sol=0.
1062 : * Returns a parametrization of the solutions with the good invariants,
1063 : * as a 3x3 matrix, where each line contains the coefficients of each of the 3
1064 : * quadratic forms. If fl!=0, the fl-th form is reduced. */
1065 : GEN
1066 7244 : qfparam(GEN G, GEN sol, long fl)
1067 : {
1068 7244 : pari_sp av = avma;
1069 : GEN U, G1, G2, a, b, c, d, e;
1070 7244 : long n, tx = typ(sol);
1071 :
1072 7244 : if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
1073 7244 : if (!is_vec_t(tx)) pari_err_TYPE("qfsolve", G);
1074 7244 : if (tx == t_VEC) sol = shallowtrans(sol);
1075 7244 : n = lg(G)-1;
1076 7244 : if (n == 0) pari_err_DOMAIN("qfsolve", "dimension", "=", gen_0, G);
1077 7244 : if (n != nbrows(G) || n != 3 || lg(sol) != 4) pari_err_DIM("qfsolve");
1078 7244 : G = Q_primpart(G); RgM_check_ZM(G,"qfsolve");
1079 7244 : check_symmetric(G);
1080 7244 : sol = Q_primpart(sol); RgV_check_ZV(sol,"qfsolve");
1081 : /* build U such that U[,3] = sol, and |det(U)| = 1 */
1082 7244 : U = completebasis(sol,1);
1083 7244 : G1 = qf_ZM_apply(G,U); /* G1 has a 0 at the bottom right corner */
1084 7244 : a = shifti(gcoeff(G1,1,2),1);
1085 7244 : b = shifti(negi(gcoeff(G1,1,3)),1);
1086 7244 : c = shifti(negi(gcoeff(G1,2,3)),1);
1087 7244 : d = gcoeff(G1,1,1);
1088 7244 : e = gcoeff(G1,2,2);
1089 7244 : G2 = mkmat3(mkcol3(b,gen_0,d), mkcol3(c,b,a), mkcol3(gen_0,c,e));
1090 7244 : sol = ZM_mul(U,G2);
1091 7244 : if (fl)
1092 : {
1093 7223 : GEN v = row(sol,fl);
1094 : int fail;
1095 7223 : a = gel(v,1);
1096 7223 : b = gmul2n(gel(v,2),-1);
1097 7223 : c = gel(v,3);
1098 7223 : U = qflllgram_indef(mkmat22(a, b, b, c), 1, &fail);
1099 7223 : U = gel(U,2);
1100 7223 : a = gcoeff(U,1,1); b = gcoeff(U,1,2);
1101 7223 : c = gcoeff(U,2,1); d = gcoeff(U,2,2);
1102 7223 : U = mkmat3(mkcol3(sqri(a),mulii(a,c),sqri(c)),
1103 : mkcol3(shifti(mulii(a,b),1), addii(mulii(a,d),mulii(b,c)),
1104 : shifti(mulii(c,d),1)),
1105 : mkcol3(sqri(b),mulii(b,d),sqri(d)));
1106 7223 : sol = ZM_mul(sol,U);
1107 : }
1108 7244 : return gc_upto(av, sol);
1109 : }
|
