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 21061-8feaff2) Lines: 2727 2890 94.4 %
Date: 2017-09-24 06:24:57 Functions: 261 272 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     1263850 : checkalg_i(GEN al)
      38             : {
      39             :   GEN mt;
      40     1263850 :   if (typ(al) != t_VEC || lg(al) != 12) return 0;
      41     1263661 :   mt = alg_get_multable(al);
      42     1263661 :   if (typ(mt) != t_VEC || lg(mt) == 1 || typ(gel(mt,1)) != t_MAT) return 0;
      43     1263640 :   if (!isintzero(alg_get_splittingfield(al)) && gequal0(alg_get_char(al))) {
      44      834988 :     if (typ(gel(al,2)) != t_VEC || lg(gel(al,2)) == 1) return 0;
      45      834981 :     checkrnf(alg_get_splittingfield(al));
      46             :   }
      47     1263612 :   return 1;
      48             : }
      49             : void
      50     1263108 : checkalg(GEN al)
      51     1263108 : { if (!checkalg_i(al)) pari_err_TYPE("checkalg [please apply alginit()]",al); }
      52             : 
      53             : static int
      54      361060 : checklat_i(GEN al, GEN lat)
      55             : {
      56             :   long N,i,j;
      57             :   GEN m,t,c;
      58      361060 :   if (typ(lat)!=t_VEC || lg(lat) != 3) return 0;
      59      361060 :   t = gel(lat,2);
      60      361060 :   if (typ(t) != t_INT && typ(t) != t_FRAC) return 0;
      61      361060 :   if (gsigne(t)<=0) return 0;
      62      361060 :   m = gel(lat,1);
      63      361060 :   if (typ(m) != t_MAT) return 0;
      64      361060 :   N = algabsdim(al);
      65      361060 :   if (lg(m)-1 != N || lg(gel(m,1))-1 != N) return 0;
      66     3249540 :   for (i=1; i<=N; i++)
      67    25996320 :     for (j=1; j<=N; j++) {
      68    23107840 :       c = gcoeff(m,i,j);
      69    23107840 :       if (typ(c) != t_INT) return 0;
      70    23107840 :       if (j<i && signe(gcoeff(m,i,j))) return 0;
      71             :     }
      72      361060 :   return 1;
      73             : }
      74      361060 : void checklat(GEN al, GEN lat)
      75      361060 : { if(!checklat_i(al,lat)) pari_err_TYPE("checklat [please apply alglathnf()]", lat); }
      76             : 
      77             : 
      78             : /**  ACCESSORS  **/
      79             : long
      80     3524878 : alg_type(GEN al)
      81             : {
      82     3524878 :   if (isintzero(alg_get_splittingfield(al)) || !gequal0(alg_get_char(al))) return al_TABLE;
      83     2114434 :   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     1956087 :     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       65373 : alg_get_dim(GEN al)
     100             : {
     101             :   long d;
     102       65373 :   switch(alg_type(al)) {
     103       30387 :     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        2163 : algdim(GEN al)
     112        2163 : { checkalg(al); return alg_get_dim(al); }
     113             : 
     114             : long
     115     1508381 : alg_get_absdim(GEN al)
     116             : {
     117     1508381 :   switch(alg_type(al)) {
     118      728014 :     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      763378 :       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      375277 : algabsdim(GEN al)
     128      375277 : { checkalg(al); return alg_get_absdim(al); }
     129             : 
     130             : /* only cyclic */
     131             : GEN
     132        9478 : alg_get_auts(GEN al)
     133             : {
     134        9478 :   if (alg_type(al) != al_CYCLIC)
     135           0 :     pari_err_TYPE("alg_get_auts [non-cyclic algebra]", al);
     136        9478 :   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        9499 : alg_get_b(GEN al)
     149             : {
     150        9499 :   if (alg_type(al) != al_CYCLIC)
     151           7 :     pari_err_TYPE("alg_get_b [non-cyclic algebra]", al);
     152        9492 :   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     7445466 : 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      941122 : alg_get_degree(GEN al)
     206             : {
     207             :   long ta;
     208      941122 :   ta = alg_type(al);
     209      941122 :   if (ta != al_CYCLIC && ta != al_CSA)
     210          21 :     pari_err_TYPE("alg_get_degree [use alginit]",al);
     211      941101 :   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       62636 : alg_get_center(GEN al)
     222             : {
     223             :   long ta;
     224       62636 :   ta = alg_type(al);
     225       62636 :   if (ta != al_CSA && ta != al_CYCLIC)
     226           7 :     pari_err_TYPE("alg_get_center [use alginit]",al);
     227       62629 :   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       20636 : alg_get_abssplitting(GEN al)
     239             : {
     240       20636 :   long ta = alg_type(al), prec;
     241       20636 :   if (ta != al_CYCLIC && ta != al_CSA)
     242           0 :     pari_err_TYPE("alg_get_abssplitting [use alginit]",al);
     243       20636 :   prec = nf_get_prec(alg_get_center(al));
     244       20636 :   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        1715 : alg_get_hasse_f(GEN al)
     259             : {
     260        1715 :   long ta = alg_type(al);
     261        1715 :   if (ta != al_CYCLIC && ta != al_CSA)
     262           7 :     pari_err_TYPE("alg_get_hasse_f [use alginit]",al);
     263        1708 :   if (ta == al_CSA) pari_err_IMPL("computation of Hasse invariants over table CSA");
     264        1701 :   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        7231 : 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     2780554 : 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     4322150 : 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      283101 : alg_get_tracebasis(GEN al) { return gel(al,11); }
     288             : 
     289             : /* lattices */
     290             : GEN
     291      487312 : alglat_get_primbasis(GEN lat) { return gel(lat,1); }
     292             : GEN
     293      578382 : 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       35021 : RgC_is_ZC(GEN x)
     311             : {
     312             :   int i;
     313      315189 :   for (i=lg(x)-1; i>0; i--)
     314      280168 :     if (typ(gel(x,i)) != t_INT) return 0;
     315       35021 :   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        1344 :     for (i=1; i<lg(E); i++)
     373             :     {
     374         833 :       long v = nfval(nf,y0,gel(P,i));
     375         833 :       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        2555 : Flm_invimage_init(GEN M, ulong p)
     468             : {
     469        2555 :   GEN v = Flm_indexrank(M, p), perm = gel(v,1);
     470        2555 :   GEN MM = rowpermute(M, perm); /* square invertible */
     471        2555 :   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      376495 : Flm_invimage_pre(GEN v, GEN y, ulong p)
     476             : {
     477      376495 :   GEN inv = gel(v,1), perm = gel(v,2);
     478      376495 :   return Flm_Flc_mul(inv, vecsmallpermute(y, perm), p);
     479             : }
     480             : 
     481             : GEN
     482        3003 : algradical(GEN al)
     483             : {
     484        3003 :   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        3003 :   checkalg(al);
     489        3003 :   P = alg_get_char(al);
     490        3003 :   mt = alg_get_multable(al);
     491        3003 :   n = alg_get_absdim(al);
     492        3003 :   dbg_printf(1)("algradical: char=%Ps, dim=%d\n", P, n);
     493        3003 :   traces = algtracematrix(al);
     494        3003 :   if (!signe(P))
     495             :   {
     496         378 :     dbg_printf(2)(" char 0, computing kernel...\n");
     497         378 :     K = ker(traces);
     498         378 :     dbg_printf(2)(" ...done.\n");
     499         378 :     ni = lg(K)-1; if (!ni) { avma = av; return gen_0; }
     500          63 :     return gerepileupto(av, K);
     501             :   }
     502        2625 :   dbg_printf(2)(" char>0, computing kernel...\n");
     503        2625 :   K = FpM_ker(traces, P);
     504        2625 :   dbg_printf(2)(" ...done.\n");
     505        2625 :   ni = lg(K)-1; if (!ni) { avma = av; return gen_0; }
     506        1862 :   if (abscmpiu(P,n)>0) return gerepileupto(av, K);
     507             : 
     508             :   /* tough case, p <= n. Ronyai's algorithm */
     509        1484 :   p = P[2]; l = 1;
     510        1484 :   expo = p; modu = p*p;
     511        1484 :   dbg_printf(2)(" char>0, hard case.\n");
     512        1484 :   while (modu<=(ulong)n) { l++; modu *= p; }
     513        1484 :   MT = ZMV_to_FlmV(mt, modu);
     514        1484 :   I = ZM_to_Flm(K,p); /* I_0 */
     515        3864 :   for (i=1; i<=l; i++) {/*compute I_i, expo = p^i, modu = p^(l+1) > n*/
     516             :     long j, lig,col;
     517        2555 :     GEN v = cgetg(ni+1, t_VECSMALL);
     518        2555 :     GEN invI = Flm_invimage_init(I, p);
     519        2555 :     dbg_printf(2)(" computing I_%d:\n", i);
     520        2555 :     traces = cgetg(ni+1,t_MAT);
     521       22981 :     for (j = 1; j <= ni; j++)
     522             :     {
     523       20426 :       GEN M = algbasismultable_Flm(MT, gel(I,j), modu);
     524       20426 :       uel(v,j) = algtracei(M, p,expo,modu);
     525             :     }
     526       22981 :     for (col=1; col<=ni; col++)
     527             :     {
     528       20426 :       GEN t = cgetg(n+1,t_VECSMALL); gel(traces,col) = t;
     529       20426 :       x = gel(I, col); /*col-th basis vector of I_{i-1}*/
     530      396921 :       for (lig=1; lig<=n; lig++)
     531             :       {
     532      376495 :         GEN y = _tablemul_ej_Fl(MT,x,lig,p);
     533      376495 :         GEN z = Flm_invimage_pre(invI, y, p);
     534      376495 :         uel(t,lig) = Flv_dotproduct(v, z, p);
     535             :       }
     536             :     }
     537        2555 :     dbg_printf(2)(" computing kernel...\n");
     538        2555 :     K = Flm_ker(traces, p);
     539        2555 :     dbg_printf(2)(" ...done.\n");
     540        2555 :     ni = lg(K)-1; if (!ni) { avma = av; return gen_0; }
     541        2380 :     I = Flm_mul(I,K,p);
     542        2380 :     expo *= p;
     543             :   }
     544        1309 :   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       19733 : alg_quotient0(GEN al, GEN S, GEN Si, long nq, GEN p, int maps)
     581             : {
     582       19733 :   GEN mt = cgetg(nq+1,t_VEC), P, Pi, d;
     583             :   long i;
     584       19733 :   dbg_printf(3)("  alg_quotient0: char=%Ps, dim=%d, dim I=%d\n", p, alg_get_absdim(al), lg(S)-1);
     585       86527 :   for (i=1; i<=nq; i++) {
     586       66794 :     GEN mti = algleftmultable(al,gel(S,i));
     587       66794 :     if (signe(p)) gel(mt,i) = FpM_mul(Si, FpM_mul(mti,S,p), p);
     588        8379 :     else          gel(mt,i) = RgM_mul(Si, RgM_mul(mti,S));
     589             :   }
     590       19733 :   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       19733 :   al = algtableinit_i(mt,p);
     601       19733 :   if (maps) al = mkvec3(al,Si,S); /*algebra, proj, lift*/
     602       19733 :   return al;
     603             : }
     604             : 
     605             : /*quotient of an algebra by a nontrivial two-sided ideal*/
     606             : GEN
     607        1414 : alg_quotient(GEN al, GEN I, int maps)
     608             : {
     609        1414 :   pari_sp av = avma;
     610             :   GEN p, IS, ISi, S, Si;
     611             :   long n, ni;
     612             : 
     613        1414 :   checkalg(al);
     614        1414 :   p = alg_get_char(al);
     615        1414 :   n = alg_get_absdim(al);
     616        1414 :   ni = lg(I)-1;
     617             : 
     618             :   /*force first vector of complement to be the identity*/
     619        1414 :   IS = shallowconcat(I, gcoeff(alg_get_multable(al),1,1));
     620        1414 :   if (signe(p)) {
     621        1393 :     IS = FpM_suppl(IS,p);
     622        1393 :     ISi = FpM_inv(IS,p);
     623             :   }
     624             :   else {
     625          21 :     IS = suppl(IS);
     626          21 :     ISi = RgM_inv(IS);
     627             :   }
     628        1414 :   S = vecslice(IS, ni+1, n);
     629        1414 :   Si = rowslice(ISi, ni+1, n);
     630        1414 :   return gerepilecopy(av, alg_quotient0(al, S, Si, n-ni, p, maps));
     631             : }
     632             : 
     633             : static GEN
     634       18333 : 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       18333 :   if (gequal0(gel(m,1))) return zeromat(nbrows(m),0);
     639             : 
     640       18319 :   if (signe(p)) ir = FpM_indexrank(m,p);
     641        1848 :   else          ir = indexrank(m);
     642             : 
     643       18319 :   icol = gel(ir,2);
     644       18319 :   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        9100 : alg_centralproj(GEN al, GEN z, int maps)
     669             : {
     670        9100 :   pari_sp av = avma;
     671             :   GEN S, U, Ui, alq, p;
     672        9100 :   long i, iu, lz = lg(z);
     673             : 
     674        9100 :   checkalg(al);
     675        9100 :   if (typ(z) != t_VEC) pari_err_TYPE("alcentralproj",z);
     676        9093 :   p = alg_get_char(al);
     677        9093 :   dbg_printf(3)("  alg_centralproj: char=%Ps, dim=%d, #z=%d\n", p, alg_get_absdim(al), lz-1);
     678        9093 :   S = cgetg(lz,t_VEC); /*S[i] = Im(z_i)*/
     679       27426 :   for (i=1; i<lz; i++)
     680             :   {
     681       18333 :     GEN mti = algleftmultable(al, gel(z,i));
     682       18333 :     gel(S,i) = image_keep_first(mti,p);
     683             :   }
     684        9093 :   U = shallowconcat1(S); /*U = [Im(z_1)|Im(z_2)|...|Im(z_nz)], n x n*/
     685        9093 :   if (lg(U)-1 < alg_get_absdim(al)) pari_err_TYPE("alcentralproj [z[i]'s not surjective]",z);
     686        9086 :   if (signe(p)) Ui = FpM_inv(U,p);
     687         924 :   else          Ui = RgM_inv(U);
     688             :   if (!Ui) pari_err_BUG("alcentralproj"); /*LCOV_EXCL_LINE*/
     689             : 
     690        9086 :   alq = cgetg(lz,t_VEC);
     691       27405 :   for (iu=0,i=1; i<lz; i++)
     692             :   {
     693       18319 :     long nq = lg(gel(S,i))-1, ju = iu + nq;
     694       18319 :     GEN Si = rowslice(Ui, iu+1, ju);
     695       18319 :     gel(alq, i) = alg_quotient0(al,gel(S,i),Si,nq,p,maps);
     696       18319 :     iu = ju;
     697             :   }
     698        9086 :   return gerepilecopy(av, alq);
     699             : }
     700             : 
     701             : /* al is an al_TABLE */
     702             : static GEN
     703       29750 : algtablecenter(GEN al)
     704             : {
     705       29750 :   pari_sp av = avma;
     706             :   long n, i, j, k, ic;
     707             :   GEN C, cij, mt, p;
     708             : 
     709       29750 :   n = alg_get_absdim(al);
     710       29750 :   mt = alg_get_multable(al);
     711       29750 :   p = alg_get_char(al);
     712       29750 :   C = cgetg(n+1,t_MAT);
     713       91084 :   for (j=1; j<=n; j++)
     714             :   {
     715       61334 :     gel(C,j) = cgetg(n*n-n+1,t_COL);
     716       61334 :     ic = 1;
     717      505120 :     for (i=2; i<=n; i++) {
     718      443786 :       if (signe(p)) cij = FpC_sub(gmael(mt,i,j),gmael(mt,j,i),p);
     719      109732 :       else          cij = RgC_sub(gmael(mt,i,j),gmael(mt,j,i));
     720      443786 :       for (k=1; k<=n; k++, ic++) gcoeff(C,ic,j) = gel(cij, k);
     721             :     }
     722             :   }
     723       29750 :   if (signe(p)) return gerepileupto(av, FpM_ker(C,p));
     724        2744 :   else          return gerepileupto(av, ker(C));
     725             : }
     726             : 
     727             : GEN
     728        1603 : algcenter(GEN al)
     729             : {
     730        1603 :   checkalg(al);
     731        1603 :   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        2940 : algprimesubalg(GEN al)
     738             : {
     739        2940 :   pari_sp av = avma;
     740             :   GEN p, Z, F, K;
     741             :   long nz, i;
     742        2940 :   checkalg(al);
     743        2940 :   p = alg_get_char(al);
     744        2940 :   if (!signe(p)) pari_err_DOMAIN("algprimesubalg","characteristic","=",gen_0,p);
     745             : 
     746        2926 :   Z = algtablecenter(al);
     747        2926 :   nz = lg(Z)-1;
     748        2926 :   if (nz==1) return Z;
     749             : 
     750        2338 :   F = cgetg(nz+1, t_MAT);
     751       14609 :   for (i=1; i<=nz; i++) {
     752       12271 :     GEN zi = gel(Z,i);
     753       12271 :     gel(F,i) = FpC_sub(algpow(al,zi,p),zi,p);
     754             :   }
     755        2338 :   K = FpM_ker(F,p);
     756        2338 :   return gerepileupto(av, FpM_mul(Z,K,p));
     757             : }
     758             : 
     759             : 
     760             : static GEN
     761        8463 : _FpX_mul(void* D, GEN x, GEN y) { return FpX_mul(x,y,(GEN)D); }
     762             : static GEN
     763       23443 : _FpX_pow(void* D, GEN x, GEN n) { return FpX_powu(x,itos(n),(GEN)D); }
     764             : static GEN
     765       14980 : FpX_factorback(GEN fa, GEN p)
     766             : {
     767       14980 :   return gen_factorback(gel(fa,1), zv_to_ZV(gel(fa,2)), &_FpX_mul, &_FpX_pow, (void*)p);
     768             : }
     769             : 
     770             : static GEN
     771       16723 : out_decompose(GEN t, GEN Z, GEN P, GEN p)
     772             : {
     773       16723 :   GEN ali = gel(t,1), projm = gel(t,2), liftm = gel(t,3), pZ;
     774       16723 :   if (signe(p)) pZ = FpM_image(FpM_mul(projm,Z,p),p);
     775        1771 :   else          pZ = image(RgM_mul(projm,Z));
     776       16723 :   return mkvec5(ali, projm, liftm, pZ, P);
     777             : }
     778             : /* fa factorization of charpol(x) */
     779             : static GEN
     780        8386 : alg_decompose_from_facto(GEN al, GEN x, GEN fa, GEN Z, int mini)
     781             : {
     782        8386 :   long k = lgcols(fa)-1, k2 = mini? 1: k/2;
     783        8386 :   GEN v1 = rowslice(fa,1,k2);
     784        8386 :   GEN v2 = rowslice(fa,k2+1,k);
     785        8386 :   GEN alq, P,Q, mx, p = alg_get_char(al);
     786        8386 :   dbg_printf(3)("  alg_decompose_from_facto\n");
     787        8386 :   if (signe(p)) {
     788        7490 :     P = FpX_factorback(v1, p);
     789        7490 :     Q = FpX_factorback(v2, p);
     790        7490 :     P = FpX_mul(P, FpXQ_inv(P,Q,p), p);
     791             :   }
     792             :   else {
     793         896 :     P = factorback(v1);
     794         896 :     Q = factorback(v2);
     795         896 :     P = RgX_mul(P, RgXQ_inv(P,Q));
     796             :   }
     797        8386 :   mx = algleftmultable(al, x);
     798        8386 :   P = algpoleval(al, P, mx);
     799        8386 :   if (signe(p)) Q = FpC_sub(col_ei(lg(P)-1,1), P, p);
     800         896 :   else          Q = gsub(gen_1, P);
     801        8386 :   if(gequal0(P) || gequal0(Q)) return NULL;
     802        8386 :   alq = alg_centralproj(al, mkvec2(P,Q), 1);
     803             : 
     804        8386 :   P = out_decompose(gel(alq,1), Z, P, p); if (mini) return P;
     805        8337 :   Q = out_decompose(gel(alq,2), Z, Q, p);
     806        8337 :   return mkvec2(P,Q);
     807             : }
     808             : 
     809             : static GEN
     810        8722 : random_pm1(long n)
     811             : {
     812        8722 :   GEN z = cgetg(n+1,t_VECSMALL);
     813             :   long i;
     814        8722 :   for (i = 1; i <= n; i++) z[i] = random_bits(5)%3 - 1;
     815        8722 :   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       10388 : try_fact(GEN al, GEN x, GEN zx, GEN Z, GEN Zal, long mini)
     823             : {
     824       10388 :   GEN z, dec0, dec1, cp = algcharpoly(Zal,zx,0), fa, p = alg_get_char(al);
     825             :   long nfa, e;
     826       10388 :   dbg_printf(3)("  try_fact: zx=%Ps\n", zx);
     827       10388 :   if (signe(p)) fa = FpX_factor(cp,p);
     828        1540 :   else          fa = factor(cp);
     829       10388 :   dbg_printf(3)("  charpoly=%Ps\n", fa);
     830       10388 :   nfa = nbrows(fa);
     831       10388 :   if (nfa == 1) {
     832        2002 :     if (signe(p)) e = gel(fa,2)[1];
     833         644 :     else          e = itos(gcoeff(fa,1,2));
     834        2002 :     return e==1 ? gen_0 : NULL;
     835             :   }
     836        8386 :   dec0 = alg_decompose_from_facto(al, x, fa, Z, mini);
     837        8386 :   if (!dec0) return NULL;
     838        8386 :   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       18648 : alg_decompose(GEN al, GEN Z, int mini)
     864             : {
     865             :   pari_sp av;
     866             :   GEN Zal, x, zx, rand, dec0, B, p;
     867       18648 :   long i, nz = lg(Z)-1;
     868             : 
     869       18648 :   if (nz==1) return gen_0;
     870        8722 :   p = alg_get_char(al);
     871        8722 :   dbg_printf(2)(" alg_decompose: char=%Ps, dim=%d, dim Z=%d\n", p, alg_get_absdim(al), nz);
     872        8722 :   Zal = alg_subalg(al,Z);
     873        8722 :   Z = gel(Zal,2);
     874        8722 :   Zal = gel(Zal,1);
     875        8722 :   av = avma;
     876        8722 :   rand = random_pm1(nz);
     877        8722 :   zx = zc_to_ZC(rand);
     878        8722 :   if (signe(p)) {
     879        7490 :     zx = FpC_red(zx,p);
     880        7490 :     x = ZM_zc_mul(Z,rand);
     881        7490 :     x = FpC_red(x,p);
     882             :   }
     883        1232 :   else x = RgM_zc_mul(Z,rand);
     884        8722 :   dec0 = try_fact(al,x,zx,Z,Zal,mini);
     885        8722 :   if (dec0) return dec0;
     886        1610 :   avma = av;
     887        1666 :   for (i=2; i<=nz; i++)
     888             :   {
     889        1659 :     dec0 = try_fact(al,gel(Z,i),col_ei(nz,i),Z,Zal,mini);
     890        1659 :     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       18536 : alg_decompose_total(GEN al, GEN Z, int maps)
     905             : {
     906             :   GEN dec, sc, p;
     907             :   long i;
     908             : 
     909       18536 :   dec = alg_decompose(al, Z, 0);
     910       18536 :   if (isintzero(dec))
     911             :   {
     912       10199 :     if (maps) {
     913        5985 :       long n = alg_get_absdim(al);
     914        5985 :       al = mkvec3(al, matid(n), matid(n));
     915             :     }
     916       10199 :     return mkvec(al);
     917             :   }
     918        8337 :   p = alg_get_char(al); if (!signe(p)) p = NULL;
     919        8337 :   sc = cgetg(lg(dec), t_VEC);
     920       25011 :   for (i=1; i<lg(sc); i++) {
     921       16674 :     GEN D = gel(dec,i), a = gel(D,1), Za = gel(D,4);
     922       16674 :     GEN S = alg_decompose_total(a, Za, maps);
     923       16674 :     gel(sc,i) = S;
     924       16674 :     if (maps)
     925             :     {
     926        9590 :       GEN projm = gel(D,2), liftm = gel(D,3);
     927        9590 :       long j, lS = lg(S);
     928       26355 :       for (j=1; j<lS; j++)
     929             :       {
     930       16765 :         GEN Sj = gel(S,j), p2 = gel(Sj,2), l2 = gel(Sj,3);
     931       16765 :         if (p) p2 = FpM_mul(p2, projm, p);
     932           0 :         else   p2 = RgM_mul(p2, projm);
     933       16765 :         if (p) l2 = FpM_mul(liftm, l2, p);
     934           0 :         else   l2 = RgM_mul(liftm, l2);
     935       16765 :         gel(Sj,2) = p2;
     936       16765 :         gel(Sj,3) = l2;
     937             :       }
     938             :     }
     939             :   }
     940        8337 :   return shallowconcat1(sc);
     941             : }
     942             : 
     943             : static GEN
     944        8757 : alg_subalg(GEN al, GEN basis)
     945             : {
     946        8757 :   GEN invbasis, mt, p = alg_get_char(al), al2;
     947        8757 :   long i, j, n = lg(basis)-1;
     948        8757 :   if (!signe(p)) p = NULL;
     949        8757 :   basis = shallowmatconcat(mkvec2(col_ei(n,1),basis));
     950             :   /* 1st column, being e1, is kept in 1st position when computing the image */
     951        8757 :   if (p)    basis = FpM_image(basis,p);
     952        1246 :   else      basis = QM_ImQ_hnf(basis);
     953        8757 :   if (p) { /*TODO change after bugfix?*/
     954        7511 :     GEN complbasis = FpM_suppl(basis,p);
     955        7511 :     invbasis = rowslice(FpM_inv(complbasis,p),1,n);
     956             :   }
     957        1246 :   else invbasis = RgM_inv(basis);
     958        8757 :   mt = cgetg(n+1,t_VEC);
     959        8757 :   gel(mt,1) = matid(n);
     960       33075 :   for(i=2; i<=n; i++) {
     961       24318 :     GEN mtx = cgetg(n+1,t_MAT), x = gel(basis,i);
     962       24318 :     gel(mtx,1) = col_ei(n,i);
     963      170198 :     for(j=2; j<=n; j++) {
     964      145880 :       GEN xy = algmul(al, x, gel(basis,j));
     965      145880 :       if (p) gel(mtx,j) = FpM_FpC_mul(invbasis, xy, p);
     966       43505 :       else   gel(mtx,j) = RgM_RgC_mul(invbasis, xy);
     967             :     }
     968       24318 :     gel(mt,i) = mtx;
     969             :   }
     970        8757 :   al2 = algtableinit_i(mt,p);
     971        8757 :   al2 = mkvec2(al2,basis);
     972        8757 :   return al2;
     973             : }
     974             : 
     975             : GEN
     976          42 : algsubalg(GEN al, GEN basis)
     977             : {
     978          42 :   pari_sp av = avma;
     979             :   GEN p;
     980          42 :   checkalg(al);
     981          42 :   if (typ(basis) != t_MAT) pari_err_TYPE("algsubalg",basis);
     982          35 :   p = alg_get_char(al);
     983          35 :   if (signe(p)) basis = RgM_to_FpM(basis,p);
     984          35 :   return gerepilecopy(av, alg_subalg(al,basis));
     985             : }
     986             : 
     987             : static int
     988       14182 : cmp_algebra(GEN x, GEN y)
     989             : {
     990       14182 :   long d = alg_get_dim(x) - alg_get_dim(y);
     991       14182 :   if (d) return d < 0? -1: 1;
     992       12439 :   d = lg(algtablecenter(x))-lg(algtablecenter(y));/* TODO precompute and store, don't compute every time when sorting */
     993       12439 :   if (d) return d < 0? -1: 1;
     994       12439 :   return cmp_universal(alg_get_multable(x), alg_get_multable(y));
     995             : }
     996             : static int
     997        7637 : cmp_algebra_maps(GEN x, GEN y)
     998        7637 : { return cmp_algebra(gel(x,1), gel(y,1)); }
     999             : 
    1000             : GEN
    1001        3066 : algsimpledec(GEN al, int maps)
    1002             : {
    1003        3066 :   pari_sp av = avma;
    1004             :   GEN Z, p, res;
    1005             :   long n;
    1006        3066 :   checkalg(al);
    1007        3066 :   p = alg_get_char(al);
    1008        3066 :   dbg_printf(1)("algsimpledec: char=%Ps, dim=%d\n", p, alg_get_absdim(al));
    1009        3066 :   if (signe(p)) Z = algprimesubalg(al);
    1010         287 :   else          Z = algtablecenter(al);
    1011             : 
    1012        3066 :   if (lg(Z) == 2) {/*dim Z = 1*/
    1013        1204 :     n = alg_get_absdim(al);
    1014        1204 :     avma = av;
    1015        1204 :     if (!maps) return mkveccopy(al);
    1016        1078 :     retmkvec(mkvec3(gcopy(al), matid(n), matid(n)));
    1017             :   }
    1018        1862 :   res = alg_decompose_total(al, Z, maps);
    1019        1862 :   gen_sort_inplace(res, (void*)(maps? &cmp_algebra_maps: &cmp_algebra),
    1020             :                    &cmp_nodata, NULL);
    1021        1862 :   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       30674 : check_mt(GEN mt, GEN p)
    1040             : {
    1041             :   long i, l;
    1042       30674 :   GEN MT = cgetg_copy(mt, &l);
    1043       30674 :   if (typ(MT) != t_VEC || l == 1) return NULL;
    1044      149919 :   for (i = 1; i < l; i++)
    1045             :   {
    1046      119315 :     GEN M = gel(mt,i);
    1047      119315 :     if (typ(M) != t_MAT || lg(M) != l || lgcols(M) != l) return NULL;
    1048      119294 :     if (p) M = RgM_to_FpM(M,p);
    1049      119294 :     if (i > 1 && ZC_is_ei(gel(M,1)) != i) return NULL; /* i = 1 checked at end*/
    1050      119266 :     gel(MT,i) = M;
    1051             :   }
    1052       30604 :   if (!ZM_isidentity(gel(MT,1))) return NULL;
    1053       30597 :   return MT;
    1054             : }
    1055             : 
    1056             : 
    1057             : int
    1058         539 : algisassociative(GEN mt0, GEN p)
    1059             : {
    1060         539 :   pari_sp av = avma;
    1061             :   long i, j, k, n;
    1062             :   GEN M, mt;
    1063             : 
    1064         539 :   if (checkalg_i(mt0)) { p = alg_get_char(mt0); mt0 = alg_get_multable(mt0); }
    1065         539 :   if (typ(p) != t_INT) pari_err_TYPE("algisassociative",p);
    1066         532 :   mt = check_mt(mt0, isintzero(p)? NULL: p);
    1067         532 :   if (!mt) pari_err_TYPE("algisassociative (mult. table)", mt0);
    1068         483 :   n = lg(mt)-1;
    1069         483 :   M = cgetg(n+1,t_MAT);
    1070         483 :   for (j=1; j<=n; j++) gel(M,j) = cgetg(n+1,t_COL);
    1071        4522 :   for (i=1; i<=n; i++)
    1072             :   {
    1073        4039 :     GEN mi = gel(mt,i);
    1074        4039 :     for (j=1; j<=n; j++) gcoeff(M,i,j) = gel(mi,j); /* ei.ej */
    1075             :   }
    1076        4025 :   for (i=2; i<=n; i++) {
    1077        3549 :     GEN mi = gel(mt,i);
    1078       43687 :     for (j=2; j<=n; j++) {
    1079      582589 :       for (k=2; k<=n; k++) {
    1080             :         GEN x, y;
    1081      542451 :         if (signe(p)) {
    1082      442939 :           x = _tablemul_ej_Fp(mt,gcoeff(M,i,j),k,p);
    1083      442939 :           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      542451 :         if (!gequal(x,y)) { avma = av; return 0; }
    1091             :       }
    1092             :     }
    1093             :   }
    1094         476 :   avma = av; return 1;
    1095             : }
    1096             : 
    1097             : int
    1098         420 : algiscommutative(GEN al) /* assumes e_1 = 1 */
    1099             : {
    1100             :   long i,j,k,N,sp;
    1101             :   GEN mt,a,b,p;
    1102         420 :   checkalg(al);
    1103         420 :   if(alg_type(al) != al_TABLE) return alg_get_degree(al)==1;
    1104         378 :   N = alg_get_absdim(al);
    1105         378 :   mt = alg_get_multable(al);
    1106         378 :   p = alg_get_char(al);
    1107         378 :   sp = signe(p);
    1108        2499 :   for(i=2; i<=N; i++)
    1109       24374 :     for(j=2; j<=N; j++)
    1110      314580 :       for(k=1; k<=N; k++) {
    1111      292383 :         a = gcoeff(gel(mt,i),k,j);
    1112      292383 :         b = gcoeff(gel(mt,j),k,i);
    1113      292383 :         if(sp) {
    1114      288253 :           if(cmpii(Fp_red(a,p), Fp_red(b,p))) return 0;
    1115             :         }
    1116        4130 :         else if(gcmp(a,b)) return 0;
    1117             :       }
    1118         322 :   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      813078 : alg_model0(GEN al, GEN x)
    1335             : {
    1336      813078 :   long t, N = alg_get_absdim(al), lx = lg(x), d, n, D, i;
    1337      813078 :   if (typ(x) == t_MAT) return al_MATRIX;
    1338      794423 :   if (typ(x) != t_COL) return al_INVALID;
    1339      794367 :   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      792029 :   switch(alg_type(al)) {
    1345             :     case al_TABLE:
    1346      638512 :       if(lx != N+1) return al_INVALID;
    1347      638491 :       return al_BASIS;
    1348             :     case al_CYCLIC:
    1349      138495 :       d = alg_get_degree(al);
    1350      138495 :       if (lx == N+1) return al_BASIS;
    1351       17388 :       if (lx == d+1) return al_ALGEBRAIC;
    1352          21 :       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      812952 : checkalgx(GEN x, long model)
    1375             : {
    1376             :   long t, i;
    1377      812952 :   switch(model) {
    1378             :     case al_BASIS:
    1379     9948183 :       for(i=1; i<lg(x); i++) {
    1380     9176426 :         t = typ(gel(x,i));
    1381     9176426 :         if (t != t_INT && t != t_FRAC)
    1382           7 :           pari_err_TYPE("checkalgx", gel(x,i));
    1383             :       }
    1384      771757 :       return;
    1385             :     case al_TRIVIAL:
    1386             :     case al_ALGEBRAIC:
    1387       90202 :       for(i=1; i<lg(x); i++) {
    1388       67676 :         t = typ(gel(x,i));
    1389       67676 :         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       22526 :       return;
    1394             :   }
    1395             : }
    1396             : 
    1397             : long
    1398      813078 : alg_model(GEN al, GEN x)
    1399             : {
    1400      813078 :   long res = alg_model0(al, x);
    1401      813078 :   if(res == al_INVALID) pari_err_TYPE("alg_model", x);
    1402      812952 :   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      122535 : _tablemul(GEN mt, GEN x, GEN y)
    1530             : {
    1531      122535 :   pari_sp av = avma;
    1532      122535 :   long D = lg(mt)-1, i;
    1533      122535 :   GEN res = NULL;
    1534     1400203 :   for (i=1; i<=D; i++) {
    1535     1277668 :     GEN c = gel(x,i);
    1536     1277668 :     if (!gequal0(c)) {
    1537      460117 :       GEN My = RgM_RgC_mul(gel(mt,i),y);
    1538      460117 :       GEN t = RgC_Rg_mul(My,c);
    1539      460117 :       res = res? RgC_add(res,t): t;
    1540             :     }
    1541             :   }
    1542      122535 :   if (!res) { avma = av; return zerocol(D); }
    1543      122528 :   return gerepileupto(av, res);
    1544             : }
    1545             : 
    1546             : static GEN
    1547      186466 : _tablemul_Fp(GEN mt, GEN x, GEN y, GEN p)
    1548             : {
    1549      186466 :   pari_sp av = avma;
    1550      186466 :   long D = lg(mt)-1, i;
    1551      186466 :   GEN res = NULL;
    1552     2457161 :   for (i=1; i<=D; i++) {
    1553     2270695 :     GEN c = gel(x,i);
    1554     2270695 :     if (signe(c)) {
    1555      366849 :       GEN My = FpM_FpC_mul(gel(mt,i),y,p);
    1556      366849 :       GEN t = FpC_Fp_mul(My,c,p);
    1557      366849 :       res = res? FpC_add(res,t,p): t;
    1558             :     }
    1559             :   }
    1560      186466 :   if (!res) { avma = av; return zerocol(D); }
    1561      185955 :   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      442939 : _tablemul_ej_Fp(GEN mt, GEN x, long j, GEN p)
    1584             : {
    1585      442939 :   pari_sp av = avma;
    1586      442939 :   long D = lg(mt)-1, i;
    1587      442939 :   GEN res = NULL;
    1588     7704053 :   for (i=1; i<=D; i++) {
    1589     7261114 :     GEN c = gel(x,i);
    1590     7261114 :     if (!gequal0(c)) {
    1591      655564 :       GEN My = gel(gel(mt,i),j);
    1592      655564 :       GEN t = FpC_Fp_mul(My,c,p);
    1593      655564 :       res = res? FpC_add(res,t,p): t;
    1594             :     }
    1595             :   }
    1596      442939 :   if (!res) { avma = av; return zerocol(D); }
    1597      442827 :   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      376495 : _tablemul_ej_Fl(GEN mt, GEN x, long j, ulong p)
    1606             : {
    1607      376495 :   pari_sp av = avma;
    1608      376495 :   long D = lg(mt)-1, i;
    1609      376495 :   GEN res = NULL;
    1610    10745616 :   for (i=1; i<=D; i++) {
    1611    10369121 :     ulong c = x[i];
    1612    10369121 :     if (c) {
    1613      727419 :       GEN My = gel(gel(mt,i),j);
    1614      727419 :       GEN t = Flv_Fl_mul(My,c, p);
    1615      727419 :       res = res? Flv_add(res,t, p): t;
    1616             :     }
    1617             :   }
    1618      376495 :   if (!res) { avma = av; return zero_Flv(D); }
    1619      376495 :   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      308581 : algbasismul(GEN al, GEN x, GEN y)
    1640             : {
    1641      308581 :   GEN mt = alg_get_multable(al), p = alg_get_char(al);
    1642      308581 :   if (signe(p)) return _tablemul_Fp(mt, x, y, p);
    1643      122115 :   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      262507 : algmul(GEN al, GEN x, GEN y)
    1685             : {
    1686      262507 :   pari_sp av = avma;
    1687             :   long tx, ty;
    1688      262507 :   checkalg(al);
    1689      262507 :   tx = alg_model(al,x);
    1690      262493 :   ty = alg_model(al,y);
    1691      262493 :   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      254072 :   if (signe(alg_get_char(al))) return algbasismul(al,x,y);
    1696      121989 :   if (tx==al_TRIVIAL) retmkcol(gmul(gel(x,1),gel(y,1)));
    1697      121744 :   if (tx==al_ALGEBRAIC && ty==al_ALGEBRAIC) return algalgmul(al,x,y);
    1698      120988 :   if (tx==al_ALGEBRAIC) x = algalgtobasis(al,x);
    1699      120988 :   if (ty==al_ALGEBRAIC) y = algalgtobasis(al,y);
    1700      120988 :   return gerepileupto(av,algbasismul(al,x,y));
    1701             : }
    1702             : 
    1703             : GEN
    1704       56924 : algsqr(GEN al, GEN x)
    1705             : {
    1706       56924 :   pari_sp av = avma;
    1707             :   long tx;
    1708       56924 :   checkalg(al);
    1709       56889 :   tx = alg_model(al,x);
    1710       56833 :   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        8078 : algalgmultable_cyc(GEN al, GEN x)
    1781             : {
    1782        8078 :   pari_sp av = avma;
    1783        8078 :   long n = alg_get_degree(al), i, j;
    1784             :   GEN res, rnf, auts, b, pol;
    1785        8078 :   rnf = alg_get_splittingfield(al);
    1786        8078 :   auts = alg_get_auts(al);
    1787        8078 :   b = alg_get_b(al);
    1788        8078 :   pol = rnf_get_pol(rnf);
    1789             : 
    1790        8078 :   res = zeromatcopy(n,n);
    1791       38269 :   for (i=0; i<n; i++)
    1792       30191 :     gcoeff(res,i+1,1) = lift_shallow(rnfbasistoalg(rnf,gel(x,i+1)));
    1793             : 
    1794       38269 :   for (i=0; i<n; i++) {
    1795       93772 :     for (j=1; j<=i; j++)
    1796       63581 :       gcoeff(res,i+1,j+1) = gmodulo(poleval(gcoeff(res,i-j+1,1),gel(auts,j)),pol);
    1797       93772 :     for (; j<n; j++)
    1798       63581 :       gcoeff(res,i+1,j+1) = gmodulo(gmul(b,poleval(gcoeff(res,n+i-j+1,1),gel(auts,j))), pol);
    1799             :   }
    1800             : 
    1801       38269 :   for (i=0; i<n; i++)
    1802       30191 :     gcoeff(res,i+1,1) = gmodulo(gcoeff(res,i+1,1),pol);
    1803             : 
    1804        8078 :   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       20426 : algbasismultable_Flm(GEN mt, GEN x, ulong m)
    1828             : {
    1829       20426 :   pari_sp av = avma;
    1830       20426 :   long D = lg(gel(mt,1))-1, i;
    1831       20426 :   GEN z = NULL;
    1832      396921 :   for (i=1; i<=D; i++)
    1833             :   {
    1834      376495 :     ulong c = x[i];
    1835      376495 :     if (c)
    1836             :     {
    1837       35098 :       GEN M = Flm_Fl_mul(gel(mt,i),c, m);
    1838       35098 :       z = z? Flm_add(z, M, m): M;
    1839             :     }
    1840             :   }
    1841       20426 :   if (!z) { avma = av; return zero_Flm(D,D); }
    1842       20426 :   return gerepileupto(av, z);
    1843             : }
    1844             : static GEN
    1845      283108 : elementabsmultable_Z(GEN mt, GEN x)
    1846             : {
    1847      283108 :   long i, l = lg(x);
    1848      283108 :   GEN z = NULL;
    1849     3271996 :   for (i = 1; i < l; i++)
    1850             :   {
    1851     2988888 :     GEN c = gel(x,i);
    1852     2988888 :     if (signe(c))
    1853             :     {
    1854     1143478 :       GEN M = ZM_Z_mul(gel(mt,i),c);
    1855     1143478 :       z = z? ZM_add(z, M): M;
    1856             :     }
    1857             :   }
    1858      283108 :   return z;
    1859             : }
    1860             : static GEN
    1861      191352 : elementabsmultable(GEN mt, GEN x)
    1862             : {
    1863      191352 :   GEN d, z = elementabsmultable_Z(mt, Q_remove_denom(x,&d));
    1864      191352 :   return (z && d)? ZM_Z_div(z, d): z;
    1865             : }
    1866             : static GEN
    1867       91756 : elementabsmultable_Fp(GEN mt, GEN x, GEN p)
    1868             : {
    1869       91756 :   GEN z = elementabsmultable_Z(mt, x);
    1870       91756 :   return z? FpM_red(z, p): z;
    1871             : }
    1872             : GEN
    1873      283108 : algbasismultable(GEN al, GEN x)
    1874             : {
    1875      283108 :   pari_sp av = avma;
    1876      283108 :   GEN z, p = alg_get_char(al), mt = alg_get_multable(al);
    1877      283108 :   z = signe(p)? elementabsmultable_Fp(mt, x, p): elementabsmultable(mt, x);
    1878      283108 :   if (!z)
    1879             :   {
    1880         476 :     long D = lg(mt)-1;
    1881         476 :     avma = av; return zeromat(D,D);
    1882             :   }
    1883      282632 :   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        8757 : algalgmultable(GEN al, GEN x)
    1895             : {
    1896        8757 :   switch(alg_type(al))
    1897             :   {
    1898        8078 :     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       66731 : algbasisrightmultable(GEN al, GEN x)
    1926             : {
    1927       66731 :   long N = alg_get_absdim(al), i,j,k;
    1928       66731 :   GEN res = zeromatcopy(N,N), c, mt = alg_get_multable(al);
    1929      600425 :   for(i=1; i<=N; i++) {
    1930      533694 :     c = gel(x,i);
    1931      533694 :     if (!gequal0(c)) {
    1932     1390865 :       for(j=1; j<=N; j++)
    1933    11155102 :       for(k=1; k<=N; k++)
    1934     9918692 :         gcoeff(res,k,j) = addii(gcoeff(res,k,j), mulii(c, gcoeff(gel(mt,j),k,i)));
    1935             :     }
    1936             :   }
    1937       66731 :   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      106400 : algleftmultable(GEN al, GEN x)
    2034             : {
    2035      106400 :   pari_sp av = avma;
    2036             :   long tx;
    2037             :   GEN res;
    2038             : 
    2039      106400 :   checkalg(al);
    2040      106400 :   tx = alg_model(al,x);
    2041      106386 :   switch(tx) {
    2042         175 :     case al_TRIVIAL : res = mkmatcopy(mkcol(gel(x,1))); break;
    2043         812 :     case al_ALGEBRAIC : res = algalgmultable(al,x); break;
    2044      105119 :     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      106379 :   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       29050 : static GEN _mul(void* data, GEN x, GEN y) { return algmul((GEN)data,x,y); }
    2207       55881 : 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       12880 : algpow(GEN al, GEN x, GEN n)
    2226             : {
    2227       12880 :   pari_sp av = avma;
    2228             :   GEN res;
    2229       12880 :   checkalg(al);
    2230       12880 :   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       12803 :     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       12873 :   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       10927 : algbasischarpoly(GEN al, GEN x, long v)
    2270             : {
    2271       10927 :   pari_sp av = avma;
    2272       10927 :   GEN p = alg_get_char(al), mx;
    2273       10927 :   if (alg_model(al,x) == al_MATRIX) mx = algleftmultable_mat(al,x);
    2274       10878 :   else                              mx = algbasismultable(al,x);
    2275       10920 :   if (signe(p)) {
    2276        8967 :     GEN res = FpM_charpoly(mx,p);
    2277        8967 :     setvarn(res,v);
    2278        8967 :     return gerepileupto(av, res);
    2279             :   }
    2280        1953 :   return gerepileupto(av, charpoly(mx,v));
    2281             : }
    2282             : 
    2283             : GEN
    2284       10948 : algcharpoly(GEN al, GEN x, long v)
    2285             : {
    2286       10948 :   checkalg(al);
    2287       10948 :   if (v<0) v=0;
    2288             : 
    2289             :   /* gneg(x[1]) left on stack */
    2290       10948 :   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       10857 :   switch(alg_type(al)) {
    2297         252 :     case al_CYCLIC: case al_CSA: return algredcharpoly(al,x,v);
    2298       10605 :     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      283164 : algabstrace(GEN al, GEN x)
    2306             : {
    2307      283164 :   pari_sp av = avma;
    2308      283164 :   GEN res = NULL, p = alg_get_char(al);
    2309      283164 :   if (signe(p)) return FpV_dotproduct(x, alg_get_tracebasis(al), p);
    2310       44667 :   switch(alg_model(al,x)) {
    2311          63 :     case al_TRIVIAL: return gcopy(gel(x,1)); break;
    2312       44604 :     case al_BASIS: res = RgV_dotproduct(x, alg_get_tracebasis(al)); break;
    2313             :   }
    2314       44604 :   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       38969 : ZM_trace(GEN x)
    2373             : {
    2374       38969 :   long i, lx = lg(x);
    2375             :   GEN t;
    2376       38969 :   if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
    2377       38087 :   t = gcoeff(x,1,1);
    2378       38087 :   for (i = 2; i < lx; i++) t = addii(t, gcoeff(x,i,i));
    2379       38087 :   return t;
    2380             : }
    2381             : static GEN
    2382       96131 : FpM_trace(GEN x, GEN p)
    2383             : {
    2384       96131 :   long i, lx = lg(x);
    2385             :   GEN t;
    2386       96131 :   if (lx < 3) return lx == 1? gen_0: gcopy(gcoeff(x,1,1));
    2387       88088 :   t = gcoeff(x,1,1);
    2388       88088 :   for (i = 2; i < lx; i++) t = Fp_add(t, gcoeff(x,i,i), p);
    2389       88088 :   return t;
    2390             : }
    2391             : 
    2392             : static GEN
    2393       31458 : algtracebasis(GEN al)
    2394             : {
    2395       31458 :   pari_sp av = avma;
    2396       31458 :   GEN mt = alg_get_multable(al), p = alg_get_char(al);
    2397       31458 :   long i, l = lg(mt);
    2398       31458 :   GEN v = cgetg(l, t_VEC);
    2399       31458 :   if (signe(p)) for (i=1; i < l; i++) gel(v,i) = FpM_trace(gel(mt,i), p);
    2400        4991 :   else          for (i=1; i < l; i++) gel(v,i) = ZM_trace(gel(mt,i));
    2401       31458 :   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       20426 : algtracei(GEN mt, ulong p, ulong expo, ulong modu)
    2408             : {
    2409       20426 :   pari_sp av = avma;
    2410       20426 :   long j, l = lg(mt);
    2411       20426 :   ulong tr = 0;
    2412       20426 :   mt = Flm_powu(mt,expo,modu);
    2413       20426 :   for (j=1; j<l; j++) tr += ucoeff(mt,j,j);
    2414       20426 :   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        6342 : algalgtonat_cyc(GEN al, GEN x)
    2450             : {
    2451        6342 :   pari_sp av = avma;
    2452        6342 :   GEN nf = alg_get_abssplitting(al), rnf = alg_get_splittingfield(al), res, c;
    2453        6342 :   long n = alg_get_degree(al), N = nf_get_degree(nf), i, i1;
    2454        6342 :   res = zerocol(N*n);
    2455       23275 :   for (i=0; i<n; i++) {
    2456       16933 :     c = gel(x,i+1);
    2457       16933 :     c = rnfeltreltoabs(rnf,c);
    2458       16933 :     if (!gequal0(c)) {
    2459        8855 :       c = algtobasis(nf,c);
    2460        8855 :       for (i1=1; i1<=N; i1++) gel(res,i*N+i1) = gel(c,i1);
    2461             :     }
    2462             :   }
    2463        6342 :   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        7182 : algalgtonat(GEN al, GEN x)
    2486             : {
    2487        7182 :   switch(alg_type(al))
    2488             :   {
    2489        6342 :     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        7217 : algalgtobasis(GEN al, GEN x)
    2555             : {
    2556             :   pari_sp av;
    2557             :   long tx;
    2558        7217 :   checkalg(al);
    2559        7217 :   if (alg_type(al) == al_TABLE) pari_err_TYPE("algalgtobasis [use alginit]", al);
    2560        7203 :   tx = alg_model(al,x);
    2561        7203 :   if (tx==al_BASIS) return gcopy(x);
    2562        7196 :   if (tx==al_MATRIX) return algalgtobasis_mat(al,x);
    2563        7182 :   av = avma;
    2564        7182 :   x = algalgtonat(al,x);
    2565        7182 :   x = RgM_RgC_mul(alg_get_invbasis(al),x);
    2566        7182 :   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       21315 : algrandom(GEN al, GEN b)
    2599             : {
    2600             :   GEN res, p, N;
    2601             :   long i, n;
    2602       21315 :   if (typ(b) != t_INT) pari_err_TYPE("algrandom",b);
    2603       21308 :   if (signe(b)<0) pari_err_DOMAIN("algrandom", "b", "<", gen_0, b);
    2604       21301 :   checkalg(al);
    2605       21294 :   n = alg_get_absdim(al);
    2606       21294 :   N = addiu(shifti(b,1), 1); /* left on stack */
    2607       21294 :   p = alg_get_char(al);
    2608       21294 :   res = cgetg(n+1,t_COL);
    2609      190638 :   for (i=1; i<= n; i++)
    2610             :   {
    2611      169344 :     pari_sp av = avma;
    2612      169344 :     gel(res,i) = gerepileuptoint(av, subii(randomi(N),b));
    2613             :   }
    2614       21294 :   if (signe(p)) res = FpC_red(res, p); /*FIXME: need garbage collection here?*/
    2615       21294 :   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        8589 : algpoleval(GEN al, GEN pol, GEN x)
    2623             : {
    2624        8589 :   pari_sp av = avma;
    2625             :   GEN p, mx, res;
    2626             :   long i;
    2627        8589 :   checkalg(al);
    2628        8589 :   p = alg_get_char(al);
    2629        8589 :   if (typ(pol) != t_POL) pari_err_TYPE("algpoleval",pol);
    2630        8582 :   mx = (typ(x) == t_MAT)? x: algleftmultable(al,x);
    2631        8582 :   res = zerocol(lg(mx)-1);
    2632        8582 :   if (signe(p)) {
    2633       34874 :     for (i=lg(pol)-1; i>1; i--)
    2634             :     {
    2635       27251 :       gel(res,1) = Fp_add(gel(res,1), gel(pol,i), p);
    2636       27251 :       if (i>2) res = FpM_FpC_mul(mx, res, p);
    2637             :     }
    2638             :   }
    2639             :   else {
    2640        6013 :     for (i=lg(pol)-1; i>1; i--)
    2641             :     {
    2642        5054 :       gel(res,1) = gadd(gel(res,1), gel(pol,i));
    2643        5054 :       if (i>2) res = RgM_RgC_mul(mx, res);
    2644             :     }
    2645             :   }
    2646        8582 :   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         238 :   for (i = 1;i < l; i++)
    3084         175 :     if (signe(gel(E,i))) {
    3085          70 :       gel(P2,j) = gel(P,i);
    3086          70 :       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         385 :     x1 = NULL;
    3114         833 :     for (j=1; j<l; j++)
    3115             :     {
    3116         448 :       GEN pr = gel(listpr,j), e;
    3117         448 :       long v = nfval(nf, x, pr);
    3118         448 :       e = stoi(v); gel(L,j) = e;
    3119         581 :       if (v) x1 = x1? idealmulpowprime(nf, x1, pr, e)
    3120         133 :                     : idealpow(nf, pr, e);
    3121             :     }
    3122         385 :     if (x1) x = idealdivexact(nf, idealhnf(nf,x), x1);
    3123             :   }
    3124             :   else
    3125             :   { /* HNF, divide as we proceed (reduce size) */
    3126         329 :     for (j=1; j<l; j++)
    3127             :     {
    3128         189 :       GEN pr = gel(listpr,j);
    3129         189 :       long v = idealval(nf, x, pr);
    3130         189 :       gel(L,j) = stoi(v);
    3131         189 :       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         413 : pr_primes(GEN v)
    3282             : {
    3283         413 :   long i, l = lg(v);
    3284         413 :   GEN w = cgetg(l,t_VEC);
    3285         413 :   for (i=1; i<l; i++) gel(w,i) = pr_get_p(gel(v,i));
    3286         413 :   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         497 : alg_maximal_primes(GEN al, GEN P)
    3896             : {
    3897         497 :   pari_sp av = avma;
    3898         497 :   long l = lg(P), i;
    3899        1274 :   for (i=1; i<l; i++)
    3900             :   {
    3901         777 :     if (i != 1) al = gerepilecopy(av, al);
    3902         777 :     al = alg_pmaximal(al,gel(P,i));
    3903             :   }
    3904         497 :   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         364 :     GEN hf = alg_get_hasse_f(al), pr = gel(hf,1);
    3938         364 :     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       29995 : algtableinit_i(GEN mt0, GEN p)
    4089             : {
    4090             :   GEN al, mt;
    4091             :   long i, n;
    4092             : 
    4093       29995 :   if (p && typ(p) != t_INT) pari_err_TYPE("algtableinit",p);
    4094       29988 :   if (p && !signe(p)) p = NULL;
    4095       29988 :   mt = check_mt(mt0,p);
    4096       29988 :   if (!mt) pari_err_TYPE("algtableinit", mt0);
    4097       29988 :   if (!p && !isint1(Q_denom(mt0)))
    4098           7 :     pari_err_DOMAIN("algtableinit", "denominator(mt)", "!=", gen_1, mt0);
    4099       29981 :   n = lg(mt)-1;
    4100       29981 :   al = cgetg(12, t_VEC);
    4101       29981 :   for (i=1; i<=6; i++) gel(al,i) = gen_0;
    4102       29981 :   gel(al,7) = matid(n);
    4103       29981 :   gel(al,8) = matid(n);
    4104       29981 :   gel(al,9) = mt;
    4105       29981 :   gel(al,10) = p? p: gen_0;
    4106       29981 :   gel(al,11)= algtracebasis(al);
    4107       29981 :   return al;
    4108             : }
    4109             : GEN
    4110         224 : algtableinit(GEN mt0, GEN p)
    4111             : {
    4112         224 :   pari_sp av = avma;
    4113         224 :   return gerepilecopy(av, algtableinit_i(mt0, p));
    4114             : }
    4115             : 
    4116             : /** REPRESENTATIONS OF GROUPS **/
    4117             : 
    4118             : static GEN
    4119         378 : list_to_regular_rep(GEN elts, long n)
    4120             : {
    4121             :   GEN reg, elts2, g;
    4122             :   long i,j;
    4123         378 :   elts = shallowcopy(elts);
    4124         378 :   gen_sort_inplace(elts, (void*)&vecsmall_lexcmp, &cmp_nodata, NULL);
    4125         378 :   reg = cgetg(n+1, t_VEC);
    4126         378 :   gel(reg,1) = identity_perm(n);
    4127        6797 :   for(i=2; i<=n; i++) {
    4128        6419 :     g = perm_inv(gel(elts,i));
    4129        6419 :     elts2 = cgetg(n+1, t_VEC);
    4130        6419 :     for(j=1; j<=n; j++) gel(elts2,j) = perm_mul(g,gel(elts,j));
    4131        6419 :     gen_sort_inplace(elts2, (void*)&vecsmall_lexcmp, &cmp_nodata, &gel(reg,i));
    4132             :   }
    4133         378 :   return reg;
    4134             : }
    4135             : 
    4136             : static GEN
    4137        6797 : matrix_perm(GEN perm, long n)
    4138             : {
    4139             :   GEN m;
    4140             :   long j;
    4141        6797 :   m = cgetg(n+1, t_MAT);
    4142      186858 :   for(j=1; j<=n; j++) {
    4143      180061 :     gel(m,j) = col_ei(n,perm[j]);
    4144             :   }
    4145        6797 :   return m;
    4146             : }
    4147             : 
    4148             : static GEN
    4149         903 : groupelts_algcenter(GEN elts, GEN p, GEN* ptr_conjclasses)
    4150             : {
    4151         903 :   pari_sp av = avma;
    4152             :   long nbcl, i, n, k, j, ci, cj, ck;
    4153             :   GEN conjclass, mt, xi, xj, xixj, repclass, al;
    4154         903 :   if(p && !signe(p)) p = NULL;
    4155         903 :   n = lg(elts)-1;
    4156         903 :   elts = gen_sort(elts,(void*)vecsmall_lexcmp,cmp_nodata);
    4157             :   /* compute conjugacy classes */
    4158         903 :   conjclass = groupelts_conjclasses(elts,&nbcl);
    4159         903 :   repclass = const_vecsmall(nbcl,0);
    4160       18907 :   for(i=1;i<=n;i++)
    4161       18004 :     if(!repclass[conjclass[i]]) repclass[conjclass[i]] = i;
    4162             : 
    4163             :   /* TODO sort conjugacy classes using a reasonable cmp */
    4164             : 
    4165             :   /* compute multiplication table of the center of the group algebra
    4166             :    * (class functions). */
    4167         903 :   mt = cgetg(nbcl+1,t_VEC);
    4168        8344 :   for(i=1;i<=nbcl;i++)
    4169        7441 :     gel(mt,i) = zeromatcopy(nbcl,nbcl);
    4170       18907 :   for(i=1;i<=n;i++) {
    4171       18004 :     xi = gel(elts,i);
    4172       18004 :     ci = conjclass[i];
    4173      803376 :     for(j=1;j<=n;j++) {
    4174      785372 :       xj = gel(elts,j);
    4175      785372 :       cj = conjclass[j];
    4176      785372 :       xixj = perm_mul(xi,xj);
    4177      785372 :       k = vecsearch(elts,xixj,NULL);
    4178      785372 :       ck = conjclass[k];
    4179      785372 :       if(repclass[ck]==k)
    4180      190946 :         gcoeff(gel(mt,ci),ck,cj) = addiu(gcoeff(gel(mt,ci),ck,cj), 1);
    4181             :     }
    4182             :   }
    4183         903 :   if(p)
    4184        7322 :     for(i=1; i<=nbcl; i++)
    4185        6531 :       gel(mt,i) = FpM_red(gel(mt,i),p);
    4186             : 
    4187         903 :   al = algtableinit_i(mt,p);
    4188         903 :   if(ptr_conjclasses) {
    4189         511 :     *ptr_conjclasses = mkvec3(elts,conjclass,repclass);
    4190         511 :     gerepileall(av,2,&al,ptr_conjclasses);
    4191             :   }
    4192         392 :   else al = gerepilecopy(av,al);
    4193         903 :   return al;
    4194             : }
    4195             : 
    4196             : GEN
    4197         917 : alggroupcenter(GEN gal, GEN p, GEN* ptr_conjclasses)
    4198             : {
    4199         917 :   GEN elts = checkgroupelts(gal);
    4200         903 :   return groupelts_algcenter(elts, p, ptr_conjclasses);
    4201             : }
    4202             : 
    4203             : static GEN
    4204         378 : groupelts_algebra(GEN elts, GEN p)
    4205             : {
    4206         378 :   pari_sp av = avma;
    4207             :   GEN mt;
    4208         378 :   long i, n = lg(elts)-1;
    4209         378 :   elts = list_to_regular_rep(elts,n);
    4210         378 :   mt = cgetg(n+1, t_VEC);
    4211         378 :   for(i=1; i<=n; i++) gel(mt,i) = matrix_perm(gel(elts,i),n);
    4212         378 :   return gerepilecopy(av, algtableinit_i(mt,p));
    4213             : }
    4214             : 
    4215             : GEN
    4216         413 : alggroup(GEN gal, GEN p)
    4217             : {
    4218         413 :   GEN elts = checkgroupelts(gal);
    4219         378 :   return groupelts_algebra(elts, p);
    4220             : }
    4221             : 
    4222             : /** MAXIMAL ORDER **/
    4223             : 
    4224             : static GEN
    4225       32550 : mattocol(GEN M, long n)
    4226             : {
    4227       32550 :   GEN C = cgetg(n*n+1, t_COL);
    4228             :   long i,j,ic;
    4229       32550 :   ic = 1;
    4230      873474 :   for (i=1; i<=n; i++)
    4231      840924 :   for (j=1; j<=n; j++, ic++) gel(C,ic) = gcoeff(M,i,j);
    4232       32550 :   return C;
    4233             : }
    4234             : 
    4235             : /*Ip is a lift of a left O/pO-ideal where O is the integral basis of al*/
    4236             : GEN
    4237        2107 : algleftordermodp(GEN al, GEN Ip, GEN p)
    4238             : {
    4239        2107 :   pari_sp av = avma;
    4240             :   GEN I, Ii, M, mt, K, imi, p2;
    4241             :   long n, i;
    4242        2107 :   n = alg_get_absdim(al);
    4243        2107 :   mt = alg_get_multable(al);
    4244        2107 :   p2 = sqri(p);
    4245             : 
    4246        2107 :   I = ZM_hnfmodid(Ip, p);
    4247        2107 :   Ii = ZM_inv(I,p);
    4248             : 
    4249        2107 :   M = cgetg(n+1, t_MAT);
    4250       34657 :   for (i=1; i<=n; i++) {
    4251       32550 :     imi = FpM_mul(Ii, FpM_mul(gel(mt,i), I, p2), p2);
    4252       32550 :     imi = ZM_Z_divexact(imi, p);
    4253       32550 :     gel(M,i) = mattocol(imi, n);
    4254             :   }
    4255             : 
    4256             :   /*TODO : FpM_invimage superbad documentation (have to read RgM_invimage) Does it really do what it claims if left matrix is not invertible ?*/
    4257        2107 :   K = FpM_ker(M, p);
    4258        2107 :   if (lg(K)==1) { avma = av; return matid(n); }
    4259         924 :   K = ZM_hnfmodid(K,p);
    4260             : 
    4261         924 :   return gerepileupto(av, ZM_Z_div(K,p));
    4262             : }
    4263             : 
    4264             : GEN
    4265        3003 : alg_ordermodp(GEN al, GEN p)
    4266             : {
    4267             :   GEN alp;
    4268        3003 :   long i, N = alg_get_absdim(al);
    4269        3003 :   alp = cgetg(12, t_VEC);
    4270        3003 :   for (i=1; i<=8; i++) gel(alp,i) = gen_0;
    4271        3003 :   gel(alp,9) = cgetg(N+1, t_VEC);
    4272        3003 :   for(i=1; i<=N; i++) gmael(alp,9,i) = FpM_red(gmael(al,9,i), p);
    4273        3003 :   gel(alp,10) = p;
    4274        3003 :   gel(alp,11) = cgetg(N+1, t_VEC);
    4275        3003 :   for(i=1; i<=N; i++) gmael(alp,11,i) = Fp_red(gmael(al,11,i), p);
    4276             : 
    4277        3003 :   return alp;
    4278             : }
    4279             : 
    4280             : static GEN
    4281        1701 : algpradical_i(GEN al, GEN p, GEN zprad, GEN projs)
    4282             : {
    4283        1701 :   pari_sp av = avma;
    4284        1701 :   GEN alp = alg_ordermodp(al, p), liftrad, projrad, alq, alrad, res, Lalp, radq;
    4285             :   long i;
    4286        1701 :   if(lg(zprad)==1) {
    4287        1407 :     liftrad = NULL;
    4288        1407 :     projrad = NULL;
    4289             :   }
    4290             :   else {
    4291         294 :     alq = alg_quotient(alp, zprad, 1);
    4292         294 :     alp = gel(alq,1);
    4293         294 :     projrad = gel(alq,2);
    4294         294 :     liftrad = gel(alq,3);
    4295             :   }
    4296             : 
    4297        1701 :   if (projs) {
    4298         343 :     if (projrad) {
    4299          21 :       projs = gcopy(projs);
    4300          63 :       for(i=1; i<lg(projs); i++)
    4301          42 :         gel(projs,i) = FpM_FpC_mul(projrad, gel(projs,i), p);
    4302             :     }
    4303         343 :     Lalp = alg_centralproj(alp,projs,1);
    4304             : 
    4305         343 :     alrad = cgetg(lg(Lalp),t_VEC);
    4306        1106 :     for(i=1; i<lg(Lalp); i++) {
    4307         763 :       alq = gel(Lalp,i);
    4308         763 :       radq = algradical(gel(alq,1));
    4309         763 :       if(gequal0(radq))
    4310         245 :         gel(alrad,i) = cgetg(1,t_MAT);
    4311             :       else {
    4312         518 :         radq = FpM_mul(gel(alq,3),radq,p);
    4313         518 :         gel(alrad,i) = radq;
    4314             :       }
    4315             :     }
    4316         343 :     alrad = shallowmatconcat(alrad);
    4317         343 :     alrad = FpM_image(alrad,p);
    4318             :   }
    4319        1358 :   else alrad = algradical(alp);
    4320             : 
    4321        1701 :   if(!gequal0(alrad)) {
    4322        1330 :     if (liftrad) alrad = FpM_mul(liftrad, alrad, p);
    4323        1330 :     res = shallowmatconcat(mkvec2(alrad, zprad));
    4324        1330 :     res = FpM_image(res,p);
    4325             :   }
    4326         371 :   else res = lg(zprad)==1 ? gen_0 : zprad;
    4327        1701 :   return gerepilecopy(av, res);
    4328             : }
    4329             : #if 0
    4330             : /* not used */
    4331             : GEN
    4332             : algpradical(GEN al, GEN p)
    4333             : {
    4334             :   GEN placeholder = cgetg(1,t_MAT); /*left on stack*/
    4335             :   return algpradical_i(al, p, placeholder, NULL);
    4336             : }
    4337             : #endif
    4338             : 
    4339             : static GEN
    4340        1302 : algpdecompose0(GEN al, GEN prad, GEN p, GEN projs)
    4341             : {
    4342        1302 :   pari_sp av = avma;
    4343        1302 :   GEN alp, quo, ss, liftm = NULL, projm = NULL, dec, res, I, Lss, deci;
    4344             :   long i, j;
    4345             : 
    4346        1302 :   alp = alg_ordermodp(al, p);
    4347        1302 :   if (!gequal0(prad)) {
    4348        1029 :     quo = alg_quotient(alp,prad,1);
    4349        1029 :     ss = gel(quo,1);
    4350        1029 :     projm = gel(quo,2);
    4351        1029 :     liftm = gel(quo,3);
    4352             :   }
    4353         273 :   else ss = alp;
    4354             : 
    4355        1302 :   if (projs) {
    4356         294 :     if (projm) {
    4357         833 :       for(i=1; i<lg(projs); i++)
    4358         574 :         gel(projs,i) = FpM_FpC_mul(projm, gel(projs,i), p);
    4359             :     }
    4360         294 :     Lss = alg_centralproj(ss, projs, 1);
    4361             : 
    4362         294 :     dec = cgetg(lg(Lss),t_VEC);
    4363         952 :     for(i=1; i<lg(Lss); i++) {
    4364         658 :       gel(dec,i) = algsimpledec(gmael(Lss,i,1), 1);
    4365         658 :       deci = gel(dec,i);
    4366        1736 :       for(j=1; j<lg(deci); j++)
    4367        1078 :        gmael(deci,j,3) = FpM_mul(gmael(Lss,i,3), gmael(deci,j,3), p);
    4368             :     }
    4369         294 :     dec = shallowconcat1(dec);
    4370             :   }
    4371        1008 :   else dec = algsimpledec(ss,1);
    4372             : 
    4373        1302 :   res = cgetg(lg(dec),t_VEC);
    4374        3864 :   for (i=1; i<lg(dec); i++) {
    4375        2562 :     I = gmael(dec,i,3);
    4376        2562 :     if (liftm) I = FpM_mul(liftm,I,p);
    4377        2562 :     I = shallowmatconcat(mkvec2(I,prad));
    4378        2562 :     gel(res,i) = I;
    4379             :   }
    4380             : 
    4381        1302 :   return gerepilecopy(av, res);
    4382             : }
    4383             : 
    4384             : /*finds a nontrivial ideal of O/prad or gen_0 if there is none.*/
    4385             : static GEN
    4386         525 : algpdecompose_i(GEN al, GEN p, GEN zprad, GEN projs)
    4387             : {
    4388         525 :   pari_sp av = avma;
    4389         525 :   GEN prad = algpradical_i(al,p,zprad,projs);
    4390         525 :   return gerepileupto(av, algpdecompose0(al, prad, p, projs));
    4391             : }
    4392             : #if 0
    4393             : /* not used */
    4394             : GEN
    4395             : algpdecompose(GEN al, GEN p)
    4396             : {
    4397             :   GEN placeholder = cgetg(1,t_MAT); /*left on stack*/
    4398             :   return algpdecompose_i(al, p, placeholder, NULL);
    4399             : }
    4400             : #endif
    4401             : 
    4402             : /* ord is assumed to be in hnf wrt the integral basis of al. */
    4403             : /* assumes that alg_get_invbasis(al) is integral. */
    4404             : GEN
    4405         924 : alg_change_overorder_shallow(GEN al, GEN ord)
    4406             : {
    4407             :   GEN al2, mt, iord, mtx, den, den2, div;
    4408             :   long i, n;
    4409         924 :   n = alg_get_absdim(al);
    4410             : 
    4411         924 :   iord = QM_inv(ord,gen_1);
    4412         924 :   al2 = shallowcopy(al);
    4413         924 :   ord = Q_remove_denom(ord,&den);
    4414             : 
    4415         924 :   gel(al2,7) = Q_remove_denom(gel(al,7), &den2);
    4416         924 :   if (den2) div = mulii(den,den2);
    4417         399 :   else      div = den;
    4418         924 :   gel(al2,7) = ZM_Z_div(ZM_mul(gel(al2,7), ord), div);
    4419             : 
    4420         924 :   gel(al2,8) = ZM_mul(iord, gel(al,8));
    4421             : 
    4422         924 :   mt = cgetg(n+1,t_VEC);
    4423         924 :   gel(mt,1) = matid(n);
    4424         924 :   div = sqri(den);
    4425       13482 :   for(i=2; i<=n; i++) {
    4426       12558 :     mtx = algbasismultable(al,gel(ord,i));
    4427       12558 :     gel(mt,i) = ZM_mul(iord, ZM_mul(mtx, ord));
    4428       12558 :     gel(mt,i) = ZM_Z_divexact(gel(mt,i), div);
    4429             :   }
    4430         924 :   gel(al2,9) = mt;
    4431             : 
    4432         924 :   gel(al2,11) = algtracebasis(al2);
    4433             : 
    4434         924 :   return al2;
    4435             : }
    4436             : 
    4437             : #if 0
    4438             : /* not used */
    4439             : /*ord is assumed to be in hnf wrt the integral basis of al.*/
    4440             : GEN
    4441             : alg_changeorder_shallow(GEN al, GEN ord)
    4442             : {
    4443             :   GEN al2, mt, iord, mtx;
    4444             :   long i, n;
    4445             :   n = alg_get_absdim(al);
    4446             : 
    4447             :   iord = RgM_inv_upper(ord);
    4448             :   al2 = shallowcopy(al);
    4449             :   gel(al2,7) = RgM_mul(gel(al,7), ord);
    4450             :   gel(al2,8) = RgM_mul(iord, gel(al,8));
    4451             : 
    4452             :   mt = cgetg(n+1,t_VEC);
    4453             :   gel(mt,1) = matid(n);
    4454             :   for(i=2; i<=n; i++) {
    4455             :     mtx = algbasismultable(al,gel(ord,i));
    4456             :     gel(mt,i) = RgM_mul(iord, RgM_mul(mtx, ord));
    4457             :   }
    4458             :   gel(al2,9) = mt;
    4459             :   gel(al2,11)= algtracebasis(al2);
    4460             : 
    4461             :   return al2;
    4462             : }
    4463             : 
    4464             : GEN
    4465             : alg_changeorder(GEN al, GEN ord)
    4466             : {
    4467             :   pari_sp av = avma;
    4468             :   GEN res = alg_changeorder_shallow(al, ord);
    4469             :   return gerepilecopy(av, res);
    4470             : }
    4471             : #endif
    4472             : 
    4473             : static GEN
    4474        5243 : algfromcenter(GEN al, GEN x)
    4475             : {
    4476        5243 :   GEN nf = alg_get_center(al);
    4477             :   long n;
    4478        5243 :   switch(alg_type(al)) {
    4479             :     case al_CYCLIC:
    4480        4641 :       n = alg_get_degree(al);
    4481        4641 :       break;
    4482             :     case al_CSA:
    4483         602 :       n = alg_get_dim(al);
    4484         602 :       break;
    4485             :     default:
    4486             :       return NULL; /*LCOV_EXCL_LINE*/
    4487             :   }
    4488        5243 :   return algalgtobasis(al, scalarcol(basistoalg(nf, x), n));
    4489             : }
    4490             : 
    4491             : /* x is an ideal of the center in hnf form */
    4492             : static GEN
    4493        1701 : algfromcenterhnf(GEN al, GEN x)
    4494             : {
    4495             :   GEN res;
    4496             :   long i;
    4497        1701 :   res = cgetg(lg(x), t_MAT);
    4498        1701 :   for(i=1; i<lg(x); i++) gel(res,i) = algfromcenter(al, gel(x,i));
    4499        1701 :   return res;
    4500             : }
    4501             : 
    4502             : /* assumes al is CSA or CYCLIC */
    4503             : static GEN
    4504         777 : algcenter_precompute(GEN al, GEN p)
    4505             : {
    4506         777 :   GEN fa, pdec, nfprad, projs, nf = alg_get_center(al);
    4507             :   long i, np;
    4508             : 
    4509         777 :   pdec = idealprimedec(nf, p);
    4510         777 :   settyp(pdec, t_COL);
    4511         777 :   np = lg(pdec)-1;
    4512         777 :   fa = mkmat2(pdec, const_col(np, gen_1));
    4513         777 :   if (dvdii(nf_get_disc(nf), p))
    4514         126 :     nfprad = idealprodprime(nf, pdec);
    4515             :   else
    4516         651 :     nfprad = scalarmat_shallow(p, nf_get_degree(nf));
    4517         777 :   fa = idealchineseinit(nf, fa);
    4518         777 :   projs = cgetg(np+1, t_VEC);
    4519         777 :   for(i=1; i<=np; i++) gel(projs, i) = idealchinese(nf, fa, vec_ei(np,i));
    4520         777 :   return mkvec2(nfprad, projs);
    4521             : }
    4522             : 
    4523             : static GEN
    4524        1701 : algcenter_prad(GEN al, GEN p, GEN pre)
    4525             : {
    4526             :   GEN nfprad, zprad, mtprad;
    4527             :   long i;
    4528        1701 :   nfprad = gel(pre,1);
    4529        1701 :   zprad = algfromcenterhnf(al, nfprad);
    4530        1701 :   zprad = FpM_image(zprad, p);
    4531        1701 :   mtprad = cgetg(lg(zprad), t_VEC);
    4532        1701 :   for(i=1; i<lg(zprad); i++) gel(mtprad, i) = algbasismultable(al, gel(zprad,i));
    4533        1701 :   mtprad = shallowmatconcat(mtprad);
    4534        1701 :   zprad = FpM_image(mtprad, p);
    4535        1701 :   return zprad;
    4536             : }
    4537             : 
    4538             : static GEN
    4539        1701 : algcenter_p_projs(GEN al, GEN p, GEN pre)
    4540             : {
    4541             :   GEN projs, zprojs;
    4542             :   long i;
    4543        1701 :   projs = gel(pre,2);
    4544        1701 :   zprojs = cgetg(lg(projs), t_VEC);
    4545        1701 :   for(i=1; i<lg(projs); i++) gel(zprojs,i) = FpC_red(algfromcenter(al, gel(projs,i)),p);
    4546        1701 :   return zprojs;
    4547             : }
    4548             : 
    4549             : /*al is assumed to be simple*/
    4550             : static GEN
    4551         777 : alg_pmaximal_i(GEN al, GEN p)
    4552             : {
    4553         777 :   GEN al2, prad, lord = gen_0, I, id, dec, zprad, projs, pre;
    4554             :   long n, i;
    4555         777 :   n = alg_get_absdim(al);
    4556         777 :   id = matid(n);
    4557         777 :   al2 = al;
    4558             : 
    4559         777 :   dbg_printf(0)("Round 2 (non-commutative) at p=%Ps, dim=%d\n", p, n);
    4560             : 
    4561         777 :   pre = algcenter_precompute(al,p);
    4562             : 
    4563             :   while (1) {
    4564        1176 :     zprad = algcenter_prad(al2, p, pre);
    4565        1176 :     projs = algcenter_p_projs(al2, p, pre);
    4566        1176 :     if (lg(projs) == 2) projs = NULL;
    4567        1176 :     prad = algpradical_i(al2,p,zprad,projs);
    4568        1176 :     if (typ(prad) == t_INT) break;
    4569        1155 :     lord = algleftordermodp(al2,prad,p);
    4570        1155 :     if (!cmp_universal(lord,id)) break;
    4571         399 :     al2 = alg_change_overorder_shallow(al2,lord);
    4572         399 :   }
    4573             : 
    4574         777 :   dec = algpdecompose0(al2,prad,p,projs);
    4575        2079 :   while (lg(dec)>2) {
    4576        1092 :     for (i=1; i<lg(dec); i++) {
    4577         952 :       I = gel(dec,i);
    4578         952 :       lord = algleftordermodp(al2,I,p);
    4579         952 :       if (cmp_universal(lord,matid(n))) break;
    4580             :     }
    4581         665 :     if (i==lg(dec)) break;
    4582         525 :     al2 = alg_change_overorder_shallow(al2,lord);
    4583         525 :     zprad = algcenter_prad(al2, p, pre);
    4584         525 :     projs = algcenter_p_projs(al2, p, pre);
    4585         525 :     if (lg(projs) == 2) projs = NULL;
    4586         525 :     dec = algpdecompose_i(al2,p,zprad,projs);
    4587             :   }
    4588         777 :   return al2;
    4589             : }
    4590             : static GEN
    4591         777 : alg_pmaximal(GEN al, GEN p)
    4592             : {
    4593         777 :   pari_sp av = avma;
    4594         777 :   return gerepilecopy(av, alg_pmaximal_i(al, p));
    4595             : }
    4596             : 
    4597             : static GEN
    4598        3094 : algtracematrix(GEN al)
    4599             : {
    4600             :   GEN M, mt;
    4601             :   long n, i, j;
    4602        3094 :   n = alg_get_absdim(al);
    4603        3094 :   mt = alg_get_multable(al);
    4604        3094 :   M = cgetg(n+1, t_MAT);
    4605       33229 :   for (i=1; i<=n; i++)
    4606             :   {
    4607       30135 :     gel(M,i) = cgetg(n+1,t_MAT);
    4608      313166 :     for (j=1; j<=i; j++)
    4609      283031 :       gcoeff(M,j,i) = gcoeff(M,i,j) = algabstrace(al,gmael(mt,i,j));
    4610             :   }
    4611        3094 :   return M;
    4612             : }
    4613             : GEN
    4614          91 : algdisc(GEN al)
    4615             : {
    4616          91 :   pari_sp av = avma;
    4617          91 :   checkalg(al);
    4618          91 :   return gerepileuptoint(av, ZM_det(algtracematrix(al)));
    4619             : }
    4620             : static GEN
    4621          84 : alg_maximal(GEN al)
    4622             : {
    4623          84 :   pari_sp av = avma;
    4624          84 :   GEN fa = absZ_factor(algdisc(al));
    4625          84 :   return gerepilecopy(av, alg_maximal_primes(al, gel(fa,1)));
    4626             : }
    4627             : 
    4628             : /** LATTICES **/
    4629             : 
    4630             : /*
    4631             :  Convention: lattice = [I,t] representing t*I, where
    4632             :  - I integral hnf over the integral basis of the algebra, and
    4633             :  - t>0 either an integer or a rational number.
    4634             : */
    4635             : 
    4636             : /* TODO use hnfmodid whenever possible using a*O <= I <= O
    4637             :  * for instance a = ZM_det_triangular(I) */
    4638             : 
    4639             : static GEN
    4640      126259 : primlat(GEN lat)
    4641             : {
    4642             :   GEN m, t, c;
    4643      126259 :   m = alglat_get_primbasis(lat);
    4644      126259 :   t = alglat_get_scalar(lat);
    4645      126259 :   m = Q_primitive_part(m,&c);
    4646      126259 :   if (c) return mkvec2(m,gmul(t,c));
    4647      109319 :   return lat;
    4648             : }
    4649             : 
    4650             : /* assumes the lattice contains d * integral basis, d=0 allowed */
    4651             : GEN
    4652      101773 : alglathnf(GEN al, GEN m, GEN d)
    4653             : {
    4654      101773 :   pari_sp av = avma;
    4655             :   long N,i,j;
    4656             :   GEN m2, c;
    4657      101773 :   checkalg(al);
    4658      101773 :   N = alg_get_absdim(al);
    4659      101773 :   if (!d) d = gen_0;
    4660      101773 :   if (typ(m) == t_VEC) m = matconcat(m);
    4661      101773 :   if (typ(m) == t_COL) m = algleftmultable(al,m);
    4662      101766 :   if (typ(m) != t_MAT) pari_err_TYPE("alglathnf",m);
    4663      101759 :   if (typ(d) != t_FRAC && typ(d) != t_INT) pari_err_TYPE("alglathnf",d);
    4664      101759 :   if (lg(m)-1 < N || lg(gel(m,1))-1 != N) pari_err_DIM("alglathnf");
    4665      915614 :   for (i=1; i<=N; i++)
    4666    13604542 :     for (j=1; j<lg(m); j++)
    4667    12790666 :       if (typ(gcoeff(m,i,j)) != t_FRAC && typ(gcoeff(m,i,j)) != t_INT)
    4668           7 :         pari_err_TYPE("alglathnf", gcoeff(m,i,j));
    4669      101731 :   m2 = Q_primitive_part(m,&c);
    4670      101731 :   if (!c) c = gen_1;
    4671      101731 :   if (!signe(d)) d = detint(m2);
    4672           2 :   else           d = gdiv(d,c); /* should be an integer */
    4673      101731 :   if (!signe(d)) pari_err_INV("alglathnf [m does not have full rank]", m2);
    4674      101717 :   m2 = ZM_hnfmodid(m2,d);
    4675      101717 :   return gerepilecopy(av, mkvec2(m2,c));
    4676             : }
    4677             : 
    4678             : static GEN
    4679       21133 : prepare_multipliers(GEN *a, GEN *b)
    4680             : {
    4681             :   GEN na, nb, da, db, d;
    4682       21133 :   na = numer(*a);
    4683       21133 :   da = denom(*a);
    4684       21133 :   nb = numer(*b);
    4685       21133 :   db = denom(*b);
    4686       21133 :   na = mulii(na,db);
    4687       21133 :   nb = mulii(nb,da);
    4688       21133 :   d = gcdii(na,nb);
    4689       21133 :   *a = diviiexact(na,d);
    4690       21133 :   *b = diviiexact(nb,d);
    4691       21133 :   return gdiv(d, mulii(da,db));
    4692             : }
    4693             : 
    4694             : static GEN
    4695       21133 : prepare_lat(GEN m1, GEN t1, GEN m2, GEN t2)
    4696             : {
    4697             :   GEN d;
    4698       21133 :   d = prepare_multipliers(&t1, &t2);
    4699       21133 :   m1 = ZM_Z_mul(m1,t1);
    4700       21133 :   m2 = ZM_Z_mul(m2,t2);
    4701       21133 :   return mkvec3(m1,m2,d);
    4702             : }
    4703             : 
    4704             : static GEN
    4705       21133 : alglataddinter(GEN al, GEN lat1, GEN lat2, GEN *sum, GEN *inter)
    4706             : {
    4707             :   GEN d, m1, m2, t1, t2, M, U, H, prep;
    4708       21133 :   checkalg(al);
    4709       21133 :   checklat(al,lat1);
    4710       21133 :   checklat(al,lat2);
    4711             : 
    4712       21133 :   m1 = alglat_get_primbasis(lat1);
    4713       21133 :   t1 = alglat_get_scalar(lat1);
    4714       21133 :   m2 = alglat_get_primbasis(lat2);
    4715       21133 :   t2 = alglat_get_scalar(lat2);
    4716       21133 :   prep = prepare_lat(m1, t1, m2, t2);
    4717       21133 :   m1 = gel(prep,1);
    4718       21133 :   m2 = gel(prep,2);
    4719       21133 :   d = gel(prep,3);
    4720       21133 :   M = matconcat(mkvec2(m1,m2));
    4721       21133 :   if (!inter) U = NULL;
    4722       21133 :   H = ZM_hnfall(M,&U,1);
    4723       21133 :   if (sum) *sum = H;
    4724       21133 :   if (inter)
    4725             :   {
    4726       14084 :     U = vecslice(U,1,lg(M)-lg(H));
    4727       14084 :     *inter = ZM_hnf(ZM_mul(m1,U));
    4728             :   }
    4729       21133 :   return d;
    4730             : }
    4731             : 
    4732             : GEN
    4733        7070 : alglatinter(GEN al, GEN lat1, GEN lat2, GEN* ptsum)
    4734             : {
    4735        7070 :   pari_sp av = avma;
    4736             :   GEN inter, d;
    4737        7070 :   d = alglataddinter(al, lat1, lat2, ptsum, &inter);
    4738        7070 :   inter = primlat(mkvec2(inter, d));
    4739        7070 :   if (ptsum)
    4740             :   {
    4741          14 :     *ptsum = primlat(mkvec2(*ptsum,d));
    4742          14 :     gerepileall(av, 2, &inter, ptsum);
    4743             :   }
    4744        7056 :   else inter = gerepilecopy(av, inter);
    4745        7070 :   return inter;
    4746             : }
    4747             : 
    4748             : GEN
    4749       14063 : alglatadd(GEN al, GEN lat1, GEN lat2, GEN* ptinter)
    4750             : {
    4751       14063 :   pari_sp av = avma;
    4752             :   GEN sum, d;
    4753       14063 :   d = alglataddinter(al, lat1, lat2, &sum, ptinter);
    4754       14063 :   sum = primlat(mkvec2(sum, d));
    4755       14063 :   if (ptinter)
    4756             :   {
    4757        7014 :     *ptinter = primlat(mkvec2(*ptinter,d));
    4758        7014 :     gerepileall(av, 2, &sum, ptinter);
    4759             :   }
    4760        7049 :   else sum = gerepilecopy(av, sum);
    4761       14063 :   return sum;
    4762             : }
    4763             : 
    4764             : int
    4765       63042 : alglatsubset(GEN al, GEN lat1, GEN lat2, GEN* ptindex)
    4766             : {
    4767             :   /*TODO version that returns the quotient as abelian group?*/
    4768             :   /* return matrices to convert coordinates from one to other? */
    4769       63042 :   pari_sp av = avma;
    4770             :   int res;
    4771             :   GEN m1, m2, m2i, m, t;
    4772       63042 :   checkalg(al);
    4773       63042 :   checklat(al,lat1);
    4774       63042 :   checklat(al,lat2);
    4775       63042 :   m1 = alglat_get_primbasis(lat1);
    4776       63042 :   m2 = alglat_get_primbasis(lat2);
    4777       63042 :   m2i = RgM_inv_upper(m2);
    4778       63042 :   t = gdiv(alglat_get_scalar(lat1), alglat_get_scalar(lat2));
    4779       63042 :   m = RgM_Rg_mul(RgM_mul(m2i,m1), t);
    4780       63042 :   res = RgM_is_ZM(m);
    4781       63042 :   if (res && ptindex)
    4782             :   {
    4783        3507 :     *ptindex = mpabs(ZM_det_triangular(m));
    4784        3507 :     gerepileall(av,1,ptindex);
    4785             :   }
    4786       59535 :   else avma = av;
    4787       63042 :   return res;
    4788             : }
    4789             : 
    4790             : GEN
    4791       10514 : alglatindex(GEN al, GEN lat1, GEN lat2)
    4792             : {
    4793       10514 :   pari_sp av = avma;
    4794             :   long N;
    4795             :   GEN res;
    4796       10514 :   checkalg(al);
    4797       10514 :   checklat(al,lat1);
    4798       10514 :   checklat(al,lat2);
    4799       10514 :   N = alg_get_absdim(al);
    4800       10514 :   res = alglat_get_scalar(lat1);
    4801       10514 :   res = gdiv(res, alglat_get_scalar(lat2));
    4802       10514 :   res = gpowgs(res, N);
    4803       10514 :   res = gmul(res,RgM_det_triangular(alglat_get_primbasis(lat1)));
    4804       10514 :   res = gdiv(res, RgM_det_triangular(alglat_get_primbasis(lat2)));
    4805       10514 :   res = gabs(res,0);
    4806       10514 :   return gerepilecopy(av, res);
    4807             : }
    4808             : 
    4809             : GEN
    4810       91077 : alglatmul(GEN al, GEN lat1, GEN lat2)
    4811             : {
    4812       91077 :   pari_sp av = avma;
    4813             :   long N,i;
    4814             :   GEN m1, m2, m, V, lat, t, d;
    4815       91077 :   checkalg(al);
    4816       91077 :   if (typ(lat1)==t_COL)
    4817             :   {
    4818       38535 :     if (typ(lat2)==t_COL)
    4819           7 :       pari_err_TYPE("alglatmul [one of lat1, lat2 has to be a lattice]", lat2);
    4820       38528 :     checklat(al,lat2);
    4821       38528 :     m = algbasismultable(al,lat1);
    4822       38528 :     m2 = alglat_get_primbasis(lat2);
    4823       38528 :     m = RgM_mul(m,m2);
    4824       38528 :     t = alglat_get_scalar(lat2);
    4825             :   }
    4826             :   else /* typ(lat1)!=t_COL */
    4827             :   {
    4828       52542 :     checklat(al,lat1);
    4829       52542 :     if (typ(lat2)==t_COL)
    4830             :     {
    4831       38528 :       lat2 = Q_remove_denom(lat2,&d);
    4832       38528 :       m = algbasisrightmultable(al,lat2);
    4833       38528 :       m1 = alglat_get_primbasis(lat1);
    4834       38528 :       m = ZM_mul(m,m1);
    4835       38528 :       t = alglat_get_scalar(lat1);
    4836       38528 :       if (d) t = gdiv(t,d);
    4837             :     }
    4838             :     else
    4839             :     {
    4840       14014 :       checklat(al,lat2);
    4841       14014 :       N = algabsdim(al);
    4842       14014 :       m1 = alglat_get_primbasis(lat1);
    4843       14014 :       m2 = alglat_get_primbasis(lat2);
    4844       14014 :       V = cgetg(N+1,t_VEC);
    4845      126126 :       for (i=1; i<=N; i++) {
    4846      112112 :         gel(V,i) = algbasismultable(al,gel(m1,i));
    4847      112112 :         gel(V,i) = ZM_mul(gel(V,i),m2);
    4848             :       }
    4849       14014 :       m = matconcat(V);
    4850       14014 :       t = gmul(alglat_get_scalar(lat1), alglat_get_scalar(lat2));
    4851             :     }
    4852             :   }
    4853             : 
    4854             :   /* TODO optimise by computing d such that d*O <= I1*I2 */
    4855       91070 :   lat = alglathnf(al,m,0);
    4856       91070 :   gel(lat,2) = gmul(alglat_get_scalar(lat), t);
    4857       91070 :   lat = primlat(lat);
    4858       91070 :   return gerepilecopy(av, lat);
    4859             : }
    4860             : 
    4861             : int
    4862       35021 : alglatcontains(GEN al, GEN lat, GEN x, GEN *ptc)
    4863             : {
    4864       35021 :   pari_sp av = avma;
    4865             :   GEN m, t, sol;
    4866       35021 :   checkalg(al);
    4867       35021 :   checklat(al,lat);
    4868       35021 :   m = alglat_get_primbasis(lat);
    4869       35021 :   t = alglat_get_scalar(lat);
    4870       35021 :   x = RgC_Rg_div(x,t);
    4871       35021 :   if (!RgC_is_ZC(x)) { avma = av; return 0; }
    4872       35021 :   sol = hnf_solve(m,x);
    4873       35021 :   if (!sol) { avma = av; return 0; }
    4874       17521 :   if (ptc)
    4875             :   {
    4876       17514 :     *ptc = sol;
    4877       17514 :     gerepileall(av,1,ptc);
    4878             :   }
    4879           7 :   else avma = av;
    4880       17521 :   return 1;
    4881             : }
    4882             : 
    4883             : GEN
    4884       17521 : alglatelement(GEN al, GEN lat, GEN c)
    4885             : {
    4886       17521 :   pari_sp av = avma;
    4887             :   GEN res;
    4888       17521 :   checkalg(al);
    4889       17521 :   checklat(al,lat);
    4890       17521 :   if (typ(c)!=t_COL) pari_err_TYPE("alglatelement", c);
    4891       17514 :   res = ZM_ZC_mul(alglat_get_primbasis(lat),c);
    4892       17514 :   res = RgC_Rg_mul(res, alglat_get_scalar(lat));
    4893       17514 :   return gerepilecopy(av,res);
    4894             : }
    4895             : 
    4896             : /* If m is injective, computes a Z-basis of the submodule of elements whose
    4897             :  * image under m is integral */
    4898             : static GEN
    4899        7042 : QM_invimZ(GEN m)
    4900             : {
    4901        7042 :   return RgM_invimage(m, QM_ImQ_hnf(m));
    4902             : }
    4903             : 
    4904             : /* An isomorphism of Z-modules M_{m,n}(Z) -> Z^{m*n} */
    4905             : static GEN
    4906       56266 : mat2col(GEN M, long m, long n)
    4907             : {
    4908             :   long i,j,k,p;
    4909             :   GEN C;
    4910       56266 :   p = m*n;
    4911       56266 :   C = cgetg(p+1,t_COL);
    4912      506198 :   for(i=1,k=1;i<=m;i++)
    4913     4048772 :     for(j=1;j<=n;j++,k++)
    4914     3598840 :       gel(C,k) = gcoeff(M,i,j);
    4915       56266 :   return C;
    4916             : }
    4917             : 
    4918             : static GEN
    4919        7028 : alglattransporter_i(GEN al, GEN lat1, GEN lat2, int right)
    4920             : {
    4921             :   GEN m1, m2, m2i, M, MT, mt, t1, t2, T, d;
    4922             :   long N, i;
    4923        7028 :   N = alg_get_absdim(al);
    4924        7028 :   m1 = alglat_get_primbasis(lat1);
    4925        7028 :   m2 = alglat_get_primbasis(lat2);
    4926        7028 :   m2i = RgM_inv(m2);
    4927        7028 :   t1 = alglat_get_scalar(lat1);
    4928        7028 :   m1 = RgM_Rg_mul(m1,t1);
    4929        7028 :   t2 = alglat_get_scalar(lat2);
    4930        7028 :   m2i = RgM_Rg_div(m2i,t2);
    4931             : 
    4932        7028 :   MT = right? NULL: alg_get_multable(al);
    4933        7028 :   M = cgetg(N+1, t_MAT);
    4934       63252 :   for (i=1; i<=N; i++) {
    4935       56224 :     if (right) mt = algbasisrightmultable(al, vec_ei(N,i));
    4936       28112 :     else       mt = gel(MT,i);
    4937       56224 :     mt = RgM_mul(m2i,mt);
    4938       56224 :     mt = RgM_mul(mt,m1);
    4939       56224 :     gel(M,i) = mat2col(mt, N, N);
    4940             :   }
    4941             : 
    4942        7028 :   T = QM_invimZ(M);
    4943        7028 :   T = Q_remove_denom(T,&d);
    4944        7028 :   T = hnf(T);
    4945        7028 :   T = primlat(mkvec2(T,gen_1));
    4946        7028 :   if (d) gel(T,2) = gdiv(gel(T,2),d);
    4947        7028 :   return T;
    4948             : }
    4949             : 
    4950             : /*
    4951             :    { x in al | x*lat1 subset lat2}
    4952             : */
    4953             : GEN
    4954        3514 : alglatlefttransporter(GEN al, GEN lat1, GEN lat2)
    4955             : {
    4956        3514 :   pari_sp av = avma;
    4957        3514 :   checkalg(al);
    4958        3514 :   checklat(al,lat1);
    4959        3514 :   checklat(al,lat2);
    4960        3514 :   return gerepilecopy(av, alglattransporter_i(al,lat1,lat2,0));
    4961             : }
    4962             : 
    4963             : /*
    4964             :    { x in al | lat1*x subset lat2}
    4965             : */
    4966             : GEN
    4967        3514 : alglatrighttransporter(GEN al, GEN lat1, GEN lat2)
    4968             : {
    4969        3514 :   pari_sp av = avma;
    4970        3514 :   checkalg(al);
    4971        3514 :   checklat(al,lat1);
    4972        3514 :   checklat(al,lat2);
    4973        3514 :   return gerepilecopy(av, alglattransporter_i(al,lat1,lat2,1));
    4974             : }
    4975             : 
    4976             : GEN
    4977          42 : algmakeintegral(GEN mt0, int maps)
    4978             : {
    4979          42 :   pari_sp av = avma;
    4980             :   long n,i;
    4981             :   GEN m,P,Pi,mt2,mt;
    4982          42 :   n = lg(mt0)-1;
    4983          42 :   mt = check_mt(mt0,NULL);
    4984          42 :   if (!mt) pari_err_TYPE("algmakeintegral", mt0);
    4985          21 :   if(isint1(Q_denom(mt0))) {
    4986           7 :     if(maps) mt = mkvec3(mt,matid(n),matid(n));
    4987           7 :     return gerepilecopy(av,mt);
    4988             :   }
    4989          14 :   dbg_printf(2)(" algmakeintegral: dim=%d, denom=%Ps\n", n, Q_denom(mt0));
    4990          14 :   m = cgetg(n+1,t_MAT);
    4991          56 :   for(i=1;i<=n;i++)
    4992          42 :     gel(m,i) = mat2col(gel(mt,i),n,n);
    4993          14 :   dbg_printf(2)(" computing order, dims m = %d x %d...\n", nbrows(m), lg(m)-1);
    4994          14 :   P = QM_invimZ(m);
    4995          14 :   dbg_printf(2)(" ...done.\n");
    4996          14 :   P = shallowmatconcat(mkvec2(col_ei(n,1),P));
    4997          14 :   P = hnf(P);
    4998          14 :   Pi = RgM_inv(P);
    4999          14 :   mt2 = change_Rgmultable(mt,P,Pi);
    5000          14 :   if(maps) mt2 = mkvec3(mt2,Pi,P); /* mt2, mt->mt2, mt2->mt */
    5001          14 :   return gerepilecopy(av,mt2);
    5002             : }
    5003             : 
    5004             : /** ORDERS **/
    5005             : 
    5006             : /** IDEALS **/
    5007             : 

Generated by: LCOV version 1.11