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 21947-4fc3047) Lines: 2718 2896 93.9 %
Date: 2018-02-24 06:16:21 Functions: 262 273 96.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11