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

Generated by: LCOV version 1.13