Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - modules - algebras.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21353-12523aa) Lines: 2730 2893 94.4 %
Date: 2017-11-24 06:20:58 Functions: 262 273 96.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11