1 /* Copyright (C) 2000 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 #include "pari.h"
15 #include "paripriv.h"
16
17 #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
18
19 /********************************************************************/
20 /** **/
21 /** ASSOCIATIVE ALGEBRAS, CENTRAL SIMPLE ALGEBRAS **/
22 /** contributed by Aurel Page (2014) **/
23 /** **/
24 /********************************************************************/
25 static GEN alg_subalg(GEN al, GEN basis);
26 static GEN alg_maximal_primes(GEN al, GEN P);
27 static GEN algnatmultable(GEN al, long D);
28 static GEN _tablemul_ej(GEN mt, GEN x, long j);
29 static GEN _tablemul_ej_Fp(GEN mt, GEN x, long j, GEN p);
30 static GEN _tablemul_ej_Fl(GEN mt, GEN x, long j, ulong p);
31 static ulong algtracei(GEN mt, ulong p, ulong expo, ulong modu);
32 static GEN alg_pmaximal(GEN al, GEN p);
33 static GEN alg_maximal(GEN al);
34 static GEN algtracematrix(GEN al);
35 static GEN algtableinit_i(GEN mt0, GEN p);
36 static GEN algbasisrightmultable(GEN al, GEN x);
37 static GEN algabstrace(GEN al, GEN x);
38 static GEN algbasismul(GEN al, GEN x, GEN y);
39 static GEN algbasismultable(GEN al, GEN x);
40 static GEN algbasismultable_Flm(GEN mt, GEN x, ulong m);
41
42 static int
checkalg_i(GEN al)43 checkalg_i(GEN al)
44 {
45 GEN mt, rnf;
46 if (typ(al) != t_VEC || lg(al) != 12) return 0;
47 mt = alg_get_multable(al);
48 if (typ(mt) != t_VEC || lg(mt) == 1 || typ(gel(mt,1)) != t_MAT) return 0;
49 rnf = alg_get_splittingfield(al);
50 if (isintzero(rnf) || !gequal0(alg_get_char(al))) return 1;
51 if (typ(gel(al,2)) != t_VEC || lg(gel(al,2)) == 1) return 0;
52 /* not checkrnf_i: beware placeholder from alg_csa_table */
53 return typ(rnf)==t_VEC && lg(rnf)==13;
54 }
55 void
checkalg(GEN al)56 checkalg(GEN al)
57 { if (!checkalg_i(al)) pari_err_TYPE("checkalg [please apply alginit()]",al); }
58
59 static int
checklat_i(GEN al,GEN lat)60 checklat_i(GEN al, GEN lat)
61 {
62 long N,i,j;
63 GEN m,t,c;
64 if (typ(lat)!=t_VEC || lg(lat) != 3) return 0;
65 t = gel(lat,2);
66 if (typ(t) != t_INT && typ(t) != t_FRAC) return 0;
67 if (gsigne(t)<=0) return 0;
68 m = gel(lat,1);
69 if (typ(m) != t_MAT) return 0;
70 N = alg_get_absdim(al);
71 if (lg(m)-1 != N || lg(gel(m,1))-1 != N) return 0;
72 for (i=1; i<=N; i++)
73 for (j=1; j<=N; j++) {
74 c = gcoeff(m,i,j);
75 if (typ(c) != t_INT) return 0;
76 if (j<i && signe(gcoeff(m,i,j))) return 0;
77 if (i==j && !signe(gcoeff(m,i,j))) return 0;
78 }
79 return 1;
80 }
checklat(GEN al,GEN lat)81 void checklat(GEN al, GEN lat)
82 { if (!checklat_i(al,lat)) pari_err_TYPE("checklat [please apply alglathnf()]", lat); }
83
84 /** ACCESSORS **/
85 long
alg_type(GEN al)86 alg_type(GEN al)
87 {
88 if (isintzero(alg_get_splittingfield(al)) || !gequal0(alg_get_char(al))) return al_TABLE;
89 switch(typ(gmael(al,2,1))) {
90 case t_MAT: return al_CSA;
91 case t_INT:
92 case t_FRAC:
93 case t_POL:
94 case t_POLMOD: return al_CYCLIC;
95 default: return al_NULL;
96 }
97 return -1; /*LCOV_EXCL_LINE*/
98 }
99 long
algtype(GEN al)100 algtype(GEN al)
101 { return checkalg_i(al)? alg_type(al): al_NULL; }
102
103 /* absdim == dim for al_TABLE. */
104 long
alg_get_dim(GEN al)105 alg_get_dim(GEN al)
106 {
107 long d;
108 switch(alg_type(al)) {
109 case al_TABLE: return lg(alg_get_multable(al))-1;
110 case al_CSA: return lg(alg_get_relmultable(al))-1;
111 case al_CYCLIC: d = alg_get_degree(al); return d*d;
112 default: pari_err_TYPE("alg_get_dim", al);
113 }
114 return -1; /*LCOV_EXCL_LINE*/
115 }
116
117 long
alg_get_absdim(GEN al)118 alg_get_absdim(GEN al)
119 {
120 switch(alg_type(al)) {
121 case al_TABLE: return lg(alg_get_multable(al))-1;
122 case al_CSA: return alg_get_dim(al)*nf_get_degree(alg_get_center(al));
123 case al_CYCLIC:
124 return rnf_get_absdegree(alg_get_splittingfield(al))*alg_get_degree(al);
125 default: pari_err_TYPE("alg_get_absdim", al);
126 }
127 return -1;/*LCOV_EXCL_LINE*/
128 }
129
130 long
algdim(GEN al,long abs)131 algdim(GEN al, long abs)
132 {
133 checkalg(al);
134 if (abs) return alg_get_absdim(al);
135 return alg_get_dim(al);
136 }
137
138 /* only cyclic */
139 GEN
alg_get_auts(GEN al)140 alg_get_auts(GEN al)
141 {
142 if (alg_type(al) != al_CYCLIC)
143 pari_err_TYPE("alg_get_auts [noncyclic algebra]", al);
144 return gel(al,2);
145 }
146 GEN
alg_get_aut(GEN al)147 alg_get_aut(GEN al)
148 {
149 if (alg_type(al) != al_CYCLIC)
150 pari_err_TYPE("alg_get_aut [noncyclic algebra]", al);
151 return gel(alg_get_auts(al),1);
152 }
153 GEN
algaut(GEN al)154 algaut(GEN al) { checkalg(al); return alg_get_aut(al); }
155 GEN
alg_get_b(GEN al)156 alg_get_b(GEN al)
157 {
158 if (alg_type(al) != al_CYCLIC)
159 pari_err_TYPE("alg_get_b [noncyclic algebra]", al);
160 return gel(al,3);
161 }
162 GEN
algb(GEN al)163 algb(GEN al) { checkalg(al); return alg_get_b(al); }
164
165 /* only CSA */
166 GEN
alg_get_relmultable(GEN al)167 alg_get_relmultable(GEN al)
168 {
169 if (alg_type(al) != al_CSA)
170 pari_err_TYPE("alg_get_relmultable [algebra not given via mult. table]", al);
171 return gel(al,2);
172 }
173 GEN
algrelmultable(GEN al)174 algrelmultable(GEN al) { checkalg(al); return alg_get_relmultable(al); }
175 GEN
alg_get_splittingdata(GEN al)176 alg_get_splittingdata(GEN al)
177 {
178 if (alg_type(al) != al_CSA)
179 pari_err_TYPE("alg_get_splittingdata [algebra not given via mult. table]",al);
180 return gel(al,3);
181 }
182 GEN
algsplittingdata(GEN al)183 algsplittingdata(GEN al) { checkalg(al); return alg_get_splittingdata(al); }
184 GEN
alg_get_splittingbasis(GEN al)185 alg_get_splittingbasis(GEN al)
186 {
187 if (alg_type(al) != al_CSA)
188 pari_err_TYPE("alg_get_splittingbasis [algebra not given via mult. table]",al);
189 return gmael(al,3,2);
190 }
191 GEN
alg_get_splittingbasisinv(GEN al)192 alg_get_splittingbasisinv(GEN al)
193 {
194 if (alg_type(al) != al_CSA)
195 pari_err_TYPE("alg_get_splittingbasisinv [algebra not given via mult. table]",al);
196 return gmael(al,3,3);
197 }
198
199 /* only cyclic and CSA */
200 GEN
alg_get_splittingfield(GEN al)201 alg_get_splittingfield(GEN al) { return gel(al,1); }
202 GEN
algsplittingfield(GEN al)203 algsplittingfield(GEN al)
204 {
205 long ta;
206 checkalg(al);
207 ta = alg_type(al);
208 if (ta != al_CYCLIC && ta != al_CSA)
209 pari_err_TYPE("alg_get_splittingfield [use alginit]",al);
210 return alg_get_splittingfield(al);
211 }
212 long
alg_get_degree(GEN al)213 alg_get_degree(GEN al)
214 {
215 long ta;
216 ta = alg_type(al);
217 if (ta != al_CYCLIC && ta != al_CSA)
218 pari_err_TYPE("alg_get_degree [use alginit]",al);
219 return rnf_get_degree(alg_get_splittingfield(al));
220 }
221 long
algdegree(GEN al)222 algdegree(GEN al)
223 {
224 checkalg(al);
225 return alg_get_degree(al);
226 }
227
228 GEN
alg_get_center(GEN al)229 alg_get_center(GEN al)
230 {
231 long ta;
232 ta = alg_type(al);
233 if (ta != al_CSA && ta != al_CYCLIC)
234 pari_err_TYPE("alg_get_center [use alginit]",al);
235 return rnf_get_nf(alg_get_splittingfield(al));
236 }
237 GEN
alg_get_splitpol(GEN al)238 alg_get_splitpol(GEN al)
239 {
240 long ta = alg_type(al);
241 if (ta != al_CYCLIC && ta != al_CSA)
242 pari_err_TYPE("alg_get_splitpol [use alginit]",al);
243 return rnf_get_pol(alg_get_splittingfield(al));
244 }
245 GEN
alg_get_abssplitting(GEN al)246 alg_get_abssplitting(GEN al)
247 {
248 long ta = alg_type(al), prec;
249 if (ta != al_CYCLIC && ta != al_CSA)
250 pari_err_TYPE("alg_get_abssplitting [use alginit]",al);
251 prec = nf_get_prec(alg_get_center(al));
252 return rnf_build_nfabs(alg_get_splittingfield(al), prec);
253 }
254 GEN
alg_get_hasse_i(GEN al)255 alg_get_hasse_i(GEN al)
256 {
257 long ta = alg_type(al);
258 if (ta != al_CYCLIC && ta != al_CSA)
259 pari_err_TYPE("alg_get_hasse_i [use alginit]",al);
260 if (ta == al_CSA) pari_err_IMPL("computation of Hasse invariants over table CSA");
261 return gel(al,4);
262 }
263 GEN
alghassei(GEN al)264 alghassei(GEN al) { checkalg(al); return alg_get_hasse_i(al); }
265 GEN
alg_get_hasse_f(GEN al)266 alg_get_hasse_f(GEN al)
267 {
268 long ta = alg_type(al);
269 if (ta != al_CYCLIC && ta != al_CSA)
270 pari_err_TYPE("alg_get_hasse_f [use alginit]",al);
271 if (ta == al_CSA) pari_err_IMPL("computation of Hasse invariants over table CSA");
272 return gel(al,5);
273 }
274 GEN
alghassef(GEN al)275 alghassef(GEN al) { checkalg(al); return alg_get_hasse_f(al); }
276
277 /* all types */
278 GEN
alg_get_basis(GEN al)279 alg_get_basis(GEN al) { return gel(al,7); }
280 GEN
algbasis(GEN al)281 algbasis(GEN al) { checkalg(al); return alg_get_basis(al); }
282 GEN
alg_get_invbasis(GEN al)283 alg_get_invbasis(GEN al) { return gel(al,8); }
284 GEN
alginvbasis(GEN al)285 alginvbasis(GEN al) { checkalg(al); return alg_get_invbasis(al); }
286 GEN
alg_get_multable(GEN al)287 alg_get_multable(GEN al) { return gel(al,9); }
288 GEN
algmultable(GEN al)289 algmultable(GEN al) { checkalg(al); return alg_get_multable(al); }
290 GEN
alg_get_char(GEN al)291 alg_get_char(GEN al) { return gel(al,10); }
292 GEN
algchar(GEN al)293 algchar(GEN al) { checkalg(al); return alg_get_char(al); }
294 GEN
alg_get_tracebasis(GEN al)295 alg_get_tracebasis(GEN al) { return gel(al,11); }
296
297 /* lattices */
298 GEN
alglat_get_primbasis(GEN lat)299 alglat_get_primbasis(GEN lat) { return gel(lat,1); }
300 GEN
alglat_get_scalar(GEN lat)301 alglat_get_scalar(GEN lat) { return gel(lat,2); }
302
303 /** ADDITIONAL **/
304
305 /* no garbage collection */
306 static GEN
backtrackfacto(GEN y0,long n,GEN red,GEN pl,GEN nf,GEN data,int (* test)(GEN,GEN),GEN * fa,GEN N,GEN I)307 backtrackfacto(GEN y0, long n, GEN red, GEN pl, GEN nf, GEN data, int (*test)(GEN,GEN), GEN* fa, GEN N, GEN I)
308 {
309 long b, i;
310 ulong lim = 1UL << 17;
311 long *v = new_chunk(n+1);
312 pari_sp av = avma;
313 for (b = 0;; b += (2*b)/(3*n) + 1)
314 {
315 GEN ny, y1, y2;
316 set_avma(av);
317 for (i = 1; i <= n; i++) v[i] = -b;
318 v[n]--;
319 for(;;)
320 {
321 i = n;
322 while (i > 0)
323 { if (v[i] == b) v[i--] = -b; else { v[i]++; break; } }
324 if (i==0) break;
325
326 y1 = y0;
327 for (i = 1; i <= n; i++) y1 = nfadd(nf, y1, ZC_z_mul(gel(red,i), v[i]));
328 if (!nfchecksigns(nf, y1, pl)) continue;
329
330 ny = absi_shallow(nfnorm(nf, y1));
331 if (!signe(ny)) continue;
332 ny = diviiexact(ny, gcdii(ny, N));
333 if (!Z_issmooth(ny, lim)) continue;
334
335 y2 = idealdivexact(nf, y1, idealadd(nf,y1,I));
336 *fa = idealfactor(nf, y2);
337 if (!data || test(data,*fa)) return y1;
338 }
339 }
340 }
341
342 /* if data == NULL, the test is skipped */
343 /* in the test, the factorization does not contain the known factors */
344 static GEN
factoredextchinesetest(GEN nf,GEN x,GEN y,GEN pl,GEN * fa,GEN data,int (* test)(GEN,GEN))345 factoredextchinesetest(GEN nf, GEN x, GEN y, GEN pl, GEN* fa, GEN data, int (*test)(GEN,GEN))
346 {
347 pari_sp av = avma;
348 long n,i;
349 GEN x1, y0, y1, red, N, I, P = gel(x,1), E = gel(x,2);
350 n = nf_get_degree(nf);
351 x = idealchineseinit(nf, mkvec2(x,pl));
352 x1 = gel(x,1);
353 red = lg(x1) == 1? matid(n): gel(x1,1);
354 y0 = idealchinese(nf, x, y);
355
356 E = shallowcopy(E);
357 if (!gequal0(y0))
358 for (i=1; i<lg(E); i++)
359 {
360 long v = nfval(nf,y0,gel(P,i));
361 if (cmpsi(v, gel(E,i)) < 0) gel(E,i) = stoi(v);
362 }
363 /* N and I : known factors */
364 I = factorbackprime(nf, P, E);
365 N = idealnorm(nf,I);
366
367 y1 = backtrackfacto(y0, n, red, pl, nf, data, test, fa, N, I);
368
369 /* restore known factors */
370 for (i=1; i<lg(E); i++) gel(E,i) = stoi(nfval(nf,y1,gel(P,i)));
371 *fa = famat_reduce(famat_mul_shallow(*fa, mkmat2(P, E)));
372
373 gerepileall(av, 2, &y1, fa);
374 return y1;
375 }
376
377 static GEN
factoredextchinese(GEN nf,GEN x,GEN y,GEN pl,GEN * fa)378 factoredextchinese(GEN nf, GEN x, GEN y, GEN pl, GEN* fa)
379 { return factoredextchinesetest(nf,x,y,pl,fa,NULL,NULL); }
380
381 /** OPERATIONS ON ASSOCIATIVE ALGEBRAS algebras.c **/
382
383 /*
384 Convention:
385 (K/F,sigma,b) = sum_{i=0..n-1} u^i*K
386 t*u = u*sigma(t)
387
388 Natural basis:
389 1<=i<=d*n^2
390 b_i = u^((i-1)/(dn))*ZKabs.((i-1)%(dn)+1)
391
392 Integral basis:
393 Basis of some order.
394
395 al:
396 1- rnf of the cyclic splitting field of degree n over the center nf of degree d
397 2- VEC of aut^i 1<=i<=n
398 3- b in nf
399 4- infinite hasse invariants (mod n) : VECSMALL of size r1, values only 0 or n/2 (if integral)
400 5- finite hasse invariants (mod n) : VEC[VEC of primes, VECSMALL of hasse inv mod n]
401 6- nf of the splitting field (absolute)
402 7* dn^2*dn^2 matrix expressing the integral basis in terms of the natural basis
403 8* dn^2*dn^2 matrix expressing the natural basis in terms of the integral basis
404 9* VEC of dn^2 matrices giving the dn^2*dn^2 left multiplication tables of the integral basis
405 10* characteristic of the base field (used only for algebras given by a multiplication table)
406 11* trace of basis elements
407
408 If al is given by a multiplication table (al_TABLE), only the * fields are present.
409 */
410
411 /* assumes same center and same variable */
412 /* currently only works for coprime degrees */
413 GEN
algtensor(GEN al1,GEN al2,long maxord)414 algtensor(GEN al1, GEN al2, long maxord) {
415 pari_sp av = avma;
416 long v, k, d1, d2;
417 GEN nf, P1, P2, aut1, aut2, b1, b2, C, rnf, aut, b, x1, x2, al;
418
419 checkalg(al1);
420 checkalg(al2);
421 if (alg_type(al1) != al_CYCLIC || alg_type(al2) != al_CYCLIC)
422 pari_err_IMPL("tensor of noncyclic algebras"); /* TODO: do it. */
423
424 nf=alg_get_center(al1);
425 if (!gequal(alg_get_center(al2),nf))
426 pari_err_OP("tensor product [not the same center]", al1, al2);
427
428 P1=alg_get_splitpol(al1); aut1=alg_get_aut(al1); b1=alg_get_b(al1);
429 P2=alg_get_splitpol(al2); aut2=alg_get_aut(al2); b2=alg_get_b(al2);
430 v=varn(P1);
431
432 d1=alg_get_degree(al1);
433 d2=alg_get_degree(al2);
434 if (ugcd(d1,d2) != 1)
435 pari_err_IMPL("tensor of cylic algebras of noncoprime degrees"); /* TODO */
436
437 if (d1==1) return gcopy(al2);
438 if (d2==1) return gcopy(al1);
439
440 C = nfcompositum(nf, P1, P2, 3);
441 rnf = rnfinit(nf,gel(C,1));
442 x1 = gel(C,2);
443 x2 = gel(C,3);
444 k = itos(gel(C,4));
445 aut = gadd(gsubst(aut2,v,x2),gmulsg(k,gsubst(aut1,v,x1)));
446 b = nfmul(nf,nfpow_u(nf,b1,d2),nfpow_u(nf,b2,d1));
447 al = alg_cyclic(rnf,aut,b,maxord);
448 return gerepilecopy(av,al);
449 }
450
451 /* M an n x d Flm of rank d, n >= d. Initialize Mx = y solver */
452 static GEN
Flm_invimage_init(GEN M,ulong p)453 Flm_invimage_init(GEN M, ulong p)
454 {
455 GEN v = Flm_indexrank(M, p), perm = gel(v,1);
456 GEN MM = rowpermute(M, perm); /* square invertible */
457 return mkvec2(Flm_inv(MM,p), perm);
458 }
459 /* assume Mx = y has a solution, v = Flm_invimage_init(M,p); return x */
460 static GEN
Flm_invimage_pre(GEN v,GEN y,ulong p)461 Flm_invimage_pre(GEN v, GEN y, ulong p)
462 {
463 GEN inv = gel(v,1), perm = gel(v,2);
464 return Flm_Flc_mul(inv, vecsmallpermute(y, perm), p);
465 }
466
467 GEN
algradical(GEN al)468 algradical(GEN al)
469 {
470 pari_sp av = avma;
471 GEN I, x, traces, K, MT, P, mt;
472 long l,i,ni, n;
473 ulong modu, expo, p;
474 checkalg(al);
475 P = alg_get_char(al);
476 mt = alg_get_multable(al);
477 n = alg_get_absdim(al);
478 dbg_printf(1)("algradical: char=%Ps, dim=%d\n", P, n);
479 traces = algtracematrix(al);
480 if (!signe(P))
481 {
482 dbg_printf(2)(" char 0, computing kernel...\n");
483 K = ker(traces);
484 dbg_printf(2)(" ...done.\n");
485 ni = lg(K)-1; if (!ni) return gc_const(av, gen_0);
486 return gerepileupto(av, K);
487 }
488 dbg_printf(2)(" char>0, computing kernel...\n");
489 K = FpM_ker(traces, P);
490 dbg_printf(2)(" ...done.\n");
491 ni = lg(K)-1; if (!ni) return gc_const(av, gen_0);
492 if (abscmpiu(P,n)>0) return gerepileupto(av, K);
493
494 /* tough case, p <= n. Ronyai's algorithm */
495 p = P[2]; l = 1;
496 expo = p; modu = p*p;
497 dbg_printf(2)(" char>0, hard case.\n");
498 while (modu<=(ulong)n) { l++; modu *= p; }
499 MT = ZMV_to_FlmV(mt, modu);
500 I = ZM_to_Flm(K,p); /* I_0 */
501 for (i=1; i<=l; i++) {/*compute I_i, expo = p^i, modu = p^(l+1) > n*/
502 long j, lig,col;
503 GEN v = cgetg(ni+1, t_VECSMALL);
504 GEN invI = Flm_invimage_init(I, p);
505 dbg_printf(2)(" computing I_%d:\n", i);
506 traces = cgetg(ni+1,t_MAT);
507 for (j = 1; j <= ni; j++)
508 {
509 GEN M = algbasismultable_Flm(MT, gel(I,j), modu);
510 uel(v,j) = algtracei(M, p,expo,modu);
511 }
512 for (col=1; col<=ni; col++)
513 {
514 GEN t = cgetg(n+1,t_VECSMALL); gel(traces,col) = t;
515 x = gel(I, col); /*col-th basis vector of I_{i-1}*/
516 for (lig=1; lig<=n; lig++)
517 {
518 GEN y = _tablemul_ej_Fl(MT,x,lig,p);
519 GEN z = Flm_invimage_pre(invI, y, p);
520 uel(t,lig) = Flv_dotproduct(v, z, p);
521 }
522 }
523 dbg_printf(2)(" computing kernel...\n");
524 K = Flm_ker(traces, p);
525 dbg_printf(2)(" ...done.\n");
526 ni = lg(K)-1; if (!ni) return gc_const(av, gen_0);
527 I = Flm_mul(I,K,p);
528 expo *= p;
529 }
530 return Flm_to_ZM(I);
531 }
532
533 /* compute the multiplication table of the element x, where mt is a
534 * multiplication table in an arbitrary ring */
535 static GEN
Rgmultable(GEN mt,GEN x)536 Rgmultable(GEN mt, GEN x)
537 {
538 long i, l = lg(x);
539 GEN z = NULL;
540 for (i = 1; i < l; i++)
541 {
542 GEN c = gel(x,i);
543 if (!gequal0(c))
544 {
545 GEN M = RgM_Rg_mul(gel(mt,i),c);
546 z = z? RgM_add(z, M): M;
547 }
548 }
549 return z;
550 }
551
552 static GEN
change_Rgmultable(GEN mt,GEN P,GEN Pi)553 change_Rgmultable(GEN mt, GEN P, GEN Pi)
554 {
555 GEN mt2;
556 long lmt = lg(mt), i;
557 mt2 = cgetg(lmt,t_VEC);
558 for (i=1;i<lmt;i++) {
559 GEN mti = Rgmultable(mt,gel(P,i));
560 gel(mt2,i) = RgM_mul(Pi, RgM_mul(mti,P));
561 }
562 return mt2;
563 }
564
565 static GEN
alg_quotient0(GEN al,GEN S,GEN Si,long nq,GEN p,long maps)566 alg_quotient0(GEN al, GEN S, GEN Si, long nq, GEN p, long maps)
567 {
568 GEN mt = cgetg(nq+1,t_VEC), P, Pi, d;
569 long i;
570 dbg_printf(3)(" alg_quotient0: char=%Ps, dim=%d, dim I=%d\n", p, alg_get_absdim(al), lg(S)-1);
571 for (i=1; i<=nq; i++) {
572 GEN mti = algbasismultable(al,gel(S,i));
573 if (signe(p)) gel(mt,i) = FpM_mul(Si, FpM_mul(mti,S,p), p);
574 else gel(mt,i) = RgM_mul(Si, RgM_mul(mti,S));
575 }
576 if (!signe(p) && !isint1(Q_denom(mt))) {
577 dbg_printf(3)(" bad case: denominator=%Ps\n", Q_denom(mt));
578 P = Q_remove_denom(Si,&d);
579 P = ZM_hnf(P);
580 P = RgM_Rg_div(P,d);
581 Pi = RgM_inv(P);
582 mt = change_Rgmultable(mt,P,Pi);
583 Si = RgM_mul(P,Si);
584 S = RgM_mul(S,Pi);
585 }
586 al = algtableinit_i(mt,p);
587 if (maps) al = mkvec3(al,Si,S); /* algebra, proj, lift */
588 return al;
589 }
590
591 /* quotient of an algebra by a nontrivial two-sided ideal */
592 GEN
alg_quotient(GEN al,GEN I,long maps)593 alg_quotient(GEN al, GEN I, long maps)
594 {
595 pari_sp av = avma;
596 GEN p, IS, ISi, S, Si;
597 long n, ni;
598
599 checkalg(al);
600 p = alg_get_char(al);
601 n = alg_get_absdim(al);
602 ni = lg(I)-1;
603
604 /* force first vector of complement to be the identity */
605 IS = shallowconcat(I, gcoeff(alg_get_multable(al),1,1));
606 if (signe(p)) {
607 IS = FpM_suppl(IS,p);
608 ISi = FpM_inv(IS,p);
609 }
610 else {
611 IS = suppl(IS);
612 ISi = RgM_inv(IS);
613 }
614 S = vecslice(IS, ni+1, n);
615 Si = rowslice(ISi, ni+1, n);
616 return gerepilecopy(av, alg_quotient0(al, S, Si, n-ni, p, maps));
617 }
618
619 static GEN
image_keep_first(GEN m,GEN p)620 image_keep_first(GEN m, GEN p) /* assume first column is nonzero or m==0, no GC */
621 {
622 GEN ir, icol, irow, M, c, x;
623 long i;
624 if (gequal0(gel(m,1))) return zeromat(nbrows(m),0);
625
626 if (signe(p)) ir = FpM_indexrank(m,p);
627 else ir = indexrank(m);
628
629 icol = gel(ir,2);
630 if (icol[1]==1) return extract0(m,icol,NULL);
631
632 irow = gel(ir,1);
633 M = extract0(m, irow, icol);
634 c = extract0(gel(m,1), irow, NULL);
635 if (signe(p)) x = FpM_FpC_invimage(M,c,p);
636 else x = inverseimage(M,c); /* TODO modulo a small prime */
637
638 for (i=1; i<lg(x); i++)
639 {
640 if (!gequal0(gel(x,i)))
641 {
642 icol[i] = 1;
643 vecsmall_sort(icol);
644 return extract0(m,icol,NULL);
645 }
646 }
647
648 return NULL; /* LCOV_EXCL_LINE */
649 }
650
651 /* z[1],...z[nz] central elements such that z[1]A + z[2]A + ... + z[nz]A = A
652 * is a direct sum. idempotents ==> first basis element is identity */
653 GEN
alg_centralproj(GEN al,GEN z,long maps)654 alg_centralproj(GEN al, GEN z, long maps)
655 {
656 pari_sp av = avma;
657 GEN S, U, Ui, alq, p;
658 long i, iu, lz = lg(z);
659
660 checkalg(al);
661 if (typ(z) != t_VEC) pari_err_TYPE("alcentralproj",z);
662 p = alg_get_char(al);
663 dbg_printf(3)(" alg_centralproj: char=%Ps, dim=%d, #z=%d\n", p, alg_get_absdim(al), lz-1);
664 S = cgetg(lz,t_VEC); /* S[i] = Im(z_i) */
665 for (i=1; i<lz; i++)
666 {
667 GEN mti = algbasismultable(al, gel(z,i));
668 gel(S,i) = image_keep_first(mti,p);
669 }
670 U = shallowconcat1(S); /* U = [Im(z_1)|Im(z_2)|...|Im(z_nz)], n x n */
671 if (lg(U)-1 < alg_get_absdim(al)) pari_err_TYPE("alcentralproj [z[i]'s not surjective]",z);
672 if (signe(p)) Ui = FpM_inv(U,p);
673 else Ui = RgM_inv(U);
674 if (!Ui) pari_err_BUG("alcentralproj"); /*LCOV_EXCL_LINE*/
675
676 alq = cgetg(lz,t_VEC);
677 for (iu=0,i=1; i<lz; i++)
678 {
679 long nq = lg(gel(S,i))-1, ju = iu + nq;
680 GEN Si = rowslice(Ui, iu+1, ju);
681 gel(alq, i) = alg_quotient0(al,gel(S,i),Si,nq,p,maps);
682 iu = ju;
683 }
684 return gerepilecopy(av, alq);
685 }
686
687 /* al is an al_TABLE */
688 static GEN
algtablecenter(GEN al)689 algtablecenter(GEN al)
690 {
691 pari_sp av = avma;
692 long n, i, j, k, ic;
693 GEN C, cij, mt, p;
694
695 n = alg_get_absdim(al);
696 mt = alg_get_multable(al);
697 p = alg_get_char(al);
698 C = cgetg(n+1,t_MAT);
699 for (j=1; j<=n; j++)
700 {
701 gel(C,j) = cgetg(n*n-n+1,t_COL);
702 ic = 1;
703 for (i=2; i<=n; i++) {
704 if (signe(p)) cij = FpC_sub(gmael(mt,i,j),gmael(mt,j,i),p);
705 else cij = RgC_sub(gmael(mt,i,j),gmael(mt,j,i));
706 for (k=1; k<=n; k++, ic++) gcoeff(C,ic,j) = gel(cij, k);
707 }
708 }
709 if (signe(p)) return gerepileupto(av, FpM_ker(C,p));
710 else return gerepileupto(av, ker(C));
711 }
712
713 GEN
algcenter(GEN al)714 algcenter(GEN al)
715 {
716 checkalg(al);
717 if (alg_type(al)==al_TABLE) return algtablecenter(al);
718 return alg_get_center(al);
719 }
720
721 /* Only in positive characteristic. Assumes that al is semisimple. */
722 GEN
algprimesubalg(GEN al)723 algprimesubalg(GEN al)
724 {
725 pari_sp av = avma;
726 GEN p, Z, F, K;
727 long nz, i;
728 checkalg(al);
729 p = alg_get_char(al);
730 if (!signe(p)) pari_err_DOMAIN("algprimesubalg","characteristic","=",gen_0,p);
731
732 Z = algtablecenter(al);
733 nz = lg(Z)-1;
734 if (nz==1) return Z;
735
736 F = cgetg(nz+1, t_MAT);
737 for (i=1; i<=nz; i++) {
738 GEN zi = gel(Z,i);
739 gel(F,i) = FpC_sub(algpow(al,zi,p),zi,p);
740 }
741 K = FpM_ker(F,p);
742 return gerepileupto(av, FpM_mul(Z,K,p));
743 }
744
745 static GEN
_FpX_mul(void * D,GEN x,GEN y)746 _FpX_mul(void* D, GEN x, GEN y) { return FpX_mul(x,y,(GEN)D); }
747 static GEN
_FpX_pow(void * D,GEN x,GEN n)748 _FpX_pow(void* D, GEN x, GEN n) { return FpX_powu(x,itos(n),(GEN)D); }
749 static GEN
FpX_factorback(GEN fa,GEN p)750 FpX_factorback(GEN fa, GEN p)
751 {
752 return gen_factorback(gel(fa,1), zv_to_ZV(gel(fa,2)), (void *)p, &_FpX_mul, &_FpX_pow);
753 }
754
755 static GEN
out_decompose(GEN t,GEN Z,GEN P,GEN p)756 out_decompose(GEN t, GEN Z, GEN P, GEN p)
757 {
758 GEN ali = gel(t,1), projm = gel(t,2), liftm = gel(t,3), pZ;
759 if (signe(p)) pZ = FpM_image(FpM_mul(projm,Z,p),p);
760 else pZ = image(RgM_mul(projm,Z));
761 return mkvec5(ali, projm, liftm, pZ, P);
762 }
763 /* fa factorization of charpol(x) */
764 static GEN
alg_decompose_from_facto(GEN al,GEN x,GEN fa,GEN Z,long mini)765 alg_decompose_from_facto(GEN al, GEN x, GEN fa, GEN Z, long mini)
766 {
767 long k = lgcols(fa)-1, k2 = mini? 1: k/2;
768 GEN v1 = rowslice(fa,1,k2);
769 GEN v2 = rowslice(fa,k2+1,k);
770 GEN alq, P, Q, p = alg_get_char(al);
771 dbg_printf(3)(" alg_decompose_from_facto\n");
772 if (signe(p)) {
773 P = FpX_factorback(v1, p);
774 Q = FpX_factorback(v2, p);
775 P = FpX_mul(P, FpXQ_inv(P,Q,p), p);
776 }
777 else {
778 P = factorback(v1);
779 Q = factorback(v2);
780 P = RgX_mul(P, RgXQ_inv(P,Q));
781 }
782 P = algpoleval(al, P, x);
783 if (signe(p)) Q = FpC_sub(col_ei(lg(P)-1,1), P, p);
784 else Q = gsub(gen_1, P);
785 if (gequal0(P) || gequal0(Q)) return NULL;
786 alq = alg_centralproj(al, mkvec2(P,Q), 1);
787
788 P = out_decompose(gel(alq,1), Z, P, p); if (mini) return P;
789 Q = out_decompose(gel(alq,2), Z, Q, p);
790 return mkvec2(P,Q);
791 }
792
793 static GEN
random_pm1(long n)794 random_pm1(long n)
795 {
796 GEN z = cgetg(n+1,t_VECSMALL);
797 long i;
798 for (i = 1; i <= n; i++) z[i] = random_bits(5)%3 - 1;
799 return z;
800 }
801
802 static GEN alg_decompose(GEN al, GEN Z, long mini, GEN* pt_primelt);
803 /* Try to split al using x's charpoly. Return gen_0 if simple, NULL if failure.
804 * And a splitting otherwise
805 * If pt_primelt!=NULL, compute a primitive element of the center when simple */
806 static GEN
try_fact(GEN al,GEN x,GEN zx,GEN Z,GEN Zal,long mini,GEN * pt_primelt)807 try_fact(GEN al, GEN x, GEN zx, GEN Z, GEN Zal, long mini, GEN* pt_primelt)
808 {
809 GEN z, dec0, dec1, cp = algcharpoly(Zal,zx,0,1), fa, p = alg_get_char(al);
810 long nfa, e;
811 dbg_printf(3)(" try_fact: zx=%Ps\n", zx);
812 if (signe(p)) fa = FpX_factor(cp,p);
813 else fa = factor(cp);
814 dbg_printf(3)(" charpoly=%Ps\n", fa);
815 nfa = nbrows(fa);
816 if (nfa == 1) {
817 if (signe(p)) e = gel(fa,2)[1];
818 else e = itos(gcoeff(fa,1,2));
819 if (e == 1) {
820 if (pt_primelt != NULL) *pt_primelt = mkvec2(x, cp);
821 return gen_0;
822 }
823 else return NULL;
824 }
825 dec0 = alg_decompose_from_facto(al, x, fa, Z, mini);
826 if (!dec0) return NULL;
827 if (!mini) return dec0;
828 dec1 = alg_decompose(gel(dec0,1), gel(dec0,4), 1, pt_primelt);
829 z = gel(dec0,5);
830 if (!isintzero(dec1)) {
831 if (signe(p)) z = FpM_FpC_mul(gel(dec0,3),dec1,p);
832 else z = RgM_RgC_mul(gel(dec0,3),dec1);
833 }
834 return z;
835 }
836 static GEN
randcol(long n,GEN b)837 randcol(long n, GEN b)
838 {
839 GEN N = addiu(shifti(b,1), 1);
840 long i;
841 GEN res = cgetg(n+1,t_COL);
842 for (i=1; i<=n; i++)
843 {
844 pari_sp av = avma;
845 gel(res,i) = gerepileuptoint(av, subii(randomi(N),b));
846 }
847 return res;
848 }
849 /* Return gen_0 if already simple. mini: only returns a central idempotent
850 * corresponding to one simple factor
851 * if pt_primelt!=NULL, sets it to a primitive element of the center when simple */
852 static GEN
alg_decompose(GEN al,GEN Z,long mini,GEN * pt_primelt)853 alg_decompose(GEN al, GEN Z, long mini, GEN* pt_primelt)
854 {
855 pari_sp av;
856 GEN Zal, x, zx, rand, dec0, B, p;
857 long i, nz = lg(Z)-1;
858
859 if (nz == 1) {
860 if (pt_primelt != 0) *pt_primelt = mkvec2(zerocol(alg_get_dim(al)), pol_x(0));
861 return gen_0;
862 }
863 p = alg_get_char(al);
864 dbg_printf(2)(" alg_decompose: char=%Ps, dim=%d, dim Z=%d\n", p, alg_get_absdim(al), nz);
865 Zal = alg_subalg(al,Z);
866 Z = gel(Zal,2);
867 Zal = gel(Zal,1);
868 av = avma;
869
870 rand = random_pm1(nz);
871 zx = zc_to_ZC(rand);
872 if (signe(p)) {
873 zx = FpC_red(zx,p);
874 x = ZM_zc_mul(Z,rand);
875 x = FpC_red(x,p);
876 }
877 else x = RgM_zc_mul(Z,rand);
878 dec0 = try_fact(al,x,zx,Z,Zal,mini,pt_primelt);
879 if (dec0) return dec0;
880 set_avma(av);
881
882 for (i=2; i<=nz; i++)
883 {
884 dec0 = try_fact(al,gel(Z,i),col_ei(nz,i),Z,Zal,mini,pt_primelt);
885 if (dec0) return dec0;
886 set_avma(av);
887 }
888 B = int2n(10);
889 for (;;)
890 {
891 GEN x = randcol(nz,B), zx = ZM_ZC_mul(Z,x);
892 dec0 = try_fact(al,x,zx,Z,Zal,mini,pt_primelt);
893 if (dec0) return dec0;
894 set_avma(av);
895 }
896 }
897
898 static GEN
alg_decompose_total(GEN al,GEN Z,long maps)899 alg_decompose_total(GEN al, GEN Z, long maps)
900 {
901 GEN dec, sc, p;
902 long i;
903
904 dec = alg_decompose(al, Z, 0, NULL);
905 if (isintzero(dec))
906 {
907 if (maps) {
908 long n = alg_get_absdim(al);
909 al = mkvec3(al, matid(n), matid(n));
910 }
911 return mkvec(al);
912 }
913 p = alg_get_char(al); if (!signe(p)) p = NULL;
914 sc = cgetg(lg(dec), t_VEC);
915 for (i=1; i<lg(sc); i++) {
916 GEN D = gel(dec,i), a = gel(D,1), Za = gel(D,4);
917 GEN S = alg_decompose_total(a, Za, maps);
918 gel(sc,i) = S;
919 if (maps)
920 {
921 GEN projm = gel(D,2), liftm = gel(D,3);
922 long j, lS = lg(S);
923 for (j=1; j<lS; j++)
924 {
925 GEN Sj = gel(S,j), p2 = gel(Sj,2), l2 = gel(Sj,3);
926 if (p) p2 = FpM_mul(p2, projm, p);
927 else p2 = RgM_mul(p2, projm);
928 if (p) l2 = FpM_mul(liftm, l2, p);
929 else l2 = RgM_mul(liftm, l2);
930 gel(Sj,2) = p2;
931 gel(Sj,3) = l2;
932 }
933 }
934 }
935 return shallowconcat1(sc);
936 }
937
938 static GEN
alg_subalg(GEN al,GEN basis)939 alg_subalg(GEN al, GEN basis)
940 {
941 GEN invbasis, mt, p = alg_get_char(al);
942 long i, j, n = lg(basis)-1;
943
944 if (!signe(p)) p = NULL;
945 basis = shallowmatconcat(mkvec2(col_ei(n,1), basis));
946 if (p)
947 {
948 basis = image_keep_first(basis,p);
949 invbasis = FpM_inv(basis,p);
950 }
951 else
952 { /* FIXME use an integral variant of image_keep_first */
953 basis = QM_ImQ_hnf(basis);
954 invbasis = RgM_inv(basis);
955 }
956 mt = cgetg(n+1,t_VEC);
957 gel(mt,1) = matid(n);
958 for (i = 2; i <= n; i++)
959 {
960 GEN mtx = cgetg(n+1,t_MAT), x = gel(basis,i);
961 gel(mtx,1) = col_ei(n,i);
962 for (j = 2; j <= n; j++)
963 {
964 GEN xy = algmul(al, x, gel(basis,j));
965 if (p) gel(mtx,j) = FpM_FpC_mul(invbasis, xy, p);
966 else gel(mtx,j) = RgM_RgC_mul(invbasis, xy);
967 }
968 gel(mt,i) = mtx;
969 }
970 return mkvec2(algtableinit_i(mt,p), basis);
971 }
972
973 GEN
algsubalg(GEN al,GEN basis)974 algsubalg(GEN al, GEN basis)
975 {
976 pari_sp av = avma;
977 GEN p;
978 checkalg(al);
979 if (typ(basis) != t_MAT) pari_err_TYPE("algsubalg",basis);
980 p = alg_get_char(al);
981 if (signe(p)) basis = RgM_to_FpM(basis,p);
982 return gerepilecopy(av, alg_subalg(al,basis));
983 }
984
985 static int
cmp_algebra(GEN x,GEN y)986 cmp_algebra(GEN x, GEN y)
987 {
988 long d;
989 d = gel(x,1)[1] - gel(y,1)[1]; if (d) return d < 0? -1: 1;
990 d = gel(x,1)[2] - gel(y,1)[2]; if (d) return d < 0? -1: 1;
991 return cmp_universal(gel(x,2), gel(y,2));
992 }
993
994 GEN
algsimpledec_ss(GEN al,long maps)995 algsimpledec_ss(GEN al, long maps)
996 {
997 pari_sp av = avma;
998 GEN Z, p, r, res, perm;
999 long i, l, n;
1000 checkalg(al);
1001 p = alg_get_char(al);
1002 dbg_printf(1)("algsimpledec_ss: char=%Ps, dim=%d\n", p, alg_get_absdim(al));
1003 if (signe(p)) Z = algprimesubalg(al);
1004 else Z = algtablecenter(al);
1005
1006 if (lg(Z) == 2) {/* dim Z = 1 */
1007 n = alg_get_absdim(al);
1008 set_avma(av);
1009 if (!maps) return mkveccopy(al);
1010 retmkvec(mkvec3(gcopy(al), matid(n), matid(n)));
1011 }
1012 res = alg_decompose_total(al, Z, maps);
1013 l = lg(res); r = cgetg(l, t_VEC);
1014 for (i = 1; i < l; i++)
1015 {
1016 GEN A = maps? gmael(res,i,1): gel(res,i);
1017 gel(r,i) = mkvec2(mkvecsmall2(alg_get_dim(A), lg(algtablecenter(A))),
1018 alg_get_multable(A));
1019 }
1020 perm = gen_indexsort(r, (void*)cmp_algebra, &cmp_nodata);
1021 return gerepilecopy(av, vecpermute(res, perm));
1022 }
1023
1024 GEN
algsimpledec(GEN al,long maps)1025 algsimpledec(GEN al, long maps)
1026 {
1027 pari_sp av = avma;
1028 int ss;
1029 GEN rad, dec, res, proj=NULL, lift=NULL;
1030 rad = algradical(al);
1031 ss = gequal0(rad);
1032 if (!ss)
1033 {
1034 al = alg_quotient(al, rad, maps);
1035 if (maps) {
1036 proj = gel(al,2);
1037 lift = gel(al,3);
1038 al = gel(al,1);
1039 }
1040 }
1041 dec = algsimpledec_ss(al, maps);
1042 if (!ss && maps) /* update maps */
1043 {
1044 GEN p = alg_get_char(al);
1045 long i;
1046 for (i=1; i<lg(dec); i++)
1047 {
1048 if (signe(p))
1049 {
1050 gmael(dec,i,2) = FpM_mul(gmael(dec,i,2), proj, p);
1051 gmael(dec,i,3) = FpM_mul(lift, gmael(dec,i,3), p);
1052 }
1053 else
1054 {
1055 gmael(dec,i,2) = RgM_mul(gmael(dec,i,2), proj);
1056 gmael(dec,i,3) = RgM_mul(lift, gmael(dec,i,3));
1057 }
1058 }
1059 }
1060 res = mkvec2(rad, dec);
1061 return gerepilecopy(av,res);
1062 }
1063
1064 static GEN alg_idempotent(GEN al, long n, long d);
1065 static GEN
try_split(GEN al,GEN x,long n,long d)1066 try_split(GEN al, GEN x, long n, long d)
1067 {
1068 GEN cp, p = alg_get_char(al), fa, e, pol, exp, P, Q, U, u, mx, mte, ire;
1069 long nfa, i, smalldim = alg_get_absdim(al)+1, dim, smalli = 0;
1070 cp = algcharpoly(al,x,0,1);
1071 fa = FpX_factor(cp,p);
1072 nfa = nbrows(fa);
1073 if (nfa == 1) return NULL;
1074 pol = gel(fa,1);
1075 exp = gel(fa,2);
1076
1077 /* charpoly is always a d-th power */
1078 for (i=1; i<lg(exp); i++) {
1079 if (exp[i]%d) pari_err(e_MISC, "the algebra must be simple (try_split 1)");
1080 exp[i] /= d;
1081 }
1082 cp = FpX_factorback(fa,p);
1083
1084 /* find smallest Fp-dimension of a characteristic space */
1085 for (i=1; i<lg(pol); i++) {
1086 dim = degree(gel(pol,i))*exp[i];
1087 if (dim < smalldim) {
1088 smalldim = dim;
1089 smalli = i;
1090 }
1091 }
1092 i = smalli;
1093 if (smalldim != n) return NULL;
1094 /* We could also compute e*al*e and try again with this smaller algebra */
1095 /* Fq-rank 1 = Fp-rank n idempotent: success */
1096
1097 /* construct idempotent */
1098 mx = algbasismultable(al,x);
1099 P = gel(pol,i);
1100 P = FpX_powu(P, exp[i], p);
1101 Q = FpX_div(cp, P, p);
1102 e = algpoleval(al, Q, mkvec2(x,mx));
1103 U = FpXQ_inv(Q, P, p);
1104 u = algpoleval(al, U, mkvec2(x,mx));
1105 e = algbasismul(al, e, u);
1106 mte = algbasisrightmultable(al,e);
1107 ire = FpM_indexrank(mte,p);
1108 if (lg(gel(ire,1))-1 != smalldim*d) pari_err(e_MISC, "the algebra must be simple (try_split 2)");
1109
1110 return mkvec3(e,mte,ire);
1111 }
1112
1113 /*
1114 * Given a simple algebra al of dimension d^2 over its center of degree n,
1115 * find an idempotent e in al with rank n (which is minimal).
1116 */
1117 static GEN
alg_idempotent(GEN al,long n,long d)1118 alg_idempotent(GEN al, long n, long d)
1119 {
1120 pari_sp av = avma;
1121 long i, N = alg_get_absdim(al);
1122 GEN e, p = alg_get_char(al), x;
1123 for(i=2; i<=N; i++) {
1124 x = col_ei(N,i);
1125 e = try_split(al, x, n, d);
1126 if (e) return e;
1127 set_avma(av);
1128 }
1129 for(;;) {
1130 x = random_FpC(N,p);
1131 e = try_split(al, x, n, d);
1132 if (e) return e;
1133 set_avma(av);
1134 }
1135 }
1136
1137 static GEN
try_descend(GEN M,GEN B,GEN p,long m,long n,long d)1138 try_descend(GEN M, GEN B, GEN p, long m, long n, long d)
1139 {
1140 GEN B2 = cgetg(m+1,t_MAT), b;
1141 long i, j, k=0;
1142 for (i=1; i<=d; i++)
1143 {
1144 k++;
1145 b = gel(B,i);
1146 gel(B2,k) = b;
1147 for (j=1; j<n; j++)
1148 {
1149 k++;
1150 b = FpM_FpC_mul(M,b,p);
1151 gel(B2,k) = b;
1152 }
1153 }
1154 if (!signe(FpM_det(B2,p))) return NULL;
1155 return FpM_inv(B2,p);
1156 }
1157
1158 /* Given an m*m matrix M with irreducible charpoly over F of degree n,
1159 * let K = F(M), which is a field, and write m=d*n.
1160 * Compute the d-dimensional K-vector space structure on V=F^m induced by M.
1161 * Return [B,C] where:
1162 * - B is m*d matrix over F giving a K-basis b_1,...,b_d of V
1163 * - C is d*m matrix over F[x] expressing the canonical F-basis of V on the b_i
1164 * Currently F = Fp TODO extend this. */
1165 static GEN
descend_i(GEN M,long n,GEN p)1166 descend_i(GEN M, long n, GEN p)
1167 {
1168 GEN B, C;
1169 long m,d,i;
1170 pari_sp av;
1171 m = lg(M)-1;
1172 d = m/n;
1173 B = cgetg(d+1,t_MAT);
1174 av = avma;
1175
1176 /* try a subset of the canonical basis */
1177 for (i=1; i<=d; i++)
1178 gel(B,i) = col_ei(m,n*(i-1)+1);
1179 C = try_descend(M,B,p,m,n,d);
1180 if (C) return mkvec2(B,C);
1181 set_avma(av);
1182
1183 /* try smallish elements */
1184 for (i=1; i<=d; i++)
1185 gel(B,i) = FpC_red(zc_to_ZC(random_pm1(m)),p);
1186 C = try_descend(M,B,p,m,n,d);
1187 if (C) return mkvec2(B,C);
1188 set_avma(av);
1189
1190 /* try random elements */
1191 for (;;)
1192 {
1193 for (i=1; i<=d; i++)
1194 gel(B,i) = random_FpC(m,p);
1195 C = try_descend(M,B,p,m,n,d);
1196 if (C) return mkvec2(B,C);
1197 set_avma(av);
1198 }
1199 }
1200 static GEN
RgC_contract(GEN C,long n,long v)1201 RgC_contract(GEN C, long n, long v) /* n>1 */
1202 {
1203 GEN C2, P;
1204 long m, d, i, j;
1205 m = lg(C)-1;
1206 d = m/n;
1207 C2 = cgetg(d+1,t_COL);
1208 for (i=1; i<=d; i++)
1209 {
1210 P = pol_xn(n-1,v);
1211 for (j=1; j<=n; j++)
1212 gel(P,j+1) = gel(C,n*(i-1)+j);
1213 P = normalizepol(P);
1214 gel(C2,i) = P;
1215 }
1216 return C2;
1217 }
1218 static GEN
RgM_contract(GEN A,long n,long v)1219 RgM_contract(GEN A, long n, long v) /* n>1 */
1220 {
1221 GEN A2 = cgetg(lg(A),t_MAT);
1222 long i;
1223 for (i=1; i<lg(A2); i++)
1224 gel(A2,i) = RgC_contract(gel(A,i),n,v);
1225 return A2;
1226 }
1227 static GEN
descend(GEN M,long n,GEN p,long v)1228 descend(GEN M, long n, GEN p, long v)
1229 {
1230 GEN res = descend_i(M,n,p);
1231 gel(res,2) = RgM_contract(gel(res,2),n,v);
1232 return res;
1233 }
1234
1235 /* isomorphism of Fp-vector spaces M_d(F_p^n) -> (F_p)^(d^2*n) */
1236 static GEN
Fq_mat2col(GEN M,long d,long n)1237 Fq_mat2col(GEN M, long d, long n)
1238 {
1239 long N = d*d*n, i, j, k;
1240 GEN C = cgetg(N+1, t_COL);
1241 for (i=1; i<=d; i++)
1242 for (j=1; j<=d; j++)
1243 for (k=0; k<n; k++)
1244 gel(C,n*(d*(i-1)+j-1)+k+1) = polcoef_i(gcoeff(M,i,j),k,-1);
1245 return C;
1246 }
1247
1248 static GEN
alg_finite_csa_split(GEN al,long v)1249 alg_finite_csa_split(GEN al, long v)
1250 {
1251 GEN Z, e, mte, ire, primelt, b, T, M, proje, lifte, extre, p, B, C, mt, mx, map, mapi, T2, ro;
1252 long n, d, N = alg_get_absdim(al), i;
1253 p = alg_get_char(al);
1254 /* compute the center */
1255 Z = algcenter(al);
1256 /* TODO option to give the center as input instead of computing it */
1257 n = lg(Z)-1;
1258
1259 /* compute a minimal rank idempotent e */
1260 if (n==N) {
1261 d = 1;
1262 e = col_ei(N,1);
1263 mte = matid(N);
1264 ire = mkvec2(identity_perm(n),identity_perm(n));
1265 }
1266 else {
1267 d = usqrt(N/n);
1268 if (d*d*n != N) pari_err(e_MISC, "the algebra must be simple (alg_finite_csa_split 1)");
1269 e = alg_idempotent(al,n,d);
1270 mte = gel(e,2);
1271 ire = gel(e,3);
1272 e = gel(e,1);
1273 }
1274
1275 /* identify the center */
1276 if (n==1)
1277 {
1278 T = pol_x(v);
1279 primelt = gen_0;
1280 }
1281 else
1282 {
1283 b = alg_decompose(al, Z, 1, &primelt);
1284 if (!gequal0(b)) pari_err(e_MISC, "the algebra must be simple (alg_finite_csa_split 2)");
1285 T = gel(primelt,2);
1286 primelt = gel(primelt,1);
1287 setvarn(T,v);
1288 }
1289
1290 /* use the ffinit polynomial */
1291 if (n>1)
1292 {
1293 T2 = init_Fq(p,n,v);
1294 setvarn(T,fetch_var_higher());
1295 ro = FpXQX_roots(T2,T,p);
1296 ro = gel(ro,1);
1297 primelt = algpoleval(al,ro,primelt);
1298 T = T2;
1299 }
1300
1301 /* descend al*e to a vector space over the center */
1302 /* lifte: al*e -> al ; proje: al*e -> al */
1303 lifte = shallowextract(mte,gel(ire,2));
1304 extre = shallowmatextract(mte,gel(ire,1),gel(ire,2));
1305 extre = FpM_inv(extre,p);
1306 proje = rowpermute(mte,gel(ire,1));
1307 proje = FpM_mul(extre,proje,p);
1308 if (n==1)
1309 {
1310 B = lifte;
1311 C = proje;
1312 }
1313 else
1314 {
1315 M = algbasismultable(al,primelt);
1316 M = FpM_mul(M,lifte,p);
1317 M = FpM_mul(proje,M,p);
1318 B = descend(M,n,p,v);
1319 C = gel(B,2);
1320 B = gel(B,1);
1321 B = FpM_mul(lifte,B,p);
1322 C = FqM_mul(C,proje,T,p);
1323 }
1324
1325 /* compute the isomorphism */
1326 mt = alg_get_multable(al);
1327 map = cgetg(N+1,t_VEC);
1328 M = cgetg(N+1,t_MAT);
1329 for (i=1; i<=N; i++)
1330 {
1331 mx = gel(mt,i);
1332 mx = FpM_mul(mx,B,p);
1333 mx = FqM_mul(C,mx,T,p);
1334 gel(map,i) = mx;
1335 gel(M,i) = Fq_mat2col(mx,d,n);
1336 }
1337 mapi = FpM_inv(M,p);
1338 if (!mapi) pari_err(e_MISC, "the algebra must be simple (alg_finite_csa_split 3)");
1339 return mkvec3(T,map,mapi);
1340 }
1341
1342 GEN
algsplit(GEN al,long v)1343 algsplit(GEN al, long v)
1344 {
1345 pari_sp av = avma;
1346 GEN res, T, map, mapi, ff, p;
1347 long i,j,k,li,lj;
1348 checkalg(al);
1349 p = alg_get_char(al);
1350 if (gequal0(p))
1351 pari_err_IMPL("splitting a characteristic 0 algebra over its center");
1352 res = alg_finite_csa_split(al, v);
1353 T = gel(res,1);
1354 map = gel(res,2);
1355 mapi = gel(res,3);
1356 ff = Tp_to_FF(T,p);
1357 for (i=1; i<lg(map); i++)
1358 {
1359 li = lg(gel(map,i));
1360 for (j=1; j<li; j++)
1361 {
1362 lj = lg(gmael(map,i,j));
1363 for (k=1; k<lj; k++)
1364 gmael3(map,i,j,k) = Fq_to_FF(gmael3(map,i,j,k),ff);
1365 }
1366 }
1367
1368 return gerepilecopy(av, mkvec2(map,mapi));
1369 }
1370
1371 /* multiplication table sanity checks */
1372 static GEN
check_mt_noid(GEN mt,GEN p)1373 check_mt_noid(GEN mt, GEN p)
1374 {
1375 long i, l;
1376 GEN MT = cgetg_copy(mt, &l);
1377 if (typ(MT) != t_VEC || l == 1) return NULL;
1378 for (i = 1; i < l; i++)
1379 {
1380 GEN M = gel(mt,i);
1381 if (typ(M) != t_MAT || lg(M) != l || lgcols(M) != l) return NULL;
1382 if (p) M = RgM_to_FpM(M,p);
1383 gel(MT,i) = M;
1384 }
1385 return MT;
1386 }
1387 static GEN
check_mt(GEN mt,GEN p)1388 check_mt(GEN mt, GEN p)
1389 {
1390 long i;
1391 GEN MT;
1392 MT = check_mt_noid(mt, p);
1393 if (!MT || !ZM_isidentity(gel(MT,1))) return NULL;
1394 for (i=2; i<lg(MT); i++)
1395 if (ZC_is_ei(gmael(MT,i,1)) != i) return NULL;
1396 return MT;
1397 }
1398
1399 static GEN
check_relmt(GEN nf,GEN mt)1400 check_relmt(GEN nf, GEN mt)
1401 {
1402 long i, l = lg(mt), j, k;
1403 GEN MT = gcopy(mt), a, b, d;
1404 if (typ(MT) != t_VEC || l == 1) return NULL;
1405 for (i = 1; i < l; i++)
1406 {
1407 GEN M = gel(MT,i);
1408 if (typ(M) != t_MAT || lg(M) != l || lgcols(M) != l) return NULL;
1409 for (k = 1; k < l; k++)
1410 for (j = 1; j < l; j++)
1411 {
1412 a = gcoeff(M,j,k);
1413 if (typ(a)==t_INT) continue;
1414 b = algtobasis(nf,a);
1415 d = Q_denom(b);
1416 if (!isint1(d))
1417 pari_err_DOMAIN("alg_csa_table", "denominator(mt)", "!=", gen_1, mt);
1418 gcoeff(M,j,k) = lift(basistoalg(nf,b));
1419 }
1420 if (i > 1 && RgC_is_ei(gel(M,1)) != i) return NULL; /* i = 1 checked at end */
1421 gel(MT,i) = M;
1422 }
1423 if (!RgM_isidentity(gel(MT,1))) return NULL;
1424 return MT;
1425 }
1426
1427 int
algisassociative(GEN mt0,GEN p)1428 algisassociative(GEN mt0, GEN p)
1429 {
1430 pari_sp av = avma;
1431 long i, j, k, n;
1432 GEN M, mt;
1433
1434 if (checkalg_i(mt0)) { p = alg_get_char(mt0); mt0 = alg_get_multable(mt0); }
1435 if (typ(p) != t_INT) pari_err_TYPE("algisassociative",p);
1436 mt = check_mt_noid(mt0, isintzero(p)? NULL: p);
1437 if (!mt) pari_err_TYPE("algisassociative (mult. table)", mt0);
1438 if (!ZM_isidentity(gel(mt,1))) return gc_bool(av,0);
1439 n = lg(mt)-1;
1440 M = cgetg(n+1,t_MAT);
1441 for (j=1; j<=n; j++) gel(M,j) = cgetg(n+1,t_COL);
1442 for (i=1; i<=n; i++)
1443 {
1444 GEN mi = gel(mt,i);
1445 for (j=1; j<=n; j++) gcoeff(M,i,j) = gel(mi,j); /* ei.ej */
1446 }
1447 for (i=2; i<=n; i++) {
1448 GEN mi = gel(mt,i);
1449 for (j=2; j<=n; j++) {
1450 for (k=2; k<=n; k++) {
1451 GEN x, y;
1452 if (signe(p)) {
1453 x = _tablemul_ej_Fp(mt,gcoeff(M,i,j),k,p);
1454 y = FpM_FpC_mul(mi,gcoeff(M,j,k),p);
1455 }
1456 else {
1457 x = _tablemul_ej(mt,gcoeff(M,i,j),k);
1458 y = RgM_RgC_mul(mi,gcoeff(M,j,k));
1459 }
1460 /* not cmp_universal: must not fail on 0 == Mod(0,2) for instance */
1461 if (!gequal(x,y)) return gc_bool(av,0);
1462 }
1463 }
1464 }
1465 return gc_bool(av,1);
1466 }
1467
1468 int
algiscommutative(GEN al)1469 algiscommutative(GEN al) /* assumes e_1 = 1 */
1470 {
1471 long i,j,k,N,sp;
1472 GEN mt,a,b,p;
1473 checkalg(al);
1474 if (alg_type(al) != al_TABLE) return alg_get_degree(al)==1;
1475 N = alg_get_absdim(al);
1476 mt = alg_get_multable(al);
1477 p = alg_get_char(al);
1478 sp = signe(p);
1479 for (i=2; i<=N; i++)
1480 for (j=2; j<=N; j++)
1481 for (k=1; k<=N; k++) {
1482 a = gcoeff(gel(mt,i),k,j);
1483 b = gcoeff(gel(mt,j),k,i);
1484 if (sp) {
1485 if (cmpii(Fp_red(a,p), Fp_red(b,p))) return 0;
1486 }
1487 else if (gcmp(a,b)) return 0;
1488 }
1489 return 1;
1490 }
1491
1492 int
algissemisimple(GEN al)1493 algissemisimple(GEN al)
1494 {
1495 pari_sp av = avma;
1496 GEN rad;
1497 checkalg(al);
1498 if (alg_type(al) != al_TABLE) return 1;
1499 rad = algradical(al);
1500 set_avma(av);
1501 return gequal0(rad);
1502 }
1503
1504 /* ss : known to be semisimple */
1505 int
algissimple(GEN al,long ss)1506 algissimple(GEN al, long ss)
1507 {
1508 pari_sp av = avma;
1509 GEN Z, dec, p;
1510 checkalg(al);
1511 if (alg_type(al) != al_TABLE) return 1;
1512 if (!ss && !algissemisimple(al)) return 0;
1513
1514 p = alg_get_char(al);
1515 if (signe(p)) Z = algprimesubalg(al);
1516 else Z = algtablecenter(al);
1517
1518 if (lg(Z) == 2) {/* dim Z = 1 */
1519 set_avma(av);
1520 return 1;
1521 }
1522 dec = alg_decompose(al, Z, 1, NULL);
1523 set_avma(av);
1524 return gequal0(dec);
1525 }
1526
1527 static long
is_place_emb(GEN nf,GEN pl)1528 is_place_emb(GEN nf, GEN pl)
1529 {
1530 long r, r1, r2;
1531 if (typ(pl) != t_INT) pari_err_TYPE("is_place_emb", pl);
1532 if (signe(pl)<=0) pari_err_DOMAIN("is_place_emb", "pl", "<=", gen_0, pl);
1533 nf_get_sign(nf,&r1,&r2); r = r1+r2;
1534 if (cmpiu(pl,r)>0) pari_err_DOMAIN("is_place_emb", "pl", ">", utoi(r), pl);
1535 return itou(pl);
1536 }
1537
1538 static long
alghasse_emb(GEN al,long emb)1539 alghasse_emb(GEN al, long emb)
1540 {
1541 GEN nf = alg_get_center(al);
1542 long r1 = nf_get_r1(nf);
1543 return (emb <= r1)? alg_get_hasse_i(al)[emb]: 0;
1544 }
1545
1546 static long
alghasse_pr(GEN al,GEN pr)1547 alghasse_pr(GEN al, GEN pr)
1548 {
1549 GEN hf = alg_get_hasse_f(al);
1550 long i = tablesearch(gel(hf,1), pr, &cmp_prime_ideal);
1551 return i? gel(hf,2)[i]: 0;
1552 }
1553
1554 static long
alghasse_0(GEN al,GEN pl)1555 alghasse_0(GEN al, GEN pl)
1556 {
1557 GEN pr, nf;
1558 if (alg_type(al)== al_CSA)
1559 pari_err_IMPL("computation of Hasse invariants over table CSA");
1560 if ((pr = get_prid(pl))) return alghasse_pr(al, pr);
1561 nf = alg_get_center(al);
1562 return alghasse_emb(al, is_place_emb(nf, pl));
1563 }
1564 GEN
alghasse(GEN al,GEN pl)1565 alghasse(GEN al, GEN pl)
1566 {
1567 long h;
1568 checkalg(al);
1569 if (alg_type(al) == al_TABLE) pari_err_TYPE("alghasse [use alginit]",al);
1570 h = alghasse_0(al,pl);
1571 return sstoQ(h, alg_get_degree(al));
1572 }
1573
1574 /* h >= 0, d >= 0 */
1575 static long
indexfromhasse(long h,long d)1576 indexfromhasse(long h, long d) { return d/ugcd(h,d); }
1577
1578 long
algindex(GEN al,GEN pl)1579 algindex(GEN al, GEN pl)
1580 {
1581 long d, res, i, l;
1582 GEN hi, hf;
1583
1584 checkalg(al);
1585 if (alg_type(al) == al_TABLE) pari_err_TYPE("algindex [use alginit]",al);
1586 d = alg_get_degree(al);
1587 if (pl) return indexfromhasse(alghasse_0(al,pl), d);
1588
1589 /* else : global index */
1590 res = 1;
1591 hi = alg_get_hasse_i(al); l = lg(hi);
1592 for (i=1; i<l && res!=d; i++) res = ulcm(res, indexfromhasse(hi[i],d));
1593 hf = gel(alg_get_hasse_f(al), 2); l = lg(hf);
1594 for (i=1; i<l && res!=d; i++) res = ulcm(res, indexfromhasse(hf[i],d));
1595 return res;
1596 }
1597
1598 int
algisdivision(GEN al,GEN pl)1599 algisdivision(GEN al, GEN pl)
1600 {
1601 checkalg(al);
1602 if (alg_type(al) == al_TABLE) {
1603 if (!algissimple(al,0)) return 0;
1604 if (algiscommutative(al)) return 1;
1605 pari_err_IMPL("algisdivision for table algebras");
1606 }
1607 return algindex(al,pl) == alg_get_degree(al);
1608 }
1609
1610 int
algissplit(GEN al,GEN pl)1611 algissplit(GEN al, GEN pl)
1612 {
1613 checkalg(al);
1614 if (alg_type(al) == al_TABLE) pari_err_TYPE("algissplit [use alginit]", al);
1615 return algindex(al,pl) == 1;
1616 }
1617
1618 int
algisramified(GEN al,GEN pl)1619 algisramified(GEN al, GEN pl)
1620 {
1621 checkalg(al);
1622 if (alg_type(al) == al_TABLE) pari_err_TYPE("algisramified [use alginit]", al);
1623 return algindex(al,pl) != 1;
1624 }
1625
1626 GEN
algramifiedplaces(GEN al)1627 algramifiedplaces(GEN al)
1628 {
1629 pari_sp av = avma;
1630 GEN ram, hf, hi, Lpr;
1631 long r1, count, i;
1632 checkalg(al);
1633 if (alg_type(al) == al_TABLE) pari_err_TYPE("algramifiedplaces [use alginit]", al);
1634 r1 = nf_get_r1(alg_get_center(al));
1635 hi = alg_get_hasse_i(al);
1636 hf = alg_get_hasse_f(al);
1637 Lpr = gel(hf,1);
1638 hf = gel(hf,2);
1639 ram = cgetg(r1+lg(Lpr), t_VEC);
1640 count = 0;
1641 for (i=1; i<=r1; i++)
1642 if (hi[i]) {
1643 count++;
1644 gel(ram,count) = stoi(i);
1645 }
1646 for (i=1; i<lg(Lpr); i++)
1647 if (hf[i]) {
1648 count++;
1649 gel(ram,count) = gel(Lpr,i);
1650 }
1651 setlg(ram, count+1);
1652 return gerepilecopy(av, ram);
1653 }
1654
1655 /** OPERATIONS ON ELEMENTS operations.c **/
1656
1657 static long
alg_model0(GEN al,GEN x)1658 alg_model0(GEN al, GEN x)
1659 {
1660 long t, N = alg_get_absdim(al), lx = lg(x), d, n, D, i;
1661 if (typ(x) == t_MAT) return al_MATRIX;
1662 if (typ(x) != t_COL) return al_INVALID;
1663 if (N == 1) {
1664 if (lx != 2) return al_INVALID;
1665 switch(typ(gel(x,1)))
1666 {
1667 case t_INT: case t_FRAC: return al_TRIVIAL; /* cannot distinguish basis and alg from size */
1668 case t_POL: case t_POLMOD: return al_ALGEBRAIC;
1669 default: return al_INVALID;
1670 }
1671 }
1672
1673 switch(alg_type(al)) {
1674 case al_TABLE:
1675 if (lx != N+1) return al_INVALID;
1676 return al_BASIS;
1677 case al_CYCLIC:
1678 d = alg_get_degree(al);
1679 if (lx == N+1) return al_BASIS;
1680 if (lx == d+1) return al_ALGEBRAIC;
1681 return al_INVALID;
1682 case al_CSA:
1683 D = alg_get_dim(al);
1684 n = nf_get_degree(alg_get_center(al));
1685 if (n == 1) {
1686 if (lx != D+1) return al_INVALID;
1687 for (i=1; i<=D; i++) {
1688 t = typ(gel(x,i));
1689 if (t == t_POL || t == t_POLMOD) return al_ALGEBRAIC;
1690 /* TODO t_COL for coefficients in basis form ? */
1691 }
1692 return al_BASIS;
1693 }
1694 else {
1695 if (lx == N+1) return al_BASIS;
1696 if (lx == D+1) return al_ALGEBRAIC;
1697 return al_INVALID;
1698 }
1699 }
1700 return al_INVALID; /* LCOV_EXCL_LINE */
1701 }
1702
1703 static void
checkalgx(GEN x,long model)1704 checkalgx(GEN x, long model)
1705 {
1706 long t, i;
1707 switch(model) {
1708 case al_BASIS:
1709 for (i=1; i<lg(x); i++) {
1710 t = typ(gel(x,i));
1711 if (t != t_INT && t != t_FRAC)
1712 pari_err_TYPE("checkalgx", gel(x,i));
1713 }
1714 return;
1715 case al_TRIVIAL:
1716 case al_ALGEBRAIC:
1717 for (i=1; i<lg(x); i++) {
1718 t = typ(gel(x,i));
1719 if (t != t_INT && t != t_FRAC && t != t_POL && t != t_POLMOD)
1720 /* TODO t_COL ? */
1721 pari_err_TYPE("checkalgx", gel(x,i));
1722 }
1723 return;
1724 }
1725 }
1726
1727 long
alg_model(GEN al,GEN x)1728 alg_model(GEN al, GEN x)
1729 {
1730 long res = alg_model0(al, x);
1731 if (res == al_INVALID) pari_err_TYPE("alg_model", x);
1732 checkalgx(x, res); return res;
1733 }
1734
1735 static GEN
alC_add_i(GEN al,GEN x,GEN y,long lx)1736 alC_add_i(GEN al, GEN x, GEN y, long lx)
1737 {
1738 GEN A = cgetg(lx, t_COL);
1739 long i;
1740 for (i=1; i<lx; i++) gel(A,i) = algadd(al, gel(x,i), gel(y,i));
1741 return A;
1742 }
1743 static GEN
alM_add(GEN al,GEN x,GEN y)1744 alM_add(GEN al, GEN x, GEN y)
1745 {
1746 long lx = lg(x), l, j;
1747 GEN z;
1748 if (lg(y) != lx) pari_err_DIM("alM_add (rows)");
1749 if (lx == 1) return cgetg(1, t_MAT);
1750 z = cgetg(lx, t_MAT); l = lgcols(x);
1751 if (lgcols(y) != l) pari_err_DIM("alM_add (columns)");
1752 for (j = 1; j < lx; j++) gel(z,j) = alC_add_i(al, gel(x,j), gel(y,j), l);
1753 return z;
1754 }
1755 GEN
algadd(GEN al,GEN x,GEN y)1756 algadd(GEN al, GEN x, GEN y)
1757 {
1758 pari_sp av = avma;
1759 long tx, ty;
1760 GEN p;
1761 checkalg(al);
1762 tx = alg_model(al,x);
1763 ty = alg_model(al,y);
1764 p = alg_get_char(al);
1765 if (signe(p)) return FpC_add(x,y,p);
1766 if (tx==ty) {
1767 if (tx!=al_MATRIX) return gadd(x,y);
1768 return gerepilecopy(av, alM_add(al,x,y));
1769 }
1770 if (tx==al_ALGEBRAIC) x = algalgtobasis(al,x);
1771 if (ty==al_ALGEBRAIC) y = algalgtobasis(al,y);
1772 return gerepileupto(av, gadd(x,y));
1773 }
1774
1775 GEN
algneg(GEN al,GEN x)1776 algneg(GEN al, GEN x) { checkalg(al); (void)alg_model(al,x); return gneg(x); }
1777
1778 static GEN
alC_sub_i(GEN al,GEN x,GEN y,long lx)1779 alC_sub_i(GEN al, GEN x, GEN y, long lx)
1780 {
1781 long i;
1782 GEN A = cgetg(lx, t_COL);
1783 for (i=1; i<lx; i++) gel(A,i) = algsub(al, gel(x,i), gel(y,i));
1784 return A;
1785 }
1786 static GEN
alM_sub(GEN al,GEN x,GEN y)1787 alM_sub(GEN al, GEN x, GEN y)
1788 {
1789 long lx = lg(x), l, j;
1790 GEN z;
1791 if (lg(y) != lx) pari_err_DIM("alM_sub (rows)");
1792 if (lx == 1) return cgetg(1, t_MAT);
1793 z = cgetg(lx, t_MAT); l = lgcols(x);
1794 if (lgcols(y) != l) pari_err_DIM("alM_sub (columns)");
1795 for (j = 1; j < lx; j++) gel(z,j) = alC_sub_i(al, gel(x,j), gel(y,j), l);
1796 return z;
1797 }
1798 GEN
algsub(GEN al,GEN x,GEN y)1799 algsub(GEN al, GEN x, GEN y)
1800 {
1801 long tx, ty;
1802 pari_sp av = avma;
1803 GEN p;
1804 checkalg(al);
1805 tx = alg_model(al,x);
1806 ty = alg_model(al,y);
1807 p = alg_get_char(al);
1808 if (signe(p)) return FpC_sub(x,y,p);
1809 if (tx==ty) {
1810 if (tx != al_MATRIX) return gsub(x,y);
1811 return gerepilecopy(av, alM_sub(al,x,y));
1812 }
1813 if (tx==al_ALGEBRAIC) x = algalgtobasis(al,x);
1814 if (ty==al_ALGEBRAIC) y = algalgtobasis(al,y);
1815 return gerepileupto(av, gsub(x,y));
1816 }
1817
1818 static GEN
algalgmul_cyc(GEN al,GEN x,GEN y)1819 algalgmul_cyc(GEN al, GEN x, GEN y)
1820 {
1821 pari_sp av = avma;
1822 long n = alg_get_degree(al), i, k;
1823 GEN xalg, yalg, res, rnf, auts, sum, b, prod, autx;
1824 rnf = alg_get_splittingfield(al);
1825 auts = alg_get_auts(al);
1826 b = alg_get_b(al);
1827
1828 xalg = cgetg(n+1, t_COL);
1829 for (i=0; i<n; i++)
1830 gel(xalg,i+1) = lift_shallow(rnfbasistoalg(rnf,gel(x,i+1)));
1831
1832 yalg = cgetg(n+1, t_COL);
1833 for (i=0; i<n; i++) gel(yalg,i+1) = rnfbasistoalg(rnf,gel(y,i+1));
1834
1835 res = cgetg(n+1,t_COL);
1836 for (k=0; k<n; k++) {
1837 gel(res,k+1) = gmul(gel(xalg,k+1),gel(yalg,1));
1838 for (i=1; i<=k; i++) {
1839 autx = poleval(gel(xalg,k-i+1),gel(auts,i));
1840 prod = gmul(autx,gel(yalg,i+1));
1841 gel(res,k+1) = gadd(gel(res,k+1), prod);
1842 }
1843
1844 sum = gen_0;
1845 for (; i<n; i++) {
1846 autx = poleval(gel(xalg,k+n-i+1),gel(auts,i));
1847 prod = gmul(autx,gel(yalg,i+1));
1848 sum = gadd(sum,prod);
1849 }
1850 sum = gmul(b,sum);
1851
1852 gel(res,k+1) = gadd(gel(res,k+1),sum);
1853 }
1854
1855 return gerepilecopy(av, res);
1856 }
1857
1858 static GEN
_tablemul(GEN mt,GEN x,GEN y)1859 _tablemul(GEN mt, GEN x, GEN y)
1860 {
1861 pari_sp av = avma;
1862 long D = lg(mt)-1, i;
1863 GEN res = NULL;
1864 for (i=1; i<=D; i++) {
1865 GEN c = gel(x,i);
1866 if (!gequal0(c)) {
1867 GEN My = RgM_RgC_mul(gel(mt,i),y);
1868 GEN t = RgC_Rg_mul(My,c);
1869 res = res? RgC_add(res,t): t;
1870 }
1871 }
1872 if (!res) { set_avma(av); return zerocol(D); }
1873 return gerepileupto(av, res);
1874 }
1875
1876 static GEN
_tablemul_Fp(GEN mt,GEN x,GEN y,GEN p)1877 _tablemul_Fp(GEN mt, GEN x, GEN y, GEN p)
1878 {
1879 pari_sp av = avma;
1880 long D = lg(mt)-1, i;
1881 GEN res = NULL;
1882 for (i=1; i<=D; i++) {
1883 GEN c = gel(x,i);
1884 if (signe(c)) {
1885 GEN My = FpM_FpC_mul(gel(mt,i),y,p);
1886 GEN t = FpC_Fp_mul(My,c,p);
1887 res = res? FpC_add(res,t,p): t;
1888 }
1889 }
1890 if (!res) { set_avma(av); return zerocol(D); }
1891 return gerepileupto(av, res);
1892 }
1893
1894 /* x*ej */
1895 static GEN
_tablemul_ej(GEN mt,GEN x,long j)1896 _tablemul_ej(GEN mt, GEN x, long j)
1897 {
1898 pari_sp av = avma;
1899 long D = lg(mt)-1, i;
1900 GEN res = NULL;
1901 for (i=1; i<=D; i++) {
1902 GEN c = gel(x,i);
1903 if (!gequal0(c)) {
1904 GEN My = gel(gel(mt,i),j);
1905 GEN t = RgC_Rg_mul(My,c);
1906 res = res? RgC_add(res,t): t;
1907 }
1908 }
1909 if (!res) { set_avma(av); return zerocol(D); }
1910 return gerepileupto(av, res);
1911 }
1912 static GEN
_tablemul_ej_Fp(GEN mt,GEN x,long j,GEN p)1913 _tablemul_ej_Fp(GEN mt, GEN x, long j, GEN p)
1914 {
1915 pari_sp av = avma;
1916 long D = lg(mt)-1, i;
1917 GEN res = NULL;
1918 for (i=1; i<=D; i++) {
1919 GEN c = gel(x,i);
1920 if (!gequal0(c)) {
1921 GEN My = gel(gel(mt,i),j);
1922 GEN t = FpC_Fp_mul(My,c,p);
1923 res = res? FpC_add(res,t,p): t;
1924 }
1925 }
1926 if (!res) { set_avma(av); return zerocol(D); }
1927 return gerepileupto(av, res);
1928 }
1929
1930 static GEN
_tablemul_ej_Fl(GEN mt,GEN x,long j,ulong p)1931 _tablemul_ej_Fl(GEN mt, GEN x, long j, ulong p)
1932 {
1933 pari_sp av = avma;
1934 long D = lg(mt)-1, i;
1935 GEN res = NULL;
1936 for (i=1; i<=D; i++) {
1937 ulong c = x[i];
1938 if (c) {
1939 GEN My = gel(gel(mt,i),j);
1940 GEN t = Flv_Fl_mul(My,c, p);
1941 res = res? Flv_add(res,t, p): t;
1942 }
1943 }
1944 if (!res) { set_avma(av); return zero_Flv(D); }
1945 return gerepileupto(av, res);
1946 }
1947
1948 static GEN
algalgmul_csa(GEN al,GEN x,GEN y)1949 algalgmul_csa(GEN al, GEN x, GEN y)
1950 {
1951 GEN z, nf = alg_get_center(al);
1952 long i;
1953 z = _tablemul(alg_get_relmultable(al), x, y);
1954 for (i=1; i<lg(z); i++)
1955 gel(z,i) = basistoalg(nf,gel(z,i));
1956 return z;
1957 }
1958
1959 /* assumes x and y in algebraic form */
1960 static GEN
algalgmul(GEN al,GEN x,GEN y)1961 algalgmul(GEN al, GEN x, GEN y)
1962 {
1963 switch(alg_type(al))
1964 {
1965 case al_CYCLIC: return algalgmul_cyc(al, x, y);
1966 case al_CSA: return algalgmul_csa(al, x, y);
1967 }
1968 return NULL; /*LCOV_EXCL_LINE*/
1969 }
1970
1971 static GEN
algbasismul(GEN al,GEN x,GEN y)1972 algbasismul(GEN al, GEN x, GEN y)
1973 {
1974 GEN mt = alg_get_multable(al), p = alg_get_char(al);
1975 if (signe(p)) return _tablemul_Fp(mt, x, y, p);
1976 return _tablemul(mt, x, y);
1977 }
1978
1979 /* x[i,]*y. Assume lg(x) > 1 and 0 < i < lgcols(x) */
1980 static GEN
alMrow_alC_mul_i(GEN al,GEN x,GEN y,long i,long lx)1981 alMrow_alC_mul_i(GEN al, GEN x, GEN y, long i, long lx)
1982 {
1983 pari_sp av = avma;
1984 GEN c = algmul(al,gcoeff(x,i,1),gel(y,1)), ZERO;
1985 long k;
1986 ZERO = zerocol(alg_get_absdim(al));
1987 for (k = 2; k < lx; k++)
1988 {
1989 GEN t = algmul(al, gcoeff(x,i,k), gel(y,k));
1990 if (!gequal(t,ZERO)) c = algadd(al, c, t);
1991 }
1992 return gerepilecopy(av, c);
1993 }
1994 /* return x * y, 1 < lx = lg(x), l = lgcols(x) */
1995 static GEN
alM_alC_mul_i(GEN al,GEN x,GEN y,long lx,long l)1996 alM_alC_mul_i(GEN al, GEN x, GEN y, long lx, long l)
1997 {
1998 GEN z = cgetg(l,t_COL);
1999 long i;
2000 for (i=1; i<l; i++) gel(z,i) = alMrow_alC_mul_i(al,x,y,i,lx);
2001 return z;
2002 }
2003 static GEN
alM_mul(GEN al,GEN x,GEN y)2004 alM_mul(GEN al, GEN x, GEN y)
2005 {
2006 long j, l, lx=lg(x), ly=lg(y);
2007 GEN z;
2008 if (ly==1) return cgetg(1,t_MAT);
2009 if (lx != lgcols(y)) pari_err_DIM("alM_mul");
2010 if (lx==1) return zeromat(0, ly-1);
2011 l = lgcols(x); z = cgetg(ly,t_MAT);
2012 for (j=1; j<ly; j++) gel(z,j) = alM_alC_mul_i(al,x,gel(y,j),lx,l);
2013 return z;
2014 }
2015
2016 GEN
algmul(GEN al,GEN x,GEN y)2017 algmul(GEN al, GEN x, GEN y)
2018 {
2019 pari_sp av = avma;
2020 long tx, ty;
2021 checkalg(al);
2022 tx = alg_model(al,x);
2023 ty = alg_model(al,y);
2024 if (tx==al_MATRIX) {
2025 if (ty==al_MATRIX) return alM_mul(al,x,y);
2026 pari_err_TYPE("algmul", y);
2027 }
2028 if (signe(alg_get_char(al))) return algbasismul(al,x,y);
2029 if (tx==al_TRIVIAL) retmkcol(gmul(gel(x,1),gel(y,1)));
2030 if (tx==al_ALGEBRAIC && ty==al_ALGEBRAIC) return algalgmul(al,x,y);
2031 if (tx==al_ALGEBRAIC) x = algalgtobasis(al,x);
2032 if (ty==al_ALGEBRAIC) y = algalgtobasis(al,y);
2033 return gerepileupto(av,algbasismul(al,x,y));
2034 }
2035
2036 GEN
algsqr(GEN al,GEN x)2037 algsqr(GEN al, GEN x)
2038 {
2039 pari_sp av = avma;
2040 long tx;
2041 checkalg(al);
2042 tx = alg_model(al,x);
2043 if (tx==al_MATRIX) return gerepilecopy(av,alM_mul(al,x,x));
2044 if (signe(alg_get_char(al))) return algbasismul(al,x,x);
2045 if (tx==al_TRIVIAL) retmkcol(gsqr(gel(x,1)));
2046 if (tx==al_ALGEBRAIC) return algalgmul(al,x,x);
2047 return gerepileupto(av,algbasismul(al,x,x));
2048 }
2049
2050 static GEN
algmtK2Z_cyc(GEN al,GEN m)2051 algmtK2Z_cyc(GEN al, GEN m)
2052 {
2053 pari_sp av = avma;
2054 GEN nf = alg_get_abssplitting(al), res, mt, rnf = alg_get_splittingfield(al), c, dc;
2055 long n = alg_get_degree(al), N = nf_get_degree(nf), Nn, i, j, i1, j1;
2056 Nn = N*n;
2057 res = zeromatcopy(Nn,Nn);
2058 for (i=0; i<n; i++)
2059 for (j=0; j<n; j++) {
2060 c = gcoeff(m,i+1,j+1);
2061 if (!gequal0(c)) {
2062 c = rnfeltreltoabs(rnf,c);
2063 c = algtobasis(nf,c);
2064 c = Q_remove_denom(c,&dc);
2065 mt = zk_multable(nf,c);
2066 if (dc) mt = ZM_Z_div(mt,dc);
2067 for (i1=1; i1<=N; i1++)
2068 for (j1=1; j1<=N; j1++)
2069 gcoeff(res,i*N+i1,j*N+j1) = gcoeff(mt,i1,j1);
2070 }
2071 }
2072 return gerepilecopy(av,res);
2073 }
2074
2075 static GEN
algmtK2Z_csa(GEN al,GEN m)2076 algmtK2Z_csa(GEN al, GEN m)
2077 {
2078 pari_sp av = avma;
2079 GEN nf = alg_get_center(al), res, mt, c, dc;
2080 long d2 = alg_get_dim(al), n = nf_get_degree(nf), D, i, j, i1, j1;
2081 D = d2*n;
2082 res = zeromatcopy(D,D);
2083 for (i=0; i<d2; i++)
2084 for (j=0; j<d2; j++) {
2085 c = gcoeff(m,i+1,j+1);
2086 if (!gequal0(c)) {
2087 c = algtobasis(nf,c);
2088 c = Q_remove_denom(c,&dc);
2089 mt = zk_multable(nf,c);
2090 if (dc) mt = ZM_Z_div(mt,dc);
2091 for (i1=1; i1<=n; i1++)
2092 for (j1=1; j1<=n; j1++)
2093 gcoeff(res,i*n+i1,j*n+j1) = gcoeff(mt,i1,j1);
2094 }
2095 }
2096 return gerepilecopy(av,res);
2097 }
2098
2099 /* assumes al is a CSA or CYCLIC */
2100 static GEN
algmtK2Z(GEN al,GEN m)2101 algmtK2Z(GEN al, GEN m)
2102 {
2103 switch(alg_type(al))
2104 {
2105 case al_CYCLIC: return algmtK2Z_cyc(al, m);
2106 case al_CSA: return algmtK2Z_csa(al, m);
2107 }
2108 return NULL; /*LCOV_EXCL_LINE*/
2109 }
2110
2111 /* left multiplication table, as a vector space of dimension n over the splitting field (by right multiplication) */
2112 static GEN
algalgmultable_cyc(GEN al,GEN x)2113 algalgmultable_cyc(GEN al, GEN x)
2114 {
2115 pari_sp av = avma;
2116 long n = alg_get_degree(al), i, j;
2117 GEN res, rnf, auts, b, pol;
2118 rnf = alg_get_splittingfield(al);
2119 auts = alg_get_auts(al);
2120 b = alg_get_b(al);
2121 pol = rnf_get_pol(rnf);
2122
2123 res = zeromatcopy(n,n);
2124 for (i=0; i<n; i++)
2125 gcoeff(res,i+1,1) = lift_shallow(rnfbasistoalg(rnf,gel(x,i+1)));
2126
2127 for (i=0; i<n; i++) {
2128 for (j=1; j<=i; j++)
2129 gcoeff(res,i+1,j+1) = gmodulo(poleval(gcoeff(res,i-j+1,1),gel(auts,j)),pol);
2130 for (; j<n; j++)
2131 gcoeff(res,i+1,j+1) = gmodulo(gmul(b,poleval(gcoeff(res,n+i-j+1,1),gel(auts,j))), pol);
2132 }
2133
2134 for (i=0; i<n; i++)
2135 gcoeff(res,i+1,1) = gmodulo(gcoeff(res,i+1,1),pol);
2136
2137 return gerepilecopy(av, res);
2138 }
2139
2140 static GEN
elementmultable(GEN mt,GEN x)2141 elementmultable(GEN mt, GEN x)
2142 {
2143 pari_sp av = avma;
2144 long D = lg(mt)-1, i;
2145 GEN z = NULL;
2146 for (i=1; i<=D; i++)
2147 {
2148 GEN c = gel(x,i);
2149 if (!gequal0(c))
2150 {
2151 GEN M = RgM_Rg_mul(gel(mt,i),c);
2152 z = z? RgM_add(z, M): M;
2153 }
2154 }
2155 if (!z) { set_avma(av); return zeromatcopy(D,D); }
2156 return gerepileupto(av, z);
2157 }
2158 /* mt a t_VEC of Flm modulo m */
2159 static GEN
algbasismultable_Flm(GEN mt,GEN x,ulong m)2160 algbasismultable_Flm(GEN mt, GEN x, ulong m)
2161 {
2162 pari_sp av = avma;
2163 long D = lg(gel(mt,1))-1, i;
2164 GEN z = NULL;
2165 for (i=1; i<=D; i++)
2166 {
2167 ulong c = x[i];
2168 if (c)
2169 {
2170 GEN M = Flm_Fl_mul(gel(mt,i),c, m);
2171 z = z? Flm_add(z, M, m): M;
2172 }
2173 }
2174 if (!z) { set_avma(av); return zero_Flm(D,D); }
2175 return gerepileupto(av, z);
2176 }
2177 static GEN
elementabsmultable_Z(GEN mt,GEN x)2178 elementabsmultable_Z(GEN mt, GEN x)
2179 {
2180 long i, l = lg(x);
2181 GEN z = NULL;
2182 for (i = 1; i < l; i++)
2183 {
2184 GEN c = gel(x,i);
2185 if (signe(c))
2186 {
2187 GEN M = ZM_Z_mul(gel(mt,i),c);
2188 z = z? ZM_add(z, M): M;
2189 }
2190 }
2191 return z;
2192 }
2193 static GEN
elementabsmultable(GEN mt,GEN x)2194 elementabsmultable(GEN mt, GEN x)
2195 {
2196 GEN d, z = elementabsmultable_Z(mt, Q_remove_denom(x,&d));
2197 return (z && d)? ZM_Z_div(z, d): z;
2198 }
2199 static GEN
elementabsmultable_Fp(GEN mt,GEN x,GEN p)2200 elementabsmultable_Fp(GEN mt, GEN x, GEN p)
2201 {
2202 GEN z = elementabsmultable_Z(mt, x);
2203 return z? FpM_red(z, p): z;
2204 }
2205 static GEN
algbasismultable(GEN al,GEN x)2206 algbasismultable(GEN al, GEN x)
2207 {
2208 pari_sp av = avma;
2209 GEN z, p = alg_get_char(al), mt = alg_get_multable(al);
2210 z = signe(p)? elementabsmultable_Fp(mt, x, p): elementabsmultable(mt, x);
2211 if (!z)
2212 {
2213 long D = lg(mt)-1;
2214 set_avma(av); return zeromat(D,D);
2215 }
2216 return gerepileupto(av, z);
2217 }
2218
2219 static GEN
algalgmultable_csa(GEN al,GEN x)2220 algalgmultable_csa(GEN al, GEN x)
2221 {
2222 GEN nf = alg_get_center(al), m;
2223 long i,j;
2224 m = elementmultable(alg_get_relmultable(al), x);
2225 for (i=1; i<lg(m); i++)
2226 for(j=1; j<lg(m); j++)
2227 gcoeff(m,i,j) = basistoalg(nf,gcoeff(m,i,j));
2228 return m;
2229 }
2230
2231 /* assumes x in algebraic form */
2232 static GEN
algalgmultable(GEN al,GEN x)2233 algalgmultable(GEN al, GEN x)
2234 {
2235 switch(alg_type(al))
2236 {
2237 case al_CYCLIC: return algalgmultable_cyc(al, x);
2238 case al_CSA: return algalgmultable_csa(al, x);
2239 }
2240 return NULL; /*LCOV_EXCL_LINE*/
2241 }
2242
2243 /* on the natural basis */
2244 /* assumes x in algebraic form */
2245 static GEN
algZmultable(GEN al,GEN x)2246 algZmultable(GEN al, GEN x) {
2247 pari_sp av = avma;
2248 GEN res = NULL, x0;
2249 long tx = alg_model(al,x);
2250 switch(tx) {
2251 case al_TRIVIAL:
2252 x0 = gel(x,1);
2253 if (typ(x0)==t_POLMOD) x0 = gel(x0,2);
2254 if (typ(x0)==t_POL) x0 = constant_coeff(x0);
2255 res = mkmatcopy(mkcol(x0));
2256 break;
2257 case al_ALGEBRAIC: res = algmtK2Z(al,algalgmultable(al,x)); break;
2258 }
2259 return gerepileupto(av,res);
2260 }
2261
2262 /* x integral */
2263 static GEN
algbasisrightmultable(GEN al,GEN x)2264 algbasisrightmultable(GEN al, GEN x)
2265 {
2266 long N = alg_get_absdim(al), i,j,k;
2267 GEN res = zeromatcopy(N,N), c, mt = alg_get_multable(al), p = alg_get_char(al);
2268 if (gequal0(p)) p = NULL;
2269 for (i=1; i<=N; i++) {
2270 c = gel(x,i);
2271 if (!gequal0(c)) {
2272 for (j=1; j<=N; j++)
2273 for(k=1; k<=N; k++) {
2274 if (p) gcoeff(res,k,j) = Fp_add(gcoeff(res,k,j), Fp_mul(c, gcoeff(gel(mt,j),k,i), p), p);
2275 else gcoeff(res,k,j) = addii(gcoeff(res,k,j), mulii(c, gcoeff(gel(mt,j),k,i)));
2276 }
2277 }
2278 }
2279 return res;
2280 }
2281
2282 /* basis for matrices : 1, E_{i,j} for (i,j)!=(1,1) */
2283 /* index : ijk = ((i-1)*N+j-1)*n + k */
2284 /* square matrices only, coefficients in basis form, shallow function */
2285 static GEN
algmat2basis(GEN al,GEN M)2286 algmat2basis(GEN al, GEN M)
2287 {
2288 long n = alg_get_absdim(al), N = lg(M)-1, i, j, k, ij, ijk;
2289 GEN res, x;
2290 res = zerocol(N*N*n);
2291 for (i=1; i<=N; i++) {
2292 for (j=1, ij=(i-1)*N+1; j<=N; j++, ij++) {
2293 x = gcoeff(M,i,j);
2294 for (k=1, ijk=(ij-1)*n+1; k<=n; k++, ijk++) {
2295 gel(res, ijk) = gel(x, k);
2296 if (i>1 && i==j) gel(res, ijk) = gsub(gel(res,ijk), gel(res,k));
2297 }
2298 }
2299 }
2300
2301 return res;
2302 }
2303
2304 static GEN
algbasis2mat(GEN al,GEN M,long N)2305 algbasis2mat(GEN al, GEN M, long N)
2306 {
2307 long n = alg_get_absdim(al), i, j, k, ij, ijk;
2308 GEN res, x;
2309 res = zeromatcopy(N,N);
2310 for (i=1; i<=N; i++)
2311 for (j=1; j<=N; j++)
2312 gcoeff(res,i,j) = zerocol(n);
2313
2314 for (i=1; i<=N; i++) {
2315 for (j=1, ij=(i-1)*N+1; j<=N; j++, ij++) {
2316 x = gcoeff(res,i,j);
2317 for (k=1, ijk=(ij-1)*n+1; k<=n; k++, ijk++) {
2318 gel(x,k) = gel(M,ijk);
2319 if (i>1 && i==j) gel(x,k) = gadd(gel(x,k), gel(M,k));
2320 }
2321 }
2322 }
2323
2324 return res;
2325 }
2326
2327 static GEN
algmatbasis_ei(GEN al,long ijk,long N)2328 algmatbasis_ei(GEN al, long ijk, long N)
2329 {
2330 long n = alg_get_absdim(al), i, j, k, ij;
2331 GEN res;
2332
2333 res = zeromatcopy(N,N);
2334 for (i=1; i<=N; i++)
2335 for (j=1; j<=N; j++)
2336 gcoeff(res,i,j) = zerocol(n);
2337
2338 k = ijk%n;
2339 if (k==0) k=n;
2340 ij = (ijk-k)/n+1;
2341
2342 if (ij==1) {
2343 for (i=1; i<=N; i++)
2344 gcoeff(res,i,i) = col_ei(n,k);
2345 return res;
2346 }
2347
2348 j = ij%N;
2349 if (j==0) j=N;
2350 i = (ij-j)/N+1;
2351
2352 gcoeff(res,i,j) = col_ei(n,k);
2353 return res;
2354 }
2355
2356 /* FIXME lazy implementation! */
2357 static GEN
algleftmultable_mat(GEN al,GEN M)2358 algleftmultable_mat(GEN al, GEN M)
2359 {
2360 long N = lg(M)-1, n = alg_get_absdim(al), D = N*N*n, j;
2361 GEN res, x, Mx;
2362 if (N == 0) return cgetg(1, t_MAT);
2363 if (N != nbrows(M)) pari_err_DIM("algleftmultable_mat (nonsquare)");
2364 res = cgetg(D+1, t_MAT);
2365 for (j=1; j<=D; j++) {
2366 x = algmatbasis_ei(al, j, N);
2367 Mx = algmul(al, M, x);
2368 gel(res, j) = algmat2basis(al, Mx);
2369 }
2370 return res;
2371 }
2372
2373 /* left multiplication table on integral basis */
2374 static GEN
algleftmultable(GEN al,GEN x)2375 algleftmultable(GEN al, GEN x)
2376 {
2377 pari_sp av = avma;
2378 long tx;
2379 GEN res;
2380
2381 checkalg(al);
2382 tx = alg_model(al,x);
2383 switch(tx) {
2384 case al_TRIVIAL : res = mkmatcopy(mkcol(gel(x,1))); break;
2385 case al_ALGEBRAIC : x = algalgtobasis(al,x);
2386 case al_BASIS : res = algbasismultable(al,x); break;
2387 case al_MATRIX : res = algleftmultable_mat(al,x); break;
2388 default : return NULL; /* LCOV_EXCL_LINE */
2389 }
2390 return gerepileupto(av,res);
2391 }
2392
2393 static GEN
algbasissplittingmatrix_csa(GEN al,GEN x)2394 algbasissplittingmatrix_csa(GEN al, GEN x)
2395 {
2396 long d = alg_get_degree(al), i, j;
2397 GEN rnf = alg_get_splittingfield(al), splba = alg_get_splittingbasis(al), splbainv = alg_get_splittingbasisinv(al), M;
2398 M = algbasismultable(al,x);
2399 M = RgM_mul(M, splba); /* TODO best order ? big matrix /Q vs small matrix /nf */
2400 M = RgM_mul(splbainv, M);
2401 for (i=1; i<=d; i++)
2402 for (j=1; j<=d; j++)
2403 gcoeff(M,i,j) = rnfeltabstorel(rnf, gcoeff(M,i,j));
2404 return M;
2405 }
2406
2407 GEN
algtomatrix(GEN al,GEN x,long abs)2408 algtomatrix(GEN al, GEN x, long abs)
2409 {
2410 pari_sp av = avma;
2411 GEN res = NULL;
2412 long ta, tx, i, j;
2413 checkalg(al);
2414 ta = alg_type(al);
2415 if (abs || ta==al_TABLE) return algleftmultable(al,x);
2416 tx = alg_model(al,x);
2417 if (tx==al_MATRIX) {
2418 if (lg(x) == 1) return cgetg(1, t_MAT);
2419 res = zeromatcopy(nbrows(x),lg(x)-1);
2420 for (j=1; j<lg(x); j++)
2421 for (i=1; i<lgcols(x); i++)
2422 gcoeff(res,i,j) = algtomatrix(al,gcoeff(x,i,j),0);
2423 res = shallowmatconcat(res);
2424 }
2425 else switch(alg_type(al))
2426 {
2427 case al_CYCLIC:
2428 if (tx==al_BASIS) x = algbasistoalg(al,x);
2429 res = algalgmultable(al,x);
2430 break;
2431 case al_CSA:
2432 if (tx==al_ALGEBRAIC) x = algalgtobasis(al,x);
2433 res = algbasissplittingmatrix_csa(al,x);
2434 break;
2435 default:
2436 pari_err_DOMAIN("algtomatrix", "alg_type(al)", "=", stoi(alg_type(al)), stoi(alg_type(al)));
2437 }
2438 return gerepilecopy(av,res);
2439 }
2440
2441 /* x^(-1)*y, NULL if no solution */
2442 static GEN
algdivl_i(GEN al,GEN x,GEN y,long tx,long ty)2443 algdivl_i(GEN al, GEN x, GEN y, long tx, long ty) {
2444 pari_sp av = avma;
2445 GEN res, p = alg_get_char(al), mtx;
2446 if (tx != ty) {
2447 if (tx==al_ALGEBRAIC) { x = algalgtobasis(al,x); tx=al_BASIS; }
2448 if (ty==al_ALGEBRAIC) { y = algalgtobasis(al,y); ty=al_BASIS; }
2449 }
2450 if (ty == al_MATRIX)
2451 {
2452 if (alg_type(al) != al_TABLE) y = algalgtobasis(al,y);
2453 y = algmat2basis(al,y);
2454 }
2455 if (signe(p)) res = FpM_FpC_invimage(algbasismultable(al,x),y,p);
2456 else
2457 {
2458 if (ty==al_ALGEBRAIC) mtx = algalgmultable(al,x);
2459 else mtx = algleftmultable(al,x);
2460 res = inverseimage(mtx,y);
2461 }
2462 if (!res || lg(res)==1) return gc_NULL(av);
2463 if (tx == al_MATRIX) {
2464 res = algbasis2mat(al, res, lg(x)-1);
2465 return gerepilecopy(av,res);
2466 }
2467 return gerepileupto(av,res);
2468 }
2469 static GEN
algdivl_i2(GEN al,GEN x,GEN y)2470 algdivl_i2(GEN al, GEN x, GEN y)
2471 {
2472 long tx, ty;
2473 checkalg(al);
2474 tx = alg_model(al,x);
2475 ty = alg_model(al,y);
2476 if (tx == al_MATRIX) {
2477 if (ty != al_MATRIX) pari_err_TYPE2("\\", x, y);
2478 if (lg(y) == 1) return cgetg(1, t_MAT);
2479 if (lg(x) == 1) return NULL;
2480 if (lgcols(x) != lgcols(y)) pari_err_DIM("algdivl");
2481 if (lg(x) != lgcols(x) || lg(y) != lgcols(y))
2482 pari_err_DIM("algdivl (nonsquare)");
2483 }
2484 return algdivl_i(al,x,y,tx,ty);
2485 }
2486
algdivl(GEN al,GEN x,GEN y)2487 GEN algdivl(GEN al, GEN x, GEN y)
2488 {
2489 GEN z;
2490 z = algdivl_i2(al,x,y);
2491 if (!z) pari_err_INV("algdivl", x);
2492 return z;
2493 }
2494
2495 int
algisdivl(GEN al,GEN x,GEN y,GEN * ptz)2496 algisdivl(GEN al, GEN x, GEN y, GEN* ptz)
2497 {
2498 pari_sp av = avma;
2499 GEN z = algdivl_i2(al,x,y);
2500 if (!z) return gc_bool(av,0);
2501 if (ptz != NULL) *ptz = z;
2502 return 1;
2503 }
2504
2505 static GEN
alginv_i(GEN al,GEN x)2506 alginv_i(GEN al, GEN x)
2507 {
2508 pari_sp av = avma;
2509 GEN res = NULL, p = alg_get_char(al);
2510 long tx = alg_model(al,x), n;
2511 switch(tx) {
2512 case al_TRIVIAL :
2513 if (signe(p)) { res = mkcol(Fp_inv(gel(x,1),p)); break; }
2514 else { res = mkcol(ginv(gel(x,1))); break; }
2515 case al_ALGEBRAIC :
2516 switch(alg_type(al)) {
2517 case al_CYCLIC: n = alg_get_degree(al); break;
2518 case al_CSA: n = alg_get_dim(al); break;
2519 default: return NULL; /* LCOV_EXCL_LINE */
2520 }
2521 res = algdivl_i(al, x, col_ei(n,1), tx, al_ALGEBRAIC); break;
2522 case al_BASIS : res = algdivl_i(al, x, col_ei(alg_get_absdim(al),1), tx,
2523 al_BASIS); break;
2524 case al_MATRIX :
2525 n = lg(x)-1;
2526 if (n==0) return cgetg(1, t_MAT);
2527 if (n != nbrows(x)) pari_err_DIM("alginv_i (nonsquare)");
2528 res = algdivl_i(al, x, col_ei(n*n*alg_get_absdim(al),1), tx, al_BASIS);
2529 /* cheat on type because wrong dimension */
2530 }
2531 if (!res) return gc_NULL(av);
2532 return gerepilecopy(av,res);
2533 }
2534 GEN
alginv(GEN al,GEN x)2535 alginv(GEN al, GEN x)
2536 {
2537 GEN z;
2538 checkalg(al);
2539 z = alginv_i(al,x);
2540 if (!z) pari_err_INV("alginv", x);
2541 return z;
2542 }
2543
2544 int
algisinv(GEN al,GEN x,GEN * ptix)2545 algisinv(GEN al, GEN x, GEN* ptix)
2546 {
2547 pari_sp av = avma;
2548 GEN ix;
2549 checkalg(al);
2550 ix = alginv_i(al,x);
2551 if (!ix) return gc_bool(av,0);
2552 if (ptix != NULL) *ptix = ix;
2553 return 1;
2554 }
2555
2556 /* x*y^(-1) */
2557 GEN
algdivr(GEN al,GEN x,GEN y)2558 algdivr(GEN al, GEN x, GEN y) { return algmul(al, x, alginv(al, y)); }
2559
_mul(void * data,GEN x,GEN y)2560 static GEN _mul(void* data, GEN x, GEN y) { return algmul((GEN)data,x,y); }
_sqr(void * data,GEN x)2561 static GEN _sqr(void* data, GEN x) { return algsqr((GEN)data,x); }
2562
2563 static GEN
algmatid(GEN al,long N)2564 algmatid(GEN al, long N)
2565 {
2566 long n = alg_get_absdim(al), i, j;
2567 GEN res, one, zero;
2568
2569 res = zeromatcopy(N,N);
2570 one = col_ei(n,1);
2571 zero = zerocol(n);
2572 for (i=1; i<=N; i++)
2573 for (j=1; j<=N; j++)
2574 gcoeff(res,i,j) = i==j ? one : zero;
2575 return res;
2576 }
2577
2578 GEN
algpow(GEN al,GEN x,GEN n)2579 algpow(GEN al, GEN x, GEN n)
2580 {
2581 pari_sp av = avma;
2582 GEN res;
2583 checkalg(al);
2584 switch(signe(n)) {
2585 case 0:
2586 if (alg_model(al,x) == al_MATRIX)
2587 res = algmatid(al,lg(x)-1);
2588 else
2589 res = col_ei(alg_get_absdim(al),1);
2590 return res;
2591 case 1:
2592 res = gen_pow_i(x, n, (void*)al, _sqr, _mul); break;
2593 default: /* -1 */
2594 res = gen_pow_i(alginv(al,x), gneg(n), (void*)al, _sqr, _mul);
2595 }
2596 return gerepilecopy(av,res);
2597 }
2598
2599 static GEN
algredcharpoly_i(GEN al,GEN x,long v)2600 algredcharpoly_i(GEN al, GEN x, long v)
2601 {
2602 GEN rnf = alg_get_splittingfield(al);
2603 GEN cp = charpoly(algtomatrix(al,x,0),v);
2604 long i, m = lg(cp);
2605 for (i=2; i<m; i++) gel(cp,i) = rnfeltdown(rnf, gel(cp,i));
2606 return cp;
2607 }
2608
2609 /* assumes al is CSA or CYCLIC */
2610 static GEN
algredcharpoly(GEN al,GEN x,long v)2611 algredcharpoly(GEN al, GEN x, long v)
2612 {
2613 pari_sp av = avma;
2614 long w = gvar(rnf_get_pol(alg_get_center(al)));
2615 if (varncmp(v,w)>=0) pari_err_PRIORITY("algredcharpoly",pol_x(v),">=",w);
2616 switch(alg_type(al))
2617 {
2618 case al_CYCLIC:
2619 case al_CSA:
2620 return gerepileupto(av, algredcharpoly_i(al, x, v));
2621 }
2622 return NULL; /*LCOV_EXCL_LINE*/
2623 }
2624
2625 static GEN
algbasischarpoly(GEN al,GEN x,long v)2626 algbasischarpoly(GEN al, GEN x, long v)
2627 {
2628 pari_sp av = avma;
2629 GEN p = alg_get_char(al), mx;
2630 if (alg_model(al,x) == al_MATRIX) mx = algleftmultable_mat(al,x);
2631 else mx = algbasismultable(al,x);
2632 if (signe(p)) {
2633 GEN res = FpM_charpoly(mx,p);
2634 setvarn(res,v);
2635 return gerepileupto(av, res);
2636 }
2637 return gerepileupto(av, charpoly(mx,v));
2638 }
2639
2640 GEN
algcharpoly(GEN al,GEN x,long v,long abs)2641 algcharpoly(GEN al, GEN x, long v, long abs)
2642 {
2643 checkalg(al);
2644 if (v<0) v=0;
2645
2646 /* gneg(x[1]) left on stack */
2647 if (alg_model(al,x) == al_TRIVIAL) {
2648 GEN p = alg_get_char(al);
2649 if (signe(p)) return deg1pol(gen_1,Fp_neg(gel(x,1),p),v);
2650 return deg1pol(gen_1,gneg(gel(x,1)),v);
2651 }
2652
2653 switch(alg_type(al)) {
2654 case al_CYCLIC: case al_CSA:
2655 if (abs)
2656 {
2657 if (alg_model(al,x)==al_ALGEBRAIC) x = algalgtobasis(al,x);
2658 }
2659 else return algredcharpoly(al,x,v);
2660 case al_TABLE: return algbasischarpoly(al,x,v);
2661 default : return NULL; /* LCOV_EXCL_LINE */
2662 }
2663 }
2664
2665 /* assumes x in basis form */
2666 static GEN
algabstrace(GEN al,GEN x)2667 algabstrace(GEN al, GEN x)
2668 {
2669 pari_sp av = avma;
2670 GEN res = NULL, p = alg_get_char(al);
2671 if (signe(p)) return FpV_dotproduct(x, alg_get_tracebasis(al), p);
2672 switch(alg_model(al,x)) {
2673 case al_TRIVIAL: return gcopy(gel(x,1)); break;
2674 case al_BASIS: res = RgV_dotproduct(x, alg_get_tracebasis(al)); break;
2675 }
2676 return gerepileupto(av,res);
2677 }
2678
2679 static GEN
algredtrace(GEN al,GEN x)2680 algredtrace(GEN al, GEN x)
2681 {
2682 pari_sp av = avma;
2683 GEN res = NULL;
2684 switch(alg_model(al,x)) {
2685 case al_TRIVIAL: return gcopy(gel(x,1)); break;
2686 case al_BASIS: return algredtrace(al, algbasistoalg(al,x));
2687 /* TODO precompute too? */
2688 case al_ALGEBRAIC:
2689 switch(alg_type(al))
2690 {
2691 case al_CYCLIC:
2692 res = rnfelttrace(alg_get_splittingfield(al),gel(x,1));
2693 break;
2694 case al_CSA:
2695 res = gtrace(algalgmultable_csa(al,x));
2696 res = gdiv(res, stoi(alg_get_degree(al)));
2697 break;
2698 default: return NULL; /* LCOV_EXCL_LINE */
2699 }
2700 }
2701 return gerepileupto(av,res);
2702 }
2703
2704 static GEN
algtrace_mat(GEN al,GEN M,long abs)2705 algtrace_mat(GEN al, GEN M, long abs) {
2706 pari_sp av = avma;
2707 long N = lg(M)-1, i;
2708 GEN res, p = alg_get_char(al);
2709 if (N == 0) return gen_0;
2710 if (N != nbrows(M)) pari_err_DIM("algtrace_mat (nonsquare)");
2711
2712 if (!signe(p)) p = NULL;
2713 res = algtrace(al, gcoeff(M,1,1), abs);
2714 for (i=2; i<=N; i++) {
2715 if (p) res = Fp_add(res, algtrace(al,gcoeff(M,i,i),abs), p);
2716 else res = gadd(res, algtrace(al,gcoeff(M,i,i),abs));
2717 }
2718 if (abs || alg_type(al) == al_TABLE) res = gmulgs(res, N); /* absolute trace */
2719 return gerepileupto(av, res);
2720 }
2721
2722 GEN
algtrace(GEN al,GEN x,long abs)2723 algtrace(GEN al, GEN x, long abs)
2724 {
2725 checkalg(al);
2726 if (alg_model(al,x) == al_MATRIX) return algtrace_mat(al,x,abs);
2727 switch(alg_type(al)) {
2728 case al_CYCLIC: case al_CSA:
2729 if (!abs) return algredtrace(al,x);
2730 if (alg_model(al,x)==al_ALGEBRAIC) x = algalgtobasis(al,x);
2731 case al_TABLE: return algabstrace(al,x);
2732 default : return NULL; /* LCOV_EXCL_LINE */
2733 }
2734 }
2735
2736 static GEN
ZM_trace(GEN x)2737 ZM_trace(GEN x)
2738 {
2739 long i, lx = lg(x);
2740 GEN t;
2741 if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
2742 t = gcoeff(x,1,1);
2743 for (i = 2; i < lx; i++) t = addii(t, gcoeff(x,i,i));
2744 return t;
2745 }
2746 static GEN
FpM_trace(GEN x,GEN p)2747 FpM_trace(GEN x, GEN p)
2748 {
2749 long i, lx = lg(x);
2750 GEN t;
2751 if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
2752 t = gcoeff(x,1,1);
2753 for (i = 2; i < lx; i++) t = Fp_add(t, gcoeff(x,i,i), p);
2754 return t;
2755 }
2756
2757 static GEN
algtracebasis(GEN al)2758 algtracebasis(GEN al)
2759 {
2760 pari_sp av = avma;
2761 GEN mt = alg_get_multable(al), p = alg_get_char(al);
2762 long i, l = lg(mt);
2763 GEN v = cgetg(l, t_VEC);
2764 if (signe(p)) for (i=1; i < l; i++) gel(v,i) = FpM_trace(gel(mt,i), p);
2765 else for (i=1; i < l; i++) gel(v,i) = ZM_trace(gel(mt,i));
2766 return gerepileupto(av,v);
2767 }
2768
2769 /* Assume: i > 0, expo := p^i <= absdim, x contained in I_{i-1} given by mult
2770 * table modulo modu=p^(i+1). Return Tr(x^(p^i)) mod modu */
2771 static ulong
algtracei(GEN mt,ulong p,ulong expo,ulong modu)2772 algtracei(GEN mt, ulong p, ulong expo, ulong modu)
2773 {
2774 pari_sp av = avma;
2775 long j, l = lg(mt);
2776 ulong tr = 0;
2777 mt = Flm_powu(mt,expo,modu);
2778 for (j=1; j<l; j++) tr += ucoeff(mt,j,j);
2779 return gc_ulong(av, (tr/expo) % p);
2780 }
2781
2782 GEN
algnorm(GEN al,GEN x,long abs)2783 algnorm(GEN al, GEN x, long abs)
2784 {
2785 pari_sp av = avma;
2786 long tx;
2787 GEN p, rnf, res, mx;
2788 checkalg(al);
2789 p = alg_get_char(al);
2790 tx = alg_model(al,x);
2791 if (signe(p)) {
2792 if (tx == al_MATRIX) mx = algleftmultable_mat(al,x);
2793 else mx = algbasismultable(al,x);
2794 return gerepileupto(av, FpM_det(mx,p));
2795 }
2796 if (tx == al_TRIVIAL) return gcopy(gel(x,1));
2797
2798 switch(alg_type(al)) {
2799 case al_CYCLIC: case al_CSA:
2800 if (abs)
2801 {
2802 if (alg_model(al,x)==al_ALGEBRAIC) x = algalgtobasis(al,x);
2803 }
2804 else
2805 {
2806 rnf = alg_get_splittingfield(al);
2807 res = rnfeltdown(rnf, det(algtomatrix(al,x,0)));
2808 break;
2809 }
2810 case al_TABLE:
2811 if (tx == al_MATRIX) mx = algleftmultable_mat(al,x);
2812 else mx = algbasismultable(al,x);
2813 res = det(mx);
2814 break;
2815 default: return NULL; /* LCOV_EXCL_LINE */
2816 }
2817 return gerepileupto(av, res);
2818 }
2819
2820 static GEN
algalgtonat_cyc(GEN al,GEN x)2821 algalgtonat_cyc(GEN al, GEN x)
2822 {
2823 pari_sp av = avma;
2824 GEN nf = alg_get_abssplitting(al), rnf = alg_get_splittingfield(al), res, c;
2825 long n = alg_get_degree(al), N = nf_get_degree(nf), i, i1;
2826 res = zerocol(N*n);
2827 for (i=0; i<n; i++) {
2828 c = gel(x,i+1);
2829 c = rnfeltreltoabs(rnf,c);
2830 if (!gequal0(c)) {
2831 c = algtobasis(nf,c);
2832 for (i1=1; i1<=N; i1++) gel(res,i*N+i1) = gel(c,i1);
2833 }
2834 }
2835 return gerepilecopy(av, res);
2836 }
2837
2838 static GEN
algalgtonat_csa(GEN al,GEN x)2839 algalgtonat_csa(GEN al, GEN x)
2840 {
2841 pari_sp av = avma;
2842 GEN nf = alg_get_center(al), res, c;
2843 long d2 = alg_get_dim(al), n = nf_get_degree(nf), i, i1;
2844 res = zerocol(d2*n);
2845 for (i=0; i<d2; i++) {
2846 c = gel(x,i+1);
2847 if (!gequal0(c)) {
2848 c = algtobasis(nf,c);
2849 for (i1=1; i1<=n; i1++) gel(res,i*n+i1) = gel(c,i1);
2850 }
2851 }
2852 return gerepilecopy(av, res);
2853 }
2854
2855 /* assumes al CSA or CYCLIC */
2856 static GEN
algalgtonat(GEN al,GEN x)2857 algalgtonat(GEN al, GEN x)
2858 {
2859 switch(alg_type(al))
2860 {
2861 case al_CYCLIC: return algalgtonat_cyc(al, x);
2862 case al_CSA: return algalgtonat_csa(al, x);
2863 }
2864 return NULL; /*LCOV_EXCL_LINE*/
2865 }
2866
2867 static GEN
algnattoalg_cyc(GEN al,GEN x)2868 algnattoalg_cyc(GEN al, GEN x)
2869 {
2870 pari_sp av = avma;
2871 GEN nf = alg_get_abssplitting(al), rnf = alg_get_splittingfield(al), res, c;
2872 long n = alg_get_degree(al), N = nf_get_degree(nf), i, i1;
2873 res = zerocol(n);
2874 c = zerocol(N);
2875 for (i=0; i<n; i++) {
2876 for (i1=1; i1<=N; i1++) gel(c,i1) = gel(x,i*N+i1);
2877 gel(res,i+1) = rnfeltabstorel(rnf,basistoalg(nf,c));
2878 }
2879 return gerepilecopy(av, res);
2880 }
2881
2882 static GEN
algnattoalg_csa(GEN al,GEN x)2883 algnattoalg_csa(GEN al, GEN x)
2884 {
2885 pari_sp av = avma;
2886 GEN nf = alg_get_center(al), res, c;
2887 long d2 = alg_get_dim(al), n = nf_get_degree(nf), i, i1;
2888 res = zerocol(d2);
2889 c = zerocol(n);
2890 for (i=0; i<d2; i++) {
2891 for (i1=1; i1<=n; i1++) gel(c,i1) = gel(x,i*n+i1);
2892 gel(res,i+1) = basistoalg(nf,c);
2893 }
2894 return gerepilecopy(av, res);
2895 }
2896
2897 /* assumes al CSA or CYCLIC */
2898 static GEN
algnattoalg(GEN al,GEN x)2899 algnattoalg(GEN al, GEN x)
2900 {
2901 switch(alg_type(al))
2902 {
2903 case al_CYCLIC: return algnattoalg_cyc(al, x);
2904 case al_CSA: return algnattoalg_csa(al, x);
2905 }
2906 return NULL; /*LCOV_EXCL_LINE*/
2907 }
2908
2909 static GEN
algalgtobasis_mat(GEN al,GEN x)2910 algalgtobasis_mat(GEN al, GEN x) /* componentwise */
2911 {
2912 pari_sp av = avma;
2913 long lx, lxj, i, j;
2914 GEN res;
2915 lx = lg(x);
2916 res = cgetg(lx, t_MAT);
2917 for (j=1; j<lx; j++) {
2918 lxj = lg(gel(x,j));
2919 gel(res,j) = cgetg(lxj, t_COL);
2920 for (i=1; i<lxj; i++)
2921 gcoeff(res,i,j) = algalgtobasis(al,gcoeff(x,i,j));
2922 }
2923 return gerepilecopy(av,res);
2924 }
2925 GEN
algalgtobasis(GEN al,GEN x)2926 algalgtobasis(GEN al, GEN x)
2927 {
2928 pari_sp av;
2929 long tx;
2930 checkalg(al);
2931 if (alg_type(al) == al_TABLE) pari_err_TYPE("algalgtobasis [use alginit]", al);
2932 tx = alg_model(al,x);
2933 if (tx==al_BASIS) return gcopy(x);
2934 if (tx==al_MATRIX) return algalgtobasis_mat(al,x);
2935 av = avma;
2936 x = algalgtonat(al,x);
2937 x = RgM_RgC_mul(alg_get_invbasis(al),x);
2938 return gerepileupto(av, x);
2939 }
2940
2941 static GEN
algbasistoalg_mat(GEN al,GEN x)2942 algbasistoalg_mat(GEN al, GEN x) /* componentwise */
2943 {
2944 long j, lx = lg(x);
2945 GEN res = cgetg(lx, t_MAT);
2946 for (j=1; j<lx; j++) {
2947 long i, lxj = lg(gel(x,j));
2948 gel(res,j) = cgetg(lxj, t_COL);
2949 for (i=1; i<lxj; i++) gcoeff(res,i,j) = algbasistoalg(al,gcoeff(x,i,j));
2950 }
2951 return res;
2952 }
2953 GEN
algbasistoalg(GEN al,GEN x)2954 algbasistoalg(GEN al, GEN x)
2955 {
2956 pari_sp av;
2957 long tx;
2958 checkalg(al);
2959 if (alg_type(al) == al_TABLE) pari_err_TYPE("algbasistoalg [use alginit]", al);
2960 tx = alg_model(al,x);
2961 if (tx==al_ALGEBRAIC) return gcopy(x);
2962 if (tx==al_MATRIX) return algbasistoalg_mat(al,x);
2963 av = avma;
2964 x = RgM_RgC_mul(alg_get_basis(al),x);
2965 x = algnattoalg(al,x);
2966 return gerepileupto(av, x);
2967 }
2968
2969 GEN
algrandom(GEN al,GEN b)2970 algrandom(GEN al, GEN b)
2971 {
2972 GEN res, p, N;
2973 long i, n;
2974 if (typ(b) != t_INT) pari_err_TYPE("algrandom",b);
2975 if (signe(b)<0) pari_err_DOMAIN("algrandom", "b", "<", gen_0, b);
2976 checkalg(al);
2977 n = alg_get_absdim(al);
2978 N = addiu(shifti(b,1), 1); /* left on stack */
2979 p = alg_get_char(al); if (!signe(p)) p = NULL;
2980 res = cgetg(n+1,t_COL);
2981 for (i=1; i<= n; i++)
2982 {
2983 pari_sp av = avma;
2984 GEN t = subii(randomi(N),b);
2985 if (p) t = modii(t, p);
2986 gel(res,i) = gerepileuptoint(av, t);
2987 }
2988 return res;
2989 }
2990
2991 /* Assumes pol has coefficients in the same ring as the COL x; x either
2992 * in basis or algebraic form or [x,mx] where mx is the mult. table of x.
2993 TODO more general version: pol with coeffs in center and x in basis form */
2994 GEN
algpoleval(GEN al,GEN pol,GEN x)2995 algpoleval(GEN al, GEN pol, GEN x)
2996 {
2997 pari_sp av = avma;
2998 GEN p, mx = NULL, res;
2999 long i;
3000 checkalg(al);
3001 p = alg_get_char(al);
3002 if (typ(pol) != t_POL) pari_err_TYPE("algpoleval", pol);
3003 if (typ(x) == t_VEC)
3004 {
3005 if (lg(x) != 3) pari_err_TYPE("algpoleval [vector must be of length 2]", x);
3006 mx = gel(x,2);
3007 x = gel(x,1);
3008 if (typ(mx)!=t_MAT || !gequal(x,gel(mx,1)))
3009 pari_err_TYPE("algpoleval [mx must be the multiplication table of x]", mx);
3010 }
3011 else
3012 {
3013 switch(alg_model(al,x))
3014 {
3015 case al_ALGEBRAIC: mx = algalgmultable(al,x); break;
3016 case al_BASIS: if (!RgX_is_QX(pol))
3017 pari_err_IMPL("algpoleval with x in basis form and pol not in Q[x]");
3018 case al_TRIVIAL: mx = algbasismultable(al,x); break;
3019 default: pari_err_TYPE("algpoleval", x);
3020 }
3021 }
3022 res = zerocol(lg(mx)-1);
3023 if (signe(p)) {
3024 for (i=lg(pol)-1; i>1; i--)
3025 {
3026 gel(res,1) = Fp_add(gel(res,1), gel(pol,i), p);
3027 if (i>2) res = FpM_FpC_mul(mx, res, p);
3028 }
3029 }
3030 else {
3031 for (i=lg(pol)-1; i>1; i--)
3032 {
3033 gel(res,1) = gadd(gel(res,1), gel(pol,i));
3034 if (i>2) res = RgM_RgC_mul(mx, res);
3035 }
3036 }
3037 return gerepileupto(av, res);
3038 }
3039
3040 /** GRUNWALD-WANG **/
3041 /*
3042 Song Wang's PhD thesis (pdf pages)
3043 p.25 definition of chi_b. K^Ker(chi_b) = K(b^(1/m))
3044 p.26 bound on the conductor (also Cohen adv. GTM 193 p.166)
3045 p.21 & p.34 description special case, also on wikipedia:
3046 http://en.wikipedia.org/wiki/Grunwald%E2%80%93Wang_theorem#Special_fields
3047 p.77 Kummer case
3048 */
3049
3050 /* n > 0. Is n = 2^k ? */
3051 static int
uispow2(ulong n)3052 uispow2(ulong n) { return !(n &(n-1)); }
3053
3054 static GEN
get_phi0(GEN bnr,GEN Lpr,GEN Ld,GEN pl,long * pr,long * pn)3055 get_phi0(GEN bnr, GEN Lpr, GEN Ld, GEN pl, long *pr, long *pn)
3056 {
3057 const long NTRY = 10; /* FIXME: magic constant */
3058 const long n = (lg(Ld)==1)? 2: vecsmall_max(Ld);
3059 GEN S = bnr_get_cyc(bnr);
3060 GEN Sst, G, globGmod, loc, X, Rglob, Rloc, H, U, Lconj;
3061 long i, j, r, nbfrob, nbloc, nz, t;
3062
3063 *pn = n;
3064 *pr = r = lg(S)-1;
3065 if (!r) return NULL;
3066 Lconj = NULL;
3067 nbloc = nbfrob = lg(Lpr)-1;
3068 if (uispow2(n))
3069 {
3070 long l = lg(pl), k = 1;
3071 GEN real = cgetg(l, t_VECSMALL);
3072 for (i=1; i<l; i++)
3073 if (pl[i]==-1) real[k++] = i;
3074 if (k > 1)
3075 {
3076 GEN nf = bnr_get_nf(bnr), I = bid_get_fact(bnr_get_bid(bnr));
3077 GEN v, y, C = idealchineseinit(bnr, I);
3078 long r1 = nf_get_r1(nf), n = nbrows(I);
3079 nbloc += k-1;
3080 Lconj = cgetg(k, t_VEC);
3081 v = const_vecsmall(r1,1);
3082 y = const_vec(n, gen_1);
3083 for (i = 1; i < k; i++)
3084 {
3085 v[i] = -1; gel(Lconj,i) = idealchinese(nf,mkvec2(C,v),y);
3086 v[i] = 1;
3087 }
3088 }
3089 }
3090
3091 /* compute Z/n-dual */
3092 Sst = cgetg(r+1, t_VECSMALL);
3093 for (i=1; i<=r; i++) Sst[i] = ugcdiu(gel(S,i), n);
3094 if (Sst[1] != n) return NULL;
3095
3096 globGmod = cgetg(r+1,t_MAT);
3097 G = cgetg(r+1,t_VECSMALL);
3098 for (i=1; i<=r; i++)
3099 {
3100 G[i] = n / Sst[i]; /* pairing between S and Sst */
3101 gel(globGmod,i) = cgetg(nbloc+1,t_VECSMALL);
3102 }
3103
3104 /* compute images of Frobenius elements (and complex conjugation) */
3105 loc = cgetg(nbloc+1,t_VECSMALL);
3106 for (i=1; i<=nbloc; i++) {
3107 long L;
3108 if (i<=nbfrob)
3109 {
3110 X = gel(Lpr,i);
3111 L = Ld[i];
3112 }
3113 else
3114 { /* X = 1 (mod f), sigma_i(x) < 0, positive at all other real places */
3115 X = gel(Lconj,i-nbfrob);
3116 L = 2;
3117 }
3118 X = ZV_to_Flv(isprincipalray(bnr,X), n);
3119 for (nz=0,j=1; j<=r; j++)
3120 {
3121 ulong c = (X[j] * G[j]) % L;
3122 ucoeff(globGmod,i,j) = c;
3123 if (c) nz = 1;
3124 }
3125 if (!nz) return NULL;
3126 loc[i] = L;
3127 }
3128
3129 /* try some random elements in the dual */
3130 Rglob = cgetg(r+1,t_VECSMALL);
3131 for (t=0; t<NTRY; t++) {
3132 for (j=1; j<=r; j++) Rglob[j] = random_Fl(Sst[j]);
3133 Rloc = zm_zc_mul(globGmod,Rglob);
3134 for (i=1; i<=nbloc; i++)
3135 if (Rloc[i] % loc[i] == 0) break;
3136 if (i > nbloc)
3137 return zv_to_ZV(Rglob);
3138 }
3139
3140 /* try to realize some random elements of the product of the local duals */
3141 H = ZM_hnfall_i(shallowconcat(zm_to_ZM(globGmod),
3142 diagonal_shallow(zv_to_ZV(loc))), &U, 2);
3143 /* H,U nbloc x nbloc */
3144 Rloc = cgetg(nbloc+1,t_COL);
3145 for (t=0; t<NTRY; t++) {
3146 /* nonzero random coordinate */ /* TODO add special case ? */
3147 for (i=1; i<=nbloc; i++) gel(Rloc,i) = stoi(1 + random_Fl(loc[i]-1));
3148 Rglob = hnf_invimage(H, Rloc);
3149 if (Rglob)
3150 {
3151 Rglob = ZM_ZC_mul(U,Rglob);
3152 return vecslice(Rglob,1,r);
3153 }
3154 }
3155 return NULL;
3156 }
3157
3158 static GEN
bnrgwsearch(GEN bnr,GEN Lpr,GEN Ld,GEN pl)3159 bnrgwsearch(GEN bnr, GEN Lpr, GEN Ld, GEN pl)
3160 {
3161 pari_sp av = avma;
3162 long n, r;
3163 GEN phi0 = get_phi0(bnr,Lpr,Ld,pl, &r,&n), gn, v, H,U;
3164 if (!phi0) return gc_const(av, gen_0);
3165 gn = stoi(n);
3166 /* compute kernel of phi0 */
3167 v = ZV_extgcd(shallowconcat(phi0, gn));
3168 U = vecslice(gel(v,2), 1,r);
3169 H = ZM_hnfmodid(rowslice(U, 1,r), gn);
3170 return gerepileupto(av, H);
3171 }
3172
3173 GEN
bnfgwgeneric(GEN bnf,GEN Lpr,GEN Ld,GEN pl,long var)3174 bnfgwgeneric(GEN bnf, GEN Lpr, GEN Ld, GEN pl, long var)
3175 {
3176 pari_sp av = avma;
3177 const long n = (lg(Ld)==1)? 2: vecsmall_max(Ld);
3178 forprime_t S;
3179 GEN bnr = NULL, ideal = gen_1, nf, dec, H = gen_0, finf, pol;
3180 ulong ell, p;
3181 long deg, i, degell;
3182 (void)uisprimepower(n, &ell);
3183 nf = bnf_get_nf(bnf);
3184 deg = nf_get_degree(nf);
3185 degell = ugcd(deg,ell-1);
3186 finf = cgetg(lg(pl),t_VEC);
3187 for (i=1; i<lg(pl); i++) gel(finf,i) = pl[i]==-1 ? gen_1 : gen_0;
3188
3189 u_forprime_init(&S, 2, ULONG_MAX);
3190 while ((p = u_forprime_next(&S))) {
3191 if (Fl_powu(p % ell, degell, ell) != 1) continue; /* ell | p^deg-1 ? */
3192 dec = idealprimedec(nf, utoipos(p));
3193 for (i=1; i<lg(dec); i++) {
3194 GEN pp = gel(dec,i);
3195 if (RgV_isin(Lpr,pp)) continue;
3196 /* TODO also accept the prime ideals at which there is a condition
3197 * (use local Artin)? */
3198 if (smodis(idealnorm(nf,pp),ell) != 1) continue; /* ell | N(pp)-1 ? */
3199 ideal = idealmul(bnf,ideal,pp);
3200 /* TODO: give factorization ? */
3201 bnr = Buchray(bnf, mkvec2(ideal,finf), nf_INIT);
3202 H = bnrgwsearch(bnr,Lpr,Ld,pl);
3203 if (H != gen_0)
3204 {
3205 pol = rnfkummer(bnr,H,nf_get_prec(nf));
3206 setvarn(pol, var);
3207 return gerepileupto(av,pol);
3208 }
3209 }
3210 }
3211 pari_err_BUG("bnfgwgeneric (no suitable p)"); /*LCOV_EXCL_LINE*/
3212 return NULL;/*LCOV_EXCL_LINE*/
3213 }
3214
3215 /* no garbage collection */
3216 static GEN
localextdeg(GEN nf,GEN pr,GEN cnd,long d,long ell,long n)3217 localextdeg(GEN nf, GEN pr, GEN cnd, long d, long ell, long n)
3218 {
3219 long g = n/d;
3220 GEN res, modpr, ppr = pr, T, p, gen, k;
3221 if (d==1) return gen_1;
3222 if (equalsi(ell,pr_get_p(pr))) { /* ell == p */
3223 res = nfadd(nf, gen_1, pr_get_gen(pr));
3224 res = nfpowmodideal(nf, res, stoi(g), cnd);
3225 }
3226 else { /* ell != p */
3227 k = powis(stoi(ell),Z_lval(subiu(pr_norm(pr),1),ell));
3228 k = divis(k,g);
3229 modpr = nf_to_Fq_init(nf, &ppr, &T, &p);
3230 (void)Fq_sqrtn(gen_1,k,T,p,&gen);
3231 res = Fq_to_nf(gen, modpr);
3232 }
3233 return res;
3234 }
3235
3236 /* Ld[i] must be nontrivial powers of the same prime ell */
3237 /* pl : -1 at real places at which the extention must ramify, 0 elsewhere */
3238 GEN
nfgwkummer(GEN nf,GEN Lpr,GEN Ld,GEN pl,long var)3239 nfgwkummer(GEN nf, GEN Lpr, GEN Ld, GEN pl, long var)
3240 {
3241 const long n = (lg(Ld)==1)? 2: vecsmall_max(Ld);
3242 pari_sp av = avma;
3243 ulong ell;
3244 long i, v;
3245 GEN cnd, y, x, pol;
3246 v = uisprimepower(n, &ell);
3247 cnd = zeromatcopy(lg(Lpr)-1,2);
3248
3249 y = vec_ei(lg(Lpr)-1,1);
3250 for (i=1; i<lg(Lpr); i++) {
3251 GEN pr = gel(Lpr,i), p = pr_get_p(pr), E;
3252 long e = pr_get_e(pr);
3253 gcoeff(cnd,i,1) = pr;
3254
3255 if (!absequalui(ell,p))
3256 E = gen_1;
3257 else
3258 E = addui(1 + v*e, divsi(e,subiu(p,1)));
3259 gcoeff(cnd,i,2) = E;
3260 gel(y,i) = localextdeg(nf, pr, idealpow(nf,pr,E), Ld[i], ell, n);
3261 }
3262
3263 /* TODO use a factoredextchinese to ease computations afterwards ? */
3264 x = idealchinese(nf, mkvec2(cnd,pl), y);
3265 x = basistoalg(nf,x);
3266 pol = gsub(gpowgs(pol_x(var),n),x);
3267
3268 return gerepileupto(av,pol);
3269 }
3270
3271 static GEN
get_vecsmall(GEN v)3272 get_vecsmall(GEN v)
3273 {
3274 switch(typ(v))
3275 {
3276 case t_VECSMALL: return v;
3277 case t_VEC: if (RgV_is_ZV(v)) return ZV_to_zv(v);
3278 }
3279 pari_err_TYPE("nfgrunwaldwang",v);
3280 return NULL;/*LCOV_EXCL_LINE*/
3281 }
3282 GEN
nfgrunwaldwang(GEN nf0,GEN Lpr,GEN Ld,GEN pl,long var)3283 nfgrunwaldwang(GEN nf0, GEN Lpr, GEN Ld, GEN pl, long var)
3284 {
3285 ulong n, ell, ell2;
3286 pari_sp av = avma;
3287 GEN nf, bnf;
3288 long t, w, i, vnf;
3289
3290 if (var < 0) var = 0;
3291 nf = get_nf(nf0,&t);
3292 if (!nf) pari_err_TYPE("nfgrunwaldwang",nf0);
3293 vnf = nf_get_varn(nf);
3294 if (varncmp(var, vnf) >= 0)
3295 pari_err_PRIORITY("nfgrunwaldwang", pol_x(var), ">=", vnf);
3296 if (typ(Lpr) != t_VEC) pari_err_TYPE("nfgrunwaldwang",Lpr);
3297 if (lg(Lpr) != lg(Ld)) pari_err_DIM("nfgrunwaldwang [#Lpr != #Ld]");
3298 if (nf_get_degree(nf)==1) Lpr = shallowcopy(Lpr);
3299 for (i=1; i<lg(Lpr); i++) {
3300 GEN pr = gel(Lpr,i);
3301 if (nf_get_degree(nf)==1 && typ(pr)==t_INT)
3302 gel(Lpr,i) = gel(idealprimedec(nf,pr), 1);
3303 else checkprid(pr);
3304 }
3305 if (lg(pl)-1 != nf_get_r1(nf))
3306 pari_err_DOMAIN("nfgrunwaldwang [pl should have r1 components]", "#pl",
3307 "!=", stoi(nf_get_r1(nf)), stoi(lg(pl)-1));
3308
3309 Ld = get_vecsmall(Ld);
3310 pl = get_vecsmall(pl);
3311 bnf = get_bnf(nf0,&t);
3312 n = (lg(Ld)==1)? 2: vecsmall_max(Ld);
3313
3314 if (!uisprimepower(n, &ell))
3315 pari_err_IMPL("nfgrunwaldwang for non prime-power local degrees (a)");
3316 for (i=1; i<lg(Ld); i++)
3317 if (Ld[i]!=1 && (!uisprimepower(Ld[i],&ell2) || ell2!=ell))
3318 pari_err_IMPL("nfgrunwaldwang for non prime-power local degrees (b)");
3319 for (i=1; i<lg(pl); i++)
3320 if (pl[i]==-1 && ell%2)
3321 pari_err_IMPL("nfgrunwaldwang for non prime-power local degrees (c)");
3322
3323 w = bnf? bnf_get_tuN(bnf): itos(gel(nfrootsof1(nf),1));
3324
3325 /* TODO choice between kummer and generic ? Let user choose between speed
3326 * and size */
3327 if (w%n==0 && lg(Ld)>1)
3328 return gerepileupto(av,nfgwkummer(nf,Lpr,Ld,pl,var));
3329 if (ell==n) {
3330 if (!bnf) bnf = Buchall(nf,0,0);
3331 return gerepileupto(av,bnfgwgeneric(bnf,Lpr,Ld,pl,var));
3332 }
3333 pari_err_IMPL("nfgrunwaldwang for nonprime degree");
3334 return NULL; /*LCOV_EXCL_LINE*/
3335 }
3336
3337 /** HASSE INVARIANTS **/
3338
3339 /* TODO long -> ulong + uel */
3340 static GEN
hasseconvert(GEN H,long n)3341 hasseconvert(GEN H, long n)
3342 {
3343 GEN h, c;
3344 long i, l;
3345 switch(typ(H)) {
3346 case t_VEC:
3347 l = lg(H); h = cgetg(l,t_VECSMALL);
3348 if (l == 1) return h;
3349 c = gel(H,1);
3350 if (typ(c) == t_VEC && l == 3)
3351 return mkvec2(gel(H,1),hasseconvert(gel(H,2),n));
3352 for (i=1; i<l; i++)
3353 {
3354 c = gel(H,i);
3355 switch(typ(c)) {
3356 case t_INT: break;
3357 case t_INTMOD:
3358 c = gel(c,2); break;
3359 case t_FRAC :
3360 c = gmulgs(c,n);
3361 if (typ(c) == t_INT) break;
3362 pari_err_DOMAIN("hasseconvert [degree should be a denominator of the invariant]", "denom(h)", "ndiv", stoi(n), Q_denom(gel(H,i)));
3363 default : pari_err_TYPE("Hasse invariant", c);
3364 }
3365 h[i] = smodis(c,n);
3366 }
3367 return h;
3368 case t_VECSMALL: return H;
3369 }
3370 pari_err_TYPE("Hasse invariant", H);
3371 return NULL;/*LCOV_EXCL_LINE*/
3372 }
3373
3374 /* assume f >= 2 */
3375 static long
cyclicrelfrob0(GEN nf,GEN aut,GEN pr,GEN q,long f,long g)3376 cyclicrelfrob0(GEN nf, GEN aut, GEN pr, GEN q, long f, long g)
3377 {
3378 GEN T, p, a, b, modpr = nf_to_Fq_init(nf,&pr,&T,&p);
3379 long s;
3380
3381 a = pol_x(nf_get_varn(nf));
3382 b = galoisapply(nf, aut, modpr_genFq(modpr));
3383 b = nf_to_Fq(nf, b, modpr);
3384 for (s = 0; !ZX_equal(a, b); s++) a = Fq_pow(a, q, T, p);
3385 return g * Fl_inv(s, f); /* < n */
3386 }
3387
3388 static long
cyclicrelfrob(GEN rnf,GEN auts,GEN pr)3389 cyclicrelfrob(GEN rnf, GEN auts, GEN pr)
3390 {
3391 pari_sp av = avma;
3392 long f,g,frob, n = rnf_get_degree(rnf);
3393 GEN P = rnfidealprimedec(rnf, pr);
3394
3395 if (pr_get_e(gel(P,1)) > pr_get_e(pr))
3396 pari_err_DOMAIN("cyclicrelfrob","e(PR/pr)",">",gen_1,pr);
3397 g = lg(P) - 1;
3398 f = n / g;
3399
3400 if (f <= 2) frob = g % n;
3401 else {
3402 GEN nf2, PR = gel(P,1);
3403 GEN autabs = rnfeltreltoabs(rnf,gel(auts,g));
3404 nf2 = obj_check(rnf,rnf_NFABS);
3405 autabs = nfadd(nf2, autabs, gmul(rnf_get_k(rnf), rnf_get_alpha(rnf)));
3406 frob = cyclicrelfrob0(nf2, autabs, PR, pr_norm(pr), f, g);
3407 }
3408 return gc_long(av, frob);
3409 }
3410
3411 static long
localhasse(GEN rnf,GEN cnd,GEN pl,GEN auts,GEN b,long k)3412 localhasse(GEN rnf, GEN cnd, GEN pl, GEN auts, GEN b, long k)
3413 {
3414 pari_sp av = avma;
3415 long v, m, h, lfa, frob, n, i;
3416 GEN previous, y, pr, nf, q, fa;
3417 nf = rnf_get_nf(rnf);
3418 n = rnf_get_degree(rnf);
3419 pr = gcoeff(cnd,k,1);
3420 v = nfval(nf, b, pr);
3421 m = lg(cnd)>1 ? nbrows(cnd) : 0;
3422
3423 /* add the valuation of b to the conductor... */
3424 previous = gcoeff(cnd,k,2);
3425 gcoeff(cnd,k,2) = addis(previous, v);
3426
3427 y = const_vec(m, gen_1);
3428 gel(y,k) = b;
3429 /* find a factored element y congruent to b mod pr^(vpr(b)+vpr(cnd)) and to 1 mod the conductor. */
3430 y = factoredextchinese(nf, cnd, y, pl, &fa);
3431 h = 0;
3432 lfa = nbrows(fa);
3433 /* sum of all Hasse invariants of (rnf/nf,aut,y) is 0, Hasse invariants at q!=pr are easy, Hasse invariant at pr is the same as for al=(rnf/nf,aut,b). */
3434 for (i=1; i<=lfa; i++) {
3435 q = gcoeff(fa,i,1);
3436 if (cmp_prime_ideal(pr,q)) {
3437 frob = cyclicrelfrob(rnf, auts, q);
3438 frob = Fl_mul(frob,umodiu(gcoeff(fa,i,2),n),n);
3439 h = Fl_add(h,frob,n);
3440 }
3441 }
3442 /* ...then restore it. */
3443 gcoeff(cnd,k,2) = previous;
3444 return gc_long(av, Fl_neg(h,n));
3445 }
3446
3447 static GEN
allauts(GEN rnf,GEN aut)3448 allauts(GEN rnf, GEN aut)
3449 {
3450 long n = rnf_get_degree(rnf), i;
3451 GEN pol = rnf_get_pol(rnf), vaut;
3452 if (n==1) n=2;
3453 vaut = cgetg(n,t_VEC);
3454 aut = lift_shallow(rnfbasistoalg(rnf,aut));
3455 gel(vaut,1) = aut;
3456 for (i=1; i<n-1; i++)
3457 gel(vaut,i+1) = RgX_rem(poleval(gel(vaut,i), aut), pol);
3458 return vaut;
3459 }
3460
3461 static GEN
clean_factor(GEN fa)3462 clean_factor(GEN fa)
3463 {
3464 GEN P2,E2, P = gel(fa,1), E = gel(fa,2);
3465 long l = lg(P), i, j = 1;
3466 P2 = cgetg(l, t_COL);
3467 E2 = cgetg(l, t_COL);
3468 for (i = 1;i < l; i++)
3469 if (signe(gel(E,i))) {
3470 gel(P2,j) = gel(P,i);
3471 gel(E2,j) = gel(E,i); j++;
3472 }
3473 setlg(P2,j);
3474 setlg(E2,j); return mkmat2(P2,E2);
3475 }
3476
3477 /* shallow concat x[1],...x[nx],y[1], ... y[ny], returning a t_COL. To be
3478 * used when we do not know whether x,y are t_VEC or t_COL */
3479 static GEN
colconcat(GEN x,GEN y)3480 colconcat(GEN x, GEN y)
3481 {
3482 long i, lx = lg(x), ly = lg(y);
3483 GEN z=cgetg(lx+ly-1, t_COL);
3484 for (i=1; i<lx; i++) z[i] = x[i];
3485 for (i=1; i<ly; i++) z[lx+i-1]= y[i];
3486 return z;
3487 }
3488
3489 /* return v(x) at all primes in listpr, replace x by cofactor */
3490 static GEN
nfmakecoprime(GEN nf,GEN * px,GEN listpr)3491 nfmakecoprime(GEN nf, GEN *px, GEN listpr)
3492 {
3493 long j, l = lg(listpr);
3494 GEN x1, x = *px, L = cgetg(l, t_COL);
3495
3496 if (typ(x) != t_MAT)
3497 { /* scalar, divide at the end (fast valuation) */
3498 x1 = NULL;
3499 for (j=1; j<l; j++)
3500 {
3501 GEN pr = gel(listpr,j), e;
3502 long v = nfval(nf, x, pr);
3503 e = stoi(v); gel(L,j) = e;
3504 if (v) x1 = x1? idealmulpowprime(nf, x1, pr, e)
3505 : idealpow(nf, pr, e);
3506 }
3507 if (x1) x = idealdivexact(nf, idealhnf(nf,x), x1);
3508 }
3509 else
3510 { /* HNF, divide as we proceed (reduce size) */
3511 for (j=1; j<l; j++)
3512 {
3513 GEN pr = gel(listpr,j);
3514 long v = idealval(nf, x, pr);
3515 gel(L,j) = stoi(v);
3516 if (v) x = idealmulpowprime(nf, x, pr, stoi(-v));
3517 }
3518 }
3519 *px = x; return L;
3520 }
3521
3522 /* Caveat: factorizations are not sorted wrt cmp_prime_ideal: Lpr comes first */
3523 static GEN
computecnd(GEN rnf,GEN Lpr)3524 computecnd(GEN rnf, GEN Lpr)
3525 {
3526 GEN id, nf, fa, Le, P,E;
3527 long n = rnf_get_degree(rnf);
3528
3529 nf = rnf_get_nf(rnf);
3530 id = rnf_get_idealdisc(rnf);
3531 Le = nfmakecoprime(nf, &id, Lpr);
3532 fa = idealfactor(nf, id); /* part of D_{L/K} coprime with Lpr */
3533 P = colconcat(Lpr,gel(fa,1));
3534 E = colconcat(Le, gel(fa,2));
3535 fa = mkmat2(P, gdiventgs(E, eulerphiu(n)));
3536 return mkvec2(fa, clean_factor(fa));
3537 }
3538
3539 /* h >= 0 */
3540 static void
nextgen(GEN gene,long h,GEN * gens,GEN * hgens,long * ngens,long * curgcd)3541 nextgen(GEN gene, long h, GEN* gens, GEN* hgens, long* ngens, long* curgcd) {
3542 long nextgcd = ugcd(h,*curgcd);
3543 if (nextgcd == *curgcd) return;
3544 (*ngens)++;
3545 gel(*gens,*ngens) = gene;
3546 gel(*hgens,*ngens) = utoi(h);
3547 *curgcd = nextgcd;
3548 return;
3549 }
3550
3551 static int
dividesmod(long d,long h,long n)3552 dividesmod(long d, long h, long n) { return !(h%cgcd(d,n)); }
3553
3554 /* ramified prime with nontrivial Hasse invariant */
3555 static GEN
localcomplete(GEN rnf,GEN pl,GEN cnd,GEN auts,long j,long n,long h,long * v)3556 localcomplete(GEN rnf, GEN pl, GEN cnd, GEN auts, long j, long n, long h, long* v)
3557 {
3558 GEN nf, gens, hgens, pr, modpr, T, p, sol, U, D, b, gene, randg, pu;
3559 long ngens, i, d, np, k, d1, d2, hg, dnf, vcnd, curgcd;
3560 nf = rnf_get_nf(rnf);
3561 pr = gcoeff(cnd,j,1);
3562 np = umodiu(pr_norm(pr), n);
3563 dnf = nf_get_degree(nf);
3564 vcnd = itos(gcoeff(cnd,j,2));
3565 ngens = 13+dnf;
3566 gens = zerovec(ngens);
3567 hgens = zerovec(ngens);
3568 *v = 0;
3569 curgcd = 0;
3570 ngens = 0;
3571
3572 if (!uisprime(n)) {
3573 gene = pr_get_gen(pr);
3574 hg = localhasse(rnf, cnd, pl, auts, gene, j);
3575 nextgen(gene, hg, &gens, &hgens, &ngens, &curgcd);
3576 }
3577
3578 if (ugcd(np,n) != 1) { /* GCD(Np,n) != 1 */
3579 pu = idealprincipalunits(nf,pr,vcnd);
3580 pu = abgrp_get_gen(pu);
3581 for (i=1; i<lg(pu) && !dividesmod(curgcd,h,n); i++) {
3582 gene = gel(pu,i);
3583 hg = localhasse(rnf, cnd, pl, auts, gene, j);
3584 nextgen(gene, hg, &gens, &hgens, &ngens, &curgcd);
3585 }
3586 }
3587
3588 d = ugcd(np-1,n);
3589 if (d != 1) { /* GCD(Np-1,n) != 1 */
3590 modpr = nf_to_Fq_init(nf, &pr, &T, &p);
3591 while (!dividesmod(curgcd,h,n)) { /* TODO gener_FpXQ_local */
3592 if (T==NULL) randg = randomi(p);
3593 else randg = random_FpX(degpol(T), varn(T),p);
3594
3595 if (!gequal0(randg) && !gequal1(randg)) {
3596 gene = Fq_to_nf(randg, modpr);
3597 hg = localhasse(rnf, cnd, pl, auts, gene, j);
3598 nextgen(gene, hg, &gens, &hgens, &ngens, &curgcd);
3599 }
3600 }
3601 }
3602
3603 setlg(gens,ngens+1);
3604 setlg(hgens,ngens+1);
3605
3606 sol = ZV_extgcd(hgens);
3607 D = gel(sol,1);
3608 U = gmael(sol,2,ngens);
3609
3610 b = gen_1;
3611 d = itou(D);
3612 d1 = ugcd(d,n);
3613 d2 = d/d1;
3614 d = ((h/d1)*Fl_inv(d2,n))%n;
3615 for (i=1; i<=ngens; i++) {
3616 k = (itos(gel(U,i))*d)%n;
3617 if (k<0) k = n-k;
3618 if (k) b = nfmul(nf, b, nfpow_u(nf, gel(gens,i),k));
3619 if (i==1) *v = k;
3620 }
3621 return b;
3622 }
3623
3624 static int
testsplits(GEN data,GEN fa)3625 testsplits(GEN data, GEN fa)
3626 {
3627 GEN rnf = gel(data,1), forbid = gel(data,2), P = gel(fa,1), E = gel(fa,2);
3628 long i, n, l = lg(P);
3629
3630 for (i = 1; i < l; i++)
3631 {
3632 GEN pr = gel(P,i);
3633 if (tablesearch(forbid, pr, &cmp_prime_ideal)) return 0;
3634 }
3635 n = rnf_get_degree(rnf);
3636 for (i = 1; i < l; i++)
3637 {
3638 long e = itos(gel(E,i)) % n;
3639 if (e)
3640 {
3641 GEN L = rnfidealprimedec(rnf, gel(P,i));
3642 long g = lg(L) - 1;
3643 if ((e * g) % n) return 0;
3644 }
3645 }
3646 return 1;
3647 }
3648
3649 /* remove entries with Hasse invariant 0 */
3650 static GEN
hassereduce(GEN hf)3651 hassereduce(GEN hf)
3652 {
3653 GEN pr,h, PR = gel(hf,1), H = gel(hf,2);
3654 long i, j, l = lg(PR);
3655
3656 pr= cgetg(l, t_VEC);
3657 h = cgetg(l, t_VECSMALL);
3658 for (i = j = 1; i < l; i++)
3659 if (H[i]) {
3660 gel(pr,j) = gel(PR,i);
3661 h[j] = H[i]; j++;
3662 }
3663 setlg(pr,j);
3664 setlg(h,j); return mkvec2(pr,h);
3665 }
3666
3667 /* v vector of prid. Return underlying list of rational primes */
3668 static GEN
pr_primes(GEN v)3669 pr_primes(GEN v)
3670 {
3671 long i, l = lg(v);
3672 GEN w = cgetg(l,t_VEC);
3673 for (i=1; i<l; i++) gel(w,i) = pr_get_p(gel(v,i));
3674 return ZV_sort_uniq(w);
3675 }
3676
3677 /* rnf complete */
3678 static GEN
alg_complete0(GEN rnf,GEN aut,GEN hf,GEN hi,long maxord)3679 alg_complete0(GEN rnf, GEN aut, GEN hf, GEN hi, long maxord)
3680 {
3681 pari_sp av = avma;
3682 GEN nf, pl, pl2, cnd, prcnd, cnds, y, Lpr, auts, b, fa, data, hfe;
3683 GEN forbid, al, ind;
3684 long D, n, d, i, j, l;
3685 nf = rnf_get_nf(rnf);
3686 n = rnf_get_degree(rnf);
3687 d = nf_get_degree(nf);
3688 D = d*n*n;
3689 checkhasse(nf,hf,hi,n);
3690 hf = hassereduce(hf);
3691 Lpr = gel(hf,1);
3692 hfe = gel(hf,2);
3693
3694 auts = allauts(rnf,aut);
3695
3696 pl = leafcopy(hi); /* conditions on the final b */
3697 pl2 = leafcopy(hi); /* conditions for computing local Hasse invariants */
3698 l = lg(pl); ind = cgetg(l, t_VECSMALL);
3699 for (i = j = 1; i < l; i++)
3700 if (hi[i]) { pl[i] = -1; pl2[i] = 1; } else ind[j++] = i;
3701 setlg(ind, j);
3702 y = nfpolsturm(nf, rnf_get_pol(rnf), ind);
3703 for (i = 1; i < j; i++)
3704 if (!signe(gel(y,i))) { pl[ind[i]] = 1; pl2[ind[i]] = 1; }
3705
3706 cnds = computecnd(rnf,Lpr);
3707 prcnd = gel(cnds,1);
3708 cnd = gel(cnds,2);
3709 y = cgetg(lgcols(prcnd),t_VEC);
3710 forbid = vectrunc_init(lg(Lpr));
3711 for (i=j=1; i<lg(Lpr); i++)
3712 {
3713 GEN pr = gcoeff(prcnd,i,1), yi;
3714 long v, e = itou( gcoeff(prcnd,i,2) );
3715 if (!e) {
3716 long frob = cyclicrelfrob(rnf,auts,pr), f1 = ugcd(frob,n);
3717 vectrunc_append(forbid, pr);
3718 yi = gen_0;
3719 v = ((hfe[i]/f1) * Fl_inv(frob/f1,n)) % n;
3720 }
3721 else
3722 yi = localcomplete(rnf, pl2, cnd, auts, j++, n, hfe[i], &v);
3723 gel(y,i) = yi;
3724 gcoeff(prcnd,i,2) = stoi(e + v);
3725 }
3726 for (; i<lgcols(prcnd); i++) gel(y,i) = gen_1;
3727 gen_sort_inplace(forbid, (void*)&cmp_prime_ideal, &cmp_nodata, NULL);
3728 data = mkvec2(rnf,forbid);
3729 b = factoredextchinesetest(nf,prcnd,y,pl,&fa,data,testsplits);
3730
3731 al = cgetg(12, t_VEC);
3732 gel(al,10)= gen_0; /* must be set first */
3733 gel(al,1) = rnf;
3734 gel(al,2) = auts;
3735 gel(al,3) = basistoalg(nf,b);
3736 gel(al,4) = hi;
3737 /* add primes | disc or b with trivial Hasse invariant to hf */
3738 Lpr = gel(prcnd,1); y = b;
3739 (void)nfmakecoprime(nf, &y, Lpr);
3740 Lpr = shallowconcat(Lpr, gel(idealfactor(nf,y), 1));
3741 settyp(Lpr,t_VEC);
3742 hf = mkvec2(Lpr, shallowconcat(hfe, const_vecsmall(lg(Lpr)-lg(hfe), 0)));
3743 gel(al,5) = hf;
3744 gel(al,6) = gen_0;
3745 gel(al,7) = matid(D);
3746 gel(al,8) = matid(D); /* TODO modify 7, 8 et 9 once LLL added */
3747 gel(al,9) = algnatmultable(al,D);
3748 gel(al,11)= algtracebasis(al);
3749 if (maxord) al = alg_maximal_primes(al, pr_primes(Lpr));
3750 return gerepilecopy(av, al);
3751 }
3752
3753 GEN
alg_complete(GEN rnf,GEN aut,GEN hf,GEN hi,long maxord)3754 alg_complete(GEN rnf, GEN aut, GEN hf, GEN hi, long maxord)
3755 {
3756 long n = rnf_get_degree(rnf);
3757 rnfcomplete(rnf);
3758 return alg_complete0(rnf,aut,hasseconvert(hf,n),hasseconvert(hi,n), maxord);
3759 }
3760
3761 void
checkhasse(GEN nf,GEN hf,GEN hi,long n)3762 checkhasse(GEN nf, GEN hf, GEN hi, long n)
3763 {
3764 GEN Lpr, Lh;
3765 long i, sum;
3766 if (typ(hf) != t_VEC || lg(hf) != 3) pari_err_TYPE("checkhasse [hf]", hf);
3767 Lpr = gel(hf,1);
3768 Lh = gel(hf,2);
3769 if (typ(Lpr) != t_VEC) pari_err_TYPE("checkhasse [Lpr]", Lpr);
3770 if (typ(Lh) != t_VECSMALL) pari_err_TYPE("checkhasse [Lh]", Lh);
3771 if (typ(hi) != t_VECSMALL) pari_err_TYPE("checkhasse [hi]", hi);
3772 if ((nf && lg(hi) != nf_get_r1(nf)+1))
3773 pari_err_DOMAIN("checkhasse [hi should have r1 components]","#hi","!=",stoi(nf_get_r1(nf)),stoi(lg(hi)-1));
3774 if (lg(Lpr) != lg(Lh))
3775 pari_err_DIM("checkhasse [Lpr and Lh should have same length]");
3776 for (i=1; i<lg(Lpr); i++) checkprid(gel(Lpr,i));
3777 if (lg(gen_sort_uniq(Lpr, (void*)cmp_prime_ideal, cmp_nodata)) < lg(Lpr))
3778 pari_err(e_MISC, "error in checkhasse [duplicate prime ideal]");
3779 sum = 0;
3780 for (i=1; i<lg(Lh); i++) sum = (sum+Lh[i])%n;
3781 for (i=1; i<lg(hi); i++) {
3782 if (hi[i] && 2*hi[i] != n) pari_err_DOMAIN("checkhasse", "Hasse invariant at real place [must be 0 or 1/2]", "!=", n%2? gen_0 : stoi(n/2), stoi(hi[i]));
3783 sum = (sum+hi[i])%n;
3784 }
3785 if (sum<0) sum = n+sum;
3786 if (sum != 0)
3787 pari_err_DOMAIN("checkhasse","sum(Hasse invariants)","!=",gen_0,Lh);
3788 }
3789
3790 static GEN
hassecoprime(GEN hf,GEN hi,long n)3791 hassecoprime(GEN hf, GEN hi, long n)
3792 {
3793 pari_sp av = avma;
3794 long l, i, j, lk, inv;
3795 GEN fa, P,E, res, hil, hfl;
3796 hi = hasseconvert(hi, n);
3797 hf = hasseconvert(hf, n);
3798 checkhasse(NULL,hf,hi,n);
3799 fa = factoru(n);
3800 P = gel(fa,1); l = lg(P);
3801 E = gel(fa,2);
3802 res = cgetg(l,t_VEC);
3803 for (i=1; i<l; i++) {
3804 lk = upowuu(P[i],E[i]);
3805 inv = Fl_invsafe((n/lk)%lk, lk);
3806 hil = gcopy(hi);
3807 hfl = gcopy(hf);
3808
3809 if (P[i] == 2)
3810 for (j=1; j<lg(hil); j++) hil[j] = hi[j]==0 ? 0 : lk/2;
3811 else
3812 for (j=1; j<lg(hil); j++) hil[j] = 0;
3813 for (j=1; j<lgcols(hfl); j++) gel(hfl,2)[j] = (gel(hf,2)[j]*inv)%lk;
3814 hfl = hassereduce(hfl);
3815 gel(res,i) = mkvec3(hfl,hil,utoi(lk));
3816 }
3817
3818 return gerepilecopy(av, res);
3819 }
3820
3821 /* no garbage collection */
3822 static GEN
genefrob(GEN nf,GEN gal,GEN r)3823 genefrob(GEN nf, GEN gal, GEN r)
3824 {
3825 long i;
3826 GEN g = identity_perm(nf_get_degree(nf)), fa = Z_factor(r), p, pr, frob;
3827 for (i=1; i<lgcols(fa); i++) {
3828 p = gcoeff(fa,i,1);
3829 pr = idealprimedec(nf, p);
3830 pr = gel(pr,1);
3831 frob = idealfrobenius(nf, gal, pr);
3832 g = perm_mul(g, perm_pow(frob, gcoeff(fa,i,2)));
3833 }
3834 return g;
3835 }
3836
3837 static GEN
rnfcycaut(GEN rnf)3838 rnfcycaut(GEN rnf)
3839 {
3840 GEN nf2 = obj_check(rnf, rnf_NFABS);
3841 GEN L, alpha, pol, salpha, s, sj, polabs, k, X, pol0, nf;
3842 long i, d, j;
3843 d = rnf_get_degree(rnf);
3844 L = galoisconj(nf2,NULL);
3845 alpha = lift_shallow(rnf_get_alpha(rnf));
3846 pol = rnf_get_pol(rnf);
3847 k = rnf_get_k(rnf);
3848 polabs = rnf_get_polabs(rnf);
3849 nf = rnf_get_nf(rnf);
3850 pol0 = nf_get_pol(nf);
3851 X = RgX_rem(pol_x(varn(pol0)), pol0);
3852
3853 /* TODO check mod prime of degree 1 */
3854 for (i=1; i<lg(L); i++) {
3855 s = gel(L,i);
3856 salpha = RgX_RgXQ_eval(alpha,s,polabs);
3857 if (!gequal(alpha,salpha)) continue;
3858
3859 s = lift_shallow(rnfeltabstorel(rnf,s));
3860 sj = s = gsub(s, gmul(k,X));
3861 for (j=1; !gequal0(gsub(sj,pol_x(varn(s)))); j++)
3862 sj = RgX_RgXQ_eval(sj,s,pol);
3863 if (j<d) continue;
3864 return s;
3865 }
3866 return NULL; /*LCOV_EXCL_LINE*/
3867 }
3868
3869 /* returns Lpr augmented with an extra, distinct prime */
3870 /* TODO be less lazy and return a small prime */
3871 static GEN
extraprime(GEN nf,GEN Lpr)3872 extraprime(GEN nf, GEN Lpr)
3873 {
3874 GEN Lpr2, p = gen_2, pr;
3875 long i;
3876 Lpr2 = cgetg(lg(Lpr)+1,t_VEC);
3877 for (i=1; i<lg(Lpr); i++)
3878 {
3879 gel(Lpr2,i) = gel(Lpr,i);
3880 p = gmax_shallow(p, pr_get_p(gel(Lpr,i)));
3881 }
3882 p = nextprime(addis(p,1));
3883 pr = gel(idealprimedec_limit_f(nf, p, 0), 1);
3884 gel(Lpr2,lg(Lpr)) = pr;
3885 return Lpr2;
3886 }
3887
3888 GEN
alg_hasse(GEN nf,long n,GEN hf,GEN hi,long var,long maxord)3889 alg_hasse(GEN nf, long n, GEN hf, GEN hi, long var, long maxord)
3890 {
3891 pari_sp av = avma;
3892 GEN primary, al = gen_0, al2, rnf, hil, hfl, Ld, pl, pol, Lpr, aut, Lpr2, Ld2;
3893 long i, lk, j, maxdeg;
3894 dbg_printf(1)("alg_hasse\n");
3895 if (n<=1) pari_err_DOMAIN("alg_hasse", "degree", "<=", gen_1, stoi(n));
3896 primary = hassecoprime(hf, hi, n);
3897 for (i=1; i<lg(primary); i++) {
3898 lk = itos(gmael(primary,i,3));
3899 hfl = gmael(primary,i,1);
3900 hil = gmael(primary,i,2);
3901 checkhasse(nf, hfl, hil, lk);
3902 dbg_printf(1)("alg_hasse: i=%d hf=%Ps hi=%Ps lk=%d\n", i, hfl, hil, lk);
3903
3904 if (lg(gel(hfl,1))>1 || lk%2==0) {
3905 maxdeg = 1;
3906 Lpr = gel(hfl,1);
3907 Ld = gcopy(gel(hfl,2));
3908 for (j=1; j<lg(Ld); j++)
3909 {
3910 Ld[j] = lk/ugcd(lk,Ld[j]);
3911 maxdeg = maxss(Ld[j],maxdeg);
3912 }
3913 pl = gcopy(hil);
3914 for (j=1; j<lg(pl); j++) if(pl[j])
3915 {
3916 pl[j] = -1;
3917 maxdeg = maxss(maxdeg,2);
3918 }
3919
3920 Lpr2 = Lpr;
3921 Ld2 = Ld;
3922 if (maxdeg<lk)
3923 {
3924 if (maxdeg==1 && lk==2 && lg(pl)>1) pl[1] = -1;
3925 else
3926 {
3927 Lpr2 = extraprime(nf,Lpr);
3928 Ld2 = cgetg(lg(Ld)+1, t_VECSMALL);
3929 for (j=1; j<lg(Ld); j++) Ld2[j] = Ld[j];
3930 Ld2[lg(Ld)] = lk;
3931 }
3932 }
3933
3934 dbg_printf(2)("alg_hasse: calling nfgrunwaldwang Lpr=%Ps Pd=%Ps pl=%Ps\n",
3935 Lpr, Ld, pl);
3936 pol = nfgrunwaldwang(nf, Lpr2, Ld2, pl, var);
3937 dbg_printf(2)("alg_hasse: calling rnfinit(%Ps)\n", pol);
3938 rnf = rnfinit0(nf,pol,1);
3939 dbg_printf(2)("alg_hasse: computing automorphism\n");
3940 aut = rnfcycaut(rnf);
3941 dbg_printf(2)("alg_hasse: calling alg_complete\n");
3942 al2 = alg_complete0(rnf,aut,hfl,hil,maxord);
3943 }
3944 else al2 = alg_matrix(nf, lk, var, cgetg(1,t_VEC), maxord);
3945
3946 if (i==1) al = al2;
3947 else al = algtensor(al,al2,maxord);
3948 }
3949 return gerepilecopy(av,al);
3950 }
3951
3952 /** CYCLIC ALGEBRA WITH GIVEN HASSE INVARIANTS **/
3953
3954 /* no garbage collection */
3955 static int
linindep(GEN pol,GEN L)3956 linindep(GEN pol, GEN L)
3957 {
3958 long i;
3959 GEN fa;
3960 for (i=1; i<lg(L); i++) {
3961 fa = nffactor(gel(L,i),pol);
3962 if (lgcols(fa)>2) return 0;
3963 }
3964 return 1;
3965 }
3966
3967 /* no garbage collection */
3968 static GEN
subcycloindep(GEN nf,long n,long v,GEN L,GEN * pr)3969 subcycloindep(GEN nf, long n, long v, GEN L, GEN *pr)
3970 {
3971 pari_sp av;
3972 forprime_t S;
3973 ulong p;
3974 u_forprime_arith_init(&S, 1, ULONG_MAX, 1, n);
3975 av = avma;
3976 while ((p = u_forprime_next(&S)))
3977 {
3978 ulong r = pgener_Fl(p);
3979 GEN pol = galoissubcyclo(utoipos(p), utoipos(Fl_powu(r,n,p)), 0, v);
3980 GEN fa = nffactor(nf, pol);
3981 if (lgcols(fa) == 2 && linindep(pol,L)) { *pr = utoipos(r); return pol; }
3982 set_avma(av);
3983 }
3984 pari_err_BUG("subcycloindep (no suitable prime = 1(mod n))"); /*LCOV_EXCL_LINE*/
3985 *pr = NULL; return NULL; /*LCOV_EXCL_LINE*/
3986 }
3987
3988 GEN
alg_matrix(GEN nf,long n,long v,GEN L,long maxord)3989 alg_matrix(GEN nf, long n, long v, GEN L, long maxord)
3990 {
3991 pari_sp av = avma;
3992 GEN pol, gal, rnf, cyclo, g, r, aut;
3993 dbg_printf(1)("alg_matrix\n");
3994 if (n<=0) pari_err_DOMAIN("alg_matrix", "n", "<=", gen_0, stoi(n));
3995 pol = subcycloindep(nf, n, v, L, &r);
3996 rnf = rnfinit(nf, pol);
3997 cyclo = nfinit(pol, nf_get_prec(nf));
3998 gal = galoisinit(cyclo, NULL);
3999 g = genefrob(cyclo,gal,r);
4000 aut = galoispermtopol(gal,g);
4001 return gerepileupto(av, alg_cyclic(rnf, aut, gen_1, maxord));
4002 }
4003
4004 GEN
alg_hilbert(GEN nf,GEN a,GEN b,long v,long maxord)4005 alg_hilbert(GEN nf, GEN a, GEN b, long v, long maxord)
4006 {
4007 pari_sp av = avma;
4008 GEN C, P, rnf, aut;
4009 dbg_printf(1)("alg_hilbert\n");
4010 checknf(nf);
4011 if (!isint1(Q_denom(a)))
4012 pari_err_DOMAIN("alg_hilbert", "denominator(a)", "!=", gen_1,a);
4013 if (!isint1(Q_denom(b)))
4014 pari_err_DOMAIN("alg_hilbert", "denominator(b)", "!=", gen_1,b);
4015
4016 if (v < 0) v = 0;
4017 C = Rg_col_ei(gneg(a), 3, 3);
4018 gel(C,1) = gen_1;
4019 P = gtopoly(C,v);
4020 rnf = rnfinit(nf, P);
4021 aut = gneg(pol_x(v));
4022 return gerepileupto(av, alg_cyclic(rnf, aut, b, maxord));
4023 }
4024
4025 GEN
alginit(GEN A,GEN B,long v,long maxord)4026 alginit(GEN A, GEN B, long v, long maxord)
4027 {
4028 long w;
4029 switch(nftyp(A))
4030 {
4031 case typ_NF:
4032 if (v<0) v=0;
4033 w = gvar(nf_get_pol(A));
4034 if (varncmp(v,w)>=0) pari_err_PRIORITY("alginit", pol_x(v), ">=", w);
4035 switch(typ(B))
4036 {
4037 long nB;
4038 case t_INT: return alg_matrix(A, itos(B), v, cgetg(1,t_VEC), maxord);
4039 case t_VEC:
4040 nB = lg(B)-1;
4041 if (nB && typ(gel(B,1)) == t_MAT) return alg_csa_table(A,B,v,maxord);
4042 switch(nB)
4043 {
4044 case 2: return alg_hilbert(A, gel(B,1), gel(B,2), v, maxord);
4045 case 3:
4046 if (typ(gel(B,1))!=t_INT)
4047 pari_err_TYPE("alginit [degree should be an integer]", gel(B,1));
4048 return alg_hasse(A, itos(gel(B,1)), gel(B,2), gel(B,3), v,
4049 maxord);
4050 }
4051 }
4052 pari_err_TYPE("alginit", B); break;
4053
4054 case typ_RNF:
4055 if (typ(B) != t_VEC || lg(B) != 3) pari_err_TYPE("alginit", B);
4056 return alg_cyclic(A, gel(B,1), gel(B,2), maxord);
4057 }
4058 pari_err_TYPE("alginit", A);
4059 return NULL;/*LCOV_EXCL_LINE*/
4060 }
4061
4062 /* assumes al CSA or CYCLIC */
4063 static GEN
algnatmultable(GEN al,long D)4064 algnatmultable(GEN al, long D)
4065 {
4066 GEN res, x;
4067 long i;
4068 res = cgetg(D+1,t_VEC);
4069 for (i=1; i<=D; i++) {
4070 x = algnattoalg(al,col_ei(D,i));
4071 gel(res,i) = algZmultable(al,x);
4072 }
4073 return res;
4074 }
4075
4076 /* no garbage collection */
4077 static void
algcomputehasse(GEN al)4078 algcomputehasse(GEN al)
4079 {
4080 long r1, k, n, m, m1, m2, m3, i, m23, m123;
4081 GEN rnf, nf, b, fab, disc2, cnd, fad, auts, pr, pl, perm, y, hi, PH, H, L;
4082
4083 rnf = alg_get_splittingfield(al);
4084 n = rnf_get_degree(rnf);
4085 nf = rnf_get_nf(rnf);
4086 b = alg_get_b(al);
4087 r1 = nf_get_r1(nf);
4088 auts = alg_get_auts(al);
4089 (void)alg_get_abssplitting(al);
4090
4091 y = nfpolsturm(nf, rnf_get_pol(rnf), NULL);
4092 pl = cgetg(r1+1, t_VECSMALL);
4093 /* real places where rnf/nf ramifies */
4094 for (k = 1; k <= r1; k++) pl[k] = !signe(gel(y,k));
4095
4096 /* infinite Hasse invariants */
4097 if (odd(n)) hi = const_vecsmall(r1, 0);
4098 else
4099 {
4100 GEN s = nfsign(nf, b);
4101 hi = cgetg(r1+1, t_VECSMALL);
4102 for (k = 1; k<=r1; k++) hi[k] = (s[k] && pl[k]) ? (n/2) : 0;
4103 }
4104
4105 fab = idealfactor(nf, b);
4106 disc2 = rnf_get_idealdisc(rnf);
4107 L = nfmakecoprime(nf, &disc2, gel(fab,1));
4108 m = lg(L)-1;
4109 /* m1 = #{pr|b: pr \nmid disc}, m3 = #{pr|b: pr | disc} */
4110 perm = cgetg(m+1, t_VECSMALL);
4111 for (i=1, m1=m, k=1; k<=m; k++)
4112 if (signe(gel(L,k))) perm[m1--] = k; else perm[i++] = k;
4113 m3 = m - m1;
4114
4115 /* disc2 : factor of disc coprime to b */
4116 fad = idealfactor(nf, disc2);
4117 /* m2 : number of prime factors of disc not dividing b */
4118 m2 = nbrows(fad);
4119 m23 = m2+m3;
4120 m123 = m1+m2+m3;
4121
4122 /* initialize the possibly ramified primes (hasse) and the factored conductor of rnf/nf (cnd) */
4123 cnd = zeromatcopy(m23,2);
4124 PH = cgetg(m123+1, t_VEC); /* ramified primes */
4125 H = cgetg(m123+1, t_VECSMALL); /* Hasse invariant */
4126 /* compute Hasse invariant at primes that are unramified in rnf/nf */
4127 for (k=1; k<=m1; k++) {/* pr | b, pr \nmid disc */
4128 long frob, e, j = perm[k];
4129 pr = gcoeff(fab,j,1);
4130 e = itos(gcoeff(fab,j,2));
4131 frob = cyclicrelfrob(rnf, auts, pr);
4132 gel(PH,k) = pr;
4133 H[k] = Fl_mul(frob, e, n);
4134 }
4135 /* compute Hasse invariant at primes that are ramified in rnf/nf */
4136 for (k=1; k<=m2; k++) {/* pr \nmid b, pr | disc */
4137 pr = gcoeff(fad,k,1);
4138 gel(PH,k+m1) = pr;
4139 gcoeff(cnd,k,1) = pr;
4140 gcoeff(cnd,k,2) = gcoeff(fad,k,2);
4141 }
4142 for (k=1; k<=m3; k++) { /* pr | (b, disc) */
4143 long j = perm[k+m1];
4144 pr = gcoeff(fab,j,1);
4145 gel(PH,k+m1+m2) = pr;
4146 gcoeff(cnd,k+m2,1) = pr;
4147 gcoeff(cnd,k+m2,2) = gel(L,j);
4148 }
4149 gel(cnd,2) = gdiventgs(gel(cnd,2), eulerphiu(n));
4150 for (k=1; k<=m23; k++) H[k+m1] = localhasse(rnf, cnd, pl, auts, b, k);
4151 gel(al,4) = hi;
4152 perm = gen_indexsort(PH, (void*)&cmp_prime_ideal, &cmp_nodata);
4153 gel(al,5) = mkvec2(vecpermute(PH,perm),vecsmallpermute(H,perm));
4154 checkhasse(nf,alg_get_hasse_f(al),alg_get_hasse_i(al),n);
4155 }
4156
4157 static GEN
alg_maximal_primes(GEN al,GEN P)4158 alg_maximal_primes(GEN al, GEN P)
4159 {
4160 pari_sp av = avma;
4161 long l = lg(P), i;
4162 for (i=1; i<l; i++)
4163 {
4164 if (i != 1) al = gerepilecopy(av, al);
4165 al = alg_pmaximal(al,gel(P,i));
4166 }
4167 return al;
4168 }
4169
4170 GEN
alg_cyclic(GEN rnf,GEN aut,GEN b,long maxord)4171 alg_cyclic(GEN rnf, GEN aut, GEN b, long maxord)
4172 {
4173 pari_sp av = avma;
4174 GEN al, nf;
4175 long D, n, d;
4176 dbg_printf(1)("alg_cyclic\n");
4177 checkrnf(rnf);
4178 if (!isint1(Q_denom(b)))
4179 pari_err_DOMAIN("alg_cyclic", "denominator(b)", "!=", gen_1,b);
4180
4181 nf = rnf_get_nf(rnf);
4182 n = rnf_get_degree(rnf);
4183 d = nf_get_degree(nf);
4184 D = d*n*n;
4185
4186 al = cgetg(12,t_VEC);
4187 gel(al,10)= gen_0; /* must be set first */
4188 gel(al,1) = rnf;
4189 gel(al,2) = allauts(rnf, aut);
4190 gel(al,3) = basistoalg(nf,b);
4191 rnf_build_nfabs(rnf, nf_get_prec(nf));
4192 gel(al,6) = gen_0;
4193 gel(al,7) = matid(D);
4194 gel(al,8) = matid(D); /* TODO modify 7, 8 et 9 once LLL added */
4195 gel(al,9) = algnatmultable(al,D);
4196 gel(al,11)= algtracebasis(al);
4197
4198 algcomputehasse(al);
4199
4200 if (maxord) {
4201 GEN hf = alg_get_hasse_f(al), pr = gel(hf,1);
4202 al = alg_maximal_primes(al, pr_primes(pr));
4203 }
4204 return gerepilecopy(av, al);
4205 }
4206
4207 static int
ismaximalsubfield(GEN al,GEN x,GEN d,long v,GEN * pt_minpol)4208 ismaximalsubfield(GEN al, GEN x, GEN d, long v, GEN *pt_minpol)
4209 {
4210 GEN cp = algbasischarpoly(al, x, v), lead;
4211 if (!ispower(cp, d, pt_minpol)) return 0;
4212 lead = leading_coeff(*pt_minpol);
4213 if (isintm1(lead)) *pt_minpol = gneg(*pt_minpol);
4214 return ZX_is_irred(*pt_minpol);
4215 }
4216
4217 static GEN
findmaximalsubfield(GEN al,GEN d,long v)4218 findmaximalsubfield(GEN al, GEN d, long v)
4219 {
4220 long count, nb=2, i, N = alg_get_absdim(al), n = nf_get_degree(alg_get_center(al));
4221 GEN x, minpol, maxc = gen_1;
4222
4223 for (i=n+1; i<=N; i+=n) {
4224 for (count=0; count<2 && i+count<=N; count++) {
4225 x = col_ei(N,i+count);
4226 if (ismaximalsubfield(al, x, d, v, &minpol)) return mkvec2(x,minpol);
4227 }
4228 }
4229
4230 while(1) {
4231 x = zerocol(N);
4232 for (count=0; count<nb; count++)
4233 {
4234 i = random_Fl(N)+1;
4235 gel(x,i) = addiu(randomi(maxc),1);
4236 if (random_bits(1)) gel(x,i) = negi(gel(x,i));
4237 }
4238 if (ismaximalsubfield(al, x, d, v, &minpol)) return mkvec2(x,minpol);
4239 if (!random_bits(3)) maxc = addiu(maxc,1);
4240 if (nb<N) nb++;
4241 }
4242
4243 return NULL; /* LCOV_EXCL_LINE */
4244 }
4245
4246 static GEN
frobeniusform(GEN al,GEN x)4247 frobeniusform(GEN al, GEN x)
4248 {
4249 GEN M, FP, P, Pi;
4250
4251 /* /!\ has to be the *right* multiplication table */
4252 M = algbasisrightmultable(al, x);
4253
4254 FP = matfrobenius(M,2,0); /* M = P^(-1)*F*P */
4255 P = gel(FP,2);
4256 Pi = RgM_inv(P);
4257 return mkvec2(P, Pi);
4258 }
4259
4260 static void
computesplitting(GEN al,long d,long v)4261 computesplitting(GEN al, long d, long v)
4262 {
4263 GEN subf, x, pol, polabs, basis, P, Pi, nf = alg_get_center(al), rnf, Lbasis, Lbasisinv, Q, pows;
4264 long i, n = nf_get_degree(nf), nd = n*d, N = alg_get_absdim(al), j, j2;
4265
4266 subf = findmaximalsubfield(al, utoipos(d), v);
4267 x = gel(subf, 1);
4268 polabs = gel(subf, 2);
4269
4270 /* Frobenius form to obtain L-vector space structure */
4271 basis = frobeniusform(al, x);
4272 P = gel(basis, 1);
4273 Pi = gel(basis, 2);
4274
4275 /* construct rnf of splitting field */
4276 pol = nffactor(nf,polabs);
4277 pol = gcoeff(pol,1,1);
4278 gel(al,1) = rnf = rnfinit(nf, pol);
4279 /* since pol is irreducible over Q, we have k=0 in rnf. */
4280 if (!gequal0(rnf_get_k(rnf)))
4281 pari_err_BUG("computesplitting (k!=0)"); /*LCOV_EXCL_LINE*/
4282 gel(al,6) = gen_0;
4283 rnf_build_nfabs(rnf, nf_get_prec(nf));
4284
4285 /* construct splitting data */
4286 Lbasis = cgetg(d+1, t_MAT);
4287 for (j=j2=1; j<=d; j++, j2+=nd)
4288 gel(Lbasis,j) = gel(Pi,j2);
4289
4290 Q = zeromatcopy(d,N);
4291 pows = pol_x_powers(nd,v);
4292 for (i=j=1; j<=N; j+=nd, i++)
4293 for (j2=0; j2<nd; j2++)
4294 gcoeff(Q,i,j+j2) = mkpolmod(gel(pows,j2+1),polabs);
4295 Lbasisinv = RgM_mul(Q,P);
4296
4297 gel(al,3) = mkvec3(x,Lbasis,Lbasisinv);
4298 }
4299
4300 /* assumes that mt defines a central simple algebra over nf */
4301 GEN
alg_csa_table(GEN nf,GEN mt0,long v,long maxord)4302 alg_csa_table(GEN nf, GEN mt0, long v, long maxord)
4303 {
4304 pari_sp av = avma;
4305 GEN al, mt;
4306 long n, D, d2 = lg(mt0)-1, d = usqrt(d2);
4307 dbg_printf(1)("alg_csa_table\n");
4308
4309 nf = checknf(nf);
4310 mt = check_relmt(nf,mt0);
4311 if (!mt) pari_err_TYPE("alg_csa_table", mt0);
4312 n = nf_get_degree(nf);
4313 D = n*d2;
4314 if (d*d != d2)
4315 pari_err_DOMAIN("alg_csa_table","(nonsquare) dimension","!=",stoi(d*d),mt);
4316
4317 al = cgetg(12, t_VEC);
4318 gel(al,10) = gen_0; /* must be set first */
4319 gel(al,1) = zerovec(12); gmael(al,1,10) = nf;
4320 gmael(al,1,1) = gpowgs(pol_x(0), d); /* placeholder before splitting field */
4321 gel(al,2) = mt;
4322 gel(al,3) = gen_0; /* placeholder */
4323 gel(al,4) = gel(al,5) = gen_0; /* TODO Hasse invariants */
4324 gel(al,5) = gel(al,6) = gen_0; /* placeholder */
4325 gel(al,7) = matid(D);
4326 gel(al,8) = matid(D);
4327 gel(al,9) = algnatmultable(al,D);
4328 gel(al,11)= algtracebasis(al);
4329 if (maxord) al = alg_maximal(al);
4330 computesplitting(al, d, v);
4331 return gerepilecopy(av, al);
4332 }
4333
4334 static GEN
algtableinit_i(GEN mt0,GEN p)4335 algtableinit_i(GEN mt0, GEN p)
4336 {
4337 GEN al, mt;
4338 long i, n;
4339
4340 if (p && !signe(p)) p = NULL;
4341 mt = check_mt(mt0,p);
4342 if (!mt) pari_err_TYPE("algtableinit", mt0);
4343 if (!p && !isint1(Q_denom(mt0)))
4344 pari_err_DOMAIN("algtableinit", "denominator(mt)", "!=", gen_1, mt0);
4345 n = lg(mt)-1;
4346 al = cgetg(12, t_VEC);
4347 for (i=1; i<=6; i++) gel(al,i) = gen_0;
4348 gel(al,7) = matid(n);
4349 gel(al,8) = matid(n);
4350 gel(al,9) = mt;
4351 gel(al,10) = p? p: gen_0;
4352 gel(al,11)= algtracebasis(al);
4353 return al;
4354 }
4355 GEN
algtableinit(GEN mt0,GEN p)4356 algtableinit(GEN mt0, GEN p)
4357 {
4358 pari_sp av = avma;
4359 if (p)
4360 {
4361 if (typ(p) != t_INT) pari_err_TYPE("algtableinit",p);
4362 if (signe(p) && !BPSW_psp(p)) pari_err_PRIME("algtableinit",p);
4363 }
4364 return gerepilecopy(av, algtableinit_i(mt0, p));
4365 }
4366
4367 /** REPRESENTATIONS OF GROUPS **/
4368
4369 static GEN
list_to_regular_rep(GEN elts,long n)4370 list_to_regular_rep(GEN elts, long n)
4371 {
4372 GEN reg, elts2, g;
4373 long i,j;
4374 elts = shallowcopy(elts);
4375 gen_sort_inplace(elts, (void*)&vecsmall_lexcmp, &cmp_nodata, NULL);
4376 reg = cgetg(n+1, t_VEC);
4377 gel(reg,1) = identity_perm(n);
4378 for (i=2; i<=n; i++) {
4379 g = perm_inv(gel(elts,i));
4380 elts2 = cgetg(n+1, t_VEC);
4381 for (j=1; j<=n; j++) gel(elts2,j) = perm_mul(g,gel(elts,j));
4382 gen_sort_inplace(elts2, (void*)&vecsmall_lexcmp, &cmp_nodata, &gel(reg,i));
4383 }
4384 return reg;
4385 }
4386
4387 static GEN
matrix_perm(GEN perm,long n)4388 matrix_perm(GEN perm, long n)
4389 {
4390 GEN m;
4391 long j;
4392 m = cgetg(n+1, t_MAT);
4393 for (j=1; j<=n; j++) {
4394 gel(m,j) = col_ei(n,perm[j]);
4395 }
4396 return m;
4397 }
4398
4399 GEN
conjclasses_algcenter(GEN cc,GEN p)4400 conjclasses_algcenter(GEN cc, GEN p)
4401 {
4402 GEN mt, elts = gel(cc,1), conjclass = gel(cc,2), rep = gel(cc,3), card;
4403 long i, nbcl = lg(rep)-1, n = lg(elts)-1;
4404 pari_sp av;
4405
4406 card = zero_Flv(nbcl);
4407 for (i=1; i<=n; i++) card[conjclass[i]]++;
4408
4409 /* multiplication table of the center of Z[G] (class functions) */
4410 mt = cgetg(nbcl+1,t_VEC);
4411 for (i=1;i<=nbcl;i++) gel(mt,i) = zero_Flm_copy(nbcl,nbcl);
4412 av = avma;
4413 for (i=1;i<=nbcl;i++)
4414 {
4415 GEN xi = gel(elts,rep[i]), mi = gel(mt,i);
4416 long j,k;
4417 for (j=1;j<=n;j++)
4418 {
4419 GEN xj = gel(elts,j);
4420 k = vecsearch(elts, perm_mul(xi,xj), NULL);
4421 ucoeff(mi, conjclass[k], conjclass[j])++;
4422 }
4423 for (k=1; k<=nbcl; k++)
4424 for (j=1; j<=nbcl; j++)
4425 {
4426 ucoeff(mi,k,j) *= card[i];
4427 ucoeff(mi,k,j) /= card[k];
4428 }
4429 set_avma(av);
4430 }
4431 for (i=1;i<=nbcl;i++) gel(mt,i) = Flm_to_ZM(gel(mt,i));
4432 return algtableinit_i(mt,p);
4433 }
4434
4435 GEN
alggroupcenter(GEN G,GEN p,GEN * pcc)4436 alggroupcenter(GEN G, GEN p, GEN *pcc)
4437 {
4438 pari_sp av = avma;
4439 GEN cc = group_to_cc(G), al = conjclasses_algcenter(cc, p);
4440 if (!pcc) al = gerepilecopy(av,al);
4441 else
4442 { *pcc = cc; gerepileall(av,2,&al,pcc); }
4443 return al;
4444 }
4445
4446 static GEN
groupelts_algebra(GEN elts,GEN p)4447 groupelts_algebra(GEN elts, GEN p)
4448 {
4449 pari_sp av = avma;
4450 GEN mt;
4451 long i, n = lg(elts)-1;
4452 elts = list_to_regular_rep(elts,n);
4453 mt = cgetg(n+1, t_VEC);
4454 for (i=1; i<=n; i++) gel(mt,i) = matrix_perm(gel(elts,i),n);
4455 return gerepilecopy(av, algtableinit_i(mt,p));
4456 }
4457
4458 GEN
alggroup(GEN gal,GEN p)4459 alggroup(GEN gal, GEN p)
4460 {
4461 GEN elts = checkgroupelts(gal);
4462 return groupelts_algebra(elts, p);
4463 }
4464
4465 /** MAXIMAL ORDER **/
4466
4467 GEN
alg_changeorder(GEN al,GEN ord)4468 alg_changeorder(GEN al, GEN ord)
4469 {
4470 GEN al2, mt, iord, mtx;
4471 long i, n;
4472 pari_sp av = avma;
4473
4474 if (!gequal0(gel(al,10)))
4475 pari_err_DOMAIN("alg_changeorder","characteristic","!=",gen_0,gel(al,10));
4476 n = alg_get_absdim(al);
4477
4478 iord = QM_inv(ord);
4479 al2 = shallowcopy(al);
4480
4481 gel(al2,7) = RgM_mul(gel(al2,7), ord);
4482
4483 gel(al2,8) = RgM_mul(iord, gel(al,8));
4484
4485 mt = cgetg(n+1,t_VEC);
4486 gel(mt,1) = matid(n);
4487 for (i=2; i<=n; i++) {
4488 mtx = algbasismultable(al,gel(ord,i));
4489 gel(mt,i) = RgM_mul(iord, RgM_mul(mtx, ord));
4490 }
4491 gel(al2,9) = mt;
4492
4493 gel(al2,11) = algtracebasis(al2);
4494
4495 return gerepilecopy(av,al2);
4496 }
4497
4498 static GEN
mattocol(GEN M,long n)4499 mattocol(GEN M, long n)
4500 {
4501 GEN C = cgetg(n*n+1, t_COL);
4502 long i,j,ic;
4503 ic = 1;
4504 for (i=1; i<=n; i++)
4505 for (j=1; j<=n; j++, ic++) gel(C,ic) = gcoeff(M,i,j);
4506 return C;
4507 }
4508
4509 /* Ip is a lift of a left O/pO-ideal where O is the integral basis of al */
4510 static GEN
algleftordermodp(GEN al,GEN Ip,GEN p)4511 algleftordermodp(GEN al, GEN Ip, GEN p)
4512 {
4513 pari_sp av = avma;
4514 GEN I, Ii, M, mt, K, imi, p2;
4515 long n, i;
4516 n = alg_get_absdim(al);
4517 mt = alg_get_multable(al);
4518 p2 = sqri(p);
4519
4520 I = ZM_hnfmodid(Ip, p);
4521 Ii = ZM_inv(I,NULL);
4522
4523 M = cgetg(n+1, t_MAT);
4524 for (i=1; i<=n; i++) {
4525 imi = FpM_mul(Ii, FpM_mul(gel(mt,i), I, p2), p2);
4526 imi = ZM_Z_divexact(imi, p);
4527 gel(M,i) = mattocol(imi, n);
4528 }
4529 K = FpM_ker(M, p);
4530 if (lg(K)==1) { set_avma(av); return matid(n); }
4531 K = ZM_hnfmodid(K,p);
4532
4533 return gerepileupto(av, ZM_Z_div(K,p));
4534 }
4535
4536 static GEN
alg_ordermodp(GEN al,GEN p)4537 alg_ordermodp(GEN al, GEN p)
4538 {
4539 GEN alp;
4540 long i, N = alg_get_absdim(al);
4541 alp = cgetg(12, t_VEC);
4542 for (i=1; i<=8; i++) gel(alp,i) = gen_0;
4543 gel(alp,9) = cgetg(N+1, t_VEC);
4544 for (i=1; i<=N; i++) gmael(alp,9,i) = FpM_red(gmael(al,9,i), p);
4545 gel(alp,10) = p;
4546 gel(alp,11) = cgetg(N+1, t_VEC);
4547 for (i=1; i<=N; i++) gmael(alp,11,i) = Fp_red(gmael(al,11,i), p);
4548
4549 return alp;
4550 }
4551
4552 static GEN
algpradical_i(GEN al,GEN p,GEN zprad,GEN projs)4553 algpradical_i(GEN al, GEN p, GEN zprad, GEN projs)
4554 {
4555 pari_sp av = avma;
4556 GEN alp = alg_ordermodp(al, p), liftrad, projrad, alq, alrad, res, Lalp, radq;
4557 long i;
4558 if (lg(zprad)==1) {
4559 liftrad = NULL;
4560 projrad = NULL;
4561 }
4562 else {
4563 alq = alg_quotient(alp, zprad, 1);
4564 alp = gel(alq,1);
4565 projrad = gel(alq,2);
4566 liftrad = gel(alq,3);
4567 }
4568
4569 if (projs) {
4570 if (projrad) {
4571 projs = gcopy(projs);
4572 for (i=1; i<lg(projs); i++)
4573 gel(projs,i) = FpM_FpC_mul(projrad, gel(projs,i), p);
4574 }
4575 Lalp = alg_centralproj(alp, projs, 1);
4576
4577 alrad = cgetg(lg(Lalp),t_VEC);
4578 for (i=1; i<lg(Lalp); i++) {
4579 alq = gel(Lalp,i);
4580 radq = algradical(gel(alq,1));
4581 if (gequal0(radq))
4582 gel(alrad,i) = cgetg(1,t_MAT);
4583 else {
4584 radq = FpM_mul(gel(alq,3),radq,p);
4585 gel(alrad,i) = radq;
4586 }
4587 }
4588 alrad = shallowmatconcat(alrad);
4589 alrad = FpM_image(alrad,p);
4590 }
4591 else alrad = algradical(alp);
4592
4593 if (!gequal0(alrad)) {
4594 if (liftrad) alrad = FpM_mul(liftrad, alrad, p);
4595 res = shallowmatconcat(mkvec2(alrad, zprad));
4596 res = FpM_image(res,p);
4597 }
4598 else res = lg(zprad)==1 ? gen_0 : zprad;
4599 return gerepilecopy(av, res);
4600 }
4601
4602 static GEN
algpdecompose0(GEN al,GEN prad,GEN p,GEN projs)4603 algpdecompose0(GEN al, GEN prad, GEN p, GEN projs)
4604 {
4605 pari_sp av = avma;
4606 GEN alp, quo, ss, liftm = NULL, projm = NULL, dec, res, I, Lss, deci;
4607 long i, j;
4608
4609 alp = alg_ordermodp(al, p);
4610 if (!gequal0(prad)) {
4611 quo = alg_quotient(alp, prad, 1);
4612 ss = gel(quo,1);
4613 projm = gel(quo,2);
4614 liftm = gel(quo,3);
4615 }
4616 else ss = alp;
4617
4618 if (projs) {
4619 if (projm) {
4620 for (i=1; i<lg(projs); i++)
4621 gel(projs,i) = FpM_FpC_mul(projm, gel(projs,i), p);
4622 }
4623 Lss = alg_centralproj(ss, projs, 1);
4624
4625 dec = cgetg(lg(Lss),t_VEC);
4626 for (i=1; i<lg(Lss); i++) {
4627 gel(dec,i) = algsimpledec_ss(gmael(Lss,i,1), 1);
4628 deci = gel(dec,i);
4629 for (j=1; j<lg(deci); j++)
4630 gmael(deci,j,3) = FpM_mul(gmael(Lss,i,3), gmael(deci,j,3), p);
4631 }
4632 dec = shallowconcat1(dec);
4633 }
4634 else dec = algsimpledec_ss(ss,1);
4635
4636 res = cgetg(lg(dec),t_VEC);
4637 for (i=1; i<lg(dec); i++) {
4638 I = gmael(dec,i,3);
4639 if (liftm) I = FpM_mul(liftm,I,p);
4640 I = shallowmatconcat(mkvec2(I,prad));
4641 gel(res,i) = I;
4642 }
4643
4644 return gerepilecopy(av, res);
4645 }
4646
4647 /* finds a nontrivial ideal of O/prad or gen_0 if there is none. */
4648 static GEN
algpdecompose_i(GEN al,GEN p,GEN zprad,GEN projs)4649 algpdecompose_i(GEN al, GEN p, GEN zprad, GEN projs)
4650 {
4651 pari_sp av = avma;
4652 GEN prad = algpradical_i(al,p,zprad,projs);
4653 return gerepileupto(av, algpdecompose0(al, prad, p, projs));
4654 }
4655
4656 /* ord is assumed to be in hnf wrt the integral basis of al. */
4657 /* assumes that alg_get_invbasis(al) is integral. */
4658 static GEN
alg_change_overorder_shallow(GEN al,GEN ord)4659 alg_change_overorder_shallow(GEN al, GEN ord)
4660 {
4661 GEN al2, mt, iord, mtx, den, den2, div;
4662 long i, n;
4663 n = alg_get_absdim(al);
4664
4665 iord = QM_inv(ord);
4666 al2 = shallowcopy(al);
4667 ord = Q_remove_denom(ord,&den);
4668
4669 gel(al2,7) = Q_remove_denom(gel(al,7), &den2);
4670 if (den2) div = mulii(den,den2);
4671 else div = den;
4672 gel(al2,7) = ZM_Z_div(ZM_mul(gel(al2,7), ord), div);
4673
4674 gel(al2,8) = ZM_mul(iord, gel(al,8));
4675
4676 mt = cgetg(n+1,t_VEC);
4677 gel(mt,1) = matid(n);
4678 div = sqri(den);
4679 for (i=2; i<=n; i++) {
4680 mtx = algbasismultable(al,gel(ord,i));
4681 gel(mt,i) = ZM_mul(iord, ZM_mul(mtx, ord));
4682 gel(mt,i) = ZM_Z_divexact(gel(mt,i), div);
4683 }
4684 gel(al2,9) = mt;
4685
4686 gel(al2,11) = algtracebasis(al2);
4687
4688 return al2;
4689 }
4690
4691 static GEN
algfromcenter(GEN al,GEN x)4692 algfromcenter(GEN al, GEN x)
4693 {
4694 GEN nf = alg_get_center(al);
4695 long n;
4696 switch(alg_type(al)) {
4697 case al_CYCLIC:
4698 n = alg_get_degree(al);
4699 break;
4700 case al_CSA:
4701 n = alg_get_dim(al);
4702 break;
4703 default:
4704 return NULL; /*LCOV_EXCL_LINE*/
4705 }
4706 return algalgtobasis(al, scalarcol(basistoalg(nf, x), n));
4707 }
4708
4709 /* x is an ideal of the center in hnf form */
4710 static GEN
algfromcenterhnf(GEN al,GEN x)4711 algfromcenterhnf(GEN al, GEN x)
4712 {
4713 GEN res;
4714 long i;
4715 res = cgetg(lg(x), t_MAT);
4716 for (i=1; i<lg(x); i++) gel(res,i) = algfromcenter(al, gel(x,i));
4717 return res;
4718 }
4719
4720 /* assumes al is CSA or CYCLIC */
4721 static GEN
algcenter_precompute(GEN al,GEN p)4722 algcenter_precompute(GEN al, GEN p)
4723 {
4724 GEN fa, pdec, nfprad, projs, nf = alg_get_center(al);
4725 long i, np;
4726
4727 pdec = idealprimedec(nf, p);
4728 settyp(pdec, t_COL);
4729 np = lg(pdec)-1;
4730 fa = mkmat2(pdec, const_col(np, gen_1));
4731 if (dvdii(nf_get_disc(nf), p))
4732 nfprad = idealprodprime(nf, pdec);
4733 else
4734 nfprad = scalarmat_shallow(p, nf_get_degree(nf));
4735 fa = idealchineseinit(nf, fa);
4736 projs = cgetg(np+1, t_VEC);
4737 for (i=1; i<=np; i++) gel(projs, i) = idealchinese(nf, fa, vec_ei(np,i));
4738 return mkvec2(nfprad, projs);
4739 }
4740
4741 static GEN
algcenter_prad(GEN al,GEN p,GEN pre)4742 algcenter_prad(GEN al, GEN p, GEN pre)
4743 {
4744 GEN nfprad, zprad, mtprad;
4745 long i;
4746 nfprad = gel(pre,1);
4747 zprad = algfromcenterhnf(al, nfprad);
4748 zprad = FpM_image(zprad, p);
4749 mtprad = cgetg(lg(zprad), t_VEC);
4750 for (i=1; i<lg(zprad); i++) gel(mtprad, i) = algbasismultable(al, gel(zprad,i));
4751 mtprad = shallowmatconcat(mtprad);
4752 zprad = FpM_image(mtprad, p);
4753 return zprad;
4754 }
4755
4756 static GEN
algcenter_p_projs(GEN al,GEN p,GEN pre)4757 algcenter_p_projs(GEN al, GEN p, GEN pre)
4758 {
4759 GEN projs, zprojs;
4760 long i;
4761 projs = gel(pre,2);
4762 zprojs = cgetg(lg(projs), t_VEC);
4763 for (i=1; i<lg(projs); i++) gel(zprojs,i) = FpC_red(algfromcenter(al, gel(projs,i)),p);
4764 return zprojs;
4765 }
4766
4767 /* al is assumed to be simple */
4768 static GEN
alg_pmaximal(GEN al,GEN p)4769 alg_pmaximal(GEN al, GEN p)
4770 {
4771 GEN al2, prad, lord = gen_0, I, id, dec, zprad, projs, pre;
4772 long n, i;
4773 n = alg_get_absdim(al);
4774 id = matid(n);
4775 al2 = al;
4776
4777 dbg_printf(0)("Round 2 (noncommutative) at p=%Ps, dim=%d\n", p, n);
4778
4779 pre = algcenter_precompute(al,p);
4780
4781 while (1) {
4782 zprad = algcenter_prad(al2, p, pre);
4783 projs = algcenter_p_projs(al2, p, pre);
4784 if (lg(projs) == 2) projs = NULL;
4785 prad = algpradical_i(al2,p,zprad,projs);
4786 if (typ(prad) == t_INT) break;
4787 lord = algleftordermodp(al2,prad,p);
4788 if (!cmp_universal(lord,id)) break;
4789 al2 = alg_change_overorder_shallow(al2,lord);
4790 }
4791
4792 dec = algpdecompose0(al2,prad,p,projs);
4793 while (lg(dec)>2) {
4794 for (i=1; i<lg(dec); i++) {
4795 I = gel(dec,i);
4796 lord = algleftordermodp(al2,I,p);
4797 if (cmp_universal(lord,matid(n))) break;
4798 }
4799 if (i==lg(dec)) break;
4800 al2 = alg_change_overorder_shallow(al2,lord);
4801 zprad = algcenter_prad(al2, p, pre);
4802 projs = algcenter_p_projs(al2, p, pre);
4803 if (lg(projs) == 2) projs = NULL;
4804 dec = algpdecompose_i(al2,p,zprad,projs);
4805 }
4806 return al2;
4807 }
4808
4809 static GEN
algtracematrix(GEN al)4810 algtracematrix(GEN al)
4811 {
4812 GEN M, mt;
4813 long n, i, j;
4814 n = alg_get_absdim(al);
4815 mt = alg_get_multable(al);
4816 M = cgetg(n+1, t_MAT);
4817 for (i=1; i<=n; i++)
4818 {
4819 gel(M,i) = cgetg(n+1,t_MAT);
4820 for (j=1; j<=i; j++)
4821 gcoeff(M,j,i) = gcoeff(M,i,j) = algabstrace(al,gmael(mt,i,j));
4822 }
4823 return M;
4824 }
4825 static GEN
algdisc_i(GEN al)4826 algdisc_i(GEN al) { return ZM_det(algtracematrix(al)); }
4827 GEN
algdisc(GEN al)4828 algdisc(GEN al)
4829 {
4830 pari_sp av = avma;
4831 checkalg(al); return gerepileuptoint(av, algdisc_i(al));
4832 }
4833 static GEN
alg_maximal(GEN al)4834 alg_maximal(GEN al)
4835 {
4836 GEN fa = absZ_factor(algdisc_i(al));
4837 return alg_maximal_primes(al, gel(fa,1));
4838 }
4839
4840 /** LATTICES **/
4841
4842 /*
4843 Convention: lattice = [I,t] representing t*I, where
4844 - I integral nonsingular upper-triangular matrix representing a lattice over
4845 the integral basis of the algebra, and
4846 - t>0 either an integer or a rational number.
4847
4848 Recommended and returned by the functions below:
4849 - I HNF and primitive
4850 */
4851
4852 /* TODO use hnfmodid whenever possible using a*O <= I <= O
4853 * for instance a = ZM_det_triangular(I) */
4854
4855 static GEN
primlat(GEN lat)4856 primlat(GEN lat)
4857 {
4858 GEN m, t, c;
4859 m = alglat_get_primbasis(lat);
4860 t = alglat_get_scalar(lat);
4861 m = Q_primitive_part(m,&c);
4862 if (c) return mkvec2(m,gmul(t,c));
4863 return lat;
4864 }
4865
4866 /* assumes the lattice contains d * integral basis, d=0 allowed */
4867 GEN
alglathnf(GEN al,GEN m,GEN d)4868 alglathnf(GEN al, GEN m, GEN d)
4869 {
4870 pari_sp av = avma;
4871 long N,i,j;
4872 GEN m2, c;
4873 checkalg(al);
4874 N = alg_get_absdim(al);
4875 if (!d) d = gen_0;
4876 if (typ(m) == t_VEC) m = matconcat(m);
4877 if (typ(m) == t_COL) m = algleftmultable(al,m);
4878 if (typ(m) != t_MAT) pari_err_TYPE("alglathnf",m);
4879 if (typ(d) != t_FRAC && typ(d) != t_INT) pari_err_TYPE("alglathnf",d);
4880 if (lg(m)-1 < N || lg(gel(m,1))-1 != N) pari_err_DIM("alglathnf");
4881 for (i=1; i<=N; i++)
4882 for (j=1; j<lg(m); j++)
4883 if (typ(gcoeff(m,i,j)) != t_FRAC && typ(gcoeff(m,i,j)) != t_INT)
4884 pari_err_TYPE("alglathnf", gcoeff(m,i,j));
4885 m2 = Q_primitive_part(m,&c);
4886 if (!c) c = gen_1;
4887 if (!signe(d)) d = detint(m2);
4888 else d = gdiv(d,c); /* should be an integer */
4889 if (!signe(d)) pari_err_INV("alglathnf [m does not have full rank]", m2);
4890 m2 = ZM_hnfmodid(m2,d);
4891 return gerepilecopy(av, mkvec2(m2,c));
4892 }
4893
4894 static GEN
prepare_multipliers(GEN * a,GEN * b)4895 prepare_multipliers(GEN *a, GEN *b)
4896 {
4897 GEN na, nb, da, db, d;
4898 na = numer_i(*a); da = denom_i(*a);
4899 nb = numer_i(*b); db = denom_i(*b);
4900 na = mulii(na,db);
4901 nb = mulii(nb,da);
4902 d = gcdii(na,nb);
4903 *a = diviiexact(na,d);
4904 *b = diviiexact(nb,d);
4905 return gdiv(d, mulii(da,db));
4906 }
4907
4908 static GEN
prepare_lat(GEN m1,GEN t1,GEN m2,GEN t2)4909 prepare_lat(GEN m1, GEN t1, GEN m2, GEN t2)
4910 {
4911 GEN d = prepare_multipliers(&t1, &t2);
4912 m1 = ZM_Z_mul(m1,t1);
4913 m2 = ZM_Z_mul(m2,t2);
4914 return mkvec3(m1,m2,d);
4915 }
4916
4917 static GEN
alglataddinter(GEN al,GEN lat1,GEN lat2,GEN * sum,GEN * inter)4918 alglataddinter(GEN al, GEN lat1, GEN lat2, GEN *sum, GEN *inter)
4919 {
4920 GEN d, m1, m2, t1, t2, M, prep, d1, d2, ds, di, K;
4921 checkalg(al);
4922 checklat(al,lat1);
4923 checklat(al,lat2);
4924
4925 m1 = alglat_get_primbasis(lat1);
4926 t1 = alglat_get_scalar(lat1);
4927 m2 = alglat_get_primbasis(lat2);
4928 t2 = alglat_get_scalar(lat2);
4929 prep = prepare_lat(m1, t1, m2, t2);
4930 m1 = gel(prep,1);
4931 m2 = gel(prep,2);
4932 d = gel(prep,3);
4933 M = matconcat(mkvec2(m1,m2));
4934 d1 = ZM_det_triangular(m1);
4935 d2 = ZM_det_triangular(m2);
4936 ds = gcdii(d1,d2);
4937 if (inter)
4938 {
4939 di = diviiexact(mulii(d1,d2),ds);
4940 K = matkermod(M,di,sum);
4941 K = rowslice(K,1,lg(m1));
4942 *inter = hnfmodid(FpM_mul(m1,K,di),di);
4943 if (sum) *sum = hnfmodid(*sum,ds);
4944 }
4945 else *sum = hnfmodid(M,ds);
4946 return d;
4947 }
4948
4949 GEN
alglatinter(GEN al,GEN lat1,GEN lat2,GEN * ptsum)4950 alglatinter(GEN al, GEN lat1, GEN lat2, GEN* ptsum)
4951 {
4952 pari_sp av = avma;
4953 GEN inter, d;
4954 d = alglataddinter(al, lat1, lat2, ptsum, &inter);
4955 inter = primlat(mkvec2(inter, d));
4956 if (ptsum)
4957 {
4958 *ptsum = primlat(mkvec2(*ptsum,d));
4959 gerepileall(av, 2, &inter, ptsum);
4960 }
4961 else inter = gerepilecopy(av, inter);
4962 return inter;
4963 }
4964
4965 GEN
alglatadd(GEN al,GEN lat1,GEN lat2,GEN * ptinter)4966 alglatadd(GEN al, GEN lat1, GEN lat2, GEN* ptinter)
4967 {
4968 pari_sp av = avma;
4969 GEN sum, d;
4970 d = alglataddinter(al, lat1, lat2, &sum, ptinter);
4971 sum = primlat(mkvec2(sum, d));
4972 if (ptinter)
4973 {
4974 *ptinter = primlat(mkvec2(*ptinter,d));
4975 gerepileall(av, 2, &sum, ptinter);
4976 }
4977 else sum = gerepilecopy(av, sum);
4978 return sum;
4979 }
4980
4981 int
alglatsubset(GEN al,GEN lat1,GEN lat2,GEN * ptindex)4982 alglatsubset(GEN al, GEN lat1, GEN lat2, GEN* ptindex)
4983 {
4984 /* TODO version that returns the quotient as abelian group? */
4985 /* return matrices to convert coordinates from one to other? */
4986 pari_sp av = avma;
4987 int res;
4988 GEN m1, m2, m2i, m, t;
4989 checkalg(al);
4990 checklat(al,lat1);
4991 checklat(al,lat2);
4992 m1 = alglat_get_primbasis(lat1);
4993 m2 = alglat_get_primbasis(lat2);
4994 m2i = RgM_inv_upper(m2);
4995 t = gdiv(alglat_get_scalar(lat1), alglat_get_scalar(lat2));
4996 m = RgM_Rg_mul(RgM_mul(m2i,m1), t);
4997 res = RgM_is_ZM(m);
4998 if (res && ptindex)
4999 {
5000 *ptindex = mpabs(ZM_det_triangular(m));
5001 gerepileall(av,1,ptindex);
5002 }
5003 else set_avma(av);
5004 return res;
5005 }
5006
5007 GEN
alglatindex(GEN al,GEN lat1,GEN lat2)5008 alglatindex(GEN al, GEN lat1, GEN lat2)
5009 {
5010 pari_sp av = avma;
5011 long N;
5012 GEN res;
5013 checkalg(al);
5014 checklat(al,lat1);
5015 checklat(al,lat2);
5016 N = alg_get_absdim(al);
5017 res = alglat_get_scalar(lat1);
5018 res = gdiv(res, alglat_get_scalar(lat2));
5019 res = gpowgs(res, N);
5020 res = gmul(res,RgM_det_triangular(alglat_get_primbasis(lat1)));
5021 res = gdiv(res, RgM_det_triangular(alglat_get_primbasis(lat2)));
5022 res = gabs(res,0);
5023 return gerepilecopy(av, res);
5024 }
5025
5026 GEN
alglatmul(GEN al,GEN lat1,GEN lat2)5027 alglatmul(GEN al, GEN lat1, GEN lat2)
5028 {
5029 pari_sp av = avma;
5030 long N,i;
5031 GEN m1, m2, m, V, lat, t, d, dp;
5032 checkalg(al);
5033 if (typ(lat1)==t_COL)
5034 {
5035 if (typ(lat2)==t_COL)
5036 pari_err_TYPE("alglatmul [one of lat1, lat2 has to be a lattice]", lat2);
5037 checklat(al,lat2);
5038 lat1 = Q_remove_denom(lat1,&d);
5039 m = algbasismultable(al,lat1);
5040 m2 = alglat_get_primbasis(lat2);
5041 dp = mulii(detint(m),ZM_det_triangular(m2));
5042 m = ZM_mul(m,m2);
5043 t = alglat_get_scalar(lat2);
5044 if (d) t = gdiv(t,d);
5045 }
5046 else /* typ(lat1)!=t_COL */
5047 {
5048 checklat(al,lat1);
5049 if (typ(lat2)==t_COL)
5050 {
5051 lat2 = Q_remove_denom(lat2,&d);
5052 m = algbasisrightmultable(al,lat2);
5053 m1 = alglat_get_primbasis(lat1);
5054 dp = mulii(detint(m),ZM_det_triangular(m1));
5055 m = ZM_mul(m,m1);
5056 t = alglat_get_scalar(lat1);
5057 if (d) t = gdiv(t,d);
5058 }
5059 else /* typ(lat2)!=t_COL */
5060 {
5061 checklat(al,lat2);
5062 N = alg_get_absdim(al);
5063 m1 = alglat_get_primbasis(lat1);
5064 m2 = alglat_get_primbasis(lat2);
5065 dp = mulii(ZM_det_triangular(m1), ZM_det_triangular(m2));
5066 V = cgetg(N+1,t_VEC);
5067 for (i=1; i<=N; i++) {
5068 gel(V,i) = algbasismultable(al,gel(m1,i));
5069 gel(V,i) = ZM_mul(gel(V,i),m2);
5070 }
5071 m = matconcat(V);
5072 t = gmul(alglat_get_scalar(lat1), alglat_get_scalar(lat2));
5073 }
5074 }
5075
5076 lat = alglathnf(al,m,dp);
5077 gel(lat,2) = gmul(alglat_get_scalar(lat), t);
5078 lat = primlat(lat);
5079 return gerepilecopy(av, lat);
5080 }
5081
5082 int
alglatcontains(GEN al,GEN lat,GEN x,GEN * ptc)5083 alglatcontains(GEN al, GEN lat, GEN x, GEN *ptc)
5084 {
5085 pari_sp av = avma;
5086 GEN m, t, sol;
5087 checkalg(al);
5088 checklat(al,lat);
5089 m = alglat_get_primbasis(lat);
5090 t = alglat_get_scalar(lat);
5091 x = RgC_Rg_div(x,t);
5092 if (!RgV_is_ZV(x)) return gc_bool(av,0);
5093 sol = hnf_solve(m,x);
5094 if (!sol) return gc_bool(av,0);
5095 if (!ptc) return gc_bool(av,1);
5096 *ptc = sol; gerepileall(av,1,ptc); return 1;
5097 }
5098
5099 GEN
alglatelement(GEN al,GEN lat,GEN c)5100 alglatelement(GEN al, GEN lat, GEN c)
5101 {
5102 pari_sp av = avma;
5103 GEN res;
5104 checkalg(al);
5105 checklat(al,lat);
5106 if (typ(c)!=t_COL) pari_err_TYPE("alglatelement", c);
5107 res = ZM_ZC_mul(alglat_get_primbasis(lat),c);
5108 res = RgC_Rg_mul(res, alglat_get_scalar(lat));
5109 return gerepilecopy(av,res);
5110 }
5111
5112 /* idem QM_invimZ, knowing result is contained in 1/c*Z^n */
5113 static GEN
QM_invimZ_mod(GEN m,GEN c)5114 QM_invimZ_mod(GEN m, GEN c)
5115 {
5116 GEN d, m0, K;
5117 m0 = Q_remove_denom(m, &d);
5118 if (d) d = mulii(d,c);
5119 else d = c;
5120 K = matkermod(m0, d, NULL);
5121 if (lg(K)==1) K = scalarmat(d, lg(m)-1);
5122 else K = hnfmodid(K, d);
5123 return RgM_Rg_div(K,c);
5124 }
5125
5126 /* If m is injective, computes a Z-basis of the submodule of elements whose
5127 * image under m is integral */
5128 static GEN
QM_invimZ(GEN m)5129 QM_invimZ(GEN m)
5130 {
5131 return RgM_invimage(m, QM_ImQ_hnf(m));
5132 }
5133
5134 /* An isomorphism of R-modules M_{m,n}(R) -> R^{m*n} */
5135 static GEN
mat2col(GEN M,long m,long n)5136 mat2col(GEN M, long m, long n)
5137 {
5138 long i,j,k,p;
5139 GEN C;
5140 p = m*n;
5141 C = cgetg(p+1,t_COL);
5142 for (i=1,k=1;i<=m;i++)
5143 for (j=1;j<=n;j++,k++)
5144 gel(C,k) = gcoeff(M,i,j);
5145 return C;
5146 }
5147
5148 static GEN
alglattransporter_i(GEN al,GEN lat1,GEN lat2,long right)5149 alglattransporter_i(GEN al, GEN lat1, GEN lat2, long right)
5150 {
5151 GEN m1, m2, m2i, M, MT, mt, t1, t2, T, c;
5152 long N, i;
5153 N = alg_get_absdim(al);
5154 m1 = alglat_get_primbasis(lat1);
5155 m2 = alglat_get_primbasis(lat2);
5156 m2i = RgM_inv_upper(m2);
5157 c = detint(m1);
5158 t1 = alglat_get_scalar(lat1);
5159 m1 = RgM_Rg_mul(m1,t1);
5160 t2 = alglat_get_scalar(lat2);
5161 m2i = RgM_Rg_div(m2i,t2);
5162
5163 MT = right? NULL: alg_get_multable(al);
5164 M = cgetg(N+1, t_MAT);
5165 for (i=1; i<=N; i++) {
5166 if (right) mt = algbasisrightmultable(al, vec_ei(N,i));
5167 else mt = gel(MT,i);
5168 mt = RgM_mul(m2i,mt);
5169 mt = RgM_mul(mt,m1);
5170 gel(M,i) = mat2col(mt, N, N);
5171 }
5172
5173 c = gdiv(t2,gmul(c,t1));
5174 c = denom_i(c);
5175 T = QM_invimZ_mod(M,c);
5176 return primlat(mkvec2(T,gen_1));
5177 }
5178
5179 /*
5180 { x in al | x*lat1 subset lat2}
5181 */
5182 GEN
alglatlefttransporter(GEN al,GEN lat1,GEN lat2)5183 alglatlefttransporter(GEN al, GEN lat1, GEN lat2)
5184 {
5185 pari_sp av = avma;
5186 checkalg(al);
5187 checklat(al,lat1);
5188 checklat(al,lat2);
5189 return gerepilecopy(av, alglattransporter_i(al,lat1,lat2,0));
5190 }
5191
5192 /*
5193 { x in al | lat1*x subset lat2}
5194 */
5195 GEN
alglatrighttransporter(GEN al,GEN lat1,GEN lat2)5196 alglatrighttransporter(GEN al, GEN lat1, GEN lat2)
5197 {
5198 pari_sp av = avma;
5199 checkalg(al);
5200 checklat(al,lat1);
5201 checklat(al,lat2);
5202 return gerepilecopy(av, alglattransporter_i(al,lat1,lat2,1));
5203 }
5204
5205 GEN
algmakeintegral(GEN mt0,long maps)5206 algmakeintegral(GEN mt0, long maps)
5207 {
5208 pari_sp av = avma;
5209 long n,i;
5210 GEN m,P,Pi,mt2,mt;
5211 n = lg(mt0)-1;
5212 mt = check_mt(mt0,NULL);
5213 if (!mt) pari_err_TYPE("algmakeintegral", mt0);
5214 if (isint1(Q_denom(mt0))) {
5215 if (maps) mt = mkvec3(mt,matid(n),matid(n));
5216 return gerepilecopy(av,mt);
5217 }
5218 dbg_printf(2)(" algmakeintegral: dim=%d, denom=%Ps\n", n, Q_denom(mt0));
5219 m = cgetg(n+1,t_MAT);
5220 for (i=1;i<=n;i++)
5221 gel(m,i) = mat2col(gel(mt,i),n,n);
5222 dbg_printf(2)(" computing order, dims m = %d x %d...\n", nbrows(m), lg(m)-1);
5223 P = QM_invimZ(m);
5224 dbg_printf(2)(" ...done.\n");
5225 P = shallowmatconcat(mkvec2(col_ei(n,1),P));
5226 P = hnf(P);
5227 Pi = RgM_inv(P);
5228 mt2 = change_Rgmultable(mt,P,Pi);
5229 if (maps) mt2 = mkvec3(mt2,Pi,P); /* mt2, mt->mt2, mt2->mt */
5230 return gerepilecopy(av,mt2);
5231 }
5232
5233 /** ORDERS **/
5234
5235 /** IDEALS **/
5236
5237