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 - basemath - base3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 20443-183d202) Lines: 1579 1683 93.8 %
Date: 2017-03-27 05:17:48 Functions: 177 187 94.7 %
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             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*                       BASIC NF OPERATIONS                       */
      17             : /*                                                                 */
      18             : /*******************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : /*******************************************************************/
      23             : /*                                                                 */
      24             : /*                OPERATIONS OVER NUMBER FIELD ELEMENTS.           */
      25             : /*     represented as column vectors over the integral basis       */
      26             : /*                                                                 */
      27             : /*******************************************************************/
      28             : static GEN
      29     5461974 : get_tab(GEN nf, long *N)
      30             : {
      31     5461974 :   GEN tab = (typ(nf) == t_MAT)? nf: gel(nf,9);
      32     5461974 :   *N = nbrows(tab); return tab;
      33             : }
      34             : 
      35             : /* x != 0, y t_INT. Return x * y (not memory clean if x = 1) */
      36             : static GEN
      37   240270678 : _mulii(GEN x, GEN y) {
      38   606856644 :   return is_pm1(x)? (signe(x) < 0)? negi(y): y
      39   366585966 :                   : mulii(x, y);
      40             : }
      41             : 
      42             : GEN
      43        1675 : tablemul_ei_ej(GEN M, long i, long j)
      44             : {
      45             :   long N;
      46        1675 :   GEN tab = get_tab(M, &N);
      47        1675 :   tab += (i-1)*N; return gel(tab,j);
      48             : }
      49             : 
      50             : /* Outputs x.ei, where ei is the i-th elt of the algebra basis.
      51             :  * x an RgV of correct length and arbitrary content (polynomials, scalars...).
      52             :  * M is the multiplication table ei ej = sum_k M_k^(i,j) ek */
      53             : GEN
      54        2310 : tablemul_ei(GEN M, GEN x, long i)
      55             : {
      56             :   long j, k, N;
      57             :   GEN v, tab;
      58             : 
      59        2310 :   if (i==1) return gcopy(x);
      60        2310 :   tab = get_tab(M, &N);
      61        2310 :   if (typ(x) != t_COL) { v = zerocol(N); gel(v,i) = gcopy(x); return v; }
      62        2310 :   tab += (i-1)*N; v = cgetg(N+1,t_COL);
      63             :   /* wi . x = [ sum_j tab[k,j] x[j] ]_k */
      64       15180 :   for (k=1; k<=N; k++)
      65             :   {
      66       12870 :     pari_sp av = avma;
      67       12870 :     GEN s = gen_0;
      68       90330 :     for (j=1; j<=N; j++)
      69             :     {
      70       77460 :       GEN c = gcoeff(tab,k,j);
      71       77460 :       if (!gequal0(c)) s = gadd(s, gmul(c, gel(x,j)));
      72             :     }
      73       12870 :     gel(v,k) = gerepileupto(av,s);
      74             :   }
      75        2310 :   return v;
      76             : }
      77             : /* as tablemul_ei, assume x a ZV of correct length */
      78             : GEN
      79     4638144 : zk_ei_mul(GEN nf, GEN x, long i)
      80             : {
      81             :   long j, k, N;
      82             :   GEN v, tab;
      83             : 
      84     4638144 :   if (i==1) return ZC_copy(x);
      85     4638134 :   tab = get_tab(nf, &N); tab += (i-1)*N;
      86     4638134 :   v = cgetg(N+1,t_COL);
      87    34842468 :   for (k=1; k<=N; k++)
      88             :   {
      89    30204334 :     pari_sp av = avma;
      90    30204334 :     GEN s = gen_0;
      91   394324410 :     for (j=1; j<=N; j++)
      92             :     {
      93   364120076 :       GEN c = gcoeff(tab,k,j);
      94   364120076 :       if (signe(c)) s = addii(s, _mulii(c, gel(x,j)));
      95             :     }
      96    30204334 :     gel(v,k) = gerepileuptoint(av, s);
      97             :   }
      98     4638134 :   return v;
      99             : }
     100             : 
     101             : /* table of multiplication by wi in R[w1,..., wN] */
     102             : GEN
     103         625 : ei_multable(GEN TAB, long i)
     104             : {
     105             :   long k,N;
     106         625 :   GEN m, tab = get_tab(TAB, &N);
     107         625 :   tab += (i-1)*N;
     108         625 :   m = cgetg(N+1,t_MAT);
     109         625 :   for (k=1; k<=N; k++) gel(m,k) = gel(tab,k);
     110         625 :   return m;
     111             : }
     112             : 
     113             : GEN
     114     1890516 : zk_multable(GEN nf, GEN x)
     115             : {
     116     1890516 :   long i, l = lg(x);
     117     1890516 :   GEN mul = cgetg(l,t_MAT);
     118     1890516 :   gel(mul,1) = x; /* assume w_1 = 1 */
     119     1890516 :   for (i=2; i<l; i++) gel(mul,i) = zk_ei_mul(nf,x,i);
     120     1890516 :   return mul;
     121             : }
     122             : GEN
     123         515 : multable(GEN M, GEN x)
     124             : {
     125             :   long i, N;
     126             :   GEN mul;
     127         515 :   if (typ(x) == t_MAT) return x;
     128           0 :   M = get_tab(M, &N);
     129           0 :   if (typ(x) != t_COL) return scalarmat(x, N);
     130           0 :   mul = cgetg(N+1,t_MAT);
     131           0 :   gel(mul,1) = x; /* assume w_1 = 1 */
     132           0 :   for (i=2; i<=N; i++) gel(mul,i) = tablemul_ei(M,x,i);
     133           0 :   return mul;
     134             : }
     135             : 
     136             : /* x integral in nf; table of multiplication by x in ZK = Z[w1,..., wN].
     137             :  * Return a t_INT if x is scalar, and a ZM otherwise */
     138             : GEN
     139     1058026 : zk_scalar_or_multable(GEN nf, GEN x)
     140             : {
     141     1058026 :   long tx = typ(x);
     142     1058026 :   if (tx == t_MAT || tx == t_INT) return x;
     143     1044576 :   x = nf_to_scalar_or_basis(nf, x);
     144     1044576 :   return (typ(x) == t_COL)? zk_multable(nf, x): x;
     145             : }
     146             : 
     147             : GEN
     148       16830 : nftrace(GEN nf, GEN x)
     149             : {
     150       16830 :   pari_sp av = avma;
     151       16830 :   nf = checknf(nf);
     152       16830 :   x = nf_to_scalar_or_basis(nf, x);
     153       50475 :   x = (typ(x) == t_COL)? RgV_dotproduct(x, gel(nf_get_Tr(nf),1))
     154       33645 :                        : gmulgs(x, nf_get_degree(nf));
     155       16830 :   return gerepileupto(av, x);
     156             : }
     157             : GEN
     158         405 : rnfelttrace(GEN rnf, GEN x)
     159             : {
     160         405 :   pari_sp av = avma;
     161         405 :   checkrnf(rnf);
     162         405 :   x = rnfeltabstorel(rnf, x);
     163         960 :   x = (typ(x) == t_POLMOD)? rnfeltdown(rnf, gtrace(x))
     164         620 :                           : gmulgs(x, rnf_get_degree(rnf));
     165         340 :   return gerepileupto(av, x);
     166             : }
     167             : 
     168             : /* assume nf is a genuine nf, fa a famat */
     169             : static GEN
     170           5 : famat_norm(GEN nf, GEN fa)
     171             : {
     172           5 :   pari_sp av = avma;
     173           5 :   GEN g = gel(fa,1), e = gel(fa,2), N = gen_1;
     174           5 :   long i, l = lg(g);
     175          15 :   for (i = 1; i < l; i++)
     176          10 :     N = gmul(N, powgi(nfnorm(nf, gel(g,i)), gel(e,i)));
     177           5 :   return gerepileupto(av, N);
     178             : }
     179             : GEN
     180       14570 : nfnorm(GEN nf, GEN x)
     181             : {
     182       14570 :   pari_sp av = avma;
     183       14570 :   nf = checknf(nf);
     184       14570 :   if (typ(x) == t_MAT) return famat_norm(nf, x);
     185       14565 :   x = nf_to_scalar_or_alg(nf, x);
     186       42505 :   x = (typ(x) == t_POL)? RgXQ_norm(x, nf_get_pol(nf))
     187       27940 :                        : gpowgs(x, nf_get_degree(nf));
     188       14565 :   return gerepileupto(av, x);
     189             : }
     190             : 
     191             : GEN
     192         165 : rnfeltnorm(GEN rnf, GEN x)
     193             : {
     194         165 :   pari_sp av = avma;
     195         165 :   checkrnf(rnf);
     196         165 :   x = rnfeltabstorel(rnf, x);
     197         270 :   x = (typ(x) == t_POLMOD)? rnfeltdown(rnf, gnorm(x))
     198         170 :                           : gpowgs(x, rnf_get_degree(rnf));
     199         100 :   return gerepileupto(av, x);
     200             : }
     201             : 
     202             : /* x + y in nf */
     203             : GEN
     204     1148270 : nfadd(GEN nf, GEN x, GEN y)
     205             : {
     206     1148270 :   pari_sp av = avma;
     207             :   GEN z;
     208             : 
     209     1148270 :   nf = checknf(nf);
     210     1148270 :   x = nf_to_scalar_or_basis(nf, x);
     211     1148270 :   y = nf_to_scalar_or_basis(nf, y);
     212     1148270 :   if (typ(x) != t_COL)
     213      805590 :   { z = (typ(y) == t_COL)? RgC_Rg_add(y, x): gadd(x,y); }
     214             :   else
     215      342680 :   { z = (typ(y) == t_COL)? RgC_add(x, y): RgC_Rg_add(x, y); }
     216     1148270 :   return gerepileupto(av, z);
     217             : }
     218             : /* x - y in nf */
     219             : GEN
     220       52625 : nfsub(GEN nf, GEN x, GEN y)
     221             : {
     222       52625 :   pari_sp av = avma;
     223             :   GEN z;
     224             : 
     225       52625 :   nf = checknf(nf);
     226       52625 :   x = nf_to_scalar_or_basis(nf, x);
     227       52625 :   y = nf_to_scalar_or_basis(nf, y);
     228       52625 :   if (typ(x) != t_COL)
     229       19050 :   { z = (typ(y) == t_COL)? Rg_RgC_sub(x,y): gsub(x,y); }
     230             :   else
     231       33575 :   { z = (typ(y) == t_COL)? RgC_sub(x,y): RgC_Rg_sub(x,y); }
     232       52625 :   return gerepileupto(av, z);
     233             : }
     234             : 
     235             : /* product of x and y in nf */
     236             : GEN
     237     2566675 : nfmul(GEN nf, GEN x, GEN y)
     238             : {
     239             :   GEN z;
     240     2566675 :   pari_sp av = avma;
     241             : 
     242     2566675 :   if (x == y) return nfsqr(nf,x);
     243             : 
     244     2237615 :   nf = checknf(nf);
     245     2237615 :   x = nf_to_scalar_or_basis(nf, x);
     246     2237615 :   y = nf_to_scalar_or_basis(nf, y);
     247     2237615 :   if (typ(x) != t_COL)
     248             :   {
     249     1758430 :     if (isintzero(x)) return gen_0;
     250     1362720 :     z = (typ(y) == t_COL)? RgC_Rg_mul(y, x): gmul(x,y); }
     251             :   else
     252             :   {
     253      479185 :     if (typ(y) != t_COL)
     254             :     {
     255      328725 :       if (isintzero(y)) return gen_0;
     256       81180 :       z = RgC_Rg_mul(x, y);
     257             :     }
     258             :     else
     259             :     {
     260             :       GEN dx, dy;
     261      150460 :       x = Q_remove_denom(x, &dx);
     262      150460 :       y = Q_remove_denom(y, &dy);
     263      150460 :       z = nfmuli(nf,x,y);
     264      150460 :       dx = mul_denom(dx,dy);
     265      150460 :       if (dx) z = RgC_Rg_div(z, dx);
     266             :     }
     267             :   }
     268     1594360 :   return gerepileupto(av, z);
     269             : }
     270             : /* square of x in nf */
     271             : GEN
     272      413070 : nfsqr(GEN nf, GEN x)
     273             : {
     274      413070 :   pari_sp av = avma;
     275             :   GEN z;
     276             : 
     277      413070 :   nf = checknf(nf);
     278      413070 :   x = nf_to_scalar_or_basis(nf, x);
     279      413070 :   if (typ(x) != t_COL) z = gsqr(x);
     280             :   else
     281             :   {
     282             :     GEN dx;
     283       50310 :     x = Q_remove_denom(x, &dx);
     284       50310 :     z = nfsqri(nf,x);
     285       50310 :     if (dx) z = RgC_Rg_div(z, sqri(dx));
     286             :   }
     287      413070 :   return gerepileupto(av, z);
     288             : }
     289             : 
     290             : /* x a ZC, v a t_COL of ZC/Z */
     291             : GEN
     292       75515 : zkC_multable_mul(GEN v, GEN x)
     293             : {
     294       75515 :   long i, l = lg(v);
     295       75515 :   GEN y = cgetg(l, t_COL);
     296      273975 :   for (i = 1; i < l; i++)
     297             :   {
     298      198460 :     GEN c = gel(v,i);
     299      198460 :     if (typ(c)!=t_COL) {
     300           0 :       if (!isintzero(c)) c = ZC_Z_mul(gel(x,1), c);
     301             :     } else {
     302      198460 :       c = ZM_ZC_mul(x,c);
     303      198460 :       if (ZV_isscalar(c)) c = gel(c,1);
     304             :     }
     305      198460 :     gel(y,i) = c;
     306             :   }
     307       75515 :   return y;
     308             : }
     309             : 
     310             : GEN
     311       16950 : nfC_multable_mul(GEN v, GEN x)
     312             : {
     313       16950 :   long i, l = lg(v);
     314       16950 :   GEN y = cgetg(l, t_COL);
     315       98030 :   for (i = 1; i < l; i++)
     316             :   {
     317       81080 :     GEN c = gel(v,i);
     318       81080 :     if (typ(c)!=t_COL) {
     319       64480 :       if (!isintzero(c)) c = RgC_Rg_mul(gel(x,1), c);
     320             :     } else {
     321       16600 :       c = RgM_RgC_mul(x,c);
     322       16600 :       if (QV_isscalar(c)) c = gel(c,1);
     323             :     }
     324       81080 :     gel(y,i) = c;
     325             :   }
     326       16950 :   return y;
     327             : }
     328             : 
     329             : GEN
     330       52450 : nfC_nf_mul(GEN nf, GEN v, GEN x)
     331             : {
     332             :   long tx;
     333             :   GEN y;
     334             : 
     335       52450 :   x = nf_to_scalar_or_basis(nf, x);
     336       52450 :   tx = typ(x);
     337       52450 :   if (tx != t_COL)
     338             :   {
     339             :     long l, i;
     340       36135 :     if (tx == t_INT)
     341             :     {
     342       35105 :       long s = signe(x);
     343       35105 :       if (!s) return zerocol(lg(v)-1);
     344       32865 :       if (is_pm1(x)) return s > 0? leafcopy(v): RgC_neg(v);
     345             :     }
     346       10040 :     l = lg(v); y = cgetg(l, t_COL);
     347       75990 :     for (i=1; i < l; i++)
     348             :     {
     349       65950 :       GEN c = gel(v,i);
     350       65950 :       if (typ(c) != t_COL) c = gmul(c, x); else c = RgC_Rg_mul(c, x);
     351       65950 :       gel(y,i) = c;
     352             :     }
     353       10040 :     return y;
     354             :   }
     355             :   else
     356             :   {
     357             :     GEN dx;
     358       16315 :     x = zk_multable(nf, Q_remove_denom(x,&dx));
     359       16315 :     y = nfC_multable_mul(v, x);
     360       16315 :     return dx? RgC_Rg_div(y, dx): y;
     361             :   }
     362             : }
     363             : static GEN
     364        2320 : mulbytab(GEN M, GEN c)
     365        2320 : { return typ(c) == t_COL? RgM_RgC_mul(M,c): RgC_Rg_mul(gel(M,1), c); }
     366             : GEN
     367         515 : tablemulvec(GEN M, GEN x, GEN v)
     368             : {
     369             :   long l, i;
     370             :   GEN y;
     371             : 
     372         515 :   if (typ(x) == t_COL && RgV_isscalar(x))
     373             :   {
     374           0 :     x = gel(x,1);
     375           0 :     return typ(v) == t_POL? RgX_Rg_mul(v,x): RgV_Rg_mul(v,x);
     376             :   }
     377         515 :   x = multable(M, x); /* multiplication table by x */
     378         515 :   y = cgetg_copy(v, &l);
     379         515 :   if (typ(v) == t_POL)
     380             :   {
     381         515 :     y[1] = v[1];
     382         515 :     for (i=2; i < l; i++) gel(y,i) = mulbytab(x, gel(v,i));
     383         515 :     y = normalizepol(y);
     384             :   }
     385             :   else
     386             :   {
     387           0 :     for (i=1; i < l; i++) gel(y,i) = mulbytab(x, gel(v,i));
     388             :   }
     389         515 :   return y;
     390             : }
     391             : 
     392             : GEN
     393      228345 : zkmultable_capZ(GEN mx) { return Q_denom(zkmultable_inv(mx)); }
     394             : 
     395             : GEN
     396      240620 : zkmultable_inv(GEN mx)
     397      240620 : { return ZM_gauss(mx, col_ei(lg(mx)-1,1)); }
     398             : 
     399             : /* nf a true nf, x a ZC */
     400             : GEN
     401       12275 : zk_inv(GEN nf, GEN x) { return zkmultable_inv(zk_multable(nf,x)); }
     402             : 
     403             : /* inverse of x in nf */
     404             : GEN
     405        8220 : nfinv(GEN nf, GEN x)
     406             : {
     407        8220 :   pari_sp av = avma;
     408             :   GEN z;
     409             : 
     410        8220 :   nf = checknf(nf);
     411        8220 :   x = nf_to_scalar_or_basis(nf, x);
     412        8220 :   if (typ(x) == t_COL)
     413             :   {
     414             :     GEN d;
     415        4170 :     x = Q_remove_denom(x, &d);
     416        4170 :     z = zk_inv(nf, x);
     417        4170 :     if (d) z = RgC_Rg_mul(z, d);
     418             :   }
     419             :   else
     420        4050 :     z = ginv(x);
     421        8220 :   return gerepileupto(av, z);
     422             : }
     423             : 
     424             : /* quotient of x and y in nf */
     425             : GEN
     426       10055 : nfdiv(GEN nf, GEN x, GEN y)
     427             : {
     428       10055 :   pari_sp av = avma;
     429             :   GEN z;
     430             : 
     431       10055 :   nf = checknf(nf);
     432       10055 :   y = nf_to_scalar_or_basis(nf, y);
     433       10055 :   if (typ(y) != t_COL)
     434             :   {
     435        4880 :     x = nf_to_scalar_or_basis(nf, x);
     436        4880 :     z = (typ(x) == t_COL)? RgC_Rg_div(x, y): gdiv(x,y);
     437             :   }
     438             :   else
     439             :   {
     440             :     GEN d;
     441        5175 :     y = Q_remove_denom(y, &d);
     442        5175 :     z = nfmul(nf, x, zk_inv(nf,y));
     443        5175 :     if (d) z = RgC_Rg_mul(z, d);
     444             :   }
     445       10055 :   return gerepileupto(av, z);
     446             : }
     447             : 
     448             : /* product of INTEGERS (t_INT or ZC) x and y in nf
     449             :  * compute xy as ( sum_i x_i sum_j y_j m^{i,j}_k )_k */
     450             : GEN
     451      425175 : nfmuli(GEN nf, GEN x, GEN y)
     452             : {
     453             :   long i, j, k, N;
     454      425175 :   GEN s, v, TAB = get_tab(nf, &N);
     455             : 
     456      425175 :   if (typ(x) == t_INT) return (typ(y) == t_COL)? ZC_Z_mul(y, x): mulii(x,y);
     457      371140 :   if (typ(y) == t_INT) return ZC_Z_mul(x, y);
     458             :   /* both x and y are ZV */
     459      352555 :   v = cgetg(N+1,t_COL);
     460     1760190 :   for (k=1; k<=N; k++)
     461             :   {
     462     1407635 :     pari_sp av = avma;
     463     1407635 :     GEN TABi = TAB;
     464     1407635 :     if (k == 1)
     465      352555 :       s = mulii(gel(x,1),gel(y,1));
     466             :     else
     467     2110160 :       s = addii(mulii(gel(x,1),gel(y,k)),
     468     2110160 :                 mulii(gel(x,k),gel(y,1)));
     469     8907985 :     for (i=2; i<=N; i++)
     470             :     {
     471     7500350 :       GEN t, xi = gel(x,i);
     472     7500350 :       TABi += N;
     473     7500350 :       if (!signe(xi)) continue;
     474             : 
     475     4384525 :       t = NULL;
     476    53886685 :       for (j=2; j<=N; j++)
     477             :       {
     478    49502160 :         GEN p1, c = gcoeff(TABi, k, j); /* m^{i,j}_k */
     479    49502160 :         if (!signe(c)) continue;
     480    24041935 :         p1 = _mulii(c, gel(y,j));
     481    24041935 :         t = t? addii(t, p1): p1;
     482             :       }
     483     4384525 :       if (t) s = addii(s, mulii(xi, t));
     484             :     }
     485     1407635 :     gel(v,k) = gerepileuptoint(av,s);
     486             :   }
     487      352555 :   return v;
     488             : }
     489             : /* square of INTEGER (t_INT or ZC) x in nf */
     490             : GEN
     491      394055 : nfsqri(GEN nf, GEN x)
     492             : {
     493             :   long i, j, k, N;
     494      394055 :   GEN s, v, TAB = get_tab(nf, &N);
     495             : 
     496      394055 :   if (typ(x) == t_INT) return sqri(x);
     497      394055 :   v = cgetg(N+1,t_COL);
     498     3022620 :   for (k=1; k<=N; k++)
     499             :   {
     500     2628565 :     pari_sp av = avma;
     501     2628565 :     GEN TABi = TAB;
     502     2628565 :     if (k == 1)
     503      394055 :       s = sqri(gel(x,1));
     504             :     else
     505     2234510 :       s = shifti(mulii(gel(x,1),gel(x,k)), 1);
     506    32392185 :     for (i=2; i<=N; i++)
     507             :     {
     508    29763620 :       GEN p1, c, t, xi = gel(x,i);
     509    29763620 :       TABi += N;
     510    29763620 :       if (!signe(xi)) continue;
     511             : 
     512     9625960 :       c = gcoeff(TABi, k, i);
     513     9625960 :       t = signe(c)? _mulii(c,xi): NULL;
     514   163902215 :       for (j=i+1; j<=N; j++)
     515             :       {
     516   154276255 :         c = gcoeff(TABi, k, j);
     517   154276255 :         if (!signe(c)) continue;
     518    83499240 :         p1 = _mulii(c, shifti(gel(x,j),1));
     519    83499240 :         t = t? addii(t, p1): p1;
     520             :       }
     521     9625960 :       if (t) s = addii(s, mulii(xi, t));
     522             :     }
     523     2628565 :     gel(v,k) = gerepileuptoint(av,s);
     524             :   }
     525      394055 :   return v;
     526             : }
     527             : 
     528             : /* both x and y are RgV */
     529             : GEN
     530           0 : tablemul(GEN TAB, GEN x, GEN y)
     531             : {
     532             :   long i, j, k, N;
     533             :   GEN s, v;
     534           0 :   if (typ(x) != t_COL) return gmul(x, y);
     535           0 :   if (typ(y) != t_COL) return gmul(y, x);
     536           0 :   N = lg(x)-1;
     537           0 :   v = cgetg(N+1,t_COL);
     538           0 :   for (k=1; k<=N; k++)
     539             :   {
     540           0 :     pari_sp av = avma;
     541           0 :     GEN TABi = TAB;
     542           0 :     if (k == 1)
     543           0 :       s = gmul(gel(x,1),gel(y,1));
     544             :     else
     545           0 :       s = gadd(gmul(gel(x,1),gel(y,k)),
     546           0 :                gmul(gel(x,k),gel(y,1)));
     547           0 :     for (i=2; i<=N; i++)
     548             :     {
     549           0 :       GEN t, xi = gel(x,i);
     550           0 :       TABi += N;
     551           0 :       if (gequal0(xi)) continue;
     552             : 
     553           0 :       t = NULL;
     554           0 :       for (j=2; j<=N; j++)
     555             :       {
     556           0 :         GEN p1, c = gcoeff(TABi, k, j); /* m^{i,j}_k */
     557           0 :         if (gequal0(c)) continue;
     558           0 :         p1 = gmul(c, gel(y,j));
     559           0 :         t = t? gadd(t, p1): p1;
     560             :       }
     561           0 :       if (t) s = gadd(s, gmul(xi, t));
     562             :     }
     563           0 :     gel(v,k) = gerepileupto(av,s);
     564             :   }
     565           0 :   return v;
     566             : }
     567             : GEN
     568        3880 : tablesqr(GEN TAB, GEN x)
     569             : {
     570             :   long i, j, k, N;
     571             :   GEN s, v;
     572             : 
     573        3880 :   if (typ(x) != t_COL) return gsqr(x);
     574        3880 :   N = lg(x)-1;
     575        3880 :   v = cgetg(N+1,t_COL);
     576             : 
     577       29570 :   for (k=1; k<=N; k++)
     578             :   {
     579       25690 :     pari_sp av = avma;
     580       25690 :     GEN TABi = TAB;
     581       25690 :     if (k == 1)
     582        3880 :       s = gsqr(gel(x,1));
     583             :     else
     584       21810 :       s = gmul2n(gmul(gel(x,1),gel(x,k)), 1);
     585      179840 :     for (i=2; i<=N; i++)
     586             :     {
     587      154150 :       GEN p1, c, t, xi = gel(x,i);
     588      154150 :       TABi += N;
     589      154150 :       if (gequal0(xi)) continue;
     590             : 
     591       51630 :       c = gcoeff(TABi, k, i);
     592       51630 :       t = !gequal0(c)? gmul(c,xi): NULL;
     593      233925 :       for (j=i+1; j<=N; j++)
     594             :       {
     595      182295 :         c = gcoeff(TABi, k, j);
     596      182295 :         if (gequal0(c)) continue;
     597       96030 :         p1 = gmul(gmul2n(c,1), gel(x,j));
     598       96030 :         t = t? gadd(t, p1): p1;
     599             :       }
     600       51630 :       if (t) s = gadd(s, gmul(xi, t));
     601             :     }
     602       25690 :     gel(v,k) = gerepileupto(av,s);
     603             :   }
     604        3880 :   return v;
     605             : }
     606             : 
     607             : static GEN
     608       18165 : _mul(void *data, GEN x, GEN y) { return nfmuli((GEN)data,x,y); }
     609             : static GEN
     610       65860 : _sqr(void *data, GEN x) { return nfsqri((GEN)data,x); }
     611             : 
     612             : /* Compute z^n in nf, left-shift binary powering */
     613             : GEN
     614       63875 : nfpow(GEN nf, GEN z, GEN n)
     615             : {
     616       63875 :   pari_sp av = avma;
     617             :   long s;
     618             :   GEN x, cx;
     619             : 
     620       63875 :   if (typ(n)!=t_INT) pari_err_TYPE("nfpow",n);
     621       63875 :   nf = checknf(nf);
     622       63875 :   s = signe(n); if (!s) return gen_1;
     623       63875 :   x = nf_to_scalar_or_basis(nf, z);
     624       63875 :   if (typ(x) != t_COL) return powgi(x,n);
     625       63535 :   if (s < 0)
     626             :   { /* simplified nfinv */
     627             :     GEN d;
     628        1000 :     x = Q_remove_denom(x, &d);
     629        1000 :     x = zk_inv(nf, x);
     630        1000 :     x = primitive_part(x, &cx);
     631        1000 :     cx = mul_content(cx, d);
     632        1000 :     n = absi(n);
     633             :   }
     634             :   else
     635       62535 :     x = primitive_part(x, &cx);
     636       63535 :   x = gen_pow(x, n, (void*)nf, _sqr, _mul);
     637       63535 :   if (cx) x = gmul(x, powgi(cx, n));
     638       63535 :   return av==avma? gcopy(x): gerepileupto(av,x);
     639             : }
     640             : /* Compute z^n in nf, left-shift binary powering */
     641             : GEN
     642       22245 : nfpow_u(GEN nf, GEN z, ulong n)
     643             : {
     644       22245 :   pari_sp av = avma;
     645             :   GEN x, cx;
     646             : 
     647       22245 :   nf = checknf(nf);
     648       22245 :   if (!n) return gen_1;
     649       22245 :   x = nf_to_scalar_or_basis(nf, z);
     650       22245 :   if (typ(x) != t_COL) return gpowgs(x,n);
     651        2245 :   x = primitive_part(x, &cx);
     652        2245 :   x = gen_powu(x, n, (void*)nf, _sqr, _mul);
     653        2245 :   if (cx) x = gmul(x, powgi(cx, utoipos(n)));
     654        2245 :   return av==avma? gcopy(x): gerepileupto(av,x);
     655             : }
     656             : 
     657             : static GEN
     658      236485 : _nf_red(void *E, GEN x) { (void)E; return x; }
     659             : 
     660             : static GEN
     661     1070475 : _nf_add(void *E, GEN x, GEN y) { return nfadd((GEN)E,x,y); }
     662             : 
     663             : static GEN
     664       60225 : _nf_neg(void *E, GEN x) { (void)E; return gneg(x); }
     665             : 
     666             : static GEN
     667     1263070 : _nf_mul(void *E, GEN x, GEN y) { return nfmul((GEN)E,x,y); }
     668             : 
     669             : static GEN
     670        3980 : _nf_inv(void *E, GEN x) { return nfinv((GEN)E,x); }
     671             : 
     672             : static GEN
     673         895 : _nf_s(void *E, long x) { (void)E; return stoi(x); }
     674             : 
     675             : static const struct bb_field nf_field={_nf_red,_nf_add,_nf_mul,_nf_neg,
     676             :                                         _nf_inv,&gequal0,_nf_s };
     677             : 
     678       16660 : const struct bb_field *get_nf_field(void **E, GEN nf)
     679       16660 : { *E = (void*)nf; return &nf_field; }
     680             : 
     681             : GEN
     682          10 : nfM_det(GEN nf, GEN M)
     683             : {
     684             :   void *E;
     685          10 :   const struct bb_field *S = get_nf_field(&E, nf);
     686          10 :   return gen_det(M, E, S);
     687             : }
     688             : GEN
     689         885 : nfM_inv(GEN nf, GEN M)
     690             : {
     691             :   void *E;
     692         885 :   const struct bb_field *S = get_nf_field(&E, nf);
     693         885 :   return gen_Gauss(M, matid(lg(M)-1), E, S);
     694             : }
     695             : GEN
     696         820 : nfM_mul(GEN nf, GEN A, GEN B)
     697             : {
     698             :   void *E;
     699         820 :   const struct bb_field *S = get_nf_field(&E, nf);
     700         820 :   return gen_matmul(A, B, E, S);
     701             : }
     702             : GEN
     703       14945 : nfM_nfC_mul(GEN nf, GEN A, GEN B)
     704             : {
     705             :   void *E;
     706       14945 :   const struct bb_field *S = get_nf_field(&E, nf);
     707       14945 :   return gen_matcolmul(A, B, E, S);
     708             : }
     709             : 
     710             : /* valuation of integral x (ZV), with resp. to prime ideal pr */
     711             : long
     712     3546035 : ZC_nfvalrem(GEN x, GEN pr, GEN *newx)
     713             : {
     714             :   long i, v, l;
     715     3546035 :   GEN r, y, p = pr_get_p(pr), mul = pr_get_tau(pr);
     716             : 
     717             :   /* p inert */
     718     3546035 :   if (typ(mul) == t_INT) return newx? ZV_pvalrem(x, p, newx):ZV_pval(x, p);
     719     3541820 :   y = cgetg_copy(x, &l); /* will hold the new x */
     720     3541820 :   x = leafcopy(x);
     721     5229290 :   for(v=0;; v++)
     722             :   {
     723    17996110 :     for (i=1; i<l; i++)
     724             :     { /* is (x.b)[i] divisible by p ? */
     725    16308640 :       gel(y,i) = dvmdii(ZMrow_ZC_mul(mul,x,i),p,&r);
     726    16308640 :       if (r != gen_0) { if (newx) *newx = x; return v; }
     727             :     }
     728     1687470 :     swap(x, y);
     729     1687470 :   }
     730             : }
     731             : long
     732     3388635 : ZC_nfval(GEN x, GEN P)
     733     3388635 : { return ZC_nfvalrem(x, P, NULL); }
     734             : 
     735             : /* v_P(x) != 0, x a ZV. Simpler version of ZC_nfvalrem */
     736             : int
     737      137960 : ZC_prdvd(GEN x, GEN P)
     738             : {
     739      137960 :   pari_sp av = avma;
     740             :   long i, l;
     741      137960 :   GEN p = pr_get_p(P), mul = pr_get_tau(P);
     742      137960 :   if (typ(mul) == t_INT) return ZV_Z_dvd(x, p);
     743      137855 :   l = lg(x);
     744      559805 :   for (i=1; i<l; i++)
     745      510305 :     if (remii(ZMrow_ZC_mul(mul,x,i), p) != gen_0) { avma = av; return 0; }
     746       49500 :   avma = av; return 1;
     747             : }
     748             : 
     749             : int
     750          20 : pr_equal(GEN P, GEN Q)
     751             : {
     752          20 :   GEN gQ, p = pr_get_p(P);
     753          20 :   long e = pr_get_e(P), f = pr_get_f(P), n;
     754          20 :   if (!equalii(p, pr_get_p(Q)) || e != pr_get_e(Q) || f != pr_get_f(Q))
     755          10 :     return 0;
     756          10 :   gQ = pr_get_gen(Q); n = lg(gQ)-1;
     757          10 :   if (2*e*f > n) return 1; /* room for only one such pr */
     758           5 :   return ZV_equal(pr_get_gen(P), gQ) || ZC_prdvd(gQ, P);
     759             : }
     760             : 
     761             : long
     762      937810 : nfval(GEN nf, GEN x, GEN pr)
     763             : {
     764      937810 :   pari_sp av = avma;
     765             :   long w, e;
     766             :   GEN cx, p;
     767             : 
     768      937810 :   if (gequal0(x)) return LONG_MAX;
     769      937445 :   nf = checknf(nf);
     770      937445 :   checkprid(pr);
     771      937445 :   p = pr_get_p(pr);
     772      937445 :   e = pr_get_e(pr);
     773      937445 :   x = nf_to_scalar_or_basis(nf, x);
     774      937445 :   if (typ(x) != t_COL) return e*Q_pval(x,p);
     775       92450 :   x = Q_primitive_part(x, &cx);
     776       92450 :   w = ZC_nfval(x,pr);
     777       92450 :   if (cx) w += e*Q_pval(cx,p);
     778       92450 :   avma = av; return w;
     779             : }
     780             : 
     781             : /* want to write p^v = uniformizer^(e*v) * z^v, z coprime to pr */
     782             : /* z := tau^e / p^(e-1), algebraic integer coprime to pr; return z^v */
     783             : static GEN
     784       14290 : powp(GEN nf, GEN pr, long v)
     785             : {
     786             :   GEN b, z;
     787             :   long e;
     788       14290 :   if (!v) return gen_1;
     789       14180 :   b = pr_get_tau(pr);
     790       14180 :   if (typ(b) == t_INT) return gen_1;
     791         820 :   e = pr_get_e(pr);
     792         820 :   z = gel(b,1);
     793         820 :   if (e != 1) z = gdiv(nfpow_u(nf, z, e), powiu(pr_get_p(pr),e-1));
     794         820 :   if (v < 0) { v = -v; z = nfinv(nf, z); }
     795         820 :   if (v != 1) z = nfpow_u(nf, z, v);
     796         820 :   return z;
     797             : }
     798             : long
     799       46160 : nfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
     800             : {
     801       46160 :   pari_sp av = avma;
     802             :   long w, e;
     803             :   GEN cx, p, t;
     804             : 
     805       46160 :   if (!py) return nfval(nf,x,pr);
     806       46075 :   if (gequal0(x)) { *py = gcopy(x); return LONG_MAX; }
     807       46035 :   nf = checknf(nf);
     808       46035 :   checkprid(pr);
     809       46035 :   p = pr_get_p(pr);
     810       46035 :   e = pr_get_e(pr);
     811       46035 :   x = nf_to_scalar_or_basis(nf, x);
     812       46035 :   if (typ(x) != t_COL) {
     813       37520 :     w = Q_pvalrem(x,p, py);
     814       37520 :     if (!w) { *py = gerepilecopy(av, x); return 0; }
     815       13455 :     *py = gerepileupto(av, gmul(powp(nf, pr, w), *py));
     816       13455 :     return e*w;
     817             :   }
     818        8515 :   x = Q_primitive_part(x, &cx);
     819        8515 :   w = ZC_nfvalrem(x,pr, py);
     820        8515 :   if (cx)
     821             :   {
     822         835 :     long v = Q_pvalrem(cx,p, &t);
     823         835 :     *py = nfmul(nf, *py, gmul(powp(nf,pr,v), t));
     824         835 :     *py = gerepileupto(av, *py);
     825         835 :     w += e*v;
     826             :   }
     827             :   else
     828        7680 :     *py = gerepilecopy(av, *py);
     829        8515 :   return w;
     830             : }
     831             : GEN
     832         105 : gpnfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
     833             : {
     834         105 :   long v = nfvalrem(nf,x,pr,py);
     835         105 :   return v == LONG_MAX? mkoo(): stoi(v);
     836             : }
     837             : 
     838             : /* true nf */
     839             : GEN
     840       52705 : coltoalg(GEN nf, GEN x)
     841             : {
     842       52705 :   return mkpolmod( nf_to_scalar_or_alg(nf, x), nf_get_pol(nf) );
     843             : }
     844             : 
     845             : GEN
     846       54345 : basistoalg(GEN nf, GEN x)
     847             : {
     848             :   GEN z, T;
     849             : 
     850       54345 :   nf = checknf(nf);
     851       54345 :   switch(typ(x))
     852             :   {
     853             :     case t_COL: {
     854       47585 :       pari_sp av = avma;
     855       47585 :       return gerepilecopy(av, coltoalg(nf, x));
     856             :     }
     857             :     case t_POLMOD:
     858          85 :       T = nf_get_pol(nf);
     859          85 :       if (!RgX_equal_var(T,gel(x,1)))
     860           0 :         pari_err_MODULUS("basistoalg", T,gel(x,1));
     861          85 :       return gcopy(x);
     862             :     case t_POL:
     863         575 :       T = nf_get_pol(nf);
     864         575 :       if (varn(T) != varn(x)) pari_err_VAR("basistoalg",x,T);
     865         575 :       z = cgetg(3,t_POLMOD);
     866         575 :       gel(z,1) = ZX_copy(T);
     867         575 :       gel(z,2) = RgX_rem(x, T); return z;
     868             :     case t_INT:
     869             :     case t_FRAC:
     870        6100 :       T = nf_get_pol(nf);
     871        6100 :       z = cgetg(3,t_POLMOD);
     872        6100 :       gel(z,1) = ZX_copy(T);
     873        6100 :       gel(z,2) = gcopy(x); return z;
     874             :     default:
     875           0 :       pari_err_TYPE("basistoalg",x);
     876             :       return NULL; /* LCOV_EXCL_LINE */
     877             :   }
     878             : }
     879             : 
     880             : /* Assume nf is a genuine nf. */
     881             : GEN
     882    11294395 : nf_to_scalar_or_basis(GEN nf, GEN x)
     883             : {
     884    11294395 :   switch(typ(x))
     885             :   {
     886             :     case t_INT: case t_FRAC:
     887     7133060 :       return x;
     888             :     case t_POLMOD:
     889      149440 :       x = checknfelt_mod(nf,x,"nf_to_scalar_or_basis");
     890      149390 :       if (typ(x) != t_POL) return x;
     891             :       /* fall through */
     892             :     case t_POL:
     893             :     {
     894      244540 :       GEN T = nf_get_pol(nf);
     895      244540 :       long l = lg(x);
     896      244540 :       if (varn(x) != varn(T)) pari_err_VAR("nf_to_scalar_or_basis", x,T);
     897      244495 :       if (l >= lg(T)) { x = RgX_rem(x, T); l = lg(x); }
     898      244495 :       if (l == 2) return gen_0;
     899      210640 :       if (l == 3) return gel(x,2);
     900      182065 :       return poltobasis(nf,x);
     901             :     }
     902             :     case t_COL:
     903     3892305 :       if (lg(x)-1 != nf_get_degree(nf)) break;
     904     3892260 :       return QV_isscalar(x)? gel(x,1): x;
     905             :   }
     906          50 :   pari_err_TYPE("nf_to_scalar_or_basis",x);
     907             :   return NULL; /* LCOV_EXCL_LINE */
     908             : }
     909             : /* Let x be a polynomial with coefficients in Q or nf. Return the same
     910             :  * polynomial with coefficients expressed as vectors (on the integral basis).
     911             :  * No consistency checks, not memory-clean. */
     912             : GEN
     913        1790 : RgX_to_nfX(GEN nf, GEN x)
     914             : {
     915             :   long i, l;
     916        1790 :   GEN y = cgetg_copy(x, &l); y[1] = x[1];
     917        1790 :   for (i=2; i<l; i++) gel(y,i) = nf_to_scalar_or_basis(nf, gel(x,i));
     918        1790 :   return y;
     919             : }
     920             : 
     921             : /* Assume nf is a genuine nf. */
     922             : GEN
     923      106675 : nf_to_scalar_or_alg(GEN nf, GEN x)
     924             : {
     925      106675 :   switch(typ(x))
     926             :   {
     927             :     case t_INT: case t_FRAC:
     928        5095 :       return x;
     929             :     case t_POLMOD:
     930         950 :       x = checknfelt_mod(nf,x,"nf_to_scalar_or_alg");
     931         950 :       if (typ(x) != t_POL) return x;
     932             :       /* fall through */
     933             :     case t_POL:
     934             :     {
     935        9655 :       GEN T = nf_get_pol(nf);
     936        9655 :       long l = lg(x);
     937        9655 :       if (varn(x) != varn(T)) pari_err_VAR("nf_to_scalar_or_alg", x,T);
     938        9655 :       if (l >= lg(T)) { x = RgX_rem(x, T); l = lg(x); }
     939        9655 :       if (l == 2) return gen_0;
     940        9655 :       if (l == 3) return gel(x,2);
     941        9500 :       return x;
     942             :     }
     943             :     case t_COL:
     944             :     {
     945             :       GEN dx;
     946       91885 :       if (lg(x)-1 != nf_get_degree(nf)) break;
     947      183770 :       if (QV_isscalar(x)) return gel(x,1);
     948       72035 :       x = Q_remove_denom(x, &dx);
     949       72035 :       x = RgV_RgC_mul(nf_get_zkprimpart(nf), x);
     950       72035 :       dx = mul_denom(dx, nf_get_zkden(nf));
     951       72035 :       return gdiv(x,dx);
     952             :     }
     953             :   }
     954          35 :   pari_err_TYPE("nf_to_scalar_or_alg",x);
     955             :   return NULL; /* LCOV_EXCL_LINE */
     956             : }
     957             : 
     958             : /* gmul(A, RgX_to_RgC(x)), A t_MAT of compatible dimensions */
     959             : GEN
     960         950 : RgM_RgX_mul(GEN A, GEN x)
     961             : {
     962         950 :   long i, l = lg(x)-1;
     963             :   GEN z;
     964         950 :   if (l == 1) return zerocol(nbrows(A));
     965         950 :   z = gmul(gel(x,2), gel(A,1));
     966        1805 :   for (i = 2; i < l; i++)
     967         855 :     if (!gequal0(gel(x,i+1))) z = gadd(z, gmul(gel(x,i+1), gel(A,i)));
     968         950 :   return z;
     969             : }
     970             : GEN
     971     1258510 : ZM_ZX_mul(GEN A, GEN x)
     972             : {
     973     1258510 :   long i, l = lg(x)-1;
     974             :   GEN z;
     975     1258510 :   if (l == 1) return zerocol(nbrows(A));
     976     1258170 :   z = ZC_Z_mul(gel(A,1), gel(x,2));
     977     5449466 :   for (i = 2; i < l ; i++)
     978     4191296 :     if (signe(gel(x,i+1))) z = ZC_add(z, ZC_Z_mul(gel(A,i), gel(x,i+1)));
     979     1258170 :   return z;
     980             : }
     981             : /* x a t_POL, nf a genuine nf. No garbage collecting. No check.  */
     982             : GEN
     983     1164715 : poltobasis(GEN nf, GEN x)
     984             : {
     985     1164715 :   GEN d, T = nf_get_pol(nf);
     986     1164715 :   if (varn(x) != varn(T)) pari_err_VAR( "poltobasis", x,T);
     987     1164675 :   if (degpol(x) >= degpol(T)) x = RgX_rem(x,T);
     988     1164675 :   x = Q_remove_denom(x, &d);
     989     1164660 :   if (!RgX_is_ZX(x)) pari_err_TYPE("poltobasis",x);
     990     1164660 :   x = ZM_ZX_mul(nf_get_invzk(nf), x);
     991     1164660 :   if (d) x = RgC_Rg_div(x, d);
     992     1164660 :   return x;
     993             : }
     994             : 
     995             : GEN
     996       54470 : algtobasis(GEN nf, GEN x)
     997             : {
     998             :   pari_sp av;
     999             : 
    1000       54470 :   nf = checknf(nf);
    1001       54470 :   switch(typ(x))
    1002             :   {
    1003             :     case t_POLMOD:
    1004       27220 :       if (!RgX_equal_var(nf_get_pol(nf),gel(x,1)))
    1005           5 :         pari_err_MODULUS("algtobasis", nf_get_pol(nf),gel(x,1));
    1006       27215 :       x = gel(x,2);
    1007       27215 :       switch(typ(x))
    1008             :       {
    1009             :         case t_INT:
    1010        2985 :         case t_FRAC: return scalarcol(x, nf_get_degree(nf));
    1011             :         case t_POL:
    1012       24230 :           av = avma;
    1013       24230 :           return gerepileupto(av,poltobasis(nf,x));
    1014             :       }
    1015           0 :       break;
    1016             : 
    1017             :     case t_POL:
    1018       16240 :       av = avma;
    1019       16240 :       return gerepileupto(av,poltobasis(nf,x));
    1020             : 
    1021             :     case t_COL:
    1022        4645 :       if (lg(x)-1 != nf_get_degree(nf)) pari_err_DIM("nfalgtobasis");
    1023        4645 :       return gcopy(x);
    1024             : 
    1025             :     case t_INT:
    1026        6365 :     case t_FRAC: return scalarcol(x, nf_get_degree(nf));
    1027             :   }
    1028           0 :   pari_err_TYPE("algtobasis",x);
    1029             :   return NULL; /* LCOV_EXCL_LINE */
    1030             : }
    1031             : 
    1032             : GEN
    1033       25500 : rnfbasistoalg(GEN rnf,GEN x)
    1034             : {
    1035       25500 :   const char *f = "rnfbasistoalg";
    1036             :   long lx, i;
    1037       25500 :   pari_sp av = avma;
    1038             :   GEN z, nf, relpol, T;
    1039             : 
    1040       25500 :   checkrnf(rnf);
    1041       25500 :   nf = rnf_get_nf(rnf);
    1042       25500 :   T = nf_get_pol(nf);
    1043       25500 :   relpol = QXQX_to_mod_shallow(rnf_get_pol(rnf), T);
    1044       25500 :   switch(typ(x))
    1045             :   {
    1046             :     case t_COL:
    1047         570 :       z = cgetg_copy(x, &lx);
    1048        1670 :       for (i=1; i<lx; i++)
    1049             :       {
    1050        1135 :         GEN c = nf_to_scalar_or_alg(nf, gel(x,i));
    1051        1100 :         if (typ(c) == t_POL) c = mkpolmod(c,T);
    1052        1100 :         gel(z,i) = c;
    1053             :       }
    1054         535 :       z = RgV_RgC_mul(gel(rnf_get_zk(rnf),1), z);
    1055         490 :       return gerepileupto(av, gmodulo(z,relpol));
    1056             : 
    1057             :     case t_POLMOD:
    1058       16960 :       x = polmod_nffix(f, rnf, x, 0);
    1059       16810 :       if (typ(x) != t_POL) break;
    1060        6835 :       retmkpolmod(RgX_copy(x), RgX_copy(relpol));
    1061             :     case t_POL:
    1062         590 :       if (varn(x) == varn(T)) { RgX_check_QX(x,f); x = gmodulo(x,T); break; }
    1063         430 :       if (varn(x) == varn(relpol))
    1064             :       {
    1065         395 :         x = RgX_nffix(f,nf_get_pol(nf),x,0);
    1066         395 :         return gmodulo(x, relpol);
    1067             :       }
    1068          35 :       pari_err_VAR(f, x,relpol);
    1069             :   }
    1070       17480 :   retmkpolmod(scalarpol(x, varn(relpol)), RgX_copy(relpol));
    1071             : }
    1072             : 
    1073             : GEN
    1074         600 : matbasistoalg(GEN nf,GEN x)
    1075             : {
    1076             :   long i, j, li, lx;
    1077         600 :   GEN z = cgetg_copy(x, &lx);
    1078             : 
    1079         600 :   if (lx == 1) return z;
    1080         595 :   switch(typ(x))
    1081             :   {
    1082             :     case t_VEC: case t_COL:
    1083          20 :       for (i=1; i<lx; i++) gel(z,i) = basistoalg(nf, gel(x,i));
    1084          20 :       return z;
    1085         575 :     case t_MAT: break;
    1086           0 :     default: pari_err_TYPE("matbasistoalg",x);
    1087             :   }
    1088         575 :   li = lgcols(x);
    1089        2110 :   for (j=1; j<lx; j++)
    1090             :   {
    1091        1535 :     GEN c = cgetg(li,t_COL), xj = gel(x,j);
    1092        1535 :     gel(z,j) = c;
    1093        1535 :     for (i=1; i<li; i++) gel(c,i) = basistoalg(nf, gel(xj,i));
    1094             :   }
    1095         575 :   return z;
    1096             : }
    1097             : 
    1098             : GEN
    1099        1770 : matalgtobasis(GEN nf,GEN x)
    1100             : {
    1101             :   long i, j, li, lx;
    1102        1770 :   GEN z = cgetg_copy(x, &lx);
    1103             : 
    1104        1770 :   if (lx == 1) return z;
    1105        1730 :   switch(typ(x))
    1106             :   {
    1107             :     case t_VEC: case t_COL:
    1108        1725 :       for (i=1; i<lx; i++) gel(z,i) = algtobasis(nf, gel(x,i));
    1109        1725 :       return z;
    1110           5 :     case t_MAT: break;
    1111           0 :     default: pari_err_TYPE("matalgtobasis",x);
    1112             :   }
    1113           5 :   li = lgcols(x);
    1114          10 :   for (j=1; j<lx; j++)
    1115             :   {
    1116           5 :     GEN c = cgetg(li,t_COL), xj = gel(x,j);
    1117           5 :     gel(z,j) = c;
    1118           5 :     for (i=1; i<li; i++) gel(c,i) = algtobasis(nf, gel(xj,i));
    1119             :   }
    1120           5 :   return z;
    1121             : }
    1122             : GEN
    1123        2245 : RgM_to_nfM(GEN nf,GEN x)
    1124             : {
    1125             :   long i, j, li, lx;
    1126        2245 :   GEN z = cgetg_copy(x, &lx);
    1127             : 
    1128        2245 :   if (lx == 1) return z;
    1129        2245 :   li = lgcols(x);
    1130       16245 :   for (j=1; j<lx; j++)
    1131             :   {
    1132       14000 :     GEN c = cgetg(li,t_COL), xj = gel(x,j);
    1133       14000 :     gel(z,j) = c;
    1134       14000 :     for (i=1; i<li; i++) gel(c,i) = nf_to_scalar_or_basis(nf, gel(xj,i));
    1135             :   }
    1136        2245 :   return z;
    1137             : }
    1138             : GEN
    1139       26150 : RgC_to_nfC(GEN nf,GEN x)
    1140             : {
    1141       26150 :   long i, lx = lg(x);
    1142       26150 :   GEN z = cgetg(lx, t_COL);
    1143       26150 :   for (i=1; i<lx; i++) gel(z,i) = nf_to_scalar_or_basis(nf, gel(x,i));
    1144       26150 :   return z;
    1145             : }
    1146             : 
    1147             : /* x a t_POLMOD, supposedly in rnf = K[z]/(T), K = Q[y]/(Tnf) */
    1148             : GEN
    1149       43355 : polmod_nffix(const char *f, GEN rnf, GEN x, int lift)
    1150       43355 : { return polmod_nffix2(f, rnf_get_nfpol(rnf), rnf_get_pol(rnf), x,lift); }
    1151             : GEN
    1152       43420 : polmod_nffix2(const char *f, GEN T, GEN relpol, GEN x, int lift)
    1153             : {
    1154       43420 :   if (RgX_equal_var(gel(x,1),relpol))
    1155             :   {
    1156       39630 :     x = gel(x,2);
    1157       39630 :     if (typ(x) == t_POL && varn(x) == varn(relpol))
    1158             :     {
    1159       20380 :       x = RgX_nffix(f, T, x, lift);
    1160       20380 :       switch(lg(x))
    1161             :       {
    1162         210 :         case 2: return gen_0;
    1163        2350 :         case 3: return gel(x,2);
    1164             :       }
    1165       17820 :       return x;
    1166             :     }
    1167             :   }
    1168       23040 :   return Rg_nffix(f, T, x, lift);
    1169             : }
    1170             : GEN
    1171         840 : rnfalgtobasis(GEN rnf,GEN x)
    1172             : {
    1173         840 :   const char *f = "rnfalgtobasis";
    1174         840 :   pari_sp av = avma;
    1175             :   GEN T, relpol;
    1176             : 
    1177         840 :   checkrnf(rnf);
    1178         840 :   relpol = rnf_get_pol(rnf);
    1179         840 :   T = rnf_get_nfpol(rnf);
    1180         840 :   switch(typ(x))
    1181             :   {
    1182             :     case t_COL:
    1183          35 :       if (lg(x)-1 != rnf_get_degree(rnf)) pari_err_DIM(f);
    1184          20 :       x = RgV_nffix(f, T, x, 0);
    1185          15 :       return gerepilecopy(av, x);
    1186             : 
    1187             :     case t_POLMOD:
    1188         745 :       x = polmod_nffix(f, rnf, x, 0);
    1189         715 :       if (typ(x) != t_POL) break;
    1190         505 :       return gerepileupto(av, RgM_RgX_mul(rnf_get_invzk(rnf), x));
    1191             :     case t_POL:
    1192          40 :       if (varn(x) == varn(T)) { RgX_check_QX(x,f); x = mkpolmod(x,T); break; }
    1193          25 :       x = RgX_nffix(f, T, x, 0);
    1194          20 :       if (degpol(x) >= degpol(relpol)) x = RgX_rem(x,relpol);
    1195          20 :       return gerepileupto(av, RgM_RgX_mul(rnf_get_invzk(rnf), x));
    1196             :   }
    1197         240 :   return gerepileupto(av, scalarcol(x, rnf_get_degree(rnf)));
    1198             : }
    1199             : 
    1200             : /* Given a and b in nf, gives an algebraic integer y in nf such that a-b.y
    1201             :  * is "small" */
    1202             : GEN
    1203         185 : nfdiveuc(GEN nf, GEN a, GEN b)
    1204             : {
    1205         185 :   pari_sp av = avma;
    1206         185 :   a = nfdiv(nf,a,b);
    1207         185 :   return gerepileupto(av, ground(a));
    1208             : }
    1209             : 
    1210             : /* Given a and b in nf, gives a "small" algebraic integer r in nf
    1211             :  * of the form a-b.y */
    1212             : GEN
    1213         185 : nfmod(GEN nf, GEN a, GEN b)
    1214             : {
    1215         185 :   pari_sp av = avma;
    1216         185 :   GEN p1 = gneg_i(nfmul(nf,b,ground(nfdiv(nf,a,b))));
    1217         185 :   return gerepileupto(av, nfadd(nf,a,p1));
    1218             : }
    1219             : 
    1220             : /* Given a and b in nf, gives a two-component vector [y,r] in nf such
    1221             :  * that r=a-b.y is "small". */
    1222             : GEN
    1223         185 : nfdivrem(GEN nf, GEN a, GEN b)
    1224             : {
    1225         185 :   pari_sp av = avma;
    1226         185 :   GEN p1,z, y = ground(nfdiv(nf,a,b));
    1227             : 
    1228         185 :   p1 = gneg_i(nfmul(nf,b,y));
    1229         185 :   z = cgetg(3,t_VEC);
    1230         185 :   gel(z,1) = gcopy(y);
    1231         185 :   gel(z,2) = nfadd(nf,a,p1); return gerepileupto(av, z);
    1232             : }
    1233             : 
    1234             : /*************************************************************************/
    1235             : /**                                                                     **/
    1236             : /**                        REAL EMBEDDINGS                              **/
    1237             : /**                                                                     **/
    1238             : /*************************************************************************/
    1239             : static GEN
    1240        8590 : sarch_get_cyc(GEN sarch) { return gel(sarch,1); }
    1241             : static GEN
    1242      146830 : sarch_get_archp(GEN sarch) { return gel(sarch,2); }
    1243             : static GEN
    1244       32135 : sarch_get_MI(GEN sarch) { return gel(sarch,3); }
    1245             : static GEN
    1246       32135 : sarch_get_lambda(GEN sarch) { return gel(sarch,4); }
    1247             : static GEN
    1248       32135 : sarch_get_F(GEN sarch) { return gel(sarch,5); }
    1249             : 
    1250             : /* true nf, return number of positive roots of char_x */
    1251             : static long
    1252        1115 : num_positive(GEN nf, GEN x)
    1253             : {
    1254        1115 :   GEN T = nf_get_pol(nf);
    1255        1115 :   GEN charx = ZXQ_charpoly(nf_to_scalar_or_alg(nf,x), T, 0);
    1256             :   long np;
    1257        1115 :   charx = ZX_radical(charx);
    1258        1115 :   np = ZX_sturmpart(charx, mkvec2(gen_0,mkoo()));
    1259        1115 :   return np * (degpol(T) / degpol(charx));
    1260             : }
    1261             : 
    1262             : /* x a QC: return sigma_k(x) where 1 <= k <= r1+r2; correct but inefficient
    1263             :  * if x in Q. M = nf_get_M(nf) */
    1264             : static GEN
    1265         420 : nfembed_i(GEN M, GEN x, long k)
    1266             : {
    1267         420 :   long i, l = lg(M);
    1268         420 :   GEN z = gel(x,1);
    1269         420 :   for (i = 2; i < l; i++) z = gadd(z, gmul(gcoeff(M,k,i), gel(x,i)));
    1270         420 :   return z;
    1271             : }
    1272             : GEN
    1273        1135 : nfembed(GEN nf, GEN x, long k)
    1274             : {
    1275        1135 :   pari_sp av = avma;
    1276        1135 :   nf = checknf(nf);
    1277        1135 :   x = nf_to_scalar_or_basis(nf,x);
    1278        1135 :   if (typ(x) != t_COL) return gerepilecopy(av, x);
    1279           0 :   return gerepileupto(av, nfembed_i(nf_get_M(nf),x,k));
    1280             : }
    1281             : 
    1282             : /* x a ZC */
    1283             : static GEN
    1284      274385 : zk_embed(GEN M, GEN x, long k)
    1285             : {
    1286      274385 :   long i, l = lg(x);
    1287      274385 :   GEN z = gel(x,1); /* times M[k,1], which is 1 */
    1288      274385 :   for (i = 2; i < l; i++) z = mpadd(z, mpmul(gcoeff(M,k,i), gel(x,i)));
    1289      274385 :   return z;
    1290             : }
    1291             : 
    1292             : /* Given floating point approximation z of sigma_k(x), decide its sign
    1293             :  * [0/+, 1/- and -1 for FAIL] */
    1294             : static long
    1295      266320 : eval_sign_embed(GEN z)
    1296             : { /* dubious, fail */
    1297      266320 :   if (typ(z) == t_REAL && realprec(z) <= LOWDEFAULTPREC) return -1;
    1298      265590 :   return (signe(z) < 1)? 1: 0;
    1299             : }
    1300             : /* return v such that (-1)^v = sign(sigma_k(x)), x primitive ZC */
    1301             : static long
    1302      217625 : eval_sign(GEN M, GEN x, long k)
    1303      217625 : { return eval_sign_embed( zk_embed(M, x, k) ); }
    1304             : 
    1305             : /* check that signs[i..#signs] == s; signs = NULL encodes "totally positive" */
    1306             : static int
    1307           0 : oksigns(long l, GEN signs, long i, long s)
    1308             : {
    1309           0 :   if (!signs) return s == 0;
    1310           0 :   for (; i < l; i++)
    1311           0 :     if (signs[i] != s) return 0;
    1312           0 :   return 1;
    1313             : }
    1314             : /* check that signs[i] = s and signs[i+1..#signs] = 1-s */
    1315             : static int
    1316           0 : oksigns2(long l, GEN signs, long i, long s)
    1317             : {
    1318           0 :   if (!signs) return s == 0 && i == l-1;
    1319           0 :   return signs[i] == s && oksigns(l, signs, i+1, 1-s);
    1320             : }
    1321             : 
    1322             : /* true nf, x a ZC (primitive for efficiency), embx its embeddings or NULL */
    1323             : static int
    1324       42595 : nfchecksigns_i(GEN nf, GEN x, GEN embx, GEN signs, GEN archp)
    1325             : {
    1326       42595 :   long l = lg(archp), i;
    1327       42595 :   GEN M = nf_get_M(nf), sarch = NULL;
    1328       42595 :   long np = -1;
    1329       62485 :   for (i = 1; i < l; i++)
    1330             :   {
    1331             :     long s;
    1332       48830 :     if (embx)
    1333       48695 :       s = eval_sign_embed(gel(embx,i));
    1334             :     else
    1335         135 :       s = eval_sign(M, x, archp[i]);
    1336             :     /* 0 / + or 1 / -; -1 for FAIL */
    1337       48830 :     if (s < 0) /* failure */
    1338             :     {
    1339           0 :       long ni, r1 = nf_get_r1(nf);
    1340             :       GEN xi;
    1341           0 :       if (np < 0)
    1342             :       {
    1343           0 :         np = num_positive(nf, x);
    1344           0 :         if (np == 0)  return oksigns(l, signs, i, 1);
    1345           0 :         if (np == r1) return oksigns(l, signs, i, 0);
    1346           0 :         sarch = nfarchstar(nf, NULL, identity_perm(r1));
    1347             :       }
    1348           0 :       xi = set_sign_mod_divisor(nf, vecsmall_ei(r1, archp[i]), gen_1, sarch);
    1349           0 :       xi = Q_primpart(xi);
    1350           0 :       ni = num_positive(nf, nfmuli(nf,x,xi));
    1351           0 :       if (ni == 0)  return oksigns2(l, signs, i, 0);
    1352           0 :       if (ni == r1) return oksigns2(l, signs, i, 1);
    1353           0 :       s = ni < np? 0: 1;
    1354             :     }
    1355       48830 :     if (s != (signs? signs[i]: 0)) return 0;
    1356             :   }
    1357       13655 :   return 1;
    1358             : }
    1359             : static void
    1360         115 : pl_convert(GEN pl, GEN *psigns, GEN *parchp)
    1361             : {
    1362         115 :   long i, j, l = lg(pl);
    1363         115 :   GEN signs = cgetg(l, t_VECSMALL);
    1364         115 :   GEN archp = cgetg(l, t_VECSMALL);
    1365         290 :   for (i = j = 1; i < l; i++)
    1366             :   {
    1367         175 :     if (!pl[i]) continue;
    1368         165 :     archp[j] = i;
    1369         165 :     signs[j] = (pl[i] < 0)? 1: 0;
    1370         165 :     j++;
    1371             :   }
    1372         115 :   setlg(archp, j); *parchp = archp;
    1373         115 :   setlg(signs, j); *psigns = signs;
    1374         115 : }
    1375             : /* pl : requested signs for real embeddings, 0 = no sign constraint */
    1376             : int
    1377         365 : nfchecksigns(GEN nf, GEN x, GEN pl)
    1378             : {
    1379         365 :   pari_sp av = avma;
    1380             :   GEN signs, archp;
    1381             :   int res;
    1382         365 :   nf = checknf(nf);
    1383         365 :   x = nf_to_scalar_or_basis(nf,x);
    1384         365 :   if (typ(x) != t_COL)
    1385             :   {
    1386         250 :     long i, l = lg(pl), s = gsigne(x);
    1387         555 :     for (i = 1; i < l; i++)
    1388         305 :       if (pl[i] && pl[i] != s) { avma = av; return 0; }
    1389         250 :     avma = av; return 1;
    1390             :   }
    1391         115 :   pl_convert(pl, &signs, &archp);
    1392         115 :   res = nfchecksigns_i(nf, x, NULL, signs, archp);
    1393         115 :   avma = av; return res;
    1394             : }
    1395             : 
    1396             : /* signs = NULL: totally positive, else sign[i] = 0 (+) or 1 (-) */
    1397             : static GEN
    1398       32135 : get_C(GEN lambda, long l, GEN signs)
    1399             : {
    1400             :   long i;
    1401             :   GEN C, mlambda;
    1402       32135 :   if (!signs) return const_vec(l-1, lambda);
    1403        3950 :   C = cgetg(l, t_COL); mlambda = gneg(lambda);
    1404        3950 :   for (i = 1; i < l; i++) gel(C,i) = signs[i]? mlambda: lambda;
    1405        3950 :   return C;
    1406             : }
    1407             : /* signs = NULL: totally positive at archp */
    1408             : static GEN
    1409       49515 : nfsetsigns(GEN nf, GEN signs, GEN x, GEN sarch)
    1410             : {
    1411       49515 :   long i, l = lg(sarch_get_archp(sarch));
    1412             :   GEN ex;
    1413             :   /* Is signature already correct ? */
    1414       49515 :   if (typ(x) != t_COL)
    1415             :   {
    1416        7035 :     long s = gsigne(x) < 0? 1: 0;
    1417        7035 :     if (!signs)
    1418        1695 :       i = (s == 1)? 1: l;
    1419             :     else
    1420             :     {
    1421        8595 :       for (i = 1; i < l; i++)
    1422        6470 :         if (signs[i] != s) break;
    1423             :     }
    1424        7035 :     ex = (i < l)? const_col(l-1, x): NULL;
    1425             :   }
    1426             :   else
    1427             :   {
    1428       42480 :     pari_sp av = avma;
    1429       42480 :     GEN cex, M = nf_get_M(nf), archp = sarch_get_archp(sarch);
    1430       42480 :     GEN xp = Q_primitive_part(x,&cex);
    1431       42480 :     ex = cgetg(l,t_COL);
    1432       42480 :     for (i = 1; i < l; i++) gel(ex,i) = zk_embed(M,xp,archp[i]);
    1433       42480 :     if (nfchecksigns_i(nf, xp, ex, signs, archp)) { ex = NULL; avma = av; }
    1434       28910 :     else if (cex) ex = RgC_Rg_mul(ex, cex); /* put back content */
    1435             :   }
    1436       49515 :   if (ex)
    1437             :   { /* If no, fix it */
    1438       32135 :     GEN lambda = sarch_get_lambda(sarch);
    1439       32135 :     GEN MI = sarch_get_MI(sarch);
    1440       32135 :     GEN F = sarch_get_F(sarch);
    1441       32135 :     GEN t = RgC_sub(get_C(lambda, l, signs), ex);
    1442             :     long e;
    1443       32135 :     t = grndtoi(RgM_RgC_mul(MI,t), &e);
    1444       32135 :     if (lg(F) != 1) t = ZM_ZC_mul(F, t);
    1445       32135 :     x = typ(x) == t_COL? RgC_add(t, x): RgC_Rg_add(t, x);
    1446             :   }
    1447       49515 :   return x;
    1448             : }
    1449             : /* - sarch = nfarchstar(nf, F);
    1450             :  * - x encodes a vector of signs at arch.archp: either a t_VECSMALL
    1451             :  *   (vector of signs as {0,1}-vector), NULL (totally positive at archp),
    1452             :  *   or a non-zero number field element (replaced by its signature at archp);
    1453             :  * - y is a non-zero number field element
    1454             :  * Return z = y (mod F) with signs(y, archp) = signs(x) (a {0,1}-vector) */
    1455             : GEN
    1456       54835 : set_sign_mod_divisor(GEN nf, GEN x, GEN y, GEN sarch)
    1457             : {
    1458       54835 :   GEN archp = sarch_get_archp(sarch);
    1459       54835 :   if (lg(archp) == 1) return y;
    1460       48605 :   nf = checknf(nf);
    1461       48605 :   if (x && typ(x) != t_VECSMALL) x = nfsign_arch(nf, x, archp);
    1462       48605 :   y = nf_to_scalar_or_basis(nf,y);
    1463       48605 :   return nfsetsigns(nf, x, y, sarch);
    1464             : }
    1465             : 
    1466             : static GEN
    1467        4100 : setsigns_init(GEN nf, GEN archp, GEN F, GEN DATA)
    1468             : {
    1469        4100 :   GEN lambda, Mr = rowpermute(nf_get_M(nf), archp), MI = F? RgM_mul(Mr,F): Mr;
    1470        4100 :   lambda = gmul2n(matrixnorm(MI,DEFAULTPREC), -1);
    1471        4100 :   if (lg(archp) < lg(MI))
    1472             :   {
    1473        2915 :     GEN perm = gel(indexrank(MI), 2);
    1474        2915 :     if (!F) F = matid(nf_get_degree(nf));
    1475        2915 :     MI = vecpermute(MI, perm);
    1476        2915 :     F = vecpermute(F, perm);
    1477             :   }
    1478        4100 :   if (!F) F = cgetg(1,t_MAT);
    1479        4100 :   MI = RgM_inv(MI);
    1480        4100 :   return mkvec5(DATA, archp, MI, lambda, F);
    1481             : }
    1482             : /* F non-0 integral ideal in HNF (or NULL: Z_K), compute elements in 1+F
    1483             :  * whose sign matrix at archp is identity; archp in 'indices' format */
    1484             : GEN
    1485        5605 : nfarchstar(GEN nf, GEN F, GEN archp)
    1486             : {
    1487        5605 :   long nba = lg(archp) - 1;
    1488        5605 :   if (!nba) return mkvec2(cgetg(1,t_VEC), archp);
    1489        3195 :   if (F && equali1(gcoeff(F,1,1))) F = NULL;
    1490        3195 :   if (F) F = idealpseudored(F, nf_get_roundG(nf));
    1491        3195 :   return setsigns_init(nf, archp, F, const_vec(nba, gen_2));
    1492             : }
    1493             : 
    1494             : /*************************************************************************/
    1495             : /**                                                                     **/
    1496             : /**                         IDEALCHINESE                                **/
    1497             : /**                                                                     **/
    1498             : /*************************************************************************/
    1499             : static int
    1500        1445 : isprfact(GEN x)
    1501             : {
    1502             :   long i, l;
    1503             :   GEN L, E;
    1504        1445 :   if (typ(x) != t_MAT || lg(x) != 3) return 0;
    1505        1445 :   L = gel(x,1); l = lg(L);
    1506        1445 :   E = gel(x,2);
    1507        3420 :   for(i=1; i<l; i++)
    1508             :   {
    1509        1975 :     checkprid(gel(L,i));
    1510        1975 :     if (typ(gel(E,i)) != t_INT) return 0;
    1511             :   }
    1512        1445 :   return 1;
    1513             : }
    1514             : 
    1515             : /* initialize projectors mod pr[i]^e[i] for idealchinese */
    1516             : static GEN
    1517        1445 : pr_init(GEN nf, GEN fa, GEN w, GEN dw)
    1518             : {
    1519        1445 :   GEN U, E, F, L = gel(fa,1), E0 = gel(fa,2);
    1520        1445 :   long i, r = lg(L);
    1521             : 
    1522        1445 :   if (w && lg(w) != r) pari_err_TYPE("idealchinese", w);
    1523        1445 :   if (r == 1 && !dw) return cgetg(1,t_VEC);
    1524        1440 :   E = leafcopy(E0); /* do not destroy fa[2] */
    1525        3415 :   for (i = 1; i < r; i++)
    1526        1975 :     if (signe(gel(E,i)) < 0) gel(E,i) = gen_0;
    1527        1440 :   F = factorbackprime(nf, L, E);
    1528        1440 :   if (dw)
    1529             :   {
    1530         490 :     F = ZM_Z_mul(F, dw);
    1531        1120 :     for (i = 1; i < r; i++)
    1532             :     {
    1533         630 :       GEN pr = gel(L,i);
    1534         630 :       long e = itos(gel(E0,i)), v = idealval(nf, dw, pr);
    1535         630 :       if (e >= 0)
    1536         625 :         gel(E,i) = addiu(gel(E,i), v);
    1537           5 :       else if (v + e <= 0)
    1538           0 :         F = idealmulpowprime(nf, F, pr, stoi(-v)); /* coprime to pr */
    1539             :       else
    1540             :       {
    1541           5 :         F = idealmulpowprime(nf, F, pr, stoi(e));
    1542           5 :         gel(E,i) = stoi(v + e);
    1543             :       }
    1544             :     }
    1545             :   }
    1546        1440 :   U = cgetg(r, t_VEC);
    1547        3415 :   for (i = 1; i < r; i++)
    1548             :   {
    1549             :     GEN u;
    1550        1975 :     if (w && gequal0(gel(w,i))) u = gen_0; /* unused */
    1551             :     else
    1552             :     {
    1553        1935 :       GEN pr = gel(L,i), e = gel(E,i), t;
    1554        1935 :       t = idealdivpowprime(nf,F, pr, e);
    1555        1935 :       u = hnfmerge_get_1(t, idealpow(nf, pr, e));
    1556        1935 :       if (!u) pari_err_COPRIME("idealchinese", t,pr);
    1557             :     }
    1558        1975 :     gel(U,i) = u;
    1559             :   }
    1560        1440 :   F = idealpseudored(F, nf_get_roundG(nf));
    1561        1440 :   return mkvec2(F, U);
    1562             : }
    1563             : 
    1564             : static GEN
    1565         905 : pl_normalize(GEN nf, GEN pl)
    1566             : {
    1567         905 :   const char *fun = "idealchinese";
    1568         905 :   if (lg(pl)-1 != nf_get_r1(nf)) pari_err_TYPE(fun,pl);
    1569         905 :   switch(typ(pl))
    1570             :   {
    1571         485 :     case t_VEC: RgV_check_ZV(pl,fun); pl = ZV_to_zv(pl);
    1572             :       /* fall through */
    1573         905 :     case t_VECSMALL: break;
    1574           0 :     default: pari_err_TYPE(fun,pl);
    1575             :   }
    1576         905 :   return pl;
    1577             : }
    1578             : 
    1579             : static int
    1580        3435 : is_chineseinit(GEN x)
    1581             : {
    1582             :   GEN fa, pl;
    1583             :   long l;
    1584        3435 :   if (typ(x) != t_VEC || lg(x)!=3) return 0;
    1585        2555 :   fa = gel(x,1);
    1586        2555 :   pl = gel(x,2);
    1587        2555 :   if (typ(fa) != t_VEC || typ(pl) != t_VEC) return 0;
    1588        1090 :   l = lg(fa);
    1589        1090 :   if (l != 1)
    1590             :   {
    1591        1075 :     if (l != 3 || typ(gel(fa,1)) != t_MAT || typ(gel(fa,2)) != t_VEC)
    1592           5 :       return 0;
    1593             :   }
    1594        1085 :   l = lg(pl);
    1595        1085 :   if (l != 1)
    1596             :   {
    1597         345 :     if (l != 6 || typ(gel(pl,3)) != t_MAT || typ(gel(pl,1)) != t_VECSMALL
    1598         345 :                                           || typ(gel(pl,2)) != t_VECSMALL)
    1599           0 :       return 0;
    1600             :   }
    1601        1085 :   return 1;
    1602             : }
    1603             : 
    1604             : /* nf a true 'nf' */
    1605             : static GEN
    1606        1490 : chineseinit_i(GEN nf, GEN fa, GEN w, GEN dw)
    1607             : {
    1608        1490 :   const char *fun = "idealchineseinit";
    1609        1490 :   GEN archp = NULL, pl = NULL;
    1610        1490 :   switch(typ(fa))
    1611             :   {
    1612             :     case t_VEC:
    1613         905 :       if (is_chineseinit(fa))
    1614             :       {
    1615           0 :         if (dw) pari_err_DOMAIN(fun, "denom(y)", "!=", gen_1, w);
    1616           0 :         return fa;
    1617             :       }
    1618         905 :       if (lg(fa) != 3) pari_err_TYPE(fun, fa);
    1619             :       /* of the form [x,s] */
    1620         905 :       pl = pl_normalize(nf, gel(fa,2));
    1621         905 :       fa = gel(fa,1);
    1622         905 :       archp = vecsmall01_to_indices(pl);
    1623             :       /* keep pr_init, reset pl */
    1624         905 :       if (is_chineseinit(fa)) { fa = gel(fa,1); break; }
    1625             :       /* fall through */
    1626             :     case t_MAT: /* factorization? */
    1627        1445 :       if (isprfact(fa)) { fa = pr_init(nf, fa, w, dw); break; }
    1628           0 :     default: pari_err_TYPE(fun,fa);
    1629             :   }
    1630             : 
    1631        1490 :   if (pl)
    1632             :   {
    1633         905 :     GEN F = (lg(fa) == 1)? NULL: gel(fa,1);
    1634         905 :     long i, r = lg(archp);
    1635         905 :     GEN signs = cgetg(r, t_VECSMALL);
    1636         905 :     for (i = 1; i < r; i++) signs[i] = (pl[archp[i]] < 0)? 1: 0;
    1637         905 :     pl = setsigns_init(nf, archp, F, signs);
    1638             :   }
    1639             :   else
    1640         585 :     pl = cgetg(1,t_VEC);
    1641        1490 :   return mkvec2(fa, pl);
    1642             : }
    1643             : 
    1644             : /* Given a prime ideal factorization x, possibly with 0 or negative exponents,
    1645             :  * and a vector w of elements of nf, gives b such that
    1646             :  * v_p(b-w_p)>=v_p(x) for all prime ideals p in the ideal factorization
    1647             :  * and v_p(b)>=0 for all other p, using the standard proof given in GTM 138. */
    1648             : GEN
    1649        2530 : idealchinese(GEN nf, GEN x, GEN w)
    1650             : {
    1651        2530 :   const char *fun = "idealchinese";
    1652        2530 :   pari_sp av = avma;
    1653             :   GEN x1, x2, s, dw, F;
    1654             : 
    1655        2530 :   nf = checknf(nf);
    1656        2530 :   if (!w) return gerepilecopy(av, chineseinit_i(nf,x,NULL,NULL));
    1657             : 
    1658        1625 :   if (typ(w) != t_VEC) pari_err_TYPE(fun,w);
    1659        1625 :   w = Q_remove_denom(matalgtobasis(nf,w), &dw);
    1660        1625 :   if (!is_chineseinit(x)) x = chineseinit_i(nf,x,w,dw);
    1661             :   /* x is a 'chineseinit' */
    1662        1625 :   x1 = gel(x,1); s = NULL;
    1663        1625 :   if (lg(x1) == 1) F = NULL;
    1664             :   else
    1665             :   {
    1666        1610 :     GEN  U = gel(x1,2);
    1667        1610 :     long i, r = lg(w);
    1668        1610 :     F = gel(x1,1);
    1669        3945 :     for (i=1; i<r; i++)
    1670        2335 :       if (!gequal0(gel(w,i)))
    1671             :       {
    1672        1960 :         GEN t = nfmuli(nf, gel(U,i), gel(w,i));
    1673        1960 :         s = s? ZC_add(s,t): t;
    1674             :       }
    1675        1610 :     if (s) s = ZC_reducemodmatrix(s, F);
    1676             :   }
    1677        1625 :   if (!s) { s = zerocol(nf_get_degree(nf)); dw = NULL; }
    1678             : 
    1679        1625 :   x2 = gel(x,2);
    1680        1625 :   if (lg(x2) != 1) s = nfsetsigns(nf, gel(x2,1), s, x2);
    1681        1625 :   if (dw) s = RgC_Rg_div(s,dw);
    1682        1625 :   return gerepileupto(av, s);
    1683             : }
    1684             : 
    1685             : /*************************************************************************/
    1686             : /**                                                                     **/
    1687             : /**                           (Z_K/I)^*                                 **/
    1688             : /**                                                                     **/
    1689             : /*************************************************************************/
    1690             : GEN
    1691         905 : vecsmall01_to_indices(GEN v)
    1692             : {
    1693         905 :   long i, k, l = lg(v);
    1694         905 :   GEN p = new_chunk(l) + l;
    1695        2630 :   for (k=1, i=l-1; i; i--)
    1696        1725 :     if (v[i]) { *--p = i; k++; }
    1697         905 :   *--p = evallg(k) | evaltyp(t_VECSMALL);
    1698         905 :   avma = (pari_sp)p; return p;
    1699             : }
    1700             : GEN
    1701      229740 : vec01_to_indices(GEN v)
    1702             : {
    1703             :   long i, k, l;
    1704             :   GEN p;
    1705             : 
    1706      229740 :   switch (typ(v))
    1707             :   {
    1708      225425 :    case t_VECSMALL: return v;
    1709        4315 :    case t_VEC: break;
    1710           0 :    default: pari_err_TYPE("vec01_to_indices",v);
    1711             :   }
    1712        4315 :   l = lg(v);
    1713        4315 :   p = new_chunk(l) + l;
    1714       11795 :   for (k=1, i=l-1; i; i--)
    1715        7480 :     if (signe(gel(v,i))) { *--p = i; k++; }
    1716        4315 :   *--p = evallg(k) | evaltyp(t_VECSMALL);
    1717        4315 :   avma = (pari_sp)p; return p;
    1718             : }
    1719             : GEN
    1720        3435 : indices_to_vec01(GEN p, long r)
    1721             : {
    1722        3435 :   long i, l = lg(p);
    1723        3435 :   GEN v = zerovec(r);
    1724        3435 :   for (i = 1; i < l; i++) gel(v, p[i]) = gen_1;
    1725        3435 :   return v;
    1726             : }
    1727             : 
    1728             : /* return (column) vector of R1 signatures of x (0 or 1) */
    1729             : GEN
    1730      225425 : nfsign_arch(GEN nf, GEN x, GEN arch)
    1731             : {
    1732      225425 :   GEN sarch, M, V, archp = vec01_to_indices(arch);
    1733      225425 :   long i, s, np, n = lg(archp)-1;
    1734             :   pari_sp av;
    1735             : 
    1736      225425 :   if (!n) return cgetg(1,t_VECSMALL);
    1737      225365 :   nf = checknf(nf);
    1738      225365 :   if (typ(x) == t_MAT)
    1739             :   { /* factorisation */
    1740       65690 :     GEN g = gel(x,1), e = gel(x,2);
    1741       65690 :     V = zero_zv(n);
    1742      196125 :     for (i=1; i<lg(g); i++)
    1743      130435 :       if (mpodd(gel(e,i)))
    1744      112490 :         Flv_add_inplace(V, nfsign_arch(nf,gel(g,i),archp), 2);
    1745       65690 :     avma = (pari_sp)V; return V;
    1746             :   }
    1747      159675 :   av = avma; V = cgetg(n+1,t_VECSMALL);
    1748      159675 :   x = nf_to_scalar_or_basis(nf, x);
    1749      159675 :   switch(typ(x))
    1750             :   {
    1751             :     case t_INT:
    1752       33080 :       s = signe(x);
    1753       33080 :       if (!s) pari_err_DOMAIN("nfsign_arch","element","=",gen_0,x);
    1754       33080 :       avma = av; return const_vecsmall(n, (s < 0)? 1: 0);
    1755             :     case t_FRAC:
    1756           0 :       s = signe(gel(x,1));
    1757           0 :       avma = av; return const_vecsmall(n, (s < 0)? 1: 0);
    1758             :   }
    1759      126595 :   x = Q_primpart(x); M = nf_get_M(nf); sarch = NULL; np = -1;
    1760      343355 :   for (i = 1; i <= n; i++)
    1761             :   {
    1762      217490 :     long s = eval_sign(M, x, archp[i]);
    1763      217490 :     if (s < 0) /* failure */
    1764             :     {
    1765         730 :       long ni, r1 = nf_get_r1(nf);
    1766             :       GEN xi;
    1767         730 :       if (np < 0)
    1768             :       {
    1769         730 :         np = num_positive(nf, x);
    1770         730 :         if (np == 0) { avma = av; return const_vecsmall(n, 1); }
    1771         615 :         if (np == r1){ avma = av; return const_vecsmall(n, 0); }
    1772         385 :         sarch = nfarchstar(nf, NULL, identity_perm(r1));
    1773             :       }
    1774         385 :       xi = set_sign_mod_divisor(nf, vecsmall_ei(r1, archp[i]), gen_1, sarch);
    1775         385 :       xi = Q_primpart(xi);
    1776         385 :       ni = num_positive(nf, nfmuli(nf,x,xi));
    1777         385 :       if (ni == 0) { avma = av; V = const_vecsmall(n, 1); V[i] = 0; return V; }
    1778         270 :       if (ni == r1){ avma = av; V = const_vecsmall(n, 0); V[i] = 1; return V; }
    1779           0 :       s = ni < np? 0: 1;
    1780             :     }
    1781      216760 :     V[i] = s;
    1782             :   }
    1783      125865 :   avma = (pari_sp)V; return V;
    1784             : }
    1785             : static void
    1786        1800 : chk_ind(const char *s, long i, long r1)
    1787             : {
    1788        1800 :   if (i <= 0)
    1789           5 :     pari_err_DOMAIN(s, "index", "<=", gen_0, stoi(i));
    1790        1795 :   if (i > r1)
    1791          15 :     pari_err_DOMAIN(s, "index", ">", stoi(r1), stoi(i));
    1792        1780 : }
    1793             : GEN
    1794         550 : nfeltsign(GEN nf, GEN x, GEN ind0)
    1795             : {
    1796         550 :   pari_sp av = avma;
    1797             :   long i, l, r1;
    1798             :   GEN v, ind;
    1799         550 :   nf = checknf(nf);
    1800         550 :   r1 = nf_get_r1(nf);
    1801         550 :   x = nf_to_scalar_or_basis(nf, x);
    1802         550 :   if (!ind0) ind0 = identity_perm(r1);
    1803         550 :   switch(typ(ind0))
    1804             :   {
    1805             :     case t_INT: case t_VEC: case t_COL:
    1806          40 :       ind = gtovecsmall(ind0); break;
    1807             :     case t_VECSMALL:
    1808         510 :       ind = ind0; break;
    1809             :     default:
    1810           0 :       pari_err_TYPE("nfeltsign",ind0);
    1811             :       return NULL; /* LCOV_EXCL_LINE */
    1812             :   }
    1813         550 :   l = lg(ind);
    1814         550 :   for (i = 1; i < l; i++) chk_ind("nfeltsign", ind[i], r1);
    1815         535 :   if (typ(x) != t_COL)
    1816             :   {
    1817             :     GEN s;
    1818          15 :     switch(gsigne(x))
    1819             :     {
    1820           5 :       case -1:s = gen_m1; break;
    1821           5 :       case 1: s = gen_1; break;
    1822           5 :       default: s = gen_0; break;
    1823             :     }
    1824          15 :     avma = av;
    1825          15 :     return typ(ind0) == t_INT? s: const_vec(l-1, s);
    1826             :   }
    1827         520 :   v = nfsign_arch(nf, x, ind);
    1828         520 :   if (typ(ind0) == t_INT) { avma = av; return v[1]? gen_m1: gen_1; }
    1829         515 :   settyp(v, t_VEC);
    1830         515 :   for (i = 1; i < l; i++) gel(v,i) = v[i]? gen_m1: gen_1;
    1831         515 :   return gerepileupto(av, v);
    1832             : 
    1833             : }
    1834             : 
    1835             : GEN
    1836         390 : nfeltembed(GEN nf, GEN x, GEN ind0)
    1837             : {
    1838         390 :   pari_sp av = avma;
    1839             :   long i, l, r1, r2;
    1840             :   GEN v, ind, cx, M;
    1841         390 :   nf = checknf(nf);
    1842         390 :   r1 = nf_get_r1(nf);
    1843         390 :   r2 = nf_get_r2(nf);
    1844         390 :   x = nf_to_scalar_or_basis(nf, x);
    1845         385 :   if (!ind0) ind0 = identity_perm(r1+r2);
    1846         385 :   switch(typ(ind0))
    1847             :   {
    1848             :     case t_INT: case t_VEC: case t_COL:
    1849          30 :       ind = gtovecsmall(ind0); break;
    1850             :     case t_VECSMALL:
    1851         355 :       ind = ind0; break;
    1852             :     default:
    1853           0 :       pari_err_TYPE("nfeltsign",ind0);
    1854             :       return NULL; /* LCOV_EXCL_LINE */
    1855             :   }
    1856         385 :   l = lg(ind);
    1857         385 :   for (i = 1; i < l; i++) chk_ind("nfeltembed", ind[i], r1+r2);
    1858         380 :   if (typ(x) != t_COL)
    1859             :   {
    1860         140 :     if (typ(ind0) != t_INT) x = const_vec(l-1, x);
    1861         140 :     return gerepilecopy(av, x);
    1862             :   }
    1863         240 :   x = Q_primitive_part(x, &cx); M = nf_get_M(nf);
    1864         240 :   v = cgetg(l, t_VEC);
    1865         660 :   for (i = 1; i < l; i++)
    1866             :   {
    1867         420 :     GEN t = nfembed_i(M, x, ind[i]);
    1868         420 :     if (cx) t = gmul(t, cx);
    1869         420 :     gel(v,i) = t;
    1870             :   }
    1871         240 :   if (typ(ind0) == t_INT) v = gel(v,1);
    1872         240 :   return gerepilecopy(av, v);
    1873             : }
    1874             : 
    1875             : /* return the vector of signs of x; the matrix of such if x is a vector
    1876             :  * of nf elements */
    1877             : GEN
    1878         395 : nfsign(GEN nf, GEN x)
    1879             : {
    1880             :   long i, l;
    1881             :   GEN archp, S;
    1882             : 
    1883         395 :   nf = checknf(nf);
    1884         395 :   archp = identity_perm( nf_get_r1(nf) );
    1885         395 :   if (typ(x) != t_VEC) return nfsign_arch(nf, x, archp);
    1886         130 :   l = lg(x); S = cgetg(l, t_MAT);
    1887         130 :   for (i=1; i<l; i++) gel(S,i) = nfsign_arch(nf, gel(x,i), archp);
    1888         130 :   return S;
    1889             : }
    1890             : 
    1891             : /* x integral elt, A integral ideal in HNF; reduce x mod A */
    1892             : static GEN
    1893      400785 : zk_modHNF(GEN x, GEN A)
    1894      400785 : { return (typ(x) == t_COL)?  ZC_hnfrem(x, A): modii(x, gcoeff(A,1,1)); }
    1895             : 
    1896             : /* given an element x in Z_K and an integral ideal y in HNF, coprime with x,
    1897             :    outputs an element inverse of x modulo y */
    1898             : GEN
    1899         105 : nfinvmodideal(GEN nf, GEN x, GEN y)
    1900             : {
    1901         105 :   pari_sp av = avma;
    1902         105 :   GEN a, yZ = gcoeff(y,1,1);
    1903             : 
    1904         105 :   if (is_pm1(yZ)) return gen_0;
    1905         105 :   x = nf_to_scalar_or_basis(nf, x);
    1906         105 :   if (typ(x) == t_INT) return gerepileupto(av, Fp_inv(x, yZ));
    1907             : 
    1908          55 :   a = hnfmerge_get_1(idealhnf_principal(nf,x), y);
    1909          55 :   if (!a) pari_err_INV("nfinvmodideal", x);
    1910          55 :   return gerepileupto(av, zk_modHNF(nfdiv(nf,a,x), y));
    1911             : }
    1912             : 
    1913             : static GEN
    1914      191980 : nfsqrmodideal(GEN nf, GEN x, GEN id)
    1915      191980 : { return zk_modHNF(nfsqri(nf,x), id); }
    1916             : static GEN
    1917      449140 : nfmulmodideal(GEN nf, GEN x, GEN y, GEN id)
    1918      449140 : { return x? zk_modHNF(nfmuli(nf,x,y), id): y; }
    1919             : /* assume x integral, k integer, A in HNF */
    1920             : GEN
    1921      278375 : nfpowmodideal(GEN nf,GEN x,GEN k,GEN A)
    1922             : {
    1923      278375 :   long s = signe(k);
    1924             :   pari_sp av;
    1925             :   GEN y;
    1926             : 
    1927      278375 :   if (!s) return gen_1;
    1928      278375 :   av = avma;
    1929      278375 :   x = nf_to_scalar_or_basis(nf, x);
    1930      278375 :   if (typ(x) != t_COL) return Fp_pow(x, k, gcoeff(A,1,1));
    1931      146385 :   if (s < 0) { x = nfinvmodideal(nf, x,A); k = absi(k); }
    1932      146385 :   for(y = NULL;;)
    1933             :   {
    1934      338365 :     if (mpodd(k)) y = nfmulmodideal(nf,y,x,A);
    1935      338365 :     k = shifti(k,-1); if (!signe(k)) break;
    1936      191980 :     x = nfsqrmodideal(nf,x,A);
    1937      191980 :   }
    1938      146385 :   return gerepileupto(av, y);
    1939             : }
    1940             : 
    1941             : /* a * g^n mod id */
    1942             : static GEN
    1943      247730 : nfmulpowmodideal(GEN nf, GEN a, GEN g, GEN n, GEN id)
    1944             : {
    1945      247730 :   return nfmulmodideal(nf, a, nfpowmodideal(nf,g,n,id), id);
    1946             : }
    1947             : 
    1948             : /* assume (num(g[i]), id) = 1 for all i. Return prod g[i]^e[i] mod id.
    1949             :  * EX = multiple of exponent of (O_K/id)^* */
    1950             : GEN
    1951       95985 : famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id, GEN EX)
    1952             : {
    1953       95985 :   GEN EXo2, plus = NULL, minus = NULL, idZ = gcoeff(id,1,1);
    1954       95985 :   long i, lx = lg(g);
    1955             : 
    1956       95985 :   if (is_pm1(idZ)) return gen_1; /* id = Z_K */
    1957       95985 :   EXo2 = (expi(EX) > 10)? shifti(EX,-1): NULL;
    1958      513095 :   for (i = 1; i < lx; i++)
    1959             :   {
    1960      417110 :     GEN h, n = centermodii(gel(e,i), EX, EXo2);
    1961      417110 :     long sn = signe(n);
    1962      417110 :     if (!sn) continue;
    1963             : 
    1964      186155 :     h = nf_to_scalar_or_basis(nf, gel(g,i));
    1965      186155 :     switch(typ(h))
    1966             :     {
    1967      109735 :       case t_INT: break;
    1968             :       case t_FRAC:
    1969           0 :         h = Fp_div(gel(h,1), gel(h,2), idZ); break;
    1970             :       default:
    1971             :       {
    1972             :         GEN dh;
    1973       76420 :         h = Q_remove_denom(h, &dh);
    1974       76420 :         if (dh) h = FpC_Fp_mul(h, Fp_inv(dh,idZ), idZ);
    1975             :       }
    1976             :     }
    1977      186155 :     if (sn > 0)
    1978      185150 :       plus = nfmulpowmodideal(nf, plus, h, n, id);
    1979             :     else /* sn < 0 */
    1980        1005 :       minus = nfmulpowmodideal(nf, minus, h, absi(n), id);
    1981             :   }
    1982       95985 :   if (minus) plus = nfmulmodideal(nf, plus, nfinvmodideal(nf,minus,id), id);
    1983       95985 :   return plus? plus: gen_1;
    1984             : }
    1985             : 
    1986             : /* given 2 integral ideals x, y in HNF s.t x | y | x^2, compute (1+x)/(1+y) in
    1987             :  * the form [[cyc],[gen], U], where U := ux^-1 as a pair [ZM, denom(U)] */
    1988             : static GEN
    1989        9265 : zidealij(GEN x, GEN y)
    1990             : {
    1991        9265 :   GEN U, G, cyc, xp = gcoeff(x,1,1), xi = hnf_invscale(x, xp);
    1992             :   long j, N;
    1993             : 
    1994             :   /* x^(-1) y = relations between the 1 + x_i (HNF) */
    1995        9265 :   cyc = ZM_snf_group(ZM_Z_divexact(ZM_mul(xi, y), xp), &U, &G);
    1996        9265 :   N = lg(cyc); G = ZM_mul(x,G); settyp(G, t_VEC); /* new generators */
    1997       43720 :   for (j=1; j<N; j++)
    1998             :   {
    1999       34455 :     GEN c = gel(G,j);
    2000       34455 :     gel(c,1) = addiu(gel(c,1), 1); /* 1 + g_j */
    2001       34455 :     if (ZV_isscalar(c)) gel(G,j) = gel(c,1);
    2002             :   }
    2003        9265 :   return mkvec4(cyc, G, ZM_mul(U,xi), xp);
    2004             : }
    2005             : 
    2006             : /* lg(x) > 1, x + 1; shallow */
    2007             : static GEN
    2008        1440 : ZC_add1(GEN x)
    2009             : {
    2010        1440 :   long i, l = lg(x);
    2011        1440 :   GEN y = cgetg(l, t_COL);
    2012        1440 :   for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
    2013        1440 :   gel(y,1) = addiu(gel(x,1), 1); return y;
    2014             : }
    2015             : /* lg(x) > 1, x - 1; shallow */
    2016             : static GEN
    2017        1100 : ZC_sub1(GEN x)
    2018             : {
    2019        1100 :   long i, l = lg(x);
    2020        1100 :   GEN y = cgetg(l, t_COL);
    2021        1100 :   for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
    2022        1100 :   gel(y,1) = subiu(gel(x,1), 1); return y;
    2023             : }
    2024             : 
    2025             : /* x,y are t_INT or ZC */
    2026             : static GEN
    2027           0 : zkadd(GEN x, GEN y)
    2028             : {
    2029           0 :   long tx = typ(x);
    2030           0 :   if (tx == typ(y))
    2031           0 :     return tx == t_INT? addii(x,y): ZC_add(x,y);
    2032             :   else
    2033           0 :     return tx == t_INT? ZC_Z_add(y,x): ZC_Z_add(x,y);
    2034             : }
    2035             : /* x a t_INT or ZC, x+1; shallow */
    2036             : static GEN
    2037        2280 : zkadd1(GEN x)
    2038             : {
    2039        2280 :   long tx = typ(x);
    2040        2280 :   return tx == t_INT? addiu(x,1): ZC_add1(x);
    2041             : }
    2042             : /* x a t_INT or ZC, x-1; shallow */
    2043             : static GEN
    2044        2280 : zksub1(GEN x)
    2045             : {
    2046        2280 :   long tx = typ(x);
    2047        2280 :   return tx == t_INT? subiu(x,1): ZC_sub1(x);
    2048             : }
    2049             : /* x,y are t_INT or ZC; x - y */
    2050             : static GEN
    2051           0 : zksub(GEN x, GEN y)
    2052             : {
    2053           0 :   long tx = typ(x), ty = typ(y);
    2054           0 :   if (tx == ty)
    2055           0 :     return tx == t_INT? subii(x,y): ZC_sub(x,y);
    2056             :   else
    2057           0 :     return tx == t_INT? Z_ZC_sub(x,y): ZC_Z_sub(x,y);
    2058             : }
    2059             : /* x is t_INT or ZM (mult. map), y is t_INT or ZC; x * y */
    2060             : static GEN
    2061        2280 : zkmul(GEN x, GEN y)
    2062             : {
    2063        2280 :   long tx = typ(x), ty = typ(y);
    2064        2280 :   if (ty == t_INT)
    2065        1180 :     return tx == t_INT? mulii(x,y): ZC_Z_mul(gel(x,1),y);
    2066             :   else
    2067        1100 :     return tx == t_INT? ZC_Z_mul(y,x): ZM_ZC_mul(x,y);
    2068             : }
    2069             : 
    2070             : /* (U,V) = 1 coprime ideals. Want z = x mod U, = y mod V; namely
    2071             :  * z =vx + uy = v(x-y) + y, where u + v = 1, u in U, v in V.
    2072             :  * zkc = [v, UV], v a t_INT or ZM (mult. by v map), UV a ZM (ideal in HNF);
    2073             :  * shallow */
    2074             : GEN
    2075           0 : zkchinese(GEN zkc, GEN x, GEN y)
    2076             : {
    2077           0 :   GEN v = gel(zkc,1), UV = gel(zkc,2), z = zkadd(zkmul(v, zksub(x,y)), y);
    2078           0 :   return zk_modHNF(z, UV);
    2079             : }
    2080             : /* special case z = x mod U, = 1 mod V; shallow */
    2081             : GEN
    2082        2280 : zkchinese1(GEN zkc, GEN x)
    2083             : {
    2084        2280 :   GEN v = gel(zkc,1), UV = gel(zkc,2), z = zkadd1(zkmul(v, zksub1(x)));
    2085        2280 :   return (typ(z) == t_INT)? z: ZC_hnfrem(z, UV);
    2086             : }
    2087             : static GEN
    2088        1755 : zkVchinese1(GEN zkc, GEN v)
    2089             : {
    2090             :   long i, ly;
    2091        1755 :   GEN y = cgetg_copy(v, &ly);
    2092        1755 :   for (i=1; i<ly; i++) gel(y,i) = zkchinese1(zkc, gel(v,i));
    2093        1755 :   return y;
    2094             : }
    2095             : 
    2096             : /* prepare to solve z = x (mod A), z = y mod (B) [zkchinese or zkchinese1] */
    2097             : GEN
    2098        1570 : zkchineseinit(GEN nf, GEN A, GEN B, GEN AB)
    2099             : {
    2100             :   GEN v;
    2101             :   long e;
    2102        1570 :   nf = checknf(nf);
    2103        1570 :   v = idealaddtoone_raw(nf, A, B);
    2104        1570 :   if ((e = gexpo(v)) > 5)
    2105             :   {
    2106         165 :     GEN b = (typ(v) == t_COL)? v: scalarcol_shallow(v, nf_get_degree(nf));
    2107         165 :     b= ZC_reducemodlll(b, AB);
    2108         165 :     if (gexpo(b) < e) v = b;
    2109             :   }
    2110        1570 :   return mkvec2(zk_scalar_or_multable(nf,v), AB);
    2111             : }
    2112             : /* prepare to solve z = x (mod A), z = 1 mod (B)
    2113             :  * and then         z = 1 (mod A), z = y mod (B) [zkchinese1 twice] */
    2114             : static GEN
    2115         185 : zkchinese1init2(GEN nf, GEN A, GEN B, GEN AB)
    2116             : {
    2117         185 :   GEN zkc = zkchineseinit(nf, A, B, AB);
    2118         185 :   GEN mv = gel(zkc,1), mu;
    2119         185 :   if (typ(mv) == t_INT) return mkvec2(zkc, mkvec2(subui(1,mv),AB));
    2120          25 :   mu = RgM_Rg_add_shallow(ZM_neg(mv), gen_1);
    2121          25 :   return mkvec2(mkvec2(mv,AB), mkvec2(mu,AB));
    2122             : }
    2123             : 
    2124             : static GEN
    2125      197485 : apply_U(GEN L, GEN a)
    2126             : {
    2127      197485 :   GEN e, U = gel(L,3), dU = gel(L,4);
    2128      197485 :   if (typ(a) == t_INT)
    2129       60410 :     e = ZC_Z_mul(gel(U,1), subiu(a, 1));
    2130             :   else
    2131             :   { /* t_COL */
    2132      137075 :     GEN t = gel(a,1);
    2133      137075 :     gel(a,1) = subiu(gel(a,1), 1); /* a -= 1 */
    2134      137075 :     e = ZM_ZC_mul(U, a);
    2135      137075 :     gel(a,1) = t; /* restore */
    2136             :   }
    2137      197485 :   return gdiv(e, dU);
    2138             : }
    2139             : 
    2140             : /* true nf; vectors of [[cyc],[g],U.X^-1]. Assume k > 1. */
    2141             : static GEN
    2142        6070 : principal_units(GEN nf, GEN pr, long k, GEN prk)
    2143             : {
    2144             :   GEN list, prb;
    2145        6070 :   ulong mask = quadratic_prec_mask(k);
    2146        6070 :   long a = 1;
    2147             : 
    2148        6070 :   if (DEBUGLEVEL>3) err_printf("treating pr^%ld, pr = %Ps\n",k,pr);
    2149        6070 :   prb = pr_hnf(nf,pr);
    2150        6070 :   list = vectrunc_init(k);
    2151       21405 :   while (mask > 1)
    2152             :   {
    2153        9265 :     GEN pra = prb;
    2154        9265 :     long b = a << 1;
    2155             : 
    2156        9265 :     if (mask & 1) b--;
    2157        9265 :     mask >>= 1;
    2158             :     /* compute 1 + pr^a / 1 + pr^b, 2a <= b */
    2159        9265 :     if(DEBUGLEVEL>3) err_printf("  treating a = %ld, b = %ld\n",a,b);
    2160        9265 :     prb = (b >= k)? prk: idealpows(nf,pr,b);
    2161        9265 :     vectrunc_append(list, zidealij(pra, prb));
    2162        9265 :     a = b;
    2163             :   }
    2164        6070 :   return list;
    2165             : }
    2166             : /* a = 1 mod (pr) return log(a) on local-gens of 1+pr/1+pr^k */
    2167             : static GEN
    2168      123310 : log_prk1(GEN nf, GEN a, long nh, GEN L2, GEN prk)
    2169             : {
    2170      123310 :   GEN y = cgetg(nh+1, t_COL);
    2171      123310 :   long j, iy, c = lg(L2)-1;
    2172      320790 :   for (j = iy = 1; j <= c; j++)
    2173             :   {
    2174      197485 :     GEN L = gel(L2,j), cyc = gel(L,1), gen = gel(L,2), E = apply_U(L,a);
    2175      197485 :     long i, nc = lg(cyc)-1;
    2176      197485 :     int last = (j == c);
    2177      767830 :     for (i = 1; i <= nc; i++, iy++)
    2178             :     {
    2179      570350 :       GEN t, e = gel(E,i);
    2180      570350 :       if (typ(e) != t_INT) pari_err_COPRIME("zlog_prk1", a, prk);
    2181      570345 :       t = Fp_neg(e, gel(cyc,i));
    2182      570345 :       gel(y,iy) = negi(t);
    2183      570345 :       if (!last && signe(t)) a = nfmulpowmodideal(nf, a, gel(gen,i), t, prk);
    2184             :     }
    2185             :   }
    2186      123305 :   return y;
    2187             : }
    2188             : /* true nf */
    2189             : static GEN
    2190        2730 : principal_units_relations(GEN nf, GEN L2, GEN prk, long nh)
    2191             : {
    2192        2730 :   GEN h = cgetg(nh+1,t_MAT);
    2193        2730 :   long ih, j, c = lg(L2)-1;
    2194        8655 :   for (j = ih = 1; j <= c; j++)
    2195             :   {
    2196        5925 :     GEN L = gel(L2,j), F = gel(L,1), G = gel(L,2);
    2197        5925 :     long k, lG = lg(G);
    2198       30565 :     for (k = 1; k < lG; k++,ih++)
    2199             :     { /* log(g^f) mod pr^e */
    2200       24640 :       GEN a = nfpowmodideal(nf,gel(G,k),gel(F,k),prk);
    2201       24640 :       gel(h,ih) = ZC_neg(log_prk1(nf, a, nh, L2, prk));
    2202       24640 :       gcoeff(h,ih,ih) = gel(F,k);
    2203             :     }
    2204             :   }
    2205        2730 :   return h;
    2206             : }
    2207             : /* true nf; e > 1; multiplicative group (1 + pr) / (1 + pr^k),
    2208             :  * prk = pr^k or NULL */
    2209             : static GEN
    2210        6070 : idealprincipalunits_i(GEN nf, GEN pr, long k, GEN *pU)
    2211             : {
    2212        6070 :   GEN cyc, gen, L2, prk = idealpows(nf, pr, k);
    2213             : 
    2214        6070 :   L2 = principal_units(nf, pr, k, prk);
    2215        6070 :   if (k == 2)
    2216             :   {
    2217        3340 :     GEN L = gel(L2,1);
    2218        3340 :     cyc = gel(L,1);
    2219        3340 :     gen = gel(L,2);
    2220        3340 :     if (pU) *pU = matid(lg(gen)-1);
    2221             :   }
    2222             :   else
    2223             :   {
    2224        2730 :     long c = lg(L2), j;
    2225        2730 :     GEN EX, h, Ui, vg = cgetg(c, t_VEC);
    2226        2730 :     for (j = 1; j < c; j++) gel(vg, j) = gmael(L2,j,2);
    2227        2730 :     vg = shallowconcat1(vg);
    2228        2730 :     h = principal_units_relations(nf, L2, prk, lg(vg)-1);
    2229        2730 :     h = ZM_hnfall_i(h, NULL, 0);
    2230        2730 :     cyc = ZM_snf_group(h, pU, &Ui);
    2231        2730 :     c = lg(Ui); gen = cgetg(c, t_VEC); EX = gel(cyc,1);
    2232       19325 :     for (j = 1; j < c; j++)
    2233       16595 :       gel(gen,j) = famat_to_nf_modideal_coprime(nf, vg, gel(Ui,j), prk, EX);
    2234             :   }
    2235        6070 :   return mkvec4(cyc, gen, prk, L2);
    2236             : }
    2237             : GEN
    2238          80 : idealprincipalunits(GEN nf, GEN pr, long k)
    2239             : {
    2240             :   pari_sp av;
    2241             :   GEN v;
    2242          80 :   nf = checknf(nf);
    2243          80 :   if (k == 1) { checkprid(pr); retmkvec3(gen_1,cgetg(1,t_VEC),cgetg(1,t_VEC)); }
    2244          75 :   av = avma; v = idealprincipalunits_i(nf, pr, k, NULL);
    2245          75 :   return gerepilecopy(av, mkvec3(powiu(pr_norm(pr), k-1), gel(v,1), gel(v,2)));
    2246             : }
    2247             : 
    2248             : /* true nf; given an ideal pr^k dividing an integral ideal x (in HNF form)
    2249             :  * compute an 'sprk', the structure of G = (Z_K/pr^k)^* [ x = NULL for x=pr^k ]
    2250             :  * Return a vector with at least 4 components [cyc],[gen],[HNF pr^k,pr,k],ff,
    2251             :  * where
    2252             :  * cyc : type of G as abelian group (SNF)
    2253             :  * gen : generators of G, coprime to x
    2254             :  * pr^k: in HNF
    2255             :  * ff  : data for log_g in (Z_K/pr)^*
    2256             :  * Two extra components are present iff k > 1: L2, U
    2257             :  * L2  : list of data structures to compute local DL in (Z_K/pr)^*,
    2258             :  *       and 1 + pr^a/ 1 + pr^b for various a < b <= min(2a, k)
    2259             :  * U   : base change matrices to convert a vector of local DL to DL wrt gen */
    2260             : static GEN
    2261        8540 : sprkinit(GEN nf, GEN pr, GEN gk, GEN x)
    2262             : {
    2263             :   GEN T, p, modpr, cyc, gen, g, g0, ord0, A, prk, U, L2;
    2264        8540 :   long k = itos(gk), f = pr_get_f(pr);
    2265             : 
    2266        8540 :   if(DEBUGLEVEL>3) err_printf("treating pr^%ld, pr = %Ps\n",k,pr);
    2267        8540 :   modpr = nf_to_Fq_init(nf, &pr,&T,&p);
    2268             :   /* (Z_K / pr)^* */
    2269        8540 :   if (f == 1)
    2270             :   {
    2271        3675 :     g0 = g = pgener_Fp(p);
    2272        3675 :     ord0 = get_arith_ZZM(subiu(p,1));
    2273             :   }
    2274             :   else
    2275             :   {
    2276        4865 :     g0 = g = gener_FpXQ(T,p, &ord0);
    2277        4865 :     g = Fq_to_nf(g, modpr);
    2278        4865 :     if (typ(g) == t_POL) g = poltobasis(nf, g);
    2279             :   }
    2280        8540 :   A = gel(ord0, 1); /* Norm(pr)-1 */
    2281        8540 :   if (k == 1)
    2282             :   {
    2283        2545 :     cyc = mkvec(A);
    2284        2545 :     gen = mkvec(g);
    2285        2545 :     prk = pr_hnf(nf,pr);
    2286        2545 :     L2 = U = NULL;
    2287             :   }
    2288             :   else
    2289             :   { /* local-gens of (1 + pr)/(1 + pr^k) = SNF-gens * U */
    2290             :     GEN AB, B, u, v, w;
    2291             :     long j, l;
    2292        5995 :     w = idealprincipalunits_i(nf, pr, k, &U);
    2293             :     /* incorporate (Z_K/pr)^*, order A coprime to B = expo(1+pr/1+pr^k)*/
    2294        5995 :     cyc = leafcopy(gel(w,1)); B = gel(cyc,1); AB = mulii(A,B);
    2295        5995 :     gen = leafcopy(gel(w,2));
    2296        5995 :     prk = gel(w,3);
    2297        5995 :     g = nfpowmodideal(nf, g, B, prk);
    2298        5995 :     g0 = Fq_pow(g0, modii(B,A), T, p); /* update primitive root */
    2299        5995 :     L2 = mkvec3(A, g, gel(w,4));
    2300        5995 :     gel(cyc,1) = AB;
    2301        5995 :     gel(gen,1) = nfmulmodideal(nf, gel(gen,1), g, prk);
    2302        5995 :     u = mulii(Fp_inv(A,B), A);
    2303        5995 :     v = subui(1, u); l = lg(U);
    2304        5995 :     for (j = 1; j < l; j++) gcoeff(U,1,j) = Fp_mul(u, gcoeff(U,1,j), AB);
    2305        5995 :     U = mkvec2(Rg_col_ei(v, lg(gen)-1, 1), U);
    2306             :   }
    2307             :   /* local-gens of (Z_K/pr^k)^* = SNF-gens * U */
    2308        8540 :   if (x)
    2309             :   {
    2310        1385 :     GEN uv = zkchineseinit(nf, idealdivpowprime(nf,x,pr,gk), prk, x);
    2311        1385 :     gen = zkVchinese1(uv, gen);
    2312             :   }
    2313        8540 :   return mkvecn(U? 6: 4, cyc, gen, prk, mkvec3(modpr,g0,ord0), L2, U);
    2314             : }
    2315             : static GEN
    2316      175945 : sprk_get_cyc(GEN s) { return gel(s,1); }
    2317             : static GEN
    2318       69410 : sprk_get_expo(GEN s)
    2319             : {
    2320       69410 :   GEN cyc = sprk_get_cyc(s);
    2321       69410 :   return lg(cyc) == 1? gen_1: gel(cyc, 1);
    2322             : }
    2323             : static GEN
    2324        4210 : sprk_get_gen(GEN s) { return gel(s,2); }
    2325             : static GEN
    2326      168080 : sprk_get_prk(GEN s) { return gel(s,3); }
    2327             : static GEN
    2328      181330 : sprk_get_ff(GEN s) { return gel(s,4); }
    2329             : static GEN
    2330       73170 : sprk_get_pr(GEN s) { GEN ff = gel(s,4); return modpr_get_pr(gel(ff,1)); }
    2331             : /* A = Npr-1, <g> = (Z_K/pr)^*, L2 to 1 + pr / 1 + pr^k */
    2332             : static void
    2333      100840 : sprk_get_L2(GEN s, GEN *A, GEN *g, GEN *L2)
    2334      100840 : { GEN v = gel(s,5); *A = gel(v,1); *g = gel(v,2); *L2 = gel(v,3); }
    2335             : static void
    2336       98670 : sprk_get_U2(GEN s, GEN *U1, GEN *U2)
    2337       98670 : { GEN v = gel(s,6); *U1 = gel(v,1); *U2 = gel(v,2); }
    2338             : static int
    2339      181330 : sprk_is_prime(GEN s) { return lg(s) == 5; }
    2340             : 
    2341             : static GEN
    2342       69410 : famat_zlog_pr(GEN nf, GEN g, GEN e, GEN sprk)
    2343             : {
    2344       69410 :   GEN pr = sprk_get_pr(sprk);
    2345       69410 :   GEN prk = sprk_get_prk(sprk);
    2346       69410 :   GEN x = famat_makecoprime(nf, g, e, pr, prk, sprk_get_expo(sprk));
    2347       69410 :   return zlog_pr(nf, x, sprk);
    2348             : }
    2349             : /* log_g(a) in (Z_K/pr)^* */
    2350             : static GEN
    2351      181330 : nf_log(GEN nf, GEN a, GEN ff)
    2352             : {
    2353      181330 :   GEN pr = gel(ff,1), g = gel(ff,2), ord = gel(ff,3);
    2354      181330 :   GEN T,p, modpr = nf_to_Fq_init(nf, &pr, &T, &p);
    2355      181330 :   return Fq_log(nf_to_Fq(nf,a,modpr), g, ord, T, p);
    2356             : }
    2357             : /* a in Z_K (t_COL or t_INT), pr prime ideal, sprk = sprkinit(nf,pr,k,x).
    2358             :  * return log(a) on SNF-generators of (Z_K/pr^k)^**/
    2359             : GEN
    2360      182070 : zlog_pr(GEN nf, GEN a, GEN sprk)
    2361             : {
    2362             :   GEN e, prk, A, g, L2, U1, U2, y;
    2363             : 
    2364      182070 :   if (typ(a) == t_MAT) return famat_zlog_pr(nf, gel(a,1), gel(a,2), sprk);
    2365             : 
    2366      181330 :   e = nf_log(nf, a, sprk_get_ff(sprk));
    2367      181330 :   if (sprk_is_prime(sprk)) return mkcol(e); /* k = 1 */
    2368       98670 :   prk = sprk_get_prk(sprk);
    2369       98670 :   sprk_get_L2(sprk, &A,&g,&L2);
    2370       98670 :   if (signe(e))
    2371             :   {
    2372       28330 :     e = Fp_neg(e, A);
    2373       28330 :     a = nfmulpowmodideal(nf, a, g, e, prk);
    2374             :   }
    2375       98670 :   sprk_get_U2(sprk, &U1,&U2);
    2376       98670 :   y = ZM_ZC_mul(U2, log_prk1(nf, a, lg(U2)-1, L2, prk));
    2377       98665 :   if (signe(e)) y = ZC_sub(y, ZC_Z_mul(U1,e));
    2378       98665 :   return vecmodii(y, sprk_get_cyc(sprk));
    2379             : }
    2380             : GEN
    2381        4330 : zlog_pr_init(GEN nf, GEN pr, long k)
    2382        4330 : { return sprkinit(checknf(nf),pr,utoipos(k),NULL);}
    2383             : GEN
    2384         270 : vzlog_pr(GEN nf, GEN v, GEN sprk)
    2385             : {
    2386         270 :   long l = lg(v), i;
    2387         270 :   GEN w = cgetg(l, t_MAT);
    2388         270 :   for (i = 1; i < l; i++) gel(w,i) = zlog_pr(nf, gel(v,i), sprk);
    2389         270 :   return w;
    2390             : }
    2391             : 
    2392             : static GEN
    2393       72040 : famat_zlog(GEN nf, GEN fa, GEN sgn, zlog_S *S)
    2394             : {
    2395       72040 :   long i, n0, n = lg(S->U)-1;
    2396             :   GEN g, e, y;
    2397       72040 :   if (lg(fa) == 1) return zerocol(n);
    2398       72040 :   g = gel(fa,1);
    2399       72040 :   e = gel(fa,2);
    2400       72040 :   y = cgetg(n+1, t_COL);
    2401       72040 :   n0 = lg(S->sprk)-1; /* n0 = n (trivial arch. part), or n-1 */
    2402       72040 :   for (i=1; i <= n0; i++) gel(y,i) = famat_zlog_pr(nf, g, e, gel(S->sprk,i));
    2403       72040 :   if (n0 != n)
    2404             :   {
    2405       65280 :     if (!sgn) sgn = nfsign_arch(nf, fa, S->archp);
    2406       65280 :     gel(y,n) = Flc_to_ZC(sgn);
    2407             :   }
    2408       72040 :   return y;
    2409             : }
    2410             : 
    2411             : /* assume that cyclic factors are normalized, in particular != [1] */
    2412             : static GEN
    2413        4060 : split_U(GEN U, GEN Sprk)
    2414             : {
    2415        4060 :   long t = 0, k, n, l = lg(Sprk);
    2416        4060 :   GEN vU = cgetg(l+1, t_VEC);
    2417        7720 :   for (k = 1; k < l; k++)
    2418             :   {
    2419        3660 :     n = lg(sprk_get_cyc(gel(Sprk,k))) - 1; /* > 0 */
    2420        3660 :     gel(vU,k) = vecslice(U, t+1, t+n);
    2421        3660 :     t += n;
    2422             :   }
    2423             :   /* t+1 .. lg(U)-1 */
    2424        4060 :   n = lg(U) - t - 1; /* can be 0 */
    2425        4060 :   if (!n) setlg(vU,l); else gel(vU,l) = vecslice(U, t+1, t+n);
    2426        4060 :   return vU;
    2427             : }
    2428             : 
    2429             : void
    2430      154895 : init_zlog(zlog_S *S, GEN bid)
    2431             : {
    2432      154895 :   GEN fa2 = bid_get_fact2(bid);
    2433      154895 :   S->U = bid_get_U(bid);
    2434      154895 :   S->hU = lg(bid_get_cyc(bid))-1;
    2435      154895 :   S->archp = bid_get_archp(bid);
    2436      154895 :   S->sprk = bid_get_sprk(bid);
    2437      154895 :   S->bid = bid;
    2438      154895 :   S->P = gel(fa2,1);
    2439      154895 :   S->k = gel(fa2,2);
    2440      154895 :   S->no2 = lg(S->P) == lg(gel(bid_get_fact(bid),1));
    2441      154895 : }
    2442             : 
    2443             : /* a a t_FRAC/t_INT, reduce mod bid */
    2444             : static GEN
    2445           5 : Q_mod_bid(GEN bid, GEN a)
    2446             : {
    2447           5 :   GEN xZ = gcoeff(bid_get_ideal(bid),1,1);
    2448           5 :   GEN b = Rg_to_Fp(a, xZ);
    2449           5 :   if (gsigne(a) < 0) b = subii(b, xZ);
    2450           5 :   return b;
    2451             : }
    2452             : /* Return decomposition of a on the CRT generators blocks attached to the
    2453             :  * S->sprk and sarch; sgn = sign(a, S->arch), NULL if unknown */
    2454             : static GEN
    2455      108955 : zlog(GEN nf, GEN a, GEN sgn, zlog_S *S)
    2456             : {
    2457             :   long k, l;
    2458             :   GEN y;
    2459      108955 :   a = nf_to_scalar_or_basis(nf, a);
    2460      108955 :   switch(typ(a))
    2461             :   {
    2462        8830 :     case t_INT: break;
    2463           5 :     case t_FRAC: a = Q_mod_bid(S->bid, a); break;
    2464             :     default: /* case t_COL: */
    2465             :     {
    2466             :       GEN den;
    2467      100120 :       check_nfelt(a, &den);
    2468      100120 :       if (den)
    2469             :       {
    2470       29280 :         a = Q_muli_to_int(a, den);
    2471       29280 :         a = mkmat2(mkcol2(a, den), mkcol2(gen_1, gen_m1));
    2472       29280 :         return famat_zlog(nf, a, sgn, S);
    2473             :       }
    2474             :     }
    2475             :   }
    2476       79675 :   if (sgn)
    2477        6470 :     sgn = (lg(sgn) == 1)? NULL: leafcopy(sgn);
    2478             :   else
    2479       73205 :     sgn = (lg(S->archp) == 1)? NULL: nfsign_arch(nf, a, S->archp);
    2480       79675 :   l = lg(S->sprk);
    2481       79675 :   y = cgetg(sgn? l+1: l, t_COL);
    2482      177555 :   for (k = 1; k < l; k++)
    2483             :   {
    2484       97885 :     GEN sprk = gel(S->sprk,k);
    2485       97885 :     gel(y,k) = zlog_pr(nf, a, sprk);
    2486             :   }
    2487       79670 :   if (sgn) gel(y,l) = Flc_to_ZC(sgn);
    2488       79670 :   return y;
    2489             : }
    2490             : 
    2491             : /* true nf */
    2492             : GEN
    2493        1770 : pr_basis_perm(GEN nf, GEN pr)
    2494             : {
    2495        1770 :   long f = pr_get_f(pr);
    2496             :   GEN perm;
    2497        1770 :   if (f == nf_get_degree(nf)) return identity_perm(f);
    2498         940 :   perm = cgetg(f+1, t_VECSMALL);
    2499         940 :   perm[1] = 1;
    2500         940 :   if (f > 1)
    2501             :   {
    2502         285 :     GEN H = pr_hnf(nf,pr);
    2503             :     long i, k;
    2504        1045 :     for (i = k = 2; k <= f; i++)
    2505             :     {
    2506         760 :       if (is_pm1(gcoeff(H,i,i))) continue;
    2507         600 :       perm[k++] = i;
    2508             :     }
    2509             :   }
    2510         940 :   return perm;
    2511             : }
    2512             : 
    2513             : /* \sum U[i]*y[i], U[i] ZM, y[i] ZC. We allow lg(y) > lg(U). */
    2514             : static GEN
    2515      151710 : ZMV_ZCV_mul(GEN U, GEN y)
    2516             : {
    2517      151710 :   long i, l = lg(U);
    2518      151710 :   GEN z = NULL;
    2519      151710 :   if (l == 1) return cgetg(1,t_COL);
    2520      433995 :   for (i = 1; i < l; i++)
    2521             :   {
    2522      282285 :     GEN u = ZM_ZC_mul(gel(U,i), gel(y,i));
    2523      282285 :     z = z? ZC_add(z, u): u;
    2524             :   }
    2525      151710 :   return z;
    2526             : }
    2527             : /* A * (U[1], ..., U[d] */
    2528             : static GEN
    2529         370 : ZM_ZMV_mul(GEN A, GEN U)
    2530             : {
    2531             :   long i, l;
    2532         370 :   GEN V = cgetg_copy(U,&l);
    2533         370 :   for (i = 1; i < l; i++) gel(V,i) = ZM_mul(A,gel(U,i));
    2534         370 :   return V;
    2535             : }
    2536             : 
    2537             : /* Log on bid.gen of generators of P_{1,I pr^{e-1}} / P_{1,I pr^e} (I,pr) = 1,
    2538             :  * defined implicitly via CRT. 'ind' is the index of pr in modulus
    2539             :  * factorization */
    2540             : GEN
    2541        7390 : log_gen_pr(zlog_S *S, long ind, GEN nf, long e)
    2542             : {
    2543        7390 :   GEN A, sprk = gel(S->sprk,ind), Uind = gel(S->U, ind);
    2544             : 
    2545        7390 :   if (e == 1) retmkmat( gel(Uind,1) );
    2546             :   else
    2547             :   {
    2548        3760 :     GEN G, pr = sprk_get_pr(sprk);
    2549             :     long i, l;
    2550        3760 :     if (e == 2)
    2551             :     {
    2552        2170 :       GEN A, g, L, L2; sprk_get_L2(sprk,&A,&g,&L2); L = gel(L2,1);
    2553        2170 :       G = gel(L,2); l = lg(G);
    2554             :     }
    2555             :     else
    2556             :     {
    2557        1590 :       GEN perm = pr_basis_perm(nf,pr), PI = nfpow_u(nf, pr_get_gen(pr), e-1);
    2558        1590 :       l = lg(perm);
    2559        1590 :       G = cgetg(l, t_VEC);
    2560        1590 :       if (typ(PI) == t_INT)
    2561             :       { /* zk_ei_mul doesn't allow t_INT */
    2562         825 :         long N = nf_get_degree(nf);
    2563         825 :         gel(G,1) = addiu(PI,1);
    2564        1275 :         for (i = 2; i < l; i++)
    2565             :         {
    2566         450 :           GEN z = col_ei(N, 1);
    2567         450 :           gel(G,i) = z; gel(z, perm[i]) = PI;
    2568             :         }
    2569             :       }
    2570             :       else
    2571             :       {
    2572         765 :         gel(G,1) = nfadd(nf, gen_1, PI);
    2573         915 :         for (i = 2; i < l; i++)
    2574         150 :           gel(G,i) = nfadd(nf, gen_1, zk_ei_mul(nf, PI, perm[i]));
    2575             :       }
    2576             :     }
    2577        3760 :     A = cgetg(l, t_MAT);
    2578        8630 :     for (i = 1; i < l; i++)
    2579        4870 :       gel(A,i) = ZM_ZC_mul(Uind, zlog_pr(nf, gel(G,i), sprk));
    2580        3760 :     return A;
    2581             :   }
    2582             : }
    2583             : /* Log on bid.gen of generator of P_{1,f} / P_{1,f v[index]}
    2584             :  * v = vector of r1 real places */
    2585             : GEN
    2586        4450 : log_gen_arch(zlog_S *S, long index)
    2587             : {
    2588        4450 :   GEN U = gel(S->U, lg(S->U)-1);
    2589        4450 :   return gel(U, index);
    2590             : }
    2591             : 
    2592             : /* compute bid.clgp: [h,cyc] or [h,cyc,gen] */
    2593             : static GEN
    2594        4815 : bid_grp(GEN nf, GEN U, GEN cyc, GEN g, GEN F, GEN sarch)
    2595             : {
    2596        4815 :   GEN G, h = ZV_prod(cyc);
    2597             :   long c;
    2598        4815 :   if (!U) return mkvec2(h,cyc);
    2599        4640 :   c = lg(U);
    2600        4640 :   G = cgetg(c,t_VEC);
    2601        4640 :   if (c > 1)
    2602             :   {
    2603        3960 :     GEN U0, Uoo, EX = gel(cyc,1); /* exponent of bid */
    2604        3960 :     long i, hU = nbrows(U), nba = lg(sarch_get_cyc(sarch))-1; /* #f_oo */
    2605        3960 :     if (!nba) { U0 = U; Uoo = NULL; }
    2606        2200 :     else if (nba == hU) { U0 = NULL; Uoo = U; }
    2607             :     else
    2608             :     {
    2609        1715 :       U0 = rowslice(U, 1, hU-nba);
    2610        1715 :       Uoo = rowslice(U, hU-nba+1, hU);
    2611             :     }
    2612       12595 :     for (i = 1; i < c; i++)
    2613             :     {
    2614        8635 :       GEN t = gen_1;
    2615        8635 :       if (U0) t = famat_to_nf_modideal_coprime(nf, g, gel(U0,i), F, EX);
    2616        8635 :       if (Uoo) t = set_sign_mod_divisor(nf, ZV_to_Flv(gel(Uoo,i),2), t, sarch);
    2617        8635 :       gel(G,i) = t;
    2618             :     }
    2619             :   }
    2620        4640 :   return mkvec3(h, cyc, G);
    2621             : }
    2622             : 
    2623             : /* remove prime ideals of norm 2 with exponent 1 from factorization */
    2624             : static GEN
    2625        4590 : famat_strip2(GEN fa)
    2626             : {
    2627        4590 :   GEN P = gel(fa,1), E = gel(fa,2), Q, F;
    2628        4590 :   long l = lg(P), i, j;
    2629        4590 :   Q = cgetg(l, t_COL);
    2630        4590 :   F = cgetg(l, t_COL);
    2631        8830 :   for (i = j = 1; i < l; i++)
    2632             :   {
    2633        4240 :     GEN pr = gel(P,i), e = gel(E,i);
    2634        4240 :     if (!absequaliu(pr_get_p(pr), 2) || itou(e) != 1 || pr_get_f(pr) != 1)
    2635             :     {
    2636        4210 :       gel(Q,j) = pr;
    2637        4210 :       gel(F,j) = e; j++;
    2638             :     }
    2639             :   }
    2640        4590 :   setlg(Q,j);
    2641        4590 :   setlg(F,j); return mkmat2(Q,F);
    2642             : }
    2643             : 
    2644             : /* Compute [[ideal,arch], [h,[cyc],[gen]], idealfact, [liste], U]
    2645             :    flag may include nf_GEN | nf_INIT */
    2646             : static GEN
    2647        4605 : Idealstar_i(GEN nf, GEN ideal, long flag)
    2648             : {
    2649             :   long i, k, nbp, R1;
    2650        4605 :   GEN y, cyc, U, u1 = NULL, fa, fa2, sprk, x, arch, archp, E, P, sarch, gen;
    2651             : 
    2652        4605 :   nf = checknf(nf);
    2653        4605 :   R1 = nf_get_r1(nf);
    2654        4605 :   if (typ(ideal) == t_VEC && lg(ideal) == 3)
    2655             :   {
    2656        3080 :     arch = gel(ideal,2);
    2657        3080 :     ideal= gel(ideal,1);
    2658        3080 :     switch(typ(arch))
    2659             :     {
    2660             :       case t_VEC:
    2661        3055 :         if (lg(arch) != R1+1)
    2662           0 :           pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
    2663        3055 :         archp = vec01_to_indices(arch);
    2664        3055 :         break;
    2665             :       case t_VECSMALL:
    2666          25 :         archp = arch;
    2667          25 :         k = lg(archp)-1;
    2668          25 :         if (k && archp[k] > R1)
    2669           5 :           pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
    2670          20 :         arch = indices_to_vec01(archp, R1);
    2671          20 :         break;
    2672             :       default:
    2673           0 :         pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
    2674           0 :         return NULL;
    2675             :     }
    2676        3075 :   }
    2677             :   else
    2678             :   {
    2679        1525 :     arch = zerovec(R1);
    2680        1525 :     archp = cgetg(1, t_VECSMALL);
    2681             :   }
    2682        4600 :   if (is_nf_factor(ideal))
    2683             :   {
    2684         245 :     fa = ideal;
    2685         245 :     x = idealfactorback(nf, gel(fa,1), gel(fa,2), 0);
    2686             :   }
    2687             :   else
    2688             :   {
    2689        4355 :     fa = idealfactor(nf, ideal);
    2690        4350 :     x = ideal;
    2691             :   }
    2692        4595 :   if (typ(x) != t_MAT)  x = idealhnf_shallow(nf, x);
    2693        4595 :   if (lg(x) == 1) pari_err_DOMAIN("Idealstar", "ideal","=",gen_0,x);
    2694        4595 :   if (typ(gcoeff(x,1,1)) != t_INT)
    2695           5 :     pari_err_DOMAIN("Idealstar","denominator(ideal)", "!=",gen_1,x);
    2696        4590 :   sarch = nfarchstar(nf, x, archp);
    2697        4590 :   fa2 = famat_strip2(fa);
    2698        4590 :   P = gel(fa2,1);
    2699        4590 :   E = gel(fa2,2);
    2700        4590 :   nbp = lg(P)-1;
    2701        4590 :   sprk = cgetg(nbp+1,t_VEC);
    2702        4590 :   if (nbp)
    2703             :   {
    2704        3415 :     GEN t = (nbp==1)? NULL: x;
    2705        3415 :     cyc = cgetg(nbp+2,t_VEC);
    2706        3415 :     gen = cgetg(nbp+1,t_VEC);
    2707        7625 :     for (i = 1; i <= nbp; i++)
    2708             :     {
    2709        4210 :       GEN L = sprkinit(nf, gel(P,i), gel(E,i), t);
    2710        4210 :       gel(sprk,i) = L;
    2711        4210 :       gel(cyc,i) = sprk_get_cyc(L);
    2712             :       /* true gens are congruent to those mod x AND positive at archp */
    2713        4210 :       gel(gen,i) = sprk_get_gen(L);
    2714             :     }
    2715        3415 :     gel(cyc,i) = sarch_get_cyc(sarch);
    2716        3415 :     cyc = shallowconcat1(cyc);
    2717        3415 :     gen = shallowconcat1(gen);
    2718        3415 :     cyc = ZV_snf_group(cyc, &U, (flag & nf_GEN)? &u1: NULL);
    2719             :   }
    2720             :   else
    2721             :   {
    2722        1175 :     cyc = sarch_get_cyc(sarch);
    2723        1175 :     gen = cgetg(1,t_VEC);
    2724        1175 :     U = matid(lg(cyc)-1);
    2725        1175 :     if (flag & nf_GEN) u1 = U;
    2726             :   }
    2727        4590 :   y = bid_grp(nf, u1, cyc, gen, x, sarch);
    2728        4590 :   if (!(flag & nf_INIT)) return y;
    2729        4020 :   U = split_U(U, sprk);
    2730        4020 :   return mkvec5(mkvec2(x, arch), y, mkvec2(fa,fa2), mkvec2(sprk, sarch), U);
    2731             : }
    2732             : GEN
    2733        4415 : Idealstar(GEN nf, GEN ideal, long flag)
    2734             : {
    2735        4415 :   pari_sp av = avma;
    2736        4415 :   if (!nf) nf = nfinit(pol_x(0), DEFAULTPREC);
    2737        4415 :   return gerepilecopy(av, Idealstar_i(nf, ideal, flag));
    2738             : }
    2739             : GEN
    2740         190 : Idealstarprk(GEN nf, GEN pr, long k, long flag)
    2741             : {
    2742         190 :   pari_sp av = avma;
    2743         190 :   GEN z = Idealstar_i(nf, mkmat2(mkcol(pr),mkcols(k)), flag);
    2744         190 :   return gerepilecopy(av, z);
    2745             : }
    2746             : 
    2747             : /* FIXME: obsolete */
    2748             : GEN
    2749           0 : zidealstarinitgen(GEN nf, GEN ideal)
    2750           0 : { return Idealstar(nf,ideal, nf_INIT|nf_GEN); }
    2751             : GEN
    2752           0 : zidealstarinit(GEN nf, GEN ideal)
    2753           0 : { return Idealstar(nf,ideal, nf_INIT); }
    2754             : GEN
    2755           0 : zidealstar(GEN nf, GEN ideal)
    2756           0 : { return Idealstar(nf,ideal, nf_GEN); }
    2757             : 
    2758             : GEN
    2759          45 : idealstar0(GEN nf, GEN ideal,long flag)
    2760             : {
    2761          45 :   switch(flag)
    2762             :   {
    2763           0 :     case 0: return Idealstar(nf,ideal, nf_GEN);
    2764          35 :     case 1: return Idealstar(nf,ideal, nf_INIT);
    2765          10 :     case 2: return Idealstar(nf,ideal, nf_INIT|nf_GEN);
    2766           0 :     default: pari_err_FLAG("idealstar");
    2767             :   }
    2768             :   return NULL; /* LCOV_EXCL_LINE */
    2769             : }
    2770             : 
    2771             : void
    2772      100120 : check_nfelt(GEN x, GEN *den)
    2773             : {
    2774      100120 :   long l = lg(x), i;
    2775      100120 :   GEN t, d = NULL;
    2776      100120 :   if (typ(x) != t_COL) pari_err_TYPE("check_nfelt", x);
    2777      384450 :   for (i=1; i<l; i++)
    2778             :   {
    2779      284330 :     t = gel(x,i);
    2780      284330 :     switch (typ(t))
    2781             :     {
    2782      222625 :       case t_INT: break;
    2783             :       case t_FRAC:
    2784       61705 :         if (!d) d = gel(t,2); else d = lcmii(d, gel(t,2));
    2785       61705 :         break;
    2786           0 :       default: pari_err_TYPE("check_nfelt", x);
    2787             :     }
    2788             :   }
    2789      100120 :   *den = d;
    2790      100120 : }
    2791             : 
    2792             : GEN
    2793      446465 : vecmodii(GEN a, GEN b)
    2794             : {
    2795             :   long i, l;
    2796      446465 :   GEN c = cgetg_copy(a, &l);
    2797      446465 :   for (i = 1; i < l; i++) gel(c,i) = modii(gel(a,i), gel(b,i));
    2798      446465 :   return c;
    2799             : }
    2800             : 
    2801             : static GEN
    2802      152840 : ideallog_i(GEN nf, GEN x, GEN sgn, zlog_S *S)
    2803             : {
    2804      152840 :   pari_sp av = avma;
    2805             :   GEN y, cyc;
    2806      152840 :   if (!S->hU) return cgetg(1, t_COL);
    2807      151720 :   cyc = bid_get_cyc(S->bid);
    2808      151720 :   if (typ(x) == t_MAT)
    2809             :   {
    2810       42765 :     if (lg(x) == 1) return zerocol(lg(cyc)-1);
    2811       42760 :     y = famat_zlog(nf, x, sgn, S);
    2812             :   }
    2813             :   else
    2814      108955 :     y = zlog(nf, x, sgn, S);
    2815      151710 :   y = ZMV_ZCV_mul(S->U, y);
    2816      151710 :   return gerepileupto(av, vecmodii(y, cyc));
    2817             : }
    2818             : 
    2819             : /* Given x (not necessarily integral), and bid as output by zidealstarinit,
    2820             :  * compute the vector of components on the generators bid[2].
    2821             :  * Assume (x,bid) = 1 and sgn is either NULL or nfsign_arch(x, bid) */
    2822             : GEN
    2823      149580 : ideallog_sgn(GEN nf, GEN x, GEN sgn, GEN bid)
    2824             : {
    2825             :   zlog_S S;
    2826      149580 :   nf = checknf(nf); checkbid(bid);
    2827      149575 :   init_zlog(&S, bid);
    2828      149575 :   if (sgn && typ(x) == t_VEC) /* vector of elements and signatures */
    2829             :   {
    2830        3205 :     long i, l = lg(x);
    2831        3205 :     GEN y = cgetg(l, t_MAT);
    2832        3205 :     for (i = 1; i < l; i++) gel(y,i) = ideallog_i(nf, gel(x,i), gel(sgn,i), &S);
    2833        3205 :     return y;
    2834             :   }
    2835      146370 :   return ideallog_i(nf, x, sgn, &S);
    2836             : }
    2837             : GEN
    2838      151140 : ideallog(GEN nf, GEN x, GEN bid)
    2839             : {
    2840      151140 :   if (!nf) return Zideallog(bid, x);
    2841      146375 :   return ideallog_sgn(nf, x, NULL, bid);
    2842             : }
    2843             : 
    2844             : /*************************************************************************/
    2845             : /**                                                                     **/
    2846             : /**               JOIN BID STRUCTURES, IDEAL LISTS                      **/
    2847             : /**                                                                     **/
    2848             : /*************************************************************************/
    2849             : /* bid1, bid2: for coprime modules m1 and m2 (without arch. part).
    2850             :  * Output: bid for m1 m2 */
    2851             : static GEN
    2852         340 : join_bid(GEN nf, GEN bid1, GEN bid2)
    2853             : {
    2854         340 :   pari_sp av = avma;
    2855             :   long nbgen, l1,l2;
    2856             :   GEN I1,I2, G1,G2, sprk1,sprk2, cyc1,cyc2, sarch;
    2857         340 :   GEN sprk, fa,fa2, U, cyc, y, u1 = NULL, x, gen;
    2858             : 
    2859         340 :   I1 = bid_get_ideal(bid1);
    2860         340 :   I2 = bid_get_ideal(bid2);
    2861         340 :   if (gequal1(gcoeff(I1,1,1))) return bid2; /* frequent trivial case */
    2862         185 :   G1 = bid_get_grp(bid1);
    2863         185 :   G2 = bid_get_grp(bid2);
    2864         185 :   x = idealmul(nf, I1,I2);
    2865         185 :   fa = famat_mul_shallow(bid_get_fact(bid1), bid_get_fact(bid2));
    2866         185 :   fa2= famat_mul_shallow(bid_get_fact2(bid1), bid_get_fact2(bid2));
    2867         185 :   sprk1 = bid_get_sprk(bid1);
    2868         185 :   sprk2 = bid_get_sprk(bid2);
    2869         185 :   sprk = shallowconcat(sprk1, sprk2);
    2870             : 
    2871         185 :   cyc1 = abgrp_get_cyc(G1); l1 = lg(cyc1);
    2872         185 :   cyc2 = abgrp_get_cyc(G2); l2 = lg(cyc2);
    2873         185 :   gen = (lg(G1)>3 && lg(G2)>3)? gen_1: NULL;
    2874         185 :   nbgen = l1+l2-2;
    2875         185 :   cyc = ZV_snf_group(shallowconcat(cyc1,cyc2), &U, gen? &u1: NULL);
    2876         185 :   if (nbgen)
    2877             :   {
    2878         185 :     GEN U1 = bid_get_U(bid1), U2 = bid_get_U(bid2);
    2879         185 :     U1 = l1==1? const_vec(lg(sprk1), cgetg(1,t_MAT))
    2880         185 :               : ZM_ZMV_mul(vecslice(U, 1, l1-1),   U1);
    2881         185 :     U2 = l2==1? const_vec(lg(sprk2), cgetg(1,t_MAT))
    2882         185 :               : ZM_ZMV_mul(vecslice(U, l1, nbgen), U2);
    2883         185 :     U = shallowconcat(U1, U2);
    2884             :   }
    2885             :   else
    2886           0 :     U = const_vec(lg(sprk), cgetg(1,t_MAT));
    2887             : 
    2888         185 :   if (gen)
    2889             :   {
    2890         185 :     GEN uv = zkchinese1init2(nf, I2, I1, x);
    2891         370 :     gen = shallowconcat(zkVchinese1(gel(uv,1), abgrp_get_gen(G1)),
    2892         185 :                         zkVchinese1(gel(uv,2), abgrp_get_gen(G2)));
    2893             :   }
    2894         185 :   sarch = bid_get_sarch(bid1); /* trivial */
    2895         185 :   y = bid_grp(nf, u1, cyc, gen, x, sarch);
    2896         185 :   x = mkvec2(x, bid_get_arch(bid1));
    2897         185 :   y = mkvec5(x, y, mkvec2(fa, fa2), mkvec2(sprk, sarch), U);
    2898         185 :   return gerepilecopy(av,y);
    2899             : }
    2900             : 
    2901             : typedef struct _ideal_data {
    2902             :   GEN nf, emb, L, pr, prL, sgnU, archp;
    2903             : } ideal_data;
    2904             : 
    2905             : /* z <- ( z | f(v[i])_{i=1..#v} ) */
    2906             : static void
    2907       31010 : concat_join(GEN *pz, GEN v, GEN (*f)(ideal_data*,GEN), ideal_data *data)
    2908             : {
    2909       31010 :   long i, nz, lv = lg(v);
    2910             :   GEN z, Z;
    2911       62020 :   if (lv == 1) return;
    2912       13530 :   z = *pz; nz = lg(z)-1;
    2913       13530 :   *pz = Z = cgetg(lv + nz, typ(z));
    2914       13530 :   for (i = 1; i <=nz; i++) gel(Z,i) = gel(z,i);
    2915       13530 :   Z += nz;
    2916       13530 :   for (i = 1; i < lv; i++) gel(Z,i) = f(data, gel(v,i));
    2917             : }
    2918             : static GEN
    2919         340 : join_idealinit(ideal_data *D, GEN x)
    2920         340 : { return join_bid(D->nf, x, D->prL); }
    2921             : static GEN
    2922       18730 : join_ideal(ideal_data *D, GEN x)
    2923       18730 : { return idealmulpowprime(D->nf, x, D->pr, D->L); }
    2924             : static GEN
    2925         325 : join_unit(ideal_data *D, GEN x)
    2926             : {
    2927         325 :   GEN bid = join_idealinit(D, gel(x,1)), u = gel(x,2), v = mkvec(D->emb);
    2928         325 :   if (lg(u) != 1) v = shallowconcat(u, v);
    2929         325 :   return mkvec2(bid, v);
    2930             : }
    2931             : 
    2932             : /*  flag & nf_GEN : generators, otherwise no
    2933             :  *  flag &2 : units, otherwise no
    2934             :  *  flag &4 : ideals in HNF, otherwise bid
    2935             :  *  flag &8 : omit ideals which cannot be conductors (pr^1 with Npr=2) */
    2936             : static GEN
    2937         250 : Ideallist(GEN bnf, ulong bound, long flag)
    2938             : {
    2939         250 :   const long cond = flag & 8;
    2940         250 :   const long do_units = flag & 2, big_id = !(flag & 4);
    2941         250 :   const long istar_flag = (flag & nf_GEN) | nf_INIT;
    2942         250 :   pari_sp av, av0 = avma;
    2943             :   long i, j;
    2944         250 :   GEN nf, z, p, fa, id, BOUND, U, empty = cgetg(1,t_VEC);
    2945             :   forprime_t S;
    2946             :   ideal_data ID;
    2947         250 :   GEN (*join_z)(ideal_data*, GEN) =
    2948             :     do_units? &join_unit
    2949         250 :               : (big_id? &join_idealinit: &join_ideal);
    2950             : 
    2951         250 :   nf = checknf(bnf);
    2952         250 :   if ((long)bound <= 0) return empty;
    2953         250 :   id = matid(nf_get_degree(nf));
    2954         250 :   if (big_id) id = Idealstar(nf,id, istar_flag);
    2955             : 
    2956             :   /* z[i] will contain all "objects" of norm i. Depending on flag, this means
    2957             :    * an ideal, a bid, or a couple [bid, log(units)]. Such objects are stored
    2958             :    * in vectors, computed one primary component at a time; join_z
    2959             :    * reconstructs the global object */
    2960         250 :   BOUND = utoipos(bound);
    2961         250 :   z = cgetg(bound+1,t_VEC);
    2962         250 :   if (do_units) {
    2963          10 :     U = bnf_build_units(bnf);
    2964          10 :     gel(z,1) = mkvec( mkvec2(id, cgetg(1,t_VEC)) );
    2965             :   } else {
    2966         240 :     U = NULL; /* -Wall */
    2967         240 :     gel(z,1) = mkvec(id);
    2968             :   }
    2969         250 :   for (i=2; i<=(long)bound; i++) gel(z,i) = empty;
    2970         250 :   ID.nf = nf;
    2971             : 
    2972         250 :   p = cgetipos(3);
    2973         250 :   u_forprime_init(&S, 2, bound);
    2974         250 :   av = avma;
    2975        4090 :   while ((p[2] = u_forprime_next(&S)))
    2976             :   {
    2977        3590 :     if (DEBUGLEVEL>1) { err_printf("%ld ",p[2]); err_flush(); }
    2978        3590 :     fa = idealprimedec_limit_norm(nf, p, BOUND);
    2979        7195 :     for (j=1; j<lg(fa); j++)
    2980             :     {
    2981        3605 :       GEN pr = gel(fa,j), z2 = leafcopy(z);
    2982        3605 :       ulong Q, q = upr_norm(pr);
    2983        3605 :       long l = (cond && q == 2)? 2: 1;
    2984             : 
    2985        3605 :       ID.pr = ID.prL = pr;
    2986        9550 :       for (Q = q; Q <= bound; l++, Q *= q) /* add pr^l */
    2987             :       {
    2988             :         ulong iQ;
    2989        5945 :         ID.L = utoipos(l);
    2990        5945 :         if (big_id) {
    2991         155 :           ID.prL = Idealstarprk(nf, pr, l, istar_flag);
    2992         155 :           if (do_units)
    2993             :           {
    2994         140 :             GEN sprk = bid_get_sprk(ID.prL);
    2995         280 :             ID.emb = lg(sprk) == 1? cgetg(1,t_VEC)
    2996         140 :                                   : vzlog_pr(nf, U, gel(sprk,1));
    2997             :           }
    2998             :         }
    2999       36955 :         for (iQ = Q,i = 1; iQ <= bound; iQ += Q,i++)
    3000       31010 :           concat_join(&gel(z,iQ), gel(z2,i), join_z, &ID);
    3001             :       }
    3002             :     }
    3003        3590 :     if (gc_needed(av,1))
    3004             :     {
    3005           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Ideallist");
    3006           0 :       z = gerepilecopy(av, z);
    3007             :     }
    3008             :   }
    3009         250 :   return gerepilecopy(av0, z);
    3010             : }
    3011             : GEN
    3012         250 : ideallist0(GEN bnf,long bound, long flag) {
    3013         250 :   if (flag<0 || flag>15) pari_err_FLAG("ideallist");
    3014         250 :   return Ideallist(bnf,bound,flag);
    3015             : }
    3016             : GEN
    3017           0 : ideallist(GEN bnf,long bound) { return Ideallist(bnf,bound,4); }
    3018             : 
    3019             : /* bid = for module m (without arch. part), arch = archimedean part.
    3020             :  * Output: bid for [m,arch] */
    3021             : static GEN
    3022          40 : join_bid_arch(GEN nf, GEN bid, GEN archp)
    3023             : {
    3024          40 :   pari_sp av = avma;
    3025             :   GEN G, U;
    3026          40 :   GEN sprk, cyc, y, u1 = NULL, x, sarch, gen;
    3027             : 
    3028          40 :   checkbid(bid);
    3029          40 :   G = bid_get_grp(bid);
    3030          40 :   x = bid_get_ideal(bid);
    3031          40 :   sarch = nfarchstar(nf, bid_get_ideal(bid), archp);
    3032          40 :   sprk = bid_get_sprk(bid);
    3033             : 
    3034          40 :   gen = (lg(G)>3)? gel(G,3): NULL;
    3035          40 :   cyc = diagonal_shallow(shallowconcat(gel(G,2), sarch_get_cyc(sarch)));
    3036          40 :   cyc = ZM_snf_group(cyc, &U, gen? &u1: NULL);
    3037          40 :   y = bid_grp(nf, u1, cyc, gen, x, sarch);
    3038          40 :   U = split_U(U, sprk);
    3039          40 :   y = mkvec5(mkvec2(x, archp), y, gel(bid,3), mkvec2(sprk, sarch), U);
    3040          40 :   return gerepilecopy(av,y);
    3041             : }
    3042             : static GEN
    3043          40 : join_arch(ideal_data *D, GEN x) {
    3044          40 :   return join_bid_arch(D->nf, x, D->archp);
    3045             : }
    3046             : static GEN
    3047          20 : join_archunit(ideal_data *D, GEN x) {
    3048          20 :   GEN bid = join_arch(D, gel(x,1)), u = gel(x,2), v = mkvec(D->emb);
    3049          20 :   if (lg(u) != 1) v = shallowconcat(u, v);
    3050          20 :   return mkvec2(bid, v);
    3051             : }
    3052             : 
    3053             : /* L from ideallist, add archimedean part */
    3054             : GEN
    3055          10 : ideallistarch(GEN bnf, GEN L, GEN arch)
    3056             : {
    3057             :   pari_sp av;
    3058          10 :   long i, j, l = lg(L), lz;
    3059             :   GEN v, z, V;
    3060             :   ideal_data ID;
    3061             :   GEN (*join_z)(ideal_data*, GEN);
    3062             : 
    3063          10 :   if (typ(L) != t_VEC) pari_err_TYPE("ideallistarch",L);
    3064          10 :   if (l == 1) return cgetg(1,t_VEC);
    3065          10 :   z = gel(L,1);
    3066          10 :   if (typ(z) != t_VEC) pari_err_TYPE("ideallistarch",z);
    3067          10 :   z = gel(z,1); /* either a bid or [bid,U] */
    3068          10 :   ID.nf = checknf(bnf);
    3069          10 :   ID.archp = vec01_to_indices(arch);
    3070          10 :   if (lg(z) == 3) { /* the latter: do units */
    3071           5 :     if (typ(z) != t_VEC) pari_err_TYPE("ideallistarch",z);
    3072           5 :     ID.emb = zm_to_ZM( rowpermute(nfsign_units(bnf,NULL,1), ID.archp) );
    3073           5 :     join_z = &join_archunit;
    3074             :   } else
    3075           5 :     join_z = &join_arch;
    3076          10 :   av = avma; V = cgetg(l, t_VEC);
    3077          50 :   for (i = 1; i < l; i++)
    3078             :   {
    3079          40 :     z = gel(L,i); lz = lg(z);
    3080          40 :     gel(V,i) = v = cgetg(lz,t_VEC);
    3081          40 :     for (j=1; j<lz; j++) gel(v,j) = join_z(&ID, gel(z,j));
    3082             :   }
    3083          10 :   return gerepilecopy(av,V);
    3084             : }

Generated by: LCOV version 1.11