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 25819-e703fe1174) Lines: 1853 1981 93.5 %
Date: 2020-09-18 06:10:04 Functions: 207 219 94.5 %
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    16243380 : get_tab(GEN nf, long *N)
      30             : {
      31    16243380 :   GEN tab = (typ(nf) == t_MAT)? nf: gel(nf,9);
      32    16243380 :   *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   675413918 : _mulii(GEN x, GEN y) {
      38  1078501795 :   return is_pm1(x)? (signe(x) < 0)? negi(y): y
      39  1078501595 :                   : mulii(x, y);
      40             : }
      41             : 
      42             : GEN
      43       17864 : tablemul_ei_ej(GEN M, long i, long j)
      44             : {
      45             :   long N;
      46       17864 :   GEN tab = get_tab(M, &N);
      47       17864 :   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    12825557 : zk_ei_mul(GEN nf, GEN x, long i)
      80             : {
      81             :   long j, k, N;
      82             :   GEN v, tab;
      83             : 
      84    12825557 :   if (i==1) return ZC_copy(x);
      85    12825543 :   tab = get_tab(nf, &N); tab += (i-1)*N;
      86    12825528 :   v = cgetg(N+1,t_COL);
      87    96247853 :   for (k=1; k<=N; k++)
      88             :   {
      89    83422457 :     pari_sp av = avma;
      90    83422457 :     GEN s = gen_0;
      91  1060469056 :     for (j=1; j<=N; j++)
      92             :     {
      93   977054304 :       GEN c = gcoeff(tab,k,j);
      94   977054304 :       if (signe(c)) s = addii(s, _mulii(c, gel(x,j)));
      95             :     }
      96    83414752 :     gel(v,k) = gerepileuptoint(av, s);
      97             :   }
      98    12825396 :   return v;
      99             : }
     100             : 
     101             : /* table of multiplication by wi in R[w1,..., wN] */
     102             : GEN
     103        2765 : ei_multable(GEN TAB, long i)
     104             : {
     105             :   long k,N;
     106        2765 :   GEN m, tab = get_tab(TAB, &N);
     107        2765 :   tab += (i-1)*N;
     108        2765 :   m = cgetg(N+1,t_MAT);
     109       26201 :   for (k=1; k<=N; k++) gel(m,k) = gel(tab,k);
     110        2765 :   return m;
     111             : }
     112             : 
     113             : GEN
     114     6048780 : zk_multable(GEN nf, GEN x)
     115             : {
     116     6048780 :   long i, l = lg(x);
     117     6048780 :   GEN mul = cgetg(l,t_MAT);
     118     6048742 :   gel(mul,1) = x; /* assume w_1 = 1 */
     119    18733818 :   for (i=2; i<l; i++) gel(mul,i) = zk_ei_mul(nf,x,i);
     120     6048575 :   return mul;
     121             : }
     122             : GEN
     123        2044 : multable(GEN M, GEN x)
     124             : {
     125             :   long i, N;
     126             :   GEN mul;
     127        2044 :   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     4452064 : zk_scalar_or_multable(GEN nf, GEN x)
     140             : {
     141     4452064 :   long tx = typ(x);
     142     4452064 :   if (tx == t_MAT || tx == t_INT) return x;
     143     4430287 :   x = nf_to_scalar_or_basis(nf, x);
     144     4430280 :   return (typ(x) == t_COL)? zk_multable(nf, x): x;
     145             : }
     146             : 
     147             : GEN
     148       15330 : nftrace(GEN nf, GEN x)
     149             : {
     150       15330 :   pari_sp av = avma;
     151       15330 :   nf = checknf(nf);
     152       15330 :   x = nf_to_scalar_or_basis(nf, x);
     153       15309 :   x = (typ(x) == t_COL)? RgV_dotproduct(x, gel(nf_get_Tr(nf),1))
     154       15330 :                        : gmulgs(x, nf_get_degree(nf));
     155       15330 :   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         616 :   x = (typ(x) == t_POLMOD)? rnfeltdown(rnf, gtrace(x))
     164         693 :                           : 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       51374 : nfnorm(GEN nf, GEN x)
     181             : {
     182       51374 :   pari_sp av = avma;
     183             :   GEN c, den;
     184             :   long n;
     185       51374 :   nf = checknf(nf);
     186       51374 :   n = nf_get_degree(nf);
     187       51374 :   if (typ(x) == t_MAT) return famat_norm(nf, x);
     188       51367 :   x = nf_to_scalar_or_basis(nf, x);
     189       51367 :   if (typ(x)!=t_COL)
     190       11480 :     return gerepileupto(av, gpowgs(x, n));
     191       39887 :   x = nf_to_scalar_or_alg(nf, Q_primitive_part(x, &c));
     192       39887 :   x = Q_remove_denom(x, &den);
     193       39887 :   x = ZX_resultant_all(nf_get_pol(nf), x, den, 0);
     194       39887 :   return gerepileupto(av, c ? gmul(x, gpowgs(c, n)): x);
     195             : }
     196             : 
     197             : static GEN
     198          70 : to_RgX(GEN P, long vx)
     199             : {
     200          70 :   return varn(P) == vx ? P: scalarpol_shallow(P, vx);
     201             : }
     202             : 
     203             : GEN
     204         231 : rnfeltnorm(GEN rnf, GEN x)
     205             : {
     206         231 :   pari_sp av = avma;
     207             :   GEN nf, pol;
     208         231 :   long v = rnf_get_varn(rnf);
     209         231 :   checkrnf(rnf);
     210         231 :   x = liftpol_shallow(rnfeltabstorel(rnf, x));
     211         140 :   nf = rnf_get_nf(rnf); pol = rnf_get_pol(rnf);
     212         280 :   x = (typ(x) == t_POL)
     213          70 :     ? rnfeltdown(rnf, nfX_resultant(nf,pol,to_RgX(x,v)))
     214         140 :     : gpowgs(x, rnf_get_degree(rnf));
     215         140 :   return gerepileupto(av, x);
     216             : }
     217             : 
     218             : /* x + y in nf */
     219             : GEN
     220    16725100 : nfadd(GEN nf, GEN x, GEN y)
     221             : {
     222    16725100 :   pari_sp av = avma;
     223             :   GEN z;
     224             : 
     225    16725100 :   nf = checknf(nf);
     226    16725100 :   x = nf_to_scalar_or_basis(nf, x);
     227    16725100 :   y = nf_to_scalar_or_basis(nf, y);
     228    16725100 :   if (typ(x) != t_COL)
     229    13396817 :   { z = (typ(y) == t_COL)? RgC_Rg_add(y, x): gadd(x,y); }
     230             :   else
     231     3328283 :   { z = (typ(y) == t_COL)? RgC_add(x, y): RgC_Rg_add(x, y); }
     232    16725100 :   return gerepileupto(av, z);
     233             : }
     234             : /* x - y in nf */
     235             : GEN
     236     1248464 : nfsub(GEN nf, GEN x, GEN y)
     237             : {
     238     1248464 :   pari_sp av = avma;
     239             :   GEN z;
     240             : 
     241     1248464 :   nf = checknf(nf);
     242     1248464 :   x = nf_to_scalar_or_basis(nf, x);
     243     1248464 :   y = nf_to_scalar_or_basis(nf, y);
     244     1248464 :   if (typ(x) != t_COL)
     245      894124 :   { z = (typ(y) == t_COL)? Rg_RgC_sub(x,y): gsub(x,y); }
     246             :   else
     247      354340 :   { z = (typ(y) == t_COL)? RgC_sub(x,y): RgC_Rg_sub(x,y); }
     248     1248464 :   return gerepileupto(av, z);
     249             : }
     250             : 
     251             : /* product of ZC x,y in nf; ( sum_i x_i sum_j y_j m^{i,j}_k )_k */
     252             : static GEN
     253     2087931 : nfmuli_ZC(GEN nf, GEN x, GEN y)
     254             : {
     255             :   long i, j, k, N;
     256     2087931 :   GEN TAB = get_tab(nf, &N), v = cgetg(N+1,t_COL);
     257             : 
     258    10797918 :   for (k = 1; k <= N; k++)
     259             :   {
     260     8709987 :     pari_sp av = avma;
     261     8709987 :     GEN s, TABi = TAB;
     262     8709987 :     if (k == 1)
     263     2087931 :       s = mulii(gel(x,1),gel(y,1));
     264             :     else
     265     6622056 :       s = addii(mulii(gel(x,1),gel(y,k)),
     266     6622056 :                 mulii(gel(x,k),gel(y,1)));
     267    74578421 :     for (i=2; i<=N; i++)
     268             :     {
     269    65868434 :       GEN t, xi = gel(x,i);
     270    65868434 :       TABi += N;
     271    65868434 :       if (!signe(xi)) continue;
     272             : 
     273    30427925 :       t = NULL;
     274   391030079 :       for (j=2; j<=N; j++)
     275             :       {
     276   360602154 :         GEN p1, c = gcoeff(TABi, k, j); /* m^{i,j}_k */
     277   360602154 :         if (!signe(c)) continue;
     278   124465518 :         p1 = _mulii(c, gel(y,j));
     279   124465518 :         t = t? addii(t, p1): p1;
     280             :       }
     281    30427925 :       if (t) s = addii(s, mulii(xi, t));
     282             :     }
     283     8709987 :     gel(v,k) = gerepileuptoint(av,s);
     284             :   }
     285     2087931 :   return v;
     286             : }
     287             : static int
     288    43161245 : is_famat(GEN x) { return typ(x) == t_MAT && lg(x) == 3; }
     289             : /* product of x and y in nf */
     290             : GEN
     291    22036368 : nfmul(GEN nf, GEN x, GEN y)
     292             : {
     293             :   GEN z;
     294    22036368 :   pari_sp av = avma;
     295             : 
     296    22036368 :   if (x == y) return nfsqr(nf,x);
     297             : 
     298    18857045 :   nf = checknf(nf);
     299    18857045 :   if (is_famat(x) || is_famat(y)) return famat_mul(x, y);
     300    18856954 :   x = nf_to_scalar_or_basis(nf, x);
     301    18856954 :   y = nf_to_scalar_or_basis(nf, y);
     302    18856954 :   if (typ(x) != t_COL)
     303             :   {
     304    14293691 :     if (isintzero(x)) return gen_0;
     305    10033309 :     z = (typ(y) == t_COL)? RgC_Rg_mul(y, x): gmul(x,y); }
     306             :   else
     307             :   {
     308     4563263 :     if (typ(y) != t_COL)
     309             :     {
     310     3181479 :       if (isintzero(y)) return gen_0;
     311      695870 :       z = RgC_Rg_mul(x, y);
     312             :     }
     313             :     else
     314             :     {
     315             :       GEN dx, dy;
     316     1381784 :       x = Q_remove_denom(x, &dx);
     317     1381784 :       y = Q_remove_denom(y, &dy);
     318     1381784 :       z = nfmuli_ZC(nf,x,y);
     319     1381784 :       dx = mul_denom(dx,dy);
     320     1381784 :       if (dx) z = ZC_Z_div(z, dx);
     321             :     }
     322             :   }
     323    12110963 :   return gerepileupto(av, z);
     324             : }
     325             : /* square of ZC x in nf */
     326             : static GEN
     327     1298551 : nfsqri_ZC(GEN nf, GEN x)
     328             : {
     329             :   long i, j, k, N;
     330     1298551 :   GEN TAB = get_tab(nf, &N), v = cgetg(N+1,t_COL);
     331             : 
     332     9896824 :   for (k = 1; k <= N; k++)
     333             :   {
     334     8598273 :     pari_sp av = avma;
     335     8598273 :     GEN s, TABi = TAB;
     336     8598273 :     if (k == 1)
     337     1298551 :       s = sqri(gel(x,1));
     338             :     else
     339     7299722 :       s = shifti(mulii(gel(x,1),gel(x,k)), 1);
     340   105674249 :     for (i=2; i<=N; i++)
     341             :     {
     342    97075976 :       GEN p1, c, t, xi = gel(x,i);
     343    97075976 :       TABi += N;
     344    97075976 :       if (!signe(xi)) continue;
     345             : 
     346    32882777 :       c = gcoeff(TABi, k, i);
     347    32882777 :       t = signe(c)? _mulii(c,xi): NULL;
     348   398871639 :       for (j=i+1; j<=N; j++)
     349             :       {
     350   365988862 :         c = gcoeff(TABi, k, j);
     351   365988862 :         if (!signe(c)) continue;
     352   166461926 :         p1 = _mulii(c, shifti(gel(x,j),1));
     353   166461926 :         t = t? addii(t, p1): p1;
     354             :       }
     355    32882777 :       if (t) s = addii(s, mulii(xi, t));
     356             :     }
     357     8598273 :     gel(v,k) = gerepileuptoint(av,s);
     358             :   }
     359     1298551 :   return v;
     360             : }
     361             : /* square of x in nf */
     362             : GEN
     363     5050843 : nfsqr(GEN nf, GEN x)
     364             : {
     365     5050843 :   pari_sp av = avma;
     366             :   GEN z;
     367             : 
     368     5050843 :   nf = checknf(nf);
     369     5050843 :   if (is_famat(x)) return famat_sqr(x);
     370     5050843 :   x = nf_to_scalar_or_basis(nf, x);
     371     5050843 :   if (typ(x) != t_COL) z = gsqr(x);
     372             :   else
     373             :   {
     374             :     GEN dx;
     375       92876 :     x = Q_remove_denom(x, &dx);
     376       92876 :     z = nfsqri_ZC(nf,x);
     377       92876 :     if (dx) z = RgC_Rg_div(z, sqri(dx));
     378             :   }
     379     5050843 :   return gerepileupto(av, z);
     380             : }
     381             : 
     382             : /* x a ZC, v a t_COL of ZC/Z */
     383             : GEN
     384      134587 : zkC_multable_mul(GEN v, GEN x)
     385             : {
     386      134587 :   long i, l = lg(v);
     387      134587 :   GEN y = cgetg(l, t_COL);
     388      472034 :   for (i = 1; i < l; i++)
     389             :   {
     390      337447 :     GEN c = gel(v,i);
     391      337447 :     if (typ(c)!=t_COL) {
     392           0 :       if (!isintzero(c)) c = ZC_Z_mul(gel(x,1), c);
     393             :     } else {
     394      337447 :       c = ZM_ZC_mul(x,c);
     395      337447 :       if (ZV_isscalar(c)) c = gel(c,1);
     396             :     }
     397      337447 :     gel(y,i) = c;
     398             :   }
     399      134587 :   return y;
     400             : }
     401             : 
     402             : GEN
     403       46585 : nfC_multable_mul(GEN v, GEN x)
     404             : {
     405       46585 :   long i, l = lg(v);
     406       46585 :   GEN y = cgetg(l, t_COL);
     407      323400 :   for (i = 1; i < l; i++)
     408             :   {
     409      276815 :     GEN c = gel(v,i);
     410      276815 :     if (typ(c)!=t_COL) {
     411      230062 :       if (!isintzero(c)) c = RgC_Rg_mul(gel(x,1), c);
     412             :     } else {
     413       46753 :       c = RgM_RgC_mul(x,c);
     414       46753 :       if (QV_isscalar(c)) c = gel(c,1);
     415             :     }
     416      276815 :     gel(y,i) = c;
     417             :   }
     418       46585 :   return y;
     419             : }
     420             : 
     421             : GEN
     422      168028 : nfC_nf_mul(GEN nf, GEN v, GEN x)
     423             : {
     424             :   long tx;
     425             :   GEN y;
     426             : 
     427      168028 :   x = nf_to_scalar_or_basis(nf, x);
     428      168028 :   tx = typ(x);
     429      168028 :   if (tx != t_COL)
     430             :   {
     431             :     long l, i;
     432      128387 :     if (tx == t_INT)
     433             :     {
     434      120001 :       long s = signe(x);
     435      120001 :       if (!s) return zerocol(lg(v)-1);
     436      113225 :       if (is_pm1(x)) return s > 0? leafcopy(v): RgC_neg(v);
     437             :     }
     438       42000 :     l = lg(v); y = cgetg(l, t_COL);
     439      307650 :     for (i=1; i < l; i++)
     440             :     {
     441      265650 :       GEN c = gel(v,i);
     442      265650 :       if (typ(c) != t_COL) c = gmul(c, x); else c = RgC_Rg_mul(c, x);
     443      265650 :       gel(y,i) = c;
     444             :     }
     445       42000 :     return y;
     446             :   }
     447             :   else
     448             :   {
     449             :     GEN dx;
     450       39641 :     x = zk_multable(nf, Q_remove_denom(x,&dx));
     451       39641 :     y = nfC_multable_mul(v, x);
     452       39641 :     return dx? RgC_Rg_div(y, dx): y;
     453             :   }
     454             : }
     455             : static GEN
     456        9296 : mulbytab(GEN M, GEN c)
     457        9296 : { return typ(c) == t_COL? RgM_RgC_mul(M,c): RgC_Rg_mul(gel(M,1), c); }
     458             : GEN
     459        2044 : tablemulvec(GEN M, GEN x, GEN v)
     460             : {
     461             :   long l, i;
     462             :   GEN y;
     463             : 
     464        2044 :   if (typ(x) == t_COL && RgV_isscalar(x))
     465             :   {
     466           0 :     x = gel(x,1);
     467           0 :     return typ(v) == t_POL? RgX_Rg_mul(v,x): RgV_Rg_mul(v,x);
     468             :   }
     469        2044 :   x = multable(M, x); /* multiplication table by x */
     470        2044 :   y = cgetg_copy(v, &l);
     471        2044 :   if (typ(v) == t_POL)
     472             :   {
     473        2044 :     y[1] = v[1];
     474       11340 :     for (i=2; i < l; i++) gel(y,i) = mulbytab(x, gel(v,i));
     475        2044 :     y = normalizepol(y);
     476             :   }
     477             :   else
     478             :   {
     479           0 :     for (i=1; i < l; i++) gel(y,i) = mulbytab(x, gel(v,i));
     480             :   }
     481        2044 :   return y;
     482             : }
     483             : 
     484             : GEN
     485      292125 : zkmultable_capZ(GEN mx) { return Q_denom(zkmultable_inv(mx)); }
     486             : GEN
     487      348966 : zkmultable_inv(GEN mx) { return ZM_gauss(mx, col_ei(lg(mx)-1,1)); }
     488             : /* nf a true nf, x a ZC */
     489             : GEN
     490       56841 : zk_inv(GEN nf, GEN x) { return zkmultable_inv(zk_multable(nf,x)); }
     491             : 
     492             : /* inverse of x in nf */
     493             : GEN
     494       68068 : nfinv(GEN nf, GEN x)
     495             : {
     496       68068 :   pari_sp av = avma;
     497             :   GEN z;
     498             : 
     499       68068 :   nf = checknf(nf);
     500       68068 :   if (is_famat(x)) return famat_inv(x);
     501       68068 :   x = nf_to_scalar_or_basis(nf, x);
     502       68068 :   if (typ(x) == t_COL)
     503             :   {
     504             :     GEN d;
     505       26831 :     x = Q_remove_denom(x, &d);
     506       26831 :     z = zk_inv(nf, x);
     507       26831 :     if (d) z = RgC_Rg_mul(z, d);
     508             :   }
     509             :   else
     510       41237 :     z = ginv(x);
     511       68068 :   return gerepileupto(av, z);
     512             : }
     513             : 
     514             : /* quotient of x and y in nf */
     515             : GEN
     516       27214 : nfdiv(GEN nf, GEN x, GEN y)
     517             : {
     518       27214 :   pari_sp av = avma;
     519             :   GEN z;
     520             : 
     521       27214 :   nf = checknf(nf);
     522       27214 :   if (is_famat(x) || is_famat(y)) return famat_div(x,y);
     523       27123 :   y = nf_to_scalar_or_basis(nf, y);
     524       27123 :   if (typ(y) != t_COL)
     525             :   {
     526       19747 :     x = nf_to_scalar_or_basis(nf, x);
     527       19747 :     z = (typ(x) == t_COL)? RgC_Rg_div(x, y): gdiv(x,y);
     528             :   }
     529             :   else
     530             :   {
     531             :     GEN d;
     532        7376 :     y = Q_remove_denom(y, &d);
     533        7376 :     z = nfmul(nf, x, zk_inv(nf,y));
     534        7376 :     if (d) z = RgC_Rg_mul(z, d);
     535             :   }
     536       27123 :   return gerepileupto(av, z);
     537             : }
     538             : 
     539             : /* product of INTEGERS (t_INT or ZC) x and y in nf */
     540             : GEN
     541      843538 : nfmuli(GEN nf, GEN x, GEN y)
     542             : {
     543      843538 :   if (typ(x) == t_INT) return (typ(y) == t_COL)? ZC_Z_mul(y, x): mulii(x,y);
     544      732295 :   if (typ(y) == t_INT) return ZC_Z_mul(x, y);
     545      706147 :   return nfmuli_ZC(nf, x, y);
     546             : }
     547             : GEN
     548     1205675 : nfsqri(GEN nf, GEN x)
     549     1205675 : { return (typ(x) == t_INT)? sqri(x): nfsqri_ZC(nf, x); }
     550             : 
     551             : /* both x and y are RgV */
     552             : GEN
     553           0 : tablemul(GEN TAB, GEN x, GEN y)
     554             : {
     555             :   long i, j, k, N;
     556             :   GEN s, v;
     557           0 :   if (typ(x) != t_COL) return gmul(x, y);
     558           0 :   if (typ(y) != t_COL) return gmul(y, x);
     559           0 :   N = lg(x)-1;
     560           0 :   v = cgetg(N+1,t_COL);
     561           0 :   for (k=1; k<=N; k++)
     562             :   {
     563           0 :     pari_sp av = avma;
     564           0 :     GEN TABi = TAB;
     565           0 :     if (k == 1)
     566           0 :       s = gmul(gel(x,1),gel(y,1));
     567             :     else
     568           0 :       s = gadd(gmul(gel(x,1),gel(y,k)),
     569           0 :                gmul(gel(x,k),gel(y,1)));
     570           0 :     for (i=2; i<=N; i++)
     571             :     {
     572           0 :       GEN t, xi = gel(x,i);
     573           0 :       TABi += N;
     574           0 :       if (gequal0(xi)) continue;
     575             : 
     576           0 :       t = NULL;
     577           0 :       for (j=2; j<=N; j++)
     578             :       {
     579           0 :         GEN p1, c = gcoeff(TABi, k, j); /* m^{i,j}_k */
     580           0 :         if (gequal0(c)) continue;
     581           0 :         p1 = gmul(c, gel(y,j));
     582           0 :         t = t? gadd(t, p1): p1;
     583             :       }
     584           0 :       if (t) s = gadd(s, gmul(xi, t));
     585             :     }
     586           0 :     gel(v,k) = gerepileupto(av,s);
     587             :   }
     588           0 :   return v;
     589             : }
     590             : GEN
     591       42091 : tablesqr(GEN TAB, GEN x)
     592             : {
     593             :   long i, j, k, N;
     594             :   GEN s, v;
     595             : 
     596       42091 :   if (typ(x) != t_COL) return gsqr(x);
     597       42091 :   N = lg(x)-1;
     598       42091 :   v = cgetg(N+1,t_COL);
     599             : 
     600      298795 :   for (k=1; k<=N; k++)
     601             :   {
     602      256704 :     pari_sp av = avma;
     603      256704 :     GEN TABi = TAB;
     604      256704 :     if (k == 1)
     605       42091 :       s = gsqr(gel(x,1));
     606             :     else
     607      214613 :       s = gmul2n(gmul(gel(x,1),gel(x,k)), 1);
     608     1611036 :     for (i=2; i<=N; i++)
     609             :     {
     610     1354332 :       GEN p1, c, t, xi = gel(x,i);
     611     1354332 :       TABi += N;
     612     1354332 :       if (gequal0(xi)) continue;
     613             : 
     614      372918 :       c = gcoeff(TABi, k, i);
     615      372918 :       t = !gequal0(c)? gmul(c,xi): NULL;
     616     1497363 :       for (j=i+1; j<=N; j++)
     617             :       {
     618     1124445 :         c = gcoeff(TABi, k, j);
     619     1124445 :         if (gequal0(c)) continue;
     620      586551 :         p1 = gmul(gmul2n(c,1), gel(x,j));
     621      586551 :         t = t? gadd(t, p1): p1;
     622             :       }
     623      372918 :       if (t) s = gadd(s, gmul(xi, t));
     624             :     }
     625      256704 :     gel(v,k) = gerepileupto(av,s);
     626             :   }
     627       42091 :   return v;
     628             : }
     629             : 
     630             : static GEN
     631      168434 : _mul(void *data, GEN x, GEN y) { return nfmuli((GEN)data,x,y); }
     632             : static GEN
     633      407485 : _sqr(void *data, GEN x) { return nfsqri((GEN)data,x); }
     634             : 
     635             : /* Compute z^n in nf, left-shift binary powering */
     636             : GEN
     637      273760 : nfpow(GEN nf, GEN z, GEN n)
     638             : {
     639      273760 :   pari_sp av = avma;
     640             :   long s;
     641             :   GEN x, cx;
     642             : 
     643      273760 :   if (typ(n)!=t_INT) pari_err_TYPE("nfpow",n);
     644      273760 :   nf = checknf(nf);
     645      273760 :   s = signe(n); if (!s) return gen_1;
     646      273760 :   if (is_famat(z)) return famat_pow(z, n);
     647      273746 :   x = nf_to_scalar_or_basis(nf, z);
     648      273746 :   if (typ(x) != t_COL) return powgi(x,n);
     649      269446 :   if (s < 0)
     650             :   { /* simplified nfinv */
     651             :     GEN d;
     652       16314 :     x = Q_remove_denom(x, &d);
     653       16314 :     x = zk_inv(nf, x);
     654       16314 :     x = primitive_part(x, &cx);
     655       16314 :     cx = mul_content(cx, d);
     656       16314 :     n = negi(n);
     657             :   }
     658             :   else
     659      253132 :     x = primitive_part(x, &cx);
     660      269446 :   x = gen_pow_i(x, n, (void*)nf, _sqr, _mul);
     661      269446 :   if (cx)
     662       48218 :     x = gerepileupto(av, gmul(x, powgi(cx, n)));
     663             :   else
     664      221228 :     x = gerepilecopy(av, x);
     665      269446 :   return x;
     666             : }
     667             : /* Compute z^n in nf, left-shift binary powering */
     668             : GEN
     669       60200 : nfpow_u(GEN nf, GEN z, ulong n)
     670             : {
     671       60200 :   pari_sp av = avma;
     672             :   GEN x, cx;
     673             : 
     674       60200 :   nf = checknf(nf);
     675       60200 :   if (!n) return gen_1;
     676       60200 :   x = nf_to_scalar_or_basis(nf, z);
     677       60200 :   if (typ(x) != t_COL) return gpowgs(x,n);
     678       31626 :   x = primitive_part(x, &cx);
     679       31626 :   x = gen_powu_i(x, n, (void*)nf, _sqr, _mul);
     680       31626 :   if (cx)
     681             :   {
     682        1736 :     x = gmul(x, powgi(cx, utoipos(n)));
     683        1736 :     return gerepileupto(av,x);
     684             :   }
     685       29890 :   return gerepilecopy(av, x);
     686             : }
     687             : 
     688             : static GEN
     689     2570687 : _nf_red(void *E, GEN x) { (void)E; return x; }
     690             : 
     691             : static GEN
     692    10367140 : _nf_add(void *E, GEN x, GEN y) { return nfadd((GEN)E,x,y); }
     693             : 
     694             : static GEN
     695      623763 : _nf_neg(void *E, GEN x) { (void)E; return gneg(x); }
     696             : 
     697             : static GEN
     698    12473391 : _nf_mul(void *E, GEN x, GEN y) { return nfmul((GEN)E,x,y); }
     699             : 
     700             : static GEN
     701       43218 : _nf_inv(void *E, GEN x) { return nfinv((GEN)E,x); }
     702             : 
     703             : static GEN
     704        8498 : _nf_s(void *E, long x) { (void)E; return stoi(x); }
     705             : 
     706             : static const struct bb_field nf_field={_nf_red,_nf_add,_nf_mul,_nf_neg,
     707             :                                         _nf_inv,&gequal0,_nf_s };
     708             : 
     709      186984 : const struct bb_field *get_nf_field(void **E, GEN nf)
     710      186984 : { *E = (void*)nf; return &nf_field; }
     711             : 
     712             : GEN
     713          14 : nfM_det(GEN nf, GEN M)
     714             : {
     715             :   void *E;
     716          14 :   const struct bb_field *S = get_nf_field(&E, nf);
     717          14 :   return gen_det(M, E, S);
     718             : }
     719             : GEN
     720        8484 : nfM_inv(GEN nf, GEN M)
     721             : {
     722             :   void *E;
     723        8484 :   const struct bb_field *S = get_nf_field(&E, nf);
     724        8484 :   return gen_Gauss(M, matid(lg(M)-1), E, S);
     725             : }
     726             : GEN
     727        8232 : nfM_mul(GEN nf, GEN A, GEN B)
     728             : {
     729             :   void *E;
     730        8232 :   const struct bb_field *S = get_nf_field(&E, nf);
     731        8232 :   return gen_matmul(A, B, E, S);
     732             : }
     733             : GEN
     734      170254 : nfM_nfC_mul(GEN nf, GEN A, GEN B)
     735             : {
     736             :   void *E;
     737      170254 :   const struct bb_field *S = get_nf_field(&E, nf);
     738      170254 :   return gen_matcolmul(A, B, E, S);
     739             : }
     740             : 
     741             : /* valuation of integral x (ZV), with resp. to prime ideal pr */
     742             : long
     743    19475371 : ZC_nfvalrem(GEN x, GEN pr, GEN *newx)
     744             : {
     745             :   long i, v, l;
     746    19475371 :   GEN r, y, p = pr_get_p(pr), mul = pr_get_tau(pr);
     747             : 
     748             :   /* p inert */
     749    19475380 :   if (typ(mul) == t_INT) return newx? ZV_pvalrem(x, p, newx):ZV_pval(x, p);
     750    19400108 :   y = cgetg_copy(x, &l); /* will hold the new x */
     751    19400115 :   x = leafcopy(x);
     752    19400120 :   for(v=0;; v++)
     753             :   {
     754    87231687 :     for (i=1; i<l; i++)
     755             :     { /* is (x.b)[i] divisible by p ? */
     756    79583095 :       gel(y,i) = dvmdii(ZMrow_ZC_mul(mul,x,i),p,&r);
     757    79583081 :       if (r != gen_0) { if (newx) *newx = x; return v; }
     758             :     }
     759     7648592 :     swap(x, y);
     760             :   }
     761             : }
     762             : long
     763    18870914 : ZC_nfval(GEN x, GEN P)
     764    18870914 : { return ZC_nfvalrem(x, P, NULL); }
     765             : 
     766             : /* v_P(x) != 0, x a ZV. Simpler version of ZC_nfvalrem */
     767             : int
     768      480519 : ZC_prdvd(GEN x, GEN P)
     769             : {
     770      480519 :   pari_sp av = avma;
     771             :   long i, l;
     772      480519 :   GEN p = pr_get_p(P), mul = pr_get_tau(P);
     773      480519 :   if (typ(mul) == t_INT) return ZV_Z_dvd(x, p);
     774      480176 :   l = lg(x);
     775     1928083 :   for (i=1; i<l; i++)
     776     1721905 :     if (!dvdii(ZMrow_ZC_mul(mul,x,i), p)) return gc_bool(av,0);
     777      206178 :   return gc_bool(av,1);
     778             : }
     779             : 
     780             : int
     781          28 : pr_equal(GEN P, GEN Q)
     782             : {
     783          28 :   GEN gQ, p = pr_get_p(P);
     784          28 :   long e = pr_get_e(P), f = pr_get_f(P), n;
     785          28 :   if (!equalii(p, pr_get_p(Q)) || e != pr_get_e(Q) || f != pr_get_f(Q))
     786          14 :     return 0;
     787          14 :   gQ = pr_get_gen(Q); n = lg(gQ)-1;
     788          14 :   if (2*e*f > n) return 1; /* room for only one such pr */
     789           7 :   return ZV_equal(pr_get_gen(P), gQ) || ZC_prdvd(gQ, P);
     790             : }
     791             : 
     792             : GEN
     793        1015 : famat_nfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
     794             : {
     795        1015 :   pari_sp av = avma;
     796        1015 :   GEN P = gel(x,1), E = gel(x,2), V = gen_0, y = NULL;
     797        1015 :   long l = lg(P), simplify = 0, i;
     798        1015 :   if (py) { *py = gen_1; y = cgetg(l, t_COL); }
     799             : 
     800      293405 :   for (i = 1; i < l; i++)
     801             :   {
     802      292390 :     GEN e = gel(E,i);
     803             :     long v;
     804      292390 :     if (!signe(e))
     805             :     {
     806           7 :       if (py) gel(y,i) = gen_1;
     807           7 :       simplify = 1; continue;
     808             :     }
     809      292383 :     v = nfvalrem(nf, gel(P,i), pr, py? &gel(y,i): NULL);
     810      292383 :     if (v == LONG_MAX) { set_avma(av); if (py) *py = gen_0; return mkoo(); }
     811      292383 :     V = addmulii(V, stoi(v), e);
     812             :   }
     813        1015 :   if (!py) V = gerepileuptoint(av, V);
     814             :   else
     815             :   {
     816          35 :     y = mkmat2(y, gel(x,2));
     817          35 :     if (simplify) y = famat_remove_trivial(y);
     818          35 :     gerepileall(av, 2, &V, &y); *py = y;
     819             :   }
     820        1015 :   return V;
     821             : }
     822             : long
     823     1603771 : nfval(GEN nf, GEN x, GEN pr)
     824             : {
     825     1603771 :   pari_sp av = avma;
     826             :   long w, e;
     827             :   GEN cx, p;
     828             : 
     829     1603771 :   if (gequal0(x)) return LONG_MAX;
     830     1602188 :   nf = checknf(nf);
     831     1602189 :   checkprid(pr);
     832     1602191 :   p = pr_get_p(pr);
     833     1602190 :   e = pr_get_e(pr);
     834     1602194 :   x = nf_to_scalar_or_basis(nf, x);
     835     1602162 :   if (typ(x) != t_COL) return e*Q_pval(x,p);
     836      374260 :   x = Q_primitive_part(x, &cx);
     837      374237 :   w = ZC_nfval(x,pr);
     838      374247 :   if (cx) w += e*Q_pval(cx,p);
     839      374256 :   return gc_long(av,w);
     840             : }
     841             : 
     842             : /* want to write p^v = uniformizer^(e*v) * z^v, z coprime to pr */
     843             : /* z := tau^e / p^(e-1), algebraic integer coprime to pr; return z^v */
     844             : static GEN
     845       22064 : powp(GEN nf, GEN pr, long v)
     846             : {
     847             :   GEN b, z;
     848             :   long e;
     849       22064 :   if (!v) return gen_1;
     850       21938 :   b = pr_get_tau(pr);
     851       21938 :   if (typ(b) == t_INT) return gen_1;
     852        3234 :   e = pr_get_e(pr);
     853        3234 :   z = gel(b,1);
     854        3234 :   if (e != 1) z = gdiv(nfpow_u(nf, z, e), powiu(pr_get_p(pr),e-1));
     855        3234 :   if (v < 0) { v = -v; z = nfinv(nf, z); }
     856        3234 :   if (v != 1) z = nfpow_u(nf, z, v);
     857        3234 :   return z;
     858             : }
     859             : long
     860      356503 : nfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
     861             : {
     862      356503 :   pari_sp av = avma;
     863             :   long w, e;
     864             :   GEN cx, p, t;
     865             : 
     866      356503 :   if (!py) return nfval(nf,x,pr);
     867       64064 :   if (gequal0(x)) { *py = gen_0; return LONG_MAX; }
     868       64008 :   nf = checknf(nf);
     869       64008 :   checkprid(pr);
     870       64008 :   p = pr_get_p(pr);
     871       64008 :   e = pr_get_e(pr);
     872       64008 :   x = nf_to_scalar_or_basis(nf, x);
     873       64008 :   if (typ(x) != t_COL) {
     874       52808 :     w = Q_pvalrem(x,p, py);
     875       52808 :     if (!w) { *py = gerepilecopy(av, x); return 0; }
     876       18963 :     *py = gerepileupto(av, gmul(powp(nf, pr, w), *py));
     877       18963 :     return e*w;
     878             :   }
     879       11200 :   x = Q_primitive_part(x, &cx);
     880       11200 :   w = ZC_nfvalrem(x,pr, py);
     881       11200 :   if (cx)
     882             :   {
     883        3101 :     long v = Q_pvalrem(cx,p, &t);
     884        3101 :     *py = nfmul(nf, *py, gmul(powp(nf,pr,v), t));
     885        3101 :     *py = gerepileupto(av, *py);
     886        3101 :     w += e*v;
     887             :   }
     888             :   else
     889        8099 :     *py = gerepilecopy(av, *py);
     890       11200 :   return w;
     891             : }
     892             : GEN
     893         154 : gpnfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
     894             : {
     895             :   long v;
     896         154 :   if (is_famat(x)) return famat_nfvalrem(nf, x, pr, py);
     897         147 :   v = nfvalrem(nf,x,pr,py);
     898         147 :   return v == LONG_MAX? mkoo(): stoi(v);
     899             : }
     900             : 
     901             : /* true nf */
     902             : GEN
     903       98453 : coltoalg(GEN nf, GEN x)
     904             : {
     905       98453 :   return mkpolmod( nf_to_scalar_or_alg(nf, x), nf_get_pol(nf) );
     906             : }
     907             : 
     908             : GEN
     909      144786 : basistoalg(GEN nf, GEN x)
     910             : {
     911             :   GEN T;
     912             : 
     913      144786 :   nf = checknf(nf);
     914      144786 :   switch(typ(x))
     915             :   {
     916       92314 :     case t_COL: {
     917       92314 :       pari_sp av = avma;
     918       92314 :       return gerepilecopy(av, coltoalg(nf, x));
     919             :     }
     920       32039 :     case t_POLMOD:
     921       32039 :       T = nf_get_pol(nf);
     922       32039 :       if (!RgX_equal_var(T,gel(x,1)))
     923           0 :         pari_err_MODULUS("basistoalg", T,gel(x,1));
     924       32039 :       return gcopy(x);
     925        1848 :     case t_POL:
     926        1848 :       T = nf_get_pol(nf);
     927        1848 :       if (varn(T) != varn(x)) pari_err_VAR("basistoalg",x,T);
     928        1848 :       retmkpolmod(RgX_rem(x, T), ZX_copy(T));
     929       18585 :     case t_INT:
     930             :     case t_FRAC:
     931       18585 :       T = nf_get_pol(nf);
     932       18585 :       retmkpolmod(gcopy(x), ZX_copy(T));
     933           0 :     default:
     934           0 :       pari_err_TYPE("basistoalg",x);
     935             :       return NULL; /* LCOV_EXCL_LINE */
     936             :   }
     937             : }
     938             : 
     939             : /* true nf, x a t_POL */
     940             : static GEN
     941     1626030 : pol_to_scalar_or_basis(GEN nf, GEN x)
     942             : {
     943     1626030 :   GEN T = nf_get_pol(nf);
     944     1626029 :   long l = lg(x);
     945     1626029 :   if (varn(x) != varn(T)) pari_err_VAR("nf_to_scalar_or_basis", x,T);
     946     1625969 :   if (l >= lg(T)) { x = RgX_rem(x, T); l = lg(x); }
     947     1625969 :   if (l == 2) return gen_0;
     948      990488 :   if (l == 3)
     949             :   {
     950      215765 :     x = gel(x,2);
     951      215765 :     if (!is_rational_t(typ(x))) pari_err_TYPE("nf_to_scalar_or_basis",x);
     952      215759 :     return x;
     953             :   }
     954      774723 :   return poltobasis(nf,x);
     955             : }
     956             : /* Assume nf is a genuine nf. */
     957             : GEN
     958    91347166 : nf_to_scalar_or_basis(GEN nf, GEN x)
     959             : {
     960    91347166 :   switch(typ(x))
     961             :   {
     962    66418514 :     case t_INT: case t_FRAC:
     963    66418514 :       return x;
     964      185445 :     case t_POLMOD:
     965      185445 :       x = checknfelt_mod(nf,x,"nf_to_scalar_or_basis");
     966      185372 :       switch(typ(x))
     967             :       {
     968       34188 :         case t_INT: case t_FRAC: return x;
     969      151184 :         case t_POL: return pol_to_scalar_or_basis(nf,x);
     970             :       }
     971           0 :       break;
     972     1474852 :     case t_POL: return pol_to_scalar_or_basis(nf,x);
     973    23268376 :     case t_COL:
     974    23268376 :       if (lg(x)-1 != nf_get_degree(nf)) break;
     975    23268311 :       return QV_isscalar(x)? gel(x,1): x;
     976             :   }
     977          54 :   pari_err_TYPE("nf_to_scalar_or_basis",x);
     978             :   return NULL; /* LCOV_EXCL_LINE */
     979             : }
     980             : /* Let x be a polynomial with coefficients in Q or nf. Return the same
     981             :  * polynomial with coefficients expressed as vectors (on the integral basis).
     982             :  * No consistency checks, not memory-clean. */
     983             : GEN
     984        8433 : RgX_to_nfX(GEN nf, GEN x)
     985             : {
     986             :   long i, l;
     987        8433 :   GEN y = cgetg_copy(x, &l); y[1] = x[1];
     988       66084 :   for (i=2; i<l; i++) gel(y,i) = nf_to_scalar_or_basis(nf, gel(x,i));
     989        8433 :   return y;
     990             : }
     991             : 
     992             : /* Assume nf is a genuine nf. */
     993             : GEN
     994     1206299 : nf_to_scalar_or_alg(GEN nf, GEN x)
     995             : {
     996     1206299 :   switch(typ(x))
     997             :   {
     998       28092 :     case t_INT: case t_FRAC:
     999       28092 :       return x;
    1000           0 :     case t_POLMOD:
    1001           0 :       x = checknfelt_mod(nf,x,"nf_to_scalar_or_alg");
    1002           0 :       if (typ(x) != t_POL) return x;
    1003             :       /* fall through */
    1004             :     case t_POL:
    1005             :     {
    1006        2044 :       GEN T = nf_get_pol(nf);
    1007        2044 :       long l = lg(x);
    1008        2044 :       if (varn(x) != varn(T)) pari_err_VAR("nf_to_scalar_or_alg", x,T);
    1009        2044 :       if (l >= lg(T)) { x = RgX_rem(x, T); l = lg(x); }
    1010        2044 :       if (l == 2) return gen_0;
    1011        2044 :       if (l == 3) return gel(x,2);
    1012        2030 :       return x;
    1013             :     }
    1014     1176116 :     case t_COL:
    1015             :     {
    1016             :       GEN dx;
    1017     1176116 :       if (lg(x)-1 != nf_get_degree(nf)) break;
    1018     2316478 :       if (QV_isscalar(x)) return gel(x,1);
    1019     1140370 :       x = Q_remove_denom(x, &dx);
    1020     1140361 :       x = RgV_RgC_mul(nf_get_zkprimpart(nf), x);
    1021     1140381 :       dx = mul_denom(dx, nf_get_zkden(nf));
    1022     1140378 :       return gdiv(x,dx);
    1023             :     }
    1024             :   }
    1025          47 :   pari_err_TYPE("nf_to_scalar_or_alg",x);
    1026             :   return NULL; /* LCOV_EXCL_LINE */
    1027             : }
    1028             : 
    1029             : /* gmul(A, RgX_to_RgC(x)), A t_MAT of compatible dimensions */
    1030             : GEN
    1031        1337 : RgM_RgX_mul(GEN A, GEN x)
    1032             : {
    1033        1337 :   long i, l = lg(x)-1;
    1034             :   GEN z;
    1035        1337 :   if (l == 1) return zerocol(nbrows(A));
    1036        1337 :   z = gmul(gel(x,2), gel(A,1));
    1037        2541 :   for (i = 2; i < l; i++)
    1038        1204 :     if (!gequal0(gel(x,i+1))) z = gadd(z, gmul(gel(x,i+1), gel(A,i)));
    1039        1337 :   return z;
    1040             : }
    1041             : GEN
    1042     3125560 : ZM_ZX_mul(GEN A, GEN x)
    1043             : {
    1044     3125560 :   long i, l = lg(x)-1;
    1045             :   GEN z;
    1046     3125560 :   if (l == 1) return zerocol(nbrows(A));
    1047     3124709 :   z = ZC_Z_mul(gel(A,1), gel(x,2));
    1048    13370814 :   for (i = 2; i < l ; i++)
    1049    10246677 :     if (signe(gel(x,i+1))) z = ZC_add(z, ZC_Z_mul(gel(A,i), gel(x,i+1)));
    1050     3124137 :   return z;
    1051             : }
    1052             : /* x a t_POL, nf a genuine nf. No garbage collecting. No check.  */
    1053             : GEN
    1054     2856974 : poltobasis(GEN nf, GEN x)
    1055             : {
    1056     2856974 :   GEN d, T = nf_get_pol(nf);
    1057     2856950 :   if (varn(x) != varn(T)) pari_err_VAR( "poltobasis", x,T);
    1058     2856817 :   if (degpol(x) >= degpol(T)) x = RgX_rem(x,T);
    1059     2856797 :   x = Q_remove_denom(x, &d);
    1060     2856790 :   if (!RgX_is_ZX(x)) pari_err_TYPE("poltobasis",x);
    1061     2856758 :   x = ZM_ZX_mul(nf_get_invzk(nf), x);
    1062     2856368 :   if (d) x = RgC_Rg_div(x, d);
    1063     2856372 :   return x;
    1064             : }
    1065             : 
    1066             : GEN
    1067      302624 : algtobasis(GEN nf, GEN x)
    1068             : {
    1069             :   pari_sp av;
    1070             : 
    1071      302624 :   nf = checknf(nf);
    1072      302624 :   switch(typ(x))
    1073             :   {
    1074      107331 :     case t_POLMOD:
    1075      107331 :       if (!RgX_equal_var(nf_get_pol(nf),gel(x,1)))
    1076           7 :         pari_err_MODULUS("algtobasis", nf_get_pol(nf),gel(x,1));
    1077      107324 :       x = gel(x,2);
    1078      107324 :       switch(typ(x))
    1079             :       {
    1080        7497 :         case t_INT:
    1081        7497 :         case t_FRAC: return scalarcol(x, nf_get_degree(nf));
    1082       99827 :         case t_POL:
    1083       99827 :           av = avma;
    1084       99827 :           return gerepileupto(av,poltobasis(nf,x));
    1085             :       }
    1086           0 :       break;
    1087             : 
    1088      100324 :     case t_POL:
    1089      100324 :       av = avma;
    1090      100324 :       return gerepileupto(av,poltobasis(nf,x));
    1091             : 
    1092       22316 :     case t_COL:
    1093       22316 :       if (!RgV_is_QV(x)) pari_err_TYPE("nfalgtobasis",x);
    1094       22309 :       if (lg(x)-1 != nf_get_degree(nf)) pari_err_DIM("nfalgtobasis");
    1095       22309 :       return gcopy(x);
    1096             : 
    1097       72653 :     case t_INT:
    1098       72653 :     case t_FRAC: return scalarcol(x, nf_get_degree(nf));
    1099             :   }
    1100           0 :   pari_err_TYPE("algtobasis",x);
    1101             :   return NULL; /* LCOV_EXCL_LINE */
    1102             : }
    1103             : 
    1104             : GEN
    1105       44212 : rnfbasistoalg(GEN rnf,GEN x)
    1106             : {
    1107       44212 :   const char *f = "rnfbasistoalg";
    1108             :   long lx, i;
    1109       44212 :   pari_sp av = avma;
    1110             :   GEN z, nf, R, T;
    1111             : 
    1112       44212 :   checkrnf(rnf);
    1113       44212 :   nf = rnf_get_nf(rnf);
    1114       44212 :   T = nf_get_pol(nf);
    1115       44212 :   R = QXQX_to_mod_shallow(rnf_get_pol(rnf), T);
    1116       44212 :   switch(typ(x))
    1117             :   {
    1118         826 :     case t_COL:
    1119         826 :       z = cgetg_copy(x, &lx);
    1120        2478 :       for (i=1; i<lx; i++)
    1121             :       {
    1122        1701 :         GEN c = nf_to_scalar_or_alg(nf, gel(x,i));
    1123        1652 :         if (typ(c) == t_POL) c = mkpolmod(c,T);
    1124        1652 :         gel(z,i) = c;
    1125             :       }
    1126         777 :       z = RgV_RgC_mul(gel(rnf_get_zk(rnf),1), z);
    1127         714 :       return gerepileupto(av, gmodulo(z,R));
    1128             : 
    1129       29757 :     case t_POLMOD:
    1130       29757 :       x = polmod_nffix(f, rnf, x, 0);
    1131       29554 :       if (typ(x) != t_POL) break;
    1132       13545 :       retmkpolmod(RgX_copy(x), RgX_copy(R));
    1133        1099 :     case t_POL:
    1134        1099 :       if (varn(x) == varn(T)) { RgX_check_QX(x,f); x = gmodulo(x,T); break; }
    1135         875 :       if (varn(x) == varn(R))
    1136             :       {
    1137         826 :         x = RgX_nffix(f,nf_get_pol(nf),x,0);
    1138         826 :         return gmodulo(x, R);
    1139             :       }
    1140          49 :       pari_err_VAR(f, x,R);
    1141             :   }
    1142       28714 :   retmkpolmod(scalarpol(x, varn(R)), RgX_copy(R));
    1143             : }
    1144             : 
    1145             : GEN
    1146        1897 : matbasistoalg(GEN nf,GEN x)
    1147             : {
    1148             :   long i, j, li, lx;
    1149        1897 :   GEN z = cgetg_copy(x, &lx);
    1150             : 
    1151        1897 :   if (lx == 1) return z;
    1152        1890 :   switch(typ(x))
    1153             :   {
    1154          56 :     case t_VEC: case t_COL:
    1155         217 :       for (i=1; i<lx; i++) gel(z,i) = basistoalg(nf, gel(x,i));
    1156          56 :       return z;
    1157        1834 :     case t_MAT: break;
    1158           0 :     default: pari_err_TYPE("matbasistoalg",x);
    1159             :   }
    1160        1834 :   li = lgcols(x);
    1161        6895 :   for (j=1; j<lx; j++)
    1162             :   {
    1163        5061 :     GEN c = cgetg(li,t_COL), xj = gel(x,j);
    1164        5061 :     gel(z,j) = c;
    1165       24451 :     for (i=1; i<li; i++) gel(c,i) = basistoalg(nf, gel(xj,i));
    1166             :   }
    1167        1834 :   return z;
    1168             : }
    1169             : 
    1170             : GEN
    1171        4725 : matalgtobasis(GEN nf,GEN x)
    1172             : {
    1173             :   long i, j, li, lx;
    1174        4725 :   GEN z = cgetg_copy(x, &lx);
    1175             : 
    1176        4725 :   if (lx == 1) return z;
    1177        4361 :   switch(typ(x))
    1178             :   {
    1179        4354 :     case t_VEC: case t_COL:
    1180       13972 :       for (i=1; i<lx; i++) gel(z,i) = algtobasis(nf, gel(x,i));
    1181        4354 :       return z;
    1182           7 :     case t_MAT: break;
    1183           0 :     default: pari_err_TYPE("matalgtobasis",x);
    1184             :   }
    1185           7 :   li = lgcols(x);
    1186          14 :   for (j=1; j<lx; j++)
    1187             :   {
    1188           7 :     GEN c = cgetg(li,t_COL), xj = gel(x,j);
    1189           7 :     gel(z,j) = c;
    1190          21 :     for (i=1; i<li; i++) gel(c,i) = algtobasis(nf, gel(xj,i));
    1191             :   }
    1192           7 :   return z;
    1193             : }
    1194             : GEN
    1195        8715 : RgM_to_nfM(GEN nf,GEN x)
    1196             : {
    1197             :   long i, j, li, lx;
    1198        8715 :   GEN z = cgetg_copy(x, &lx);
    1199             : 
    1200        8715 :   if (lx == 1) return z;
    1201        8715 :   li = lgcols(x);
    1202       67011 :   for (j=1; j<lx; j++)
    1203             :   {
    1204       58296 :     GEN c = cgetg(li,t_COL), xj = gel(x,j);
    1205       58296 :     gel(z,j) = c;
    1206      391720 :     for (i=1; i<li; i++) gel(c,i) = nf_to_scalar_or_basis(nf, gel(xj,i));
    1207             :   }
    1208        8715 :   return z;
    1209             : }
    1210             : GEN
    1211       78428 : RgC_to_nfC(GEN nf, GEN x)
    1212      576660 : { pari_APPLY_type(t_COL, nf_to_scalar_or_basis(nf, gel(x,i))) }
    1213             : 
    1214             : /* x a t_POLMOD, supposedly in rnf = K[z]/(T), K = Q[y]/(Tnf) */
    1215             : GEN
    1216      134491 : polmod_nffix(const char *f, GEN rnf, GEN x, int lift)
    1217      134491 : { return polmod_nffix2(f, rnf_get_nfpol(rnf), rnf_get_pol(rnf), x,lift); }
    1218             : GEN
    1219      134582 : polmod_nffix2(const char *f, GEN T, GEN R, GEN x, int lift)
    1220             : {
    1221      134582 :   if (RgX_equal_var(gel(x,1), R))
    1222             :   {
    1223      124432 :     x = gel(x,2);
    1224      124432 :     if (typ(x) == t_POL && varn(x) == varn(R))
    1225             :     {
    1226       94808 :       x = RgX_nffix(f, T, x, lift);
    1227       94808 :       switch(lg(x))
    1228             :       {
    1229        5782 :         case 2: return gen_0;
    1230       11347 :         case 3: return gel(x,2);
    1231             :       }
    1232       77679 :       return x;
    1233             :     }
    1234             :   }
    1235       39774 :   return Rg_nffix(f, T, x, lift);
    1236             : }
    1237             : GEN
    1238        1204 : rnfalgtobasis(GEN rnf,GEN x)
    1239             : {
    1240        1204 :   const char *f = "rnfalgtobasis";
    1241        1204 :   pari_sp av = avma;
    1242             :   GEN T, R;
    1243             : 
    1244        1204 :   checkrnf(rnf);
    1245        1204 :   R = rnf_get_pol(rnf);
    1246        1204 :   T = rnf_get_nfpol(rnf);
    1247        1204 :   switch(typ(x))
    1248             :   {
    1249          49 :     case t_COL:
    1250          49 :       if (lg(x)-1 != rnf_get_degree(rnf)) pari_err_DIM(f);
    1251          28 :       x = RgV_nffix(f, T, x, 0);
    1252          21 :       return gerepilecopy(av, x);
    1253             : 
    1254        1071 :     case t_POLMOD:
    1255        1071 :       x = polmod_nffix(f, rnf, x, 0);
    1256        1036 :       if (typ(x) != t_POL) break;
    1257         714 :       return gerepileupto(av, RgM_RgX_mul(rnf_get_invzk(rnf), x));
    1258          56 :     case t_POL:
    1259          56 :       if (varn(x) == varn(T))
    1260             :       {
    1261          21 :         RgX_check_QX(x,f);
    1262          14 :         if (degpol(x) >= degpol(T)) x = RgX_rem(x,T);
    1263          14 :         x = mkpolmod(x,T); break;
    1264             :       }
    1265          35 :       x = RgX_nffix(f, T, x, 0);
    1266          28 :       if (degpol(x) >= degpol(R)) x = RgX_rem(x, R);
    1267          28 :       return gerepileupto(av, RgM_RgX_mul(rnf_get_invzk(rnf), x));
    1268             :   }
    1269         364 :   return gerepileupto(av, scalarcol(x, rnf_get_degree(rnf)));
    1270             : }
    1271             : 
    1272             : /* Given a and b in nf, gives an algebraic integer y in nf such that a-b.y
    1273             :  * is "small" */
    1274             : GEN
    1275         259 : nfdiveuc(GEN nf, GEN a, GEN b)
    1276             : {
    1277         259 :   pari_sp av = avma;
    1278         259 :   a = nfdiv(nf,a,b);
    1279         259 :   return gerepileupto(av, ground(a));
    1280             : }
    1281             : 
    1282             : /* Given a and b in nf, gives a "small" algebraic integer r in nf
    1283             :  * of the form a-b.y */
    1284             : GEN
    1285         259 : nfmod(GEN nf, GEN a, GEN b)
    1286             : {
    1287         259 :   pari_sp av = avma;
    1288         259 :   GEN p1 = gneg_i(nfmul(nf,b,ground(nfdiv(nf,a,b))));
    1289         259 :   return gerepileupto(av, nfadd(nf,a,p1));
    1290             : }
    1291             : 
    1292             : /* Given a and b in nf, gives a two-component vector [y,r] in nf such
    1293             :  * that r=a-b.y is "small". */
    1294             : GEN
    1295         259 : nfdivrem(GEN nf, GEN a, GEN b)
    1296             : {
    1297         259 :   pari_sp av = avma;
    1298         259 :   GEN p1,z, y = ground(nfdiv(nf,a,b));
    1299             : 
    1300         259 :   p1 = gneg_i(nfmul(nf,b,y));
    1301         259 :   z = cgetg(3,t_VEC);
    1302         259 :   gel(z,1) = gcopy(y);
    1303         259 :   gel(z,2) = nfadd(nf,a,p1); return gerepileupto(av, z);
    1304             : }
    1305             : 
    1306             : /*************************************************************************/
    1307             : /**                                                                     **/
    1308             : /**                   LOGARITHMIC EMBEDDINGS                            **/
    1309             : /**                                                                     **/
    1310             : /*************************************************************************/
    1311             : 
    1312             : static int
    1313      692558 : low_prec(GEN x)
    1314             : {
    1315      692558 :   switch(typ(x))
    1316             :   {
    1317           0 :     case t_INT: return !signe(x);
    1318      692558 :     case t_REAL: return !signe(x) || realprec(x) <= DEFAULTPREC;
    1319           0 :     default: return 0;
    1320             :   }
    1321             : }
    1322             : 
    1323             : static GEN
    1324       11459 : cxlog_1(GEN nf) { return zerocol(lg(nf_get_roots(nf))-1); }
    1325             : static GEN
    1326           0 : cxlog_m1(GEN nf, long prec)
    1327             : {
    1328           0 :   long i, l = lg(nf_get_roots(nf)), r1 = nf_get_r1(nf);
    1329           0 :   GEN v = cgetg(l, t_COL), p,  P;
    1330           0 :   p = mppi(prec); P = mkcomplex(gen_0, p);
    1331           0 :   for (i = 1; i <= r1; i++) gel(v,i) = P; /* IPi*/
    1332           0 :   if (i < l) P = gmul2n(P,1);
    1333           0 :   for (     ; i < l; i++) gel(v,i) = P; /* 2IPi */
    1334           0 :   return v;
    1335             : }
    1336             : static GEN
    1337        9940 : famat_cxlog(GEN nf, GEN fa, long prec)
    1338             : {
    1339        9940 :   GEN g, e, y = NULL;
    1340             :   long i, l;
    1341             : 
    1342        9940 :   if (typ(fa) != t_MAT) pari_err_TYPE("famat_cxlog",fa);
    1343        9940 :   if (lg(fa) == 1) return cxlog_1(nf);
    1344        9940 :   g = gel(fa,1);
    1345        9940 :   e = gel(fa,2); l = lg(e);
    1346       19742 :   for (i = 1; i < l; i++)
    1347             :   {
    1348        9802 :     GEN t, x = nf_to_scalar_or_basis(nf, gel(g,i));
    1349             :     /* multiplicative arch would be better (save logs), but exponents overflow
    1350             :      * [ could keep track of expo separately, but not worth it ] */
    1351        9802 :     t = nf_cxlog(nf,x,prec); if (!t) return NULL;
    1352        9802 :     if (gel(t,1) == gen_0) continue; /* positive rational */
    1353        6193 :     t = RgC_Rg_mul(t, gel(e,i));
    1354        6193 :     y = y? RgV_add(y,t): t;
    1355             :   }
    1356        9940 :   return y ? y: cxlog_1(nf);
    1357             : }
    1358             : /* Archimedean components: [e_i Log( sigma_i(X) )], where X = primpart(x),
    1359             :  * and e_i = 1 (resp 2.) for i <= R1 (resp. > R1) */
    1360             : GEN
    1361      234795 : nf_cxlog(GEN nf, GEN x, long prec)
    1362             : {
    1363             :   long i, l, r1;
    1364             :   GEN v;
    1365      234795 :   if (typ(x) == t_MAT) return famat_cxlog(nf,x,prec);
    1366      224855 :   x = nf_to_scalar_or_basis(nf,x);
    1367      224855 :   if (typ(x) != t_COL) return gsigne(x) > 0? cxlog_1(nf): cxlog_m1(nf, prec);
    1368      220602 :   x = RgM_RgC_mul(nf_get_M(nf), Q_primpart(x));
    1369      220602 :   l = lg(x); r1 = nf_get_r1(nf);
    1370      651653 :   for (i = 1; i <= r1; i++)
    1371      431051 :     if (low_prec(gel(x,i))) return NULL;
    1372      321992 :   for (     ; i <  l;  i++)
    1373      101390 :     if (low_prec(gnorm(gel(x,i)))) return NULL;
    1374      220602 :   v = cgetg(l,t_COL);
    1375      651653 :   for (i = 1; i <= r1; i++) gel(v,i) = glog(gel(x,i),prec);
    1376      321992 :   for (     ; i <  l;  i++) gel(v,i) = gmul2n(glog(gel(x,i),prec),1);
    1377      220602 :   return v;
    1378             : }
    1379             : GEN
    1380          91 : nfV_cxlog(GEN nf, GEN x, long prec)
    1381             : {
    1382             :   long i, l;
    1383          91 :   GEN v = cgetg_copy(x, &l);
    1384         161 :   for (i = 1; i < l; i++)
    1385          70 :     if (!(gel(v,i) = nf_cxlog(nf, gel(x,i), prec))) return NULL;
    1386          91 :   return v;
    1387             : }
    1388             : 
    1389             : static GEN
    1390        7420 : scalar_logembed(GEN nf, GEN u, GEN *emb)
    1391             : {
    1392             :   GEN v, logu;
    1393        7420 :   long i, s = signe(u), RU = lg(nf_get_roots(nf))-1, R1 = nf_get_r1(nf);
    1394             : 
    1395        7420 :   if (!s) pari_err_DOMAIN("nflogembed","argument","=",gen_0,u);
    1396        7420 :   v = cgetg(RU+1, t_COL); logu = logr_abs(u);
    1397        9191 :   for (i = 1; i <= R1; i++) gel(v,i) = logu;
    1398        7420 :   if (i <= RU)
    1399             :   {
    1400        6545 :     GEN logu2 = shiftr(logu,1);
    1401       24969 :     for (   ; i <= RU; i++) gel(v,i) = logu2;
    1402             :   }
    1403        7420 :   if (emb) *emb = const_col(RU, u);
    1404        7420 :   return v;
    1405             : }
    1406             : 
    1407             : static GEN
    1408        1302 : famat_logembed(GEN nf,GEN x,GEN *emb,long prec)
    1409             : {
    1410        1302 :   GEN A, M, T, a, t, g = gel(x,1), e = gel(x,2);
    1411        1302 :   long i, l = lg(e);
    1412             : 
    1413        1302 :   if (l == 1) return scalar_logembed(nf, real_1(prec), emb);
    1414        1302 :   A = NULL; T = emb? cgetg(l, t_COL): NULL;
    1415        1302 :   if (emb) *emb = M = mkmat2(T, e);
    1416       60886 :   for (i = 1; i < l; i++)
    1417             :   {
    1418       59584 :     a = nflogembed(nf, gel(g,i), &t, prec);
    1419       59584 :     if (!a) return NULL;
    1420       59584 :     a = RgC_Rg_mul(a, gel(e,i));
    1421       59584 :     A = A? RgC_add(A, a): a;
    1422       59584 :     if (emb) gel(T,i) = t;
    1423             :   }
    1424        1302 :   return A;
    1425             : }
    1426             : 
    1427             : /* Get archimedean components: [e_i log( | sigma_i(x) | )], with e_i = 1
    1428             :  * (resp 2.) for i <= R1 (resp. > R1) and set emb to the embeddings of x.
    1429             :  * Return NULL if precision problem */
    1430             : GEN
    1431       62510 : nflogembed(GEN nf, GEN x, GEN *emb, long prec)
    1432             : {
    1433             :   long i, l, r1;
    1434             :   GEN v, t;
    1435             : 
    1436       62510 :   if (typ(x) == t_MAT) return famat_logembed(nf,x,emb,prec);
    1437       61208 :   x = nf_to_scalar_or_basis(nf,x);
    1438       61208 :   if (typ(x) != t_COL) return scalar_logembed(nf, gtofp(x,prec), emb);
    1439       53788 :   x = RgM_RgC_mul(nf_get_M(nf), x);
    1440       53788 :   l = lg(x); r1 = nf_get_r1(nf); v = cgetg(l,t_COL);
    1441       58512 :   for (i = 1; i <= r1; i++)
    1442             :   {
    1443        4731 :     t = gabs(gel(x,i),prec); if (low_prec(t)) return NULL;
    1444        4724 :     gel(v,i) = glog(t,prec);
    1445             :   }
    1446      209167 :   for (   ; i < l; i++)
    1447             :   {
    1448      155386 :     t = gnorm(gel(x,i)); if (low_prec(t)) return NULL;
    1449      155386 :     gel(v,i) = glog(t,prec);
    1450             :   }
    1451       53781 :   if (emb) *emb = x;
    1452       53781 :   return v;
    1453             : }
    1454             : 
    1455             : /*************************************************************************/
    1456             : /**                                                                     **/
    1457             : /**                        REAL EMBEDDINGS                              **/
    1458             : /**                                                                     **/
    1459             : /*************************************************************************/
    1460             : static GEN
    1461       83860 : sarch_get_cyc(GEN sarch) { return gel(sarch,1); }
    1462             : static GEN
    1463      339882 : sarch_get_archp(GEN sarch) { return gel(sarch,2); }
    1464             : static GEN
    1465       67146 : sarch_get_MI(GEN sarch) { return gel(sarch,3); }
    1466             : static GEN
    1467       67146 : sarch_get_lambda(GEN sarch) { return gel(sarch,4); }
    1468             : static GEN
    1469       67146 : sarch_get_F(GEN sarch) { return gel(sarch,5); }
    1470             : 
    1471             : /* x not a scalar, true nf, return number of positive roots of char_x */
    1472             : static long
    1473        2168 : num_positive(GEN nf, GEN x)
    1474             : {
    1475        2168 :   GEN T = nf_get_pol(nf), B, charx;
    1476             :   long dnf, vnf, N;
    1477        2168 :   x = nf_to_scalar_or_alg(nf, x); /* not a scalar */
    1478        2168 :   charx = ZXQ_charpoly(x, T, 0);
    1479        2168 :   charx = ZX_radical(charx);
    1480        2168 :   N = degpol(T) / degpol(charx);
    1481             :   /* real places are unramified ? */
    1482        2168 :   if (N == 1 || ZX_sturm(charx) * N == nf_get_r1(nf))
    1483        2161 :     return ZX_sturmpart(charx, mkvec2(gen_0,mkoo())) * N;
    1484             :   /* painful case, multiply by random square until primitive */
    1485           7 :   dnf = nf_get_degree(nf);
    1486           7 :   vnf = varn(T);
    1487           7 :   B = int2n(10);
    1488             :   for(;;)
    1489           0 :   {
    1490           7 :     GEN y = RgXQ_sqr(random_FpX(dnf, vnf, B), T);
    1491           7 :     y = RgXQ_mul(x, y, T);
    1492           7 :     charx = ZXQ_charpoly(y, T, 0);
    1493           7 :     if (ZX_is_squarefree(charx))
    1494           7 :       return ZX_sturmpart(charx, mkvec2(gen_0,mkoo())) * N;
    1495             :   }
    1496             : }
    1497             : 
    1498             : /* x a QC: return sigma_k(x) where 1 <= k <= r1+r2; correct but inefficient
    1499             :  * if x in Q. M = nf_get_M(nf) */
    1500             : static GEN
    1501          91 : nfembed_i(GEN M, GEN x, long k)
    1502             : {
    1503          91 :   long i, l = lg(M);
    1504          91 :   GEN z = gel(x,1);
    1505         364 :   for (i = 2; i < l; i++) z = gadd(z, gmul(gcoeff(M,k,i), gel(x,i)));
    1506          91 :   return z;
    1507             : }
    1508             : GEN
    1509           0 : nfembed(GEN nf, GEN x, long k)
    1510             : {
    1511           0 :   pari_sp av = avma;
    1512           0 :   nf = checknf(nf);
    1513           0 :   x = nf_to_scalar_or_basis(nf,x);
    1514           0 :   if (typ(x) != t_COL) return gerepilecopy(av, x);
    1515           0 :   return gerepileupto(av, nfembed_i(nf_get_M(nf),x,k));
    1516             : }
    1517             : 
    1518             : /* x a ZC */
    1519             : static GEN
    1520      435058 : zk_embed(GEN M, GEN x, long k)
    1521             : {
    1522      435058 :   long i, l = lg(x);
    1523      435058 :   GEN z = gel(x,1); /* times M[k,1], which is 1 */
    1524     1554809 :   for (i = 2; i < l; i++) z = mpadd(z, mpmul(gcoeff(M,k,i), gel(x,i)));
    1525      435058 :   return z;
    1526             : }
    1527             : 
    1528             : /* Given floating point approximation z of sigma_k(x), decide its sign
    1529             :  * [0/+, 1/- and -1 for FAIL] */
    1530             : static long
    1531      423076 : eval_sign_embed(GEN z)
    1532             : { /* dubious, fail */
    1533      423076 :   if (typ(z) == t_REAL && realprec(z) <= LOWDEFAULTPREC) return -1;
    1534      421638 :   return (signe(z) < 1)? 1: 0;
    1535             : }
    1536             : /* return v such that (-1)^v = sign(sigma_k(x)), x primitive ZC */
    1537             : static long
    1538      345262 : eval_sign(GEN M, GEN x, long k)
    1539      345262 : { return eval_sign_embed( zk_embed(M, x, k) ); }
    1540             : 
    1541             : /* check that signs[i..#signs] == s; signs = NULL encodes "totally positive" */
    1542             : static int
    1543           0 : oksigns(long l, GEN signs, long i, long s)
    1544             : {
    1545           0 :   if (!signs) return s == 0;
    1546           0 :   for (; i < l; i++)
    1547           0 :     if (signs[i] != s) return 0;
    1548           0 :   return 1;
    1549             : }
    1550             : /* check that signs[i] = s and signs[i+1..#signs] = 1-s */
    1551             : static int
    1552           0 : oksigns2(long l, GEN signs, long i, long s)
    1553             : {
    1554           0 :   if (!signs) return s == 0 && i == l-1;
    1555           0 :   return signs[i] == s && oksigns(l, signs, i+1, 1-s);
    1556             : }
    1557             : 
    1558             : /* true nf, x a ZC (primitive for efficiency), embx its embeddings or NULL */
    1559             : static int
    1560       67872 : nfchecksigns_i(GEN nf, GEN x, GEN embx, GEN signs, GEN archp)
    1561             : {
    1562       67872 :   long l = lg(archp), i;
    1563       67872 :   GEN M = nf_get_M(nf), sarch = NULL;
    1564       67872 :   long np = -1;
    1565      102720 :   for (i = 1; i < l; i++)
    1566             :   {
    1567             :     long s;
    1568       78696 :     if (embx)
    1569       77814 :       s = eval_sign_embed(gel(embx,i));
    1570             :     else
    1571         882 :       s = eval_sign(M, x, archp[i]);
    1572             :     /* 0 / + or 1 / -; -1 for FAIL */
    1573       78696 :     if (s < 0) /* failure */
    1574             :     {
    1575           0 :       long ni, r1 = nf_get_r1(nf);
    1576             :       GEN xi;
    1577           0 :       if (np < 0)
    1578             :       {
    1579           0 :         np = num_positive(nf, x);
    1580           0 :         if (np == 0)  return oksigns(l, signs, i, 1);
    1581           0 :         if (np == r1) return oksigns(l, signs, i, 0);
    1582           0 :         sarch = nfarchstar(nf, NULL, identity_perm(r1));
    1583             :       }
    1584           0 :       xi = set_sign_mod_divisor(nf, vecsmall_ei(r1, archp[i]), gen_1, sarch);
    1585           0 :       xi = Q_primpart(xi);
    1586           0 :       ni = num_positive(nf, nfmuli(nf,x,xi));
    1587           0 :       if (ni == 0)  return oksigns2(l, signs, i, 0);
    1588           0 :       if (ni == r1) return oksigns2(l, signs, i, 1);
    1589           0 :       s = ni < np? 0: 1;
    1590             :     }
    1591       78696 :     if (s != (signs? signs[i]: 0)) return 0;
    1592             :   }
    1593       24024 :   return 1;
    1594             : }
    1595             : static void
    1596         490 : pl_convert(GEN pl, GEN *psigns, GEN *parchp)
    1597             : {
    1598         490 :   long i, j, l = lg(pl);
    1599         490 :   GEN signs = cgetg(l, t_VECSMALL);
    1600         490 :   GEN archp = cgetg(l, t_VECSMALL);
    1601        2128 :   for (i = j = 1; i < l; i++)
    1602             :   {
    1603        1638 :     if (!pl[i]) continue;
    1604        1050 :     archp[j] = i;
    1605        1050 :     signs[j] = (pl[i] < 0)? 1: 0;
    1606        1050 :     j++;
    1607             :   }
    1608         490 :   setlg(archp, j); *parchp = archp;
    1609         490 :   setlg(signs, j); *psigns = signs;
    1610         490 : }
    1611             : /* pl : requested signs for real embeddings, 0 = no sign constraint */
    1612             : int
    1613        1015 : nfchecksigns(GEN nf, GEN x, GEN pl)
    1614             : {
    1615        1015 :   pari_sp av = avma;
    1616             :   GEN signs, archp;
    1617        1015 :   nf = checknf(nf);
    1618        1015 :   x = nf_to_scalar_or_basis(nf,x);
    1619        1015 :   if (typ(x) != t_COL)
    1620             :   {
    1621         525 :     long i, l = lg(pl), s = gsigne(x);
    1622        1064 :     for (i = 1; i < l; i++)
    1623         539 :       if (pl[i] && pl[i] != s) return gc_bool(av,0);
    1624         525 :     return gc_bool(av,1);
    1625             :   }
    1626         490 :   pl_convert(pl, &signs, &archp);
    1627         490 :   return gc_bool(av, nfchecksigns_i(nf, x, NULL, signs, archp));
    1628             : }
    1629             : 
    1630             : /* signs = NULL: totally positive, else sign[i] = 0 (+) or 1 (-) */
    1631             : static GEN
    1632       67146 : get_C(GEN lambda, long l, GEN signs)
    1633             : {
    1634             :   long i;
    1635             :   GEN C, mlambda;
    1636       67146 :   if (!signs) return const_vec(l-1, lambda);
    1637       28450 :   C = cgetg(l, t_COL); mlambda = gneg(lambda);
    1638       67123 :   for (i = 1; i < l; i++) gel(C,i) = signs[i]? mlambda: lambda;
    1639       28450 :   return C;
    1640             : }
    1641             : /* signs = NULL: totally positive at archp */
    1642             : static GEN
    1643      120668 : nfsetsigns(GEN nf, GEN signs, GEN x, GEN sarch)
    1644             : {
    1645      120668 :   long i, l = lg(sarch_get_archp(sarch));
    1646             :   GEN ex;
    1647             :   /* Is signature already correct ? */
    1648      120668 :   if (typ(x) != t_COL)
    1649             :   {
    1650       53286 :     long s = gsigne(x);
    1651       53286 :     if (!s) i = 1;
    1652       53272 :     else if (!signs)
    1653        4942 :       i = (s < 0)? 1: l;
    1654             :     else
    1655             :     {
    1656       48330 :       s = s < 0? 1: 0;
    1657       78385 :       for (i = 1; i < l; i++)
    1658       53325 :         if (signs[i] != s) break;
    1659             :     }
    1660       53286 :     ex = (i < l)? const_col(l-1, x): NULL;
    1661             :   }
    1662             :   else
    1663             :   {
    1664       67382 :     pari_sp av = avma;
    1665       67382 :     GEN cex, M = nf_get_M(nf), archp = sarch_get_archp(sarch);
    1666       67382 :     GEN xp = Q_primitive_part(x,&cex);
    1667       67382 :     ex = cgetg(l,t_COL);
    1668      157178 :     for (i = 1; i < l; i++) gel(ex,i) = zk_embed(M,xp,archp[i]);
    1669       67382 :     if (nfchecksigns_i(nf, xp, ex, signs, archp)) { ex = NULL; set_avma(av); }
    1670       43729 :     else if (cex) ex = RgC_Rg_mul(ex, cex); /* put back content */
    1671             :   }
    1672      120668 :   if (ex)
    1673             :   { /* If no, fix it */
    1674       67146 :     GEN MI = sarch_get_MI(sarch), F = sarch_get_F(sarch);
    1675       67146 :     GEN lambda = sarch_get_lambda(sarch);
    1676       67146 :     GEN t = RgC_sub(get_C(lambda, l, signs), ex);
    1677             :     long e;
    1678       67146 :     t = grndtoi(RgM_RgC_mul(MI,t), &e);
    1679       67146 :     if (lg(F) != 1) t = ZM_ZC_mul(F, t);
    1680       67146 :     x = typ(x) == t_COL? RgC_add(t, x): RgC_Rg_add(t, x);
    1681             :   }
    1682      120668 :   return x;
    1683             : }
    1684             : /* - sarch = nfarchstar(nf, F);
    1685             :  * - x encodes a vector of signs at arch.archp: either a t_VECSMALL
    1686             :  *   (vector of signs as {0,1}-vector), NULL (totally positive at archp),
    1687             :  *   or a non-zero number field element (replaced by its signature at archp);
    1688             :  * - y is a non-zero number field element
    1689             :  * Return z = y (mod F) with signs(y, archp) = signs(x) (a {0,1}-vector) */
    1690             : GEN
    1691      151832 : set_sign_mod_divisor(GEN nf, GEN x, GEN y, GEN sarch)
    1692             : {
    1693      151832 :   GEN archp = sarch_get_archp(sarch);
    1694      151832 :   if (lg(archp) == 1) return y;
    1695      119296 :   nf = checknf(nf);
    1696      119296 :   if (x && typ(x) != t_VECSMALL) x = nfsign_arch(nf, x, archp);
    1697      119296 :   y = nf_to_scalar_or_basis(nf,y);
    1698      119296 :   return nfsetsigns(nf, x, y, sarch);
    1699             : }
    1700             : 
    1701             : static GEN
    1702       23641 : setsigns_init(GEN nf, GEN archp, GEN F, GEN DATA)
    1703             : {
    1704       23641 :   GEN lambda, Mr = rowpermute(nf_get_M(nf), archp), MI = F? RgM_mul(Mr,F): Mr;
    1705       23641 :   lambda = gmul2n(matrixnorm(MI,DEFAULTPREC), -1);
    1706       23641 :   if (typ(lambda) != t_REAL) lambda = gmul(lambda, sstoQ(1001,1000));
    1707       23641 :   if (lg(archp) < lg(MI))
    1708             :   {
    1709       21175 :     GEN perm = gel(indexrank(MI), 2);
    1710       21175 :     if (!F) F = matid(nf_get_degree(nf));
    1711       21175 :     MI = vecpermute(MI, perm);
    1712       21175 :     F = vecpermute(F, perm);
    1713             :   }
    1714       23641 :   if (!F) F = cgetg(1,t_MAT);
    1715       23641 :   MI = RgM_inv(MI);
    1716       23641 :   return mkvec5(DATA, archp, MI, lambda, F);
    1717             : }
    1718             : /* F non-0 integral ideal in HNF (or NULL: Z_K), compute elements in 1+F
    1719             :  * whose sign matrix at archp is identity; archp in 'indices' format */
    1720             : GEN
    1721       47098 : nfarchstar(GEN nf, GEN F, GEN archp)
    1722             : {
    1723       47098 :   long nba = lg(archp) - 1;
    1724       47098 :   if (!nba) return mkvec2(cgetg(1,t_VEC), archp);
    1725       22276 :   if (F && equali1(gcoeff(F,1,1))) F = NULL;
    1726       22276 :   if (F) F = idealpseudored(F, nf_get_roundG(nf));
    1727       22276 :   return setsigns_init(nf, archp, F, const_vec(nba, gen_2));
    1728             : }
    1729             : 
    1730             : /*************************************************************************/
    1731             : /**                                                                     **/
    1732             : /**                         IDEALCHINESE                                **/
    1733             : /**                                                                     **/
    1734             : /*************************************************************************/
    1735             : static int
    1736        2996 : isprfact(GEN x)
    1737             : {
    1738             :   long i, l;
    1739             :   GEN L, E;
    1740        2996 :   if (typ(x) != t_MAT || lg(x) != 3) return 0;
    1741        2996 :   L = gel(x,1); l = lg(L);
    1742        2996 :   E = gel(x,2);
    1743        7217 :   for(i=1; i<l; i++)
    1744             :   {
    1745        4221 :     checkprid(gel(L,i));
    1746        4221 :     if (typ(gel(E,i)) != t_INT) return 0;
    1747             :   }
    1748        2996 :   return 1;
    1749             : }
    1750             : 
    1751             : /* initialize projectors mod pr[i]^e[i] for idealchinese */
    1752             : static GEN
    1753        2996 : pr_init(GEN nf, GEN fa, GEN w, GEN dw)
    1754             : {
    1755        2996 :   GEN U, E, F, L = gel(fa,1), E0 = gel(fa,2);
    1756        2996 :   long i, r = lg(L);
    1757             : 
    1758        2996 :   if (w && lg(w) != r) pari_err_TYPE("idealchinese", w);
    1759        2996 :   if (r == 1 && !dw) return cgetg(1,t_VEC);
    1760        2989 :   E = leafcopy(E0); /* do not destroy fa[2] */
    1761        7210 :   for (i = 1; i < r; i++)
    1762        4221 :     if (signe(gel(E,i)) < 0) gel(E,i) = gen_0;
    1763        2989 :   F = factorbackprime(nf, L, E);
    1764        2989 :   if (dw)
    1765             :   {
    1766         700 :     F = ZM_Z_mul(F, dw);
    1767        1603 :     for (i = 1; i < r; i++)
    1768             :     {
    1769         903 :       GEN pr = gel(L,i);
    1770         903 :       long e = itos(gel(E0,i)), v = idealval(nf, dw, pr);
    1771         903 :       if (e >= 0)
    1772         896 :         gel(E,i) = addiu(gel(E,i), v);
    1773           7 :       else if (v + e <= 0)
    1774           0 :         F = idealmulpowprime(nf, F, pr, stoi(-v)); /* coprime to pr */
    1775             :       else
    1776             :       {
    1777           7 :         F = idealmulpowprime(nf, F, pr, stoi(e));
    1778           7 :         gel(E,i) = stoi(v + e);
    1779             :       }
    1780             :     }
    1781             :   }
    1782        2989 :   U = cgetg(r, t_VEC);
    1783        7210 :   for (i = 1; i < r; i++)
    1784             :   {
    1785             :     GEN u;
    1786        4221 :     if (w && gequal0(gel(w,i))) u = gen_0; /* unused */
    1787             :     else
    1788             :     {
    1789        4144 :       GEN pr = gel(L,i), e = gel(E,i), t;
    1790        4144 :       t = idealdivpowprime(nf,F, pr, e);
    1791        4144 :       u = hnfmerge_get_1(t, idealpow(nf, pr, e));
    1792        4144 :       if (!u) pari_err_COPRIME("idealchinese", t,pr);
    1793             :     }
    1794        4221 :     gel(U,i) = u;
    1795             :   }
    1796        2989 :   F = idealpseudored(F, nf_get_roundG(nf));
    1797        2989 :   return mkvec2(F, U);
    1798             : }
    1799             : 
    1800             : static GEN
    1801        1785 : pl_normalize(GEN nf, GEN pl)
    1802             : {
    1803        1785 :   const char *fun = "idealchinese";
    1804        1785 :   if (lg(pl)-1 != nf_get_r1(nf)) pari_err_TYPE(fun,pl);
    1805        1785 :   switch(typ(pl))
    1806             :   {
    1807         714 :     case t_VEC: RgV_check_ZV(pl,fun); pl = ZV_to_zv(pl);
    1808             :       /* fall through */
    1809        1785 :     case t_VECSMALL: break;
    1810           0 :     default: pari_err_TYPE(fun,pl);
    1811             :   }
    1812        1785 :   return pl;
    1813             : }
    1814             : 
    1815             : static int
    1816        7147 : is_chineseinit(GEN x)
    1817             : {
    1818             :   GEN fa, pl;
    1819             :   long l;
    1820        7147 :   if (typ(x) != t_VEC || lg(x)!=3) return 0;
    1821        5467 :   fa = gel(x,1);
    1822        5467 :   pl = gel(x,2);
    1823        5467 :   if (typ(fa) != t_VEC || typ(pl) != t_VEC) return 0;
    1824        2695 :   l = lg(fa);
    1825        2695 :   if (l != 1)
    1826             :   {
    1827        2674 :     if (l != 3 || typ(gel(fa,1)) != t_MAT || typ(gel(fa,2)) != t_VEC)
    1828           7 :       return 0;
    1829             :   }
    1830        2688 :   l = lg(pl);
    1831        2688 :   if (l != 1)
    1832             :   {
    1833         511 :     if (l != 6 || typ(gel(pl,3)) != t_MAT || typ(gel(pl,1)) != t_VECSMALL
    1834         511 :                                           || typ(gel(pl,2)) != t_VECSMALL)
    1835           0 :       return 0;
    1836             :   }
    1837        2688 :   return 1;
    1838             : }
    1839             : 
    1840             : /* nf a true 'nf' */
    1841             : static GEN
    1842        3129 : chineseinit_i(GEN nf, GEN fa, GEN w, GEN dw)
    1843             : {
    1844        3129 :   const char *fun = "idealchineseinit";
    1845        3129 :   GEN archp = NULL, pl = NULL;
    1846        3129 :   switch(typ(fa))
    1847             :   {
    1848        1785 :     case t_VEC:
    1849        1785 :       if (is_chineseinit(fa))
    1850             :       {
    1851           0 :         if (dw) pari_err_DOMAIN(fun, "denom(y)", "!=", gen_1, w);
    1852           0 :         return fa;
    1853             :       }
    1854        1785 :       if (lg(fa) != 3) pari_err_TYPE(fun, fa);
    1855             :       /* of the form [x,s] */
    1856        1785 :       pl = pl_normalize(nf, gel(fa,2));
    1857        1785 :       fa = gel(fa,1);
    1858        1785 :       archp = vecsmall01_to_indices(pl);
    1859             :       /* keep pr_init, reset pl */
    1860        1785 :       if (is_chineseinit(fa)) { fa = gel(fa,1); break; }
    1861             :       /* fall through */
    1862             :     case t_MAT: /* factorization? */
    1863        2996 :       if (isprfact(fa)) { fa = pr_init(nf, fa, w, dw); break; }
    1864           0 :     default: pari_err_TYPE(fun,fa);
    1865             :   }
    1866             : 
    1867        3129 :   if (!pl) pl = cgetg(1,t_VEC);
    1868             :   else
    1869             :   {
    1870        1785 :     long r = lg(archp);
    1871        1785 :     if (r == 1) pl = cgetg(1, t_VEC);
    1872             :     else
    1873             :     {
    1874        1365 :       GEN F = (lg(fa) == 1)? NULL: gel(fa,1), signs = cgetg(r, t_VECSMALL);
    1875             :       long i;
    1876        3381 :       for (i = 1; i < r; i++) signs[i] = (pl[archp[i]] < 0)? 1: 0;
    1877        1365 :       pl = setsigns_init(nf, archp, F, signs);
    1878             :     }
    1879             :   }
    1880        3129 :   return mkvec2(fa, pl);
    1881             : }
    1882             : 
    1883             : /* Given a prime ideal factorization x, possibly with 0 or negative exponents,
    1884             :  * and a vector w of elements of nf, gives b such that
    1885             :  * v_p(b-w_p)>=v_p(x) for all prime ideals p in the ideal factorization
    1886             :  * and v_p(b)>=0 for all other p, using the standard proof given in GTM 138. */
    1887             : GEN
    1888        5684 : idealchinese(GEN nf, GEN x, GEN w)
    1889             : {
    1890        5684 :   const char *fun = "idealchinese";
    1891        5684 :   pari_sp av = avma;
    1892             :   GEN x1, x2, s, dw, F;
    1893             : 
    1894        5684 :   nf = checknf(nf);
    1895        5684 :   if (!w) return gerepilecopy(av, chineseinit_i(nf,x,NULL,NULL));
    1896             : 
    1897        3577 :   if (typ(w) != t_VEC) pari_err_TYPE(fun,w);
    1898        3577 :   w = Q_remove_denom(matalgtobasis(nf,w), &dw);
    1899        3577 :   if (!is_chineseinit(x)) x = chineseinit_i(nf,x,w,dw);
    1900             :   /* x is a 'chineseinit' */
    1901        3577 :   x1 = gel(x,1); s = NULL;
    1902        3577 :   x2 = gel(x,2);
    1903        3577 :   if (lg(x1) == 1) F = NULL;
    1904             :   else
    1905             :   {
    1906        3556 :     GEN  U = gel(x1,2);
    1907        3556 :     long i, r = lg(w);
    1908        3556 :     F = gel(x1,1);
    1909       10479 :     for (i=1; i<r; i++)
    1910        6923 :       if (!gequal0(gel(w,i)))
    1911             :       {
    1912        4165 :         GEN t = nfmuli(nf, gel(U,i), gel(w,i));
    1913        4165 :         s = s? ZC_add(s,t): t;
    1914             :       }
    1915        3556 :     if (s) s = ZC_reducemodmatrix(s, F);
    1916             :   }
    1917        3577 :   if (lg(x2) != 1) s = nfsetsigns(nf, gel(x2,1), s? s: gen_0, x2);
    1918        3577 :   if (!s) { s = zerocol(nf_get_degree(nf)); dw = NULL; }
    1919             : 
    1920        3577 :   if (dw) s = RgC_Rg_div(s,dw);
    1921        3577 :   return gerepileupto(av, s);
    1922             : }
    1923             : 
    1924             : /*************************************************************************/
    1925             : /**                                                                     **/
    1926             : /**                           (Z_K/I)^*                                 **/
    1927             : /**                                                                     **/
    1928             : /*************************************************************************/
    1929             : GEN
    1930        1785 : vecsmall01_to_indices(GEN v)
    1931             : {
    1932        1785 :   long i, k, l = lg(v);
    1933        1785 :   GEN p = new_chunk(l) + l;
    1934        4788 :   for (k=1, i=l-1; i; i--)
    1935        3003 :     if (v[i]) { *--p = i; k++; }
    1936        1785 :   *--p = evallg(k) | evaltyp(t_VECSMALL);
    1937        1785 :   set_avma((pari_sp)p); return p;
    1938             : }
    1939             : GEN
    1940      423500 : vec01_to_indices(GEN v)
    1941             : {
    1942             :   long i, k, l;
    1943             :   GEN p;
    1944             : 
    1945      423500 :   switch (typ(v))
    1946             :   {
    1947      399728 :    case t_VECSMALL: return v;
    1948       23772 :    case t_VEC: break;
    1949           0 :    default: pari_err_TYPE("vec01_to_indices",v);
    1950             :   }
    1951       23772 :   l = lg(v);
    1952       23772 :   p = new_chunk(l) + l;
    1953       70238 :   for (k=1, i=l-1; i; i--)
    1954       46466 :     if (signe(gel(v,i))) { *--p = i; k++; }
    1955       23772 :   *--p = evallg(k) | evaltyp(t_VECSMALL);
    1956       23772 :   set_avma((pari_sp)p); return p;
    1957             : }
    1958             : GEN
    1959        2996 : indices_to_vec01(GEN p, long r)
    1960             : {
    1961        2996 :   long i, l = lg(p);
    1962        2996 :   GEN v = zerovec(r);
    1963        4011 :   for (i = 1; i < l; i++) gel(v, p[i]) = gen_1;
    1964        2996 :   return v;
    1965             : }
    1966             : 
    1967             : /* return (column) vector of R1 signatures of x (0 or 1) */
    1968             : GEN
    1969      399728 : nfsign_arch(GEN nf, GEN x, GEN arch)
    1970             : {
    1971      399728 :   GEN sarch, M, V, archp = vec01_to_indices(arch);
    1972      399728 :   long i, s, np, n = lg(archp)-1;
    1973             :   pari_sp av;
    1974             : 
    1975      399728 :   if (!n) return cgetg(1,t_VECSMALL);
    1976      397495 :   nf = checknf(nf);
    1977      397495 :   if (typ(x) == t_MAT)
    1978             :   { /* factorisation */
    1979      184329 :     GEN g = gel(x,1), e = gel(x,2);
    1980      184329 :     long l = lg(g);
    1981      184329 :     V = zero_zv(n);
    1982      414837 :     for (i = 1; i < l; i++)
    1983      230508 :       if (mpodd(gel(e,i)))
    1984      206208 :         Flv_add_inplace(V, nfsign_arch(nf,gel(g,i),archp), 2);
    1985      184329 :     set_avma((pari_sp)V); return V;
    1986             :   }
    1987      213166 :   av = avma; V = cgetg(n+1,t_VECSMALL);
    1988      213166 :   x = nf_to_scalar_or_basis(nf, x);
    1989      213166 :   switch(typ(x))
    1990             :   {
    1991       29470 :     case t_INT:
    1992       29470 :       s = signe(x);
    1993       29470 :       if (!s) pari_err_DOMAIN("nfsign_arch","element","=",gen_0,x);
    1994       29470 :       set_avma(av); return const_vecsmall(n, (s < 0)? 1: 0);
    1995         105 :     case t_FRAC:
    1996         105 :       s = signe(gel(x,1));
    1997         105 :       set_avma(av); return const_vecsmall(n, (s < 0)? 1: 0);
    1998             :   }
    1999      183591 :   x = Q_primpart(x); M = nf_get_M(nf); sarch = NULL; np = -1;
    2000      526535 :   for (i = 1; i <= n; i++)
    2001             :   {
    2002      344380 :     long s = eval_sign(M, x, archp[i]);
    2003      344380 :     if (s < 0) /* failure */
    2004             :     {
    2005        1438 :       long ni, r1 = nf_get_r1(nf);
    2006             :       GEN xi;
    2007        1438 :       if (np < 0)
    2008             :       {
    2009        1438 :         np = num_positive(nf, x);
    2010        1438 :         if (np == 0) { set_avma(av); return const_vecsmall(n, 1); }
    2011        1229 :         if (np == r1){ set_avma(av); return const_vecsmall(n, 0); }
    2012         730 :         sarch = nfarchstar(nf, NULL, identity_perm(r1));
    2013             :       }
    2014         730 :       xi = set_sign_mod_divisor(nf, vecsmall_ei(r1, archp[i]), gen_1, sarch);
    2015         730 :       xi = Q_primpart(xi);
    2016         730 :       ni = num_positive(nf, nfmuli(nf,x,xi));
    2017         730 :       if (ni == 0) { set_avma(av); V = const_vecsmall(n, 1); V[i] = 0; return V; }
    2018         464 :       if (ni == r1){ set_avma(av); V = const_vecsmall(n, 0); V[i] = 1; return V; }
    2019           2 :       s = ni < np? 0: 1;
    2020             :     }
    2021      342944 :     V[i] = s;
    2022             :   }
    2023      182155 :   set_avma((pari_sp)V); return V;
    2024             : }
    2025             : static void
    2026        6384 : chk_ind(const char *s, long i, long r1)
    2027             : {
    2028        6384 :   if (i <= 0) pari_err_DOMAIN(s, "index", "<=", gen_0, stoi(i));
    2029        6370 :   if (i > r1) pari_err_DOMAIN(s, "index", ">", utoi(r1), utoi(i));
    2030        6335 : }
    2031             : static GEN
    2032        6307 : parse_embed(GEN ind, long r, const char *f)
    2033             : {
    2034             :   long l, i;
    2035        6307 :   if (!ind) return identity_perm(r);
    2036        4354 :   switch(typ(ind))
    2037             :   {
    2038         154 :     case t_INT: case t_VEC: case t_COL: ind = gtovecsmall(ind); break;
    2039        4200 :     case t_VECSMALL: break;
    2040           0 :     default: pari_err_TYPE(f, ind);
    2041             :   }
    2042        4354 :   l = lg(ind);
    2043       10689 :   for (i = 1; i < l; i++) chk_ind(f, ind[i], r);
    2044        4305 :   return ind;
    2045             : }
    2046             : GEN
    2047        4774 : nfeltsign(GEN nf, GEN x, GEN ind0)
    2048             : {
    2049        4774 :   pari_sp av = avma;
    2050             :   long i, l;
    2051             :   GEN v, ind;
    2052        4774 :   nf = checknf(nf);
    2053        4774 :   ind = parse_embed(ind0, nf_get_r1(nf), "nfeltsign");
    2054        4753 :   l = lg(ind);
    2055        4753 :   if (is_rational_t(typ(x)))
    2056             :   { /* nfsign_arch would test this, but avoid converting t_VECSMALL -> t_VEC */
    2057             :     GEN s;
    2058        2086 :     switch(gsigne(x))
    2059             :     {
    2060         525 :       case -1:s = gen_m1; break;
    2061        1554 :       case 1: s = gen_1; break;
    2062           7 :       default: s = gen_0; break;
    2063             :     }
    2064        2086 :     set_avma(av);
    2065        2086 :     return (ind0 && typ(ind0) == t_INT)? s: const_vec(l-1, s);
    2066             :   }
    2067        2667 :   v = nfsign_arch(nf, x, ind);
    2068        2667 :   if (ind0 && typ(ind0) == t_INT) { set_avma(av); return v[1]? gen_m1: gen_1; }
    2069        2653 :   settyp(v, t_VEC);
    2070        7238 :   for (i = 1; i < l; i++) gel(v,i) = v[i]? gen_m1: gen_1;
    2071        2653 :   return gerepileupto(av, v);
    2072             : }
    2073             : 
    2074             : GEN
    2075          63 : nfeltembed(GEN nf, GEN x, GEN ind0, long prec0)
    2076             : {
    2077          63 :   pari_sp av = avma;
    2078             :   long i, e, l, r1, r2, prec, prec1;
    2079             :   GEN v, ind, cx;
    2080          63 :   nf = checknf(nf); nf_get_sign(nf,&r1,&r2);
    2081          63 :   x = nf_to_scalar_or_basis(nf, x);
    2082          56 :   ind = parse_embed(ind0, r1+r2, "nfeltembed");
    2083          49 :   l = lg(ind);
    2084          49 :   if (typ(x) != t_COL)
    2085             :   {
    2086           0 :     if (!(ind0 && typ(ind0) == t_INT)) x = const_vec(l-1, x);
    2087           0 :     return gerepilecopy(av, x);
    2088             :   }
    2089          49 :   x = Q_primitive_part(x, &cx);
    2090          49 :   prec1 = prec0; e = gexpo(x);
    2091          49 :   if (e > 8) prec1 += nbits2extraprec(e);
    2092          49 :   prec = prec1;
    2093          49 :   if (nf_get_prec(nf) < prec) nf = nfnewprec_shallow(nf, prec);
    2094          49 :   v = cgetg(l, t_VEC);
    2095             :   for(;;)
    2096           7 :   {
    2097          56 :     GEN M = nf_get_M(nf);
    2098         140 :     for (i = 1; i < l; i++)
    2099             :     {
    2100          91 :       GEN t = nfembed_i(M, x, ind[i]);
    2101          91 :       long e = gexpo(t);
    2102          91 :       if (gequal0(t) || precision(t) < prec0
    2103          91 :                      || (e < 0 && prec < prec1 + nbits2extraprec(-e)) ) break;
    2104          84 :       if (cx) t = gmul(t, cx);
    2105          84 :       gel(v,i) = t;
    2106             :     }
    2107          56 :     if (i == l) break;
    2108           7 :     prec = precdbl(prec);
    2109           7 :     if (DEBUGLEVEL>1) pari_warn(warnprec,"eltnfembed", prec);
    2110           7 :     nf = nfnewprec_shallow(nf, prec);
    2111             :   }
    2112          49 :   if (ind0 && typ(ind0) == t_INT) v = gel(v,1);
    2113          49 :   return gerepilecopy(av, v);
    2114             : }
    2115             : 
    2116             : /* number of distinct roots of sigma(f) */
    2117             : GEN
    2118        1477 : nfpolsturm(GEN nf, GEN f, GEN ind0)
    2119             : {
    2120        1477 :   pari_sp av = avma;
    2121             :   long d, l, r1, single;
    2122             :   GEN ind, u, v, vr1, T, s, t;
    2123             : 
    2124        1477 :   nf = checknf(nf); T = nf_get_pol(nf); r1 = nf_get_r1(nf);
    2125        1477 :   ind = parse_embed(ind0, r1, "nfpolsturm");
    2126        1456 :   single = ind0 && typ(ind0) == t_INT;
    2127        1456 :   l = lg(ind);
    2128             : 
    2129        1456 :   if (gequal0(f)) pari_err_ROOTS0("nfpolsturm");
    2130        1449 :   if (typ(f) == t_POL && varn(f) != varn(T))
    2131             :   {
    2132        1428 :     f = RgX_nffix("nfpolsturm", T, f,1);
    2133        1428 :     if (lg(f) == 3) f = NULL;
    2134             :   }
    2135             :   else
    2136             :   {
    2137          21 :     (void)Rg_nffix("nfpolsturm", T, f, 0);
    2138          21 :     f = NULL;
    2139             :   }
    2140        1449 :   if (!f) { set_avma(av); return single? gen_0: zerovec(l-1); }
    2141        1428 :   d = degpol(f);
    2142        1428 :   if (d == 1) { set_avma(av); return single? gen_1: const_vec(l-1,gen_1); }
    2143             : 
    2144        1393 :   vr1 = const_vecsmall(l-1, 1);
    2145        1393 :   u = Q_primpart(f); s = ZV_to_zv(nfeltsign(nf, gel(u,d+2), ind));
    2146        1393 :   v = RgX_deriv(u); t = odd(d)? leafcopy(s): zv_neg(s);
    2147             :   for(;;)
    2148         182 :   {
    2149        1575 :     GEN r = RgX_neg( Q_primpart(RgX_pseudorem(u, v)) ), sr;
    2150        1575 :     long i, dr = degpol(r);
    2151        1575 :     if (dr < 0) break;
    2152        1575 :     sr = ZV_to_zv(nfeltsign(nf, gel(r,dr+2), ind));
    2153        3934 :     for (i = 1; i < l; i++)
    2154        2359 :       if (sr[i] != s[i]) { s[i] = sr[i], vr1[i]--; }
    2155        1575 :     if (odd(dr)) sr = zv_neg(sr);
    2156        3934 :     for (i = 1; i < l; i++)
    2157        2359 :       if (sr[i] != t[i]) { t[i] = sr[i], vr1[i]++; }
    2158        1575 :     if (!dr) break;
    2159         182 :     u = v; v = r;
    2160             :   }
    2161        1393 :   if (single) { set_avma(av); return stoi(vr1[1]); }
    2162        1386 :   return gerepileupto(av, zv_to_ZV(vr1));
    2163             : }
    2164             : 
    2165             : /* return the vector of signs of x; the matrix of such if x is a vector
    2166             :  * of nf elements */
    2167             : GEN
    2168        2156 : nfsign(GEN nf, GEN x)
    2169             : {
    2170             :   long i, l;
    2171             :   GEN archp, S;
    2172             : 
    2173        2156 :   nf = checknf(nf);
    2174        2156 :   archp = identity_perm( nf_get_r1(nf) );
    2175        2156 :   if (typ(x) != t_VEC) return nfsign_arch(nf, x, archp);
    2176         924 :   l = lg(x); S = cgetg(l, t_MAT);
    2177        4557 :   for (i=1; i<l; i++) gel(S,i) = nfsign_arch(nf, gel(x,i), archp);
    2178         924 :   return S;
    2179             : }
    2180             : 
    2181             : /* x integral elt, A integral ideal in HNF; reduce x mod A */
    2182             : static GEN
    2183     1130527 : zk_modHNF(GEN x, GEN A)
    2184     1130527 : { return (typ(x) == t_COL)?  ZC_hnfrem(x, A): modii(x, gcoeff(A,1,1)); }
    2185             : 
    2186             : /* given an element x in Z_K and an integral ideal y in HNF, coprime with x,
    2187             :    outputs an element inverse of x modulo y */
    2188             : GEN
    2189          91 : nfinvmodideal(GEN nf, GEN x, GEN y)
    2190             : {
    2191          91 :   pari_sp av = avma;
    2192          91 :   GEN a, yZ = gcoeff(y,1,1);
    2193             : 
    2194          91 :   if (equali1(yZ)) return gen_0;
    2195          91 :   x = nf_to_scalar_or_basis(nf, x);
    2196          91 :   if (typ(x) == t_INT) return gerepileupto(av, Fp_inv(x, yZ));
    2197             : 
    2198          49 :   a = hnfmerge_get_1(idealhnf_principal(nf,x), y);
    2199          49 :   if (!a) pari_err_INV("nfinvmodideal", x);
    2200          49 :   return gerepileupto(av, zk_modHNF(nfdiv(nf,a,x), y));
    2201             : }
    2202             : 
    2203             : static GEN
    2204      534704 : nfsqrmodideal(GEN nf, GEN x, GEN id)
    2205      534704 : { return zk_modHNF(nfsqri(nf,x), id); }
    2206             : static GEN
    2207     1478986 : nfmulmodideal(GEN nf, GEN x, GEN y, GEN id)
    2208     1478986 : { return x? zk_modHNF(nfmuli(nf,x,y), id): y; }
    2209             : /* assume x integral, k integer, A in HNF */
    2210             : GEN
    2211      994773 : nfpowmodideal(GEN nf,GEN x,GEN k,GEN A)
    2212             : {
    2213      994773 :   long s = signe(k);
    2214             :   pari_sp av;
    2215             :   GEN y;
    2216             : 
    2217      994773 :   if (!s) return gen_1;
    2218      994773 :   av = avma;
    2219      994773 :   x = nf_to_scalar_or_basis(nf, x);
    2220      994773 :   if (typ(x) != t_COL) return Fp_pow(x, k, gcoeff(A,1,1));
    2221      441759 :   if (s < 0) { x = nfinvmodideal(nf, x,A); k = negi(k); }
    2222      441759 :   for(y = NULL;;)
    2223             :   {
    2224      976463 :     if (mpodd(k)) y = nfmulmodideal(nf,y,x,A);
    2225      976463 :     k = shifti(k,-1); if (!signe(k)) break;
    2226      534704 :     x = nfsqrmodideal(nf,x,A);
    2227             :   }
    2228      441759 :   return gerepileupto(av, y);
    2229             : }
    2230             : 
    2231             : /* a * g^n mod id */
    2232             : static GEN
    2233      785388 : nfmulpowmodideal(GEN nf, GEN a, GEN g, GEN n, GEN id)
    2234             : {
    2235      785388 :   return nfmulmodideal(nf, a, nfpowmodideal(nf,g,n,id), id);
    2236             : }
    2237             : 
    2238             : /* assume (num(g[i]), id) = 1 for all i. Return prod g[i]^e[i] mod id.
    2239             :  * EX = multiple of exponent of (O_K/id)^* */
    2240             : GEN
    2241      464707 : famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id, GEN EX)
    2242             : {
    2243      464707 :   GEN EXo2, plus = NULL, minus = NULL, idZ = gcoeff(id,1,1);
    2244      464707 :   long i, lx = lg(g);
    2245             : 
    2246      464707 :   if (equali1(idZ)) return gen_1; /* id = Z_K */
    2247      464252 :   EXo2 = (expi(EX) > 10)? shifti(EX,-1): NULL;
    2248     1440767 :   for (i = 1; i < lx; i++)
    2249             :   {
    2250      976515 :     GEN h, n = centermodii(gel(e,i), EX, EXo2);
    2251      976515 :     long sn = signe(n);
    2252      976515 :     if (!sn) continue;
    2253             : 
    2254      725931 :     h = nf_to_scalar_or_basis(nf, gel(g,i));
    2255      725931 :     switch(typ(h))
    2256             :     {
    2257      403717 :       case t_INT: break;
    2258           0 :       case t_FRAC:
    2259           0 :         h = Fp_div(gel(h,1), gel(h,2), idZ); break;
    2260      322214 :       default:
    2261             :       {
    2262             :         GEN dh;
    2263      322214 :         h = Q_remove_denom(h, &dh);
    2264      322214 :         if (dh) h = FpC_Fp_mul(h, Fp_inv(dh,idZ), idZ);
    2265             :       }
    2266             :     }
    2267      725931 :     if (sn > 0)
    2268      724384 :       plus = nfmulpowmodideal(nf, plus, h, n, id);
    2269             :     else /* sn < 0 */
    2270        1547 :       minus = nfmulpowmodideal(nf, minus, h, negi(n), id);
    2271             :   }
    2272      464252 :   if (minus) plus = nfmulmodideal(nf, plus, nfinvmodideal(nf,minus,id), id);
    2273      464252 :   return plus? plus: gen_1;
    2274             : }
    2275             : 
    2276             : /* given 2 integral ideals x, y in HNF s.t x | y | x^2, compute (1+x)/(1+y) in
    2277             :  * the form [[cyc],[gen], U], where U := ux^-1 as a pair [ZM, denom(U)] */
    2278             : static GEN
    2279       22820 : zidealij(GEN x, GEN y)
    2280             : {
    2281       22820 :   GEN U, G, cyc, xp = gcoeff(x,1,1), xi = hnf_invscale(x, xp);
    2282             :   long j, N;
    2283             : 
    2284             :   /* x^(-1) y = relations between the 1 + x_i (HNF) */
    2285       22820 :   cyc = ZM_snf_group(ZM_Z_divexact(ZM_mul(xi, y), xp), &U, &G);
    2286       22820 :   N = lg(cyc); G = ZM_mul(x,G); settyp(G, t_VEC); /* new generators */
    2287       54271 :   for (j=1; j<N; j++)
    2288             :   {
    2289       31451 :     GEN c = gel(G,j);
    2290       31451 :     gel(c,1) = addiu(gel(c,1), 1); /* 1 + g_j */
    2291       31451 :     if (ZV_isscalar(c)) gel(G,j) = gel(c,1);
    2292             :   }
    2293       22820 :   return mkvec4(cyc, G, ZM_mul(U,xi), xp);
    2294             : }
    2295             : 
    2296             : /* lg(x) > 1, x + 1; shallow */
    2297             : static GEN
    2298       13055 : ZC_add1(GEN x)
    2299             : {
    2300       13055 :   long i, l = lg(x);
    2301       13055 :   GEN y = cgetg(l, t_COL);
    2302       32466 :   for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
    2303       13055 :   gel(y,1) = addiu(gel(x,1), 1); return y;
    2304             : }
    2305             : /* lg(x) > 1, x - 1; shallow */
    2306             : static GEN
    2307        6349 : ZC_sub1(GEN x)
    2308             : {
    2309        6349 :   long i, l = lg(x);
    2310        6349 :   GEN y = cgetg(l, t_COL);
    2311       17073 :   for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
    2312        6349 :   gel(y,1) = subiu(gel(x,1), 1); return y;
    2313             : }
    2314             : 
    2315             : /* x,y are t_INT or ZC */
    2316             : static GEN
    2317           0 : zkadd(GEN x, GEN y)
    2318             : {
    2319           0 :   long tx = typ(x);
    2320           0 :   if (tx == typ(y))
    2321           0 :     return tx == t_INT? addii(x,y): ZC_add(x,y);
    2322             :   else
    2323           0 :     return tx == t_INT? ZC_Z_add(y,x): ZC_Z_add(x,y);
    2324             : }
    2325             : /* x a t_INT or ZC, x+1; shallow */
    2326             : static GEN
    2327       27643 : zkadd1(GEN x)
    2328             : {
    2329       27643 :   long tx = typ(x);
    2330       27643 :   return tx == t_INT? addiu(x,1): ZC_add1(x);
    2331             : }
    2332             : /* x a t_INT or ZC, x-1; shallow */
    2333             : static GEN
    2334       27643 : zksub1(GEN x)
    2335             : {
    2336       27643 :   long tx = typ(x);
    2337       27643 :   return tx == t_INT? subiu(x,1): ZC_sub1(x);
    2338             : }
    2339             : /* x,y are t_INT or ZC; x - y */
    2340             : static GEN
    2341           0 : zksub(GEN x, GEN y)
    2342             : {
    2343           0 :   long tx = typ(x), ty = typ(y);
    2344           0 :   if (tx == ty)
    2345           0 :     return tx == t_INT? subii(x,y): ZC_sub(x,y);
    2346             :   else
    2347           0 :     return tx == t_INT? Z_ZC_sub(x,y): ZC_Z_sub(x,y);
    2348             : }
    2349             : /* x is t_INT or ZM (mult. map), y is t_INT or ZC; x * y */
    2350             : static GEN
    2351       27643 : zkmul(GEN x, GEN y)
    2352             : {
    2353       27643 :   long tx = typ(x), ty = typ(y);
    2354       27643 :   if (ty == t_INT)
    2355       21294 :     return tx == t_INT? mulii(x,y): ZC_Z_mul(gel(x,1),y);
    2356             :   else
    2357        6349 :     return tx == t_INT? ZC_Z_mul(y,x): ZM_ZC_mul(x,y);
    2358             : }
    2359             : 
    2360             : /* (U,V) = 1 coprime ideals. Want z = x mod U, = y mod V; namely
    2361             :  * z =vx + uy = v(x-y) + y, where u + v = 1, u in U, v in V.
    2362             :  * zkc = [v, UV], v a t_INT or ZM (mult. by v map), UV a ZM (ideal in HNF);
    2363             :  * shallow */
    2364             : GEN
    2365           0 : zkchinese(GEN zkc, GEN x, GEN y)
    2366             : {
    2367           0 :   GEN v = gel(zkc,1), UV = gel(zkc,2), z = zkadd(zkmul(v, zksub(x,y)), y);
    2368           0 :   return zk_modHNF(z, UV);
    2369             : }
    2370             : /* special case z = x mod U, = 1 mod V; shallow */
    2371             : GEN
    2372       27643 : zkchinese1(GEN zkc, GEN x)
    2373             : {
    2374       27643 :   GEN v = gel(zkc,1), UV = gel(zkc,2), z = zkadd1(zkmul(v, zksub1(x)));
    2375       27643 :   return (typ(z) == t_INT)? z: ZC_hnfrem(z, UV);
    2376             : }
    2377             : static GEN
    2378       25536 : zkVchinese1(GEN zkc, GEN v)
    2379             : {
    2380             :   long i, ly;
    2381       25536 :   GEN y = cgetg_copy(v, &ly);
    2382       53179 :   for (i=1; i<ly; i++) gel(y,i) = zkchinese1(zkc, gel(v,i));
    2383       25536 :   return y;
    2384             : }
    2385             : 
    2386             : /* prepare to solve z = x (mod A), z = y mod (B) [zkchinese or zkchinese1] */
    2387             : GEN
    2388       25277 : zkchineseinit(GEN nf, GEN A, GEN B, GEN AB)
    2389             : {
    2390             :   GEN v;
    2391             :   long e;
    2392       25277 :   nf = checknf(nf);
    2393       25277 :   v = idealaddtoone_raw(nf, A, B);
    2394       25277 :   if ((e = gexpo(v)) > 5)
    2395             :   {
    2396        1519 :     GEN b = (typ(v) == t_COL)? v: scalarcol_shallow(v, nf_get_degree(nf));
    2397        1519 :     b= ZC_reducemodlll(b, AB);
    2398        1519 :     if (gexpo(b) < e) v = b;
    2399             :   }
    2400       25277 :   return mkvec2(zk_scalar_or_multable(nf,v), AB);
    2401             : }
    2402             : /* prepare to solve z = x (mod A), z = 1 mod (B)
    2403             :  * and then         z = 1 (mod A), z = y mod (B) [zkchinese1 twice] */
    2404             : static GEN
    2405         259 : zkchinese1init2(GEN nf, GEN A, GEN B, GEN AB)
    2406             : {
    2407         259 :   GEN zkc = zkchineseinit(nf, A, B, AB);
    2408         259 :   GEN mv = gel(zkc,1), mu;
    2409         259 :   if (typ(mv) == t_INT) return mkvec2(zkc, mkvec2(subui(1,mv),AB));
    2410          35 :   mu = RgM_Rg_add_shallow(ZM_neg(mv), gen_1);
    2411          35 :   return mkvec2(mkvec2(mv,AB), mkvec2(mu,AB));
    2412             : }
    2413             : 
    2414             : static GEN
    2415      332850 : apply_U(GEN L, GEN a)
    2416             : {
    2417      332850 :   GEN e, U = gel(L,3), dU = gel(L,4);
    2418      332850 :   if (typ(a) == t_INT)
    2419      142077 :     e = ZC_Z_mul(gel(U,1), subiu(a, 1));
    2420             :   else
    2421             :   { /* t_COL */
    2422      190773 :     GEN t = shallowcopy(a);
    2423      190773 :     gel(t,1) = subiu(gel(t,1), 1); /* t = a - 1 */
    2424      190773 :     e = ZM_ZC_mul(U, t);
    2425             :   }
    2426      332850 :   return gdiv(e, dU);
    2427             : }
    2428             : 
    2429             : /* true nf; vectors of [[cyc],[g],U.X^-1]. Assume k > 1. */
    2430             : static GEN
    2431       15568 : principal_units(GEN nf, GEN pr, long k, GEN prk)
    2432             : {
    2433             :   GEN list, prb;
    2434       15568 :   ulong mask = quadratic_prec_mask(k);
    2435       15568 :   long a = 1;
    2436             : 
    2437       15568 :   prb = pr_hnf(nf,pr);
    2438       15568 :   list = vectrunc_init(k);
    2439       38388 :   while (mask > 1)
    2440             :   {
    2441       22820 :     GEN pra = prb;
    2442       22820 :     long b = a << 1;
    2443             : 
    2444       22820 :     if (mask & 1) b--;
    2445       22820 :     mask >>= 1;
    2446             :     /* compute 1 + pr^a / 1 + pr^b, 2a <= b */
    2447       22820 :     prb = (b >= k)? prk: idealpows(nf,pr,b);
    2448       22820 :     vectrunc_append(list, zidealij(pra, prb));
    2449       22820 :     a = b;
    2450             :   }
    2451       15568 :   return list;
    2452             : }
    2453             : /* a = 1 mod (pr) return log(a) on local-gens of 1+pr/1+pr^k */
    2454             : static GEN
    2455      209608 : log_prk1(GEN nf, GEN a, long nh, GEN L2, GEN prk)
    2456             : {
    2457      209608 :   GEN y = cgetg(nh+1, t_COL);
    2458      209608 :   long j, iy, c = lg(L2)-1;
    2459      542451 :   for (j = iy = 1; j <= c; j++)
    2460             :   {
    2461      332850 :     GEN L = gel(L2,j), cyc = gel(L,1), gen = gel(L,2), E = apply_U(L,a);
    2462      332850 :     long i, nc = lg(cyc)-1;
    2463      332850 :     int last = (j == c);
    2464     1029098 :     for (i = 1; i <= nc; i++, iy++)
    2465             :     {
    2466      696255 :       GEN t, e = gel(E,i);
    2467      696255 :       if (typ(e) != t_INT) pari_err_COPRIME("zlog_prk1", a, prk);
    2468      696248 :       t = Fp_neg(e, gel(cyc,i));
    2469      696248 :       gel(y,iy) = negi(t);
    2470      696248 :       if (!last && signe(t)) a = nfmulpowmodideal(nf, a, gel(gen,i), t, prk);
    2471             :     }
    2472             :   }
    2473      209601 :   return y;
    2474             : }
    2475             : /* true nf */
    2476             : static GEN
    2477        5999 : principal_units_relations(GEN nf, GEN L2, GEN prk, long nh)
    2478             : {
    2479        5999 :   GEN h = cgetg(nh+1,t_MAT);
    2480        5999 :   long ih, j, c = lg(L2)-1;
    2481       19250 :   for (j = ih = 1; j <= c; j++)
    2482             :   {
    2483       13251 :     GEN L = gel(L2,j), F = gel(L,1), G = gel(L,2);
    2484       13251 :     long k, lG = lg(G);
    2485       33194 :     for (k = 1; k < lG; k++,ih++)
    2486             :     { /* log(g^f) mod pr^e */
    2487       19943 :       GEN a = nfpowmodideal(nf,gel(G,k),gel(F,k),prk);
    2488       19943 :       gel(h,ih) = ZC_neg(log_prk1(nf, a, nh, L2, prk));
    2489       19943 :       gcoeff(h,ih,ih) = gel(F,k);
    2490             :     }
    2491             :   }
    2492        5999 :   return h;
    2493             : }
    2494             : /* true nf; k > 1; multiplicative group (1 + pr) / (1 + pr^k) */
    2495             : static GEN
    2496       15568 : idealprincipalunits_i(GEN nf, GEN pr, long k, GEN *pU)
    2497             : {
    2498       15568 :   GEN cyc, gen, L2, prk = idealpows(nf, pr, k);
    2499             : 
    2500       15568 :   L2 = principal_units(nf, pr, k, prk);
    2501       15568 :   if (k == 2)
    2502             :   {
    2503        9569 :     GEN L = gel(L2,1);
    2504        9569 :     cyc = gel(L,1);
    2505        9569 :     gen = gel(L,2);
    2506        9569 :     if (pU) *pU = matid(lg(gen)-1);
    2507             :   }
    2508             :   else
    2509             :   {
    2510        5999 :     long c = lg(L2), j;
    2511        5999 :     GEN EX, h, Ui, vg = cgetg(c, t_VEC);
    2512       19250 :     for (j = 1; j < c; j++) gel(vg, j) = gmael(L2,j,2);
    2513        5999 :     vg = shallowconcat1(vg);
    2514        5999 :     h = principal_units_relations(nf, L2, prk, lg(vg)-1);
    2515        5999 :     h = ZM_hnfall_i(h, NULL, 0);
    2516        5999 :     cyc = ZM_snf_group(h, pU, &Ui);
    2517        5999 :     c = lg(Ui); gen = cgetg(c, t_VEC); EX = cyc_get_expo(cyc);
    2518       20979 :     for (j = 1; j < c; j++)
    2519       14980 :       gel(gen,j) = famat_to_nf_modideal_coprime(nf, vg, gel(Ui,j), prk, EX);
    2520             :   }
    2521       15568 :   return mkvec4(cyc, gen, prk, L2);
    2522             : }
    2523             : GEN
    2524         112 : idealprincipalunits(GEN nf, GEN pr, long k)
    2525             : {
    2526             :   pari_sp av;
    2527             :   GEN v;
    2528         112 :   nf = checknf(nf);
    2529         112 :   if (k == 1) { checkprid(pr); retmkvec3(gen_1,cgetg(1,t_VEC),cgetg(1,t_VEC)); }
    2530         105 :   av = avma; v = idealprincipalunits_i(nf, pr, k, NULL);
    2531         105 :   return gerepilecopy(av, mkvec3(powiu(pr_norm(pr), k-1), gel(v,1), gel(v,2)));
    2532             : }
    2533             : 
    2534             : /* true nf; given an ideal pr^k dividing an integral ideal x (in HNF form)
    2535             :  * compute an 'sprk', the structure of G = (Z_K/pr^k)^* [ x = NULL for x=pr^k ]
    2536             :  * Return a vector with at least 4 components [cyc],[gen],[HNF pr^k,pr,k],ff,
    2537             :  * where
    2538             :  * cyc : type of G as abelian group (SNF)
    2539             :  * gen : generators of G, coprime to x
    2540             :  * pr^k: in HNF
    2541             :  * ff  : data for log_g in (Z_K/pr)^*
    2542             :  * Two extra components are present iff k > 1: L2, U
    2543             :  * L2  : list of data structures to compute local DL in (Z_K/pr)^*,
    2544             :  *       and 1 + pr^a/ 1 + pr^b for various a < b <= min(2a, k)
    2545             :  * U   : base change matrices to convert a vector of local DL to DL wrt gen
    2546             :  * If MOD is not NULL, initialize G / G^MOD instead */
    2547             : static GEN
    2548       47838 : sprkinit(GEN nf, GEN pr, long k, GEN x, GEN MOD)
    2549             : {
    2550       47838 :   GEN T, p, Ld, modpr, cyc, gen, g, g0, A, prk, U, L2, ord0 = NULL;
    2551       47838 :   long f = pr_get_f(pr);
    2552             : 
    2553       47838 :   if(DEBUGLEVEL>3) err_printf("treating pr^%ld, pr = %Ps\n",k,pr);
    2554       47838 :   modpr = nf_to_Fq_init(nf, &pr,&T,&p);
    2555       47838 :   if (MOD)
    2556             :   {
    2557        4823 :     GEN A = subiu(powiu(p,f), 1), d = gcdii(A, MOD), fa = Z_factor(d);
    2558        4823 :     ord0 = mkvec2(A, fa); /* true order, factorization of order in G/G^MOD */
    2559        4823 :     Ld = gel(fa,1);
    2560        4823 :     if (lg(Ld) > 1 && equaliu(gel(Ld,1),2)) Ld = vecslice(Ld,2,lg(Ld)-1);
    2561             :   }
    2562             :   /* (Z_K / pr)^* */
    2563       47838 :   if (f == 1)
    2564             :   {
    2565       41160 :     g0 = g = MOD? pgener_Fp_local(p, Ld): pgener_Fp(p);
    2566       41160 :     if (!ord0) ord0 = get_arith_ZZM(subiu(p,1));
    2567             :   }
    2568             :   else
    2569             :   {
    2570        6678 :     g0 = g = MOD? gener_FpXQ_local(T, p, Ld): gener_FpXQ(T,p, &ord0);
    2571        6678 :     g = Fq_to_nf(g, modpr);
    2572        6678 :     if (typ(g) == t_POL) g = poltobasis(nf, g);
    2573             :   }
    2574       47838 :   A = gel(ord0, 1); /* Norm(pr)-1 */
    2575             :   /* If MOD != NULL, d = gcd(A, MOD): g^(A/d) has order d */
    2576       47838 :   if (k == 1)
    2577             :   {
    2578       32375 :     cyc = mkvec(A);
    2579       32375 :     gen = mkvec(g);
    2580       32375 :     prk = pr_hnf(nf,pr);
    2581       32375 :     L2 = U = NULL;
    2582             :   }
    2583             :   else
    2584             :   { /* local-gens of (1 + pr)/(1 + pr^k) = SNF-gens * U */
    2585             :     GEN AB, B, u, v, w;
    2586             :     long j, l;
    2587       15463 :     w = idealprincipalunits_i(nf, pr, k, &U);
    2588             :     /* incorporate (Z_K/pr)^*, order A coprime to B = expo(1+pr/1+pr^k)*/
    2589       15463 :     cyc = leafcopy(gel(w,1)); B = cyc_get_expo(cyc); AB = mulii(A,B);
    2590       15463 :     gen = leafcopy(gel(w,2));
    2591       15463 :     prk = gel(w,3);
    2592       15463 :     g = nfpowmodideal(nf, g, B, prk);
    2593       15463 :     g0 = Fq_pow(g0, modii(B,A), T, p); /* update primitive root */
    2594       15463 :     L2 = mkvec3(A, g, gel(w,4));
    2595       15463 :     gel(cyc,1) = AB;
    2596       15463 :     gel(gen,1) = nfmulmodideal(nf, gel(gen,1), g, prk);
    2597       15463 :     u = mulii(Fp_inv(A,B), A);
    2598       15463 :     v = subui(1, u); l = lg(U);
    2599       46536 :     for (j = 1; j < l; j++) gcoeff(U,1,j) = Fp_mul(u, gcoeff(U,1,j), AB);
    2600       15463 :     U = mkvec2(Rg_col_ei(v, lg(gen)-1, 1), U);
    2601             :   }
    2602             :   /* local-gens of (Z_K/pr^k)^* = SNF-gens * U */
    2603       47838 :   if (x)
    2604             :   {
    2605       25018 :     GEN uv = zkchineseinit(nf, idealmulpowprime(nf,x,pr,utoineg(k)), prk, x);
    2606       25018 :     gen = zkVchinese1(uv, gen);
    2607             :   }
    2608       47838 :   return mkvecn(U? 6: 4, cyc, gen, prk, mkvec3(modpr,g0,ord0), L2, U);
    2609             : }
    2610             : static GEN
    2611      674714 : sprk_get_cyc(GEN s) { return gel(s,1); }
    2612             : static GEN
    2613      349207 : sprk_get_expo(GEN s) { return cyc_get_expo(sprk_get_cyc(s)); }
    2614             : static GEN
    2615       45101 : sprk_get_gen(GEN s) { return gel(s,2); }
    2616             : static GEN
    2617      829972 : sprk_get_prk(GEN s) { return gel(s,3); }
    2618             : static GEN
    2619      431086 : sprk_get_ff(GEN s) { return gel(s,4); }
    2620             : static GEN
    2621      383556 : sprk_get_pr(GEN s) { GEN ff = gel(s,4); return modpr_get_pr(gel(ff,1)); }
    2622             : /* L2 to 1 + pr / 1 + pr^k */
    2623             : static GEN
    2624      207942 : sprk_get_L2(GEN s) { return gmael(s,5,3); }
    2625             : /* lift to nf of primitive root of k(pr) */
    2626             : static GEN
    2627       10500 : sprk_get_gnf(GEN s) { return gmael(s,5,2); }
    2628             : static void
    2629      190575 : sprk_get_U2(GEN s, GEN *U1, GEN *U2)
    2630      190575 : { GEN v = gel(s,6); *U1 = gel(v,1); *U2 = gel(v,2); }
    2631             : static int
    2632      431086 : sprk_is_prime(GEN s) { return lg(s) == 5; }
    2633             : 
    2634             : static GEN
    2635      349144 : famat_zlog_pr(GEN nf, GEN g, GEN e, GEN sprk, GEN mod)
    2636             : {
    2637      349144 :   GEN x, expo = sprk_get_expo(sprk);
    2638      349144 :   if (mod) expo = gcdii(expo,mod);
    2639      349144 :   x = famat_makecoprime(nf, g, e, sprk_get_pr(sprk), sprk_get_prk(sprk), expo);
    2640      349144 :   return log_prk(nf, x, sprk, mod);
    2641             : }
    2642             : /* famat_zlog_pr assuming (g,sprk.pr) = 1 */
    2643             : static GEN
    2644          63 : famat_zlog_pr_coprime(GEN nf, GEN g, GEN e, GEN sprk, GEN MOD)
    2645             : {
    2646          63 :   GEN x = famat_to_nf_modideal_coprime(nf, g, e, sprk_get_prk(sprk),
    2647             :                                        sprk_get_expo(sprk));
    2648          63 :   return log_prk(nf, x, sprk, MOD);
    2649             : }
    2650             : 
    2651             : /* o t_INT, O = [ord,fa] format for multiple of o (for Fq_log);
    2652             :  * return o in [ord,fa] format */
    2653             : static GEN
    2654      123208 : order_update(GEN o, GEN O)
    2655             : {
    2656      123208 :   GEN p = gmael(O,2,1), z = o, P, E;
    2657      123208 :   long i, j, l = lg(p);
    2658      123208 :   P = cgetg(l, t_COL);
    2659      123208 :   E = cgetg(l, t_COL);
    2660      178103 :   for (i = j = 1; i < l; i++)
    2661             :   {
    2662      178103 :     long v = Z_pvalrem(z, gel(p,i), &z);
    2663      178103 :     if (v)
    2664             :     {
    2665      159896 :       gel(P,j) = gel(p,i);
    2666      159896 :       gel(E,j) = utoipos(v); j++;
    2667      159896 :       if (is_pm1(z)) break;
    2668             :     }
    2669             :   }
    2670      123208 :   setlg(P, j);
    2671      123208 :   setlg(E, j); return mkvec2(o, mkmat2(P,E));
    2672             : }
    2673             : 
    2674             : /* a in Z_K (t_COL or t_INT), pr prime ideal, sprk = sprkinit(nf,pr,k,x),
    2675             :  * mod positive t_INT or NULL (meaning mod=0).
    2676             :  * return log(a) modulo mod on SNF-generators of (Z_K/pr^k)^* */
    2677             : GEN
    2678      436630 : log_prk(GEN nf, GEN a, GEN sprk, GEN mod)
    2679             : {
    2680             :   GEN e, prk, g, U1, U2, y, ff, O, o, oN, gN,  N, T, p, modpr, pr, cyc;
    2681             : 
    2682      436630 :   if (typ(a) == t_MAT) return famat_zlog_pr(nf, gel(a,1), gel(a,2), sprk, mod);
    2683      431086 :   N = NULL;
    2684      431086 :   ff = sprk_get_ff(sprk);
    2685      431086 :   pr = gel(ff,1); /* modpr */
    2686      431086 :   g = gN = gel(ff,2);
    2687      431086 :   O = gel(ff,3); /* order of g = |Fq^*|, in [ord, fa] format */
    2688      431086 :   o = oN = gel(O,1); /* order as a t_INT */
    2689      431086 :   prk = sprk_get_prk(sprk);
    2690      431086 :   modpr = nf_to_Fq_init(nf, &pr, &T, &p);
    2691      431086 :   if (mod)
    2692             :   {
    2693      339904 :     GEN d = gcdii(o,mod);
    2694      339904 :     if (!equalii(o, d))
    2695             :     {
    2696      168302 :       N = diviiexact(o,d); /* > 1, coprime to p */
    2697      168302 :       a = nfpowmodideal(nf, a, N, prk);
    2698      168302 :       oN = d; /* order of g^N mod pr */
    2699             :     }
    2700             :   }
    2701      431086 :   if (equali1(oN))
    2702      121954 :     e = gen_0;
    2703             :   else
    2704             :   {
    2705      309132 :     if (N) { O = order_update(oN, O); gN = Fq_pow(g, N, T, p); }
    2706      309132 :     e = Fq_log(nf_to_Fq(nf,a,modpr), gN, O, T, p);
    2707             :   }
    2708             :   /* 0 <= e < oN is correct modulo oN */
    2709      431086 :   if (sprk_is_prime(sprk)) return mkcol(e); /* k = 1 */
    2710             : 
    2711      140938 :   sprk_get_U2(sprk, &U1,&U2);
    2712      140938 :   cyc = sprk_get_cyc(sprk);
    2713      140938 :   if (mod)
    2714             :   {
    2715      118258 :     cyc = ZV_snf_gcd(cyc, mod);
    2716      118258 :     if (signe(remii(mod,p))) return vecmodii(ZC_Z_mul(U1,e), cyc);
    2717             :   }
    2718      140028 :   if (signe(e))
    2719             :   {
    2720       10500 :     GEN E = N? mulii(e, N): e;
    2721       10500 :     a = nfmulpowmodideal(nf, a, sprk_get_gnf(sprk), Fp_neg(E, o), prk);
    2722             :   }
    2723             :   /* a = 1 mod pr */
    2724      140028 :   y = log_prk1(nf, a, lg(U2)-1, sprk_get_L2(sprk), prk);
    2725      140021 :   if (N)
    2726             :   { /* from DL(a^N) to DL(a) */
    2727       45535 :     GEN E = gel(sprk_get_cyc(sprk), 1), q = powiu(p, Z_pval(E, p));
    2728       45535 :     y = ZC_Z_mul(y, Fp_inv(N, q));
    2729             :   }
    2730      140021 :   y = ZC_lincomb(gen_1, e, ZM_ZC_mul(U2,y), U1);
    2731      140021 :   return vecmodii(y, cyc);
    2732             : }
    2733             : GEN
    2734        2737 : log_prk_init(GEN nf, GEN pr, long k, GEN MOD)
    2735        2737 : { return sprkinit(checknf(nf),pr,k,NULL,MOD);}
    2736             : GEN
    2737         497 : veclog_prk(GEN nf, GEN v, GEN sprk)
    2738             : {
    2739         497 :   long l = lg(v), i;
    2740         497 :   GEN w = cgetg(l, t_MAT);
    2741        1232 :   for (i = 1; i < l; i++) gel(w,i) = log_prk(nf, gel(v,i), sprk, NULL);
    2742         497 :   return w;
    2743             : }
    2744             : 
    2745             : static GEN
    2746      292353 : famat_zlog(GEN nf, GEN fa, GEN sgn, zlog_S *S)
    2747             : {
    2748      292353 :   long i, l0, l = lg(S->U);
    2749      292353 :   GEN g = gel(fa,1), e = gel(fa,2), y = cgetg(l, t_COL);
    2750      292353 :   l0 = lg(S->sprk); /* = l (trivial arch. part), or l-1 */
    2751      635953 :   for (i=1; i < l0; i++) gel(y,i) = famat_zlog_pr(nf, g, e, gel(S->sprk,i), S->mod);
    2752      292353 :   if (l0 != l)
    2753             :   {
    2754      171323 :     if (!sgn) sgn = nfsign_arch(nf, fa, S->archp);
    2755      171323 :     gel(y,l0) = Flc_to_ZC(sgn);
    2756             :   }
    2757      292353 :   return y;
    2758             : }
    2759             : 
    2760             : /* assume that cyclic factors are normalized, in particular != [1] */
    2761             : static GEN
    2762       44716 : split_U(GEN U, GEN Sprk)
    2763             : {
    2764       44716 :   long t = 0, k, n, l = lg(Sprk);
    2765       44716 :   GEN vU = cgetg(l+1, t_VEC);
    2766       89012 :   for (k = 1; k < l; k++)
    2767             :   {
    2768       44296 :     n = lg(sprk_get_cyc(gel(Sprk,k))) - 1; /* > 0 */
    2769       44296 :     gel(vU,k) = vecslice(U, t+1, t+n);
    2770       44296 :     t += n;
    2771             :   }
    2772             :   /* t+1 .. lg(U)-1 */
    2773       44716 :   n = lg(U) - t - 1; /* can be 0 */
    2774       44716 :   if (!n) setlg(vU,l); else gel(vU,l) = vecslice(U, t+1, t+n);
    2775       44716 :   return vU;
    2776             : }
    2777             : 
    2778             : static void
    2779      437078 : init_zlog_mod(zlog_S *S, GEN bid, GEN mod)
    2780             : {
    2781      437078 :   GEN fa2 = bid_get_fact2(bid);
    2782      437078 :   S->U = bid_get_U(bid);
    2783      437078 :   S->hU = lg(bid_get_cyc(bid))-1;
    2784      437078 :   S->archp = bid_get_archp(bid);
    2785      437078 :   S->sprk = bid_get_sprk(bid);
    2786      437078 :   S->bid = bid;
    2787      437078 :   S->mod = mod;
    2788      437078 :   S->P = gel(fa2,1);
    2789      437078 :   S->k = gel(fa2,2);
    2790      437078 :   S->no2 = lg(S->P) == lg(gel(bid_get_fact(bid),1));
    2791      437078 : }
    2792             : void
    2793       99708 : init_zlog(zlog_S *S, GEN bid)
    2794             : {
    2795       99708 :   return init_zlog_mod(S, bid, NULL);
    2796             : }
    2797             : 
    2798             : /* a a t_FRAC/t_INT, reduce mod bid */
    2799             : static GEN
    2800          14 : Q_mod_bid(GEN bid, GEN a)
    2801             : {
    2802          14 :   GEN xZ = gcoeff(bid_get_ideal(bid),1,1);
    2803          14 :   GEN b = Rg_to_Fp(a, xZ);
    2804          14 :   if (gsigne(a) < 0) b = subii(b, xZ);
    2805          14 :   return signe(b)? b: xZ;
    2806             : }
    2807             : /* Return decomposition of a on the CRT generators blocks attached to the
    2808             :  * S->sprk and sarch; sgn = sign(a, S->arch), NULL if unknown */
    2809             : static GEN
    2810       66220 : zlog(GEN nf, GEN a, GEN sgn, zlog_S *S)
    2811             : {
    2812             :   long k, l;
    2813             :   GEN y;
    2814       66220 :   a = nf_to_scalar_or_basis(nf, a);
    2815       66220 :   switch(typ(a))
    2816             :   {
    2817       35441 :     case t_INT: break;
    2818          14 :     case t_FRAC: a = Q_mod_bid(S->bid, a); break;
    2819       30765 :     default: /* case t_COL: */
    2820             :     {
    2821             :       GEN den;
    2822       30765 :       check_nfelt(a, &den);
    2823       30765 :       if (den)
    2824             :       {
    2825          98 :         a = Q_muli_to_int(a, den);
    2826          98 :         a = mkmat2(mkcol2(a, den), mkcol2(gen_1, gen_m1));
    2827          98 :         return famat_zlog(nf, a, sgn, S);
    2828             :       }
    2829             :     }
    2830             :   }
    2831       66122 :   if (sgn)
    2832       60417 :     sgn = (lg(sgn) == 1)? NULL: leafcopy(sgn);
    2833             :   else
    2834        5705 :     sgn = (lg(S->archp) == 1)? NULL: nfsign_arch(nf, a, S->archp);
    2835       66122 :   l = lg(S->sprk);
    2836       66122 :   y = cgetg(sgn? l+1: l, t_COL);
    2837      144550 :   for (k = 1; k < l; k++)
    2838             :   {
    2839       78435 :     GEN sprk = gel(S->sprk,k);
    2840       78435 :     gel(y,k) = log_prk(nf, a, sprk, S->mod);
    2841             :   }
    2842       66115 :   if (sgn) gel(y,l) = Flc_to_ZC(sgn);
    2843       66115 :   return y;
    2844             : }
    2845             : 
    2846             : /* true nf */
    2847             : GEN
    2848       16345 : pr_basis_perm(GEN nf, GEN pr)
    2849             : {
    2850       16345 :   long f = pr_get_f(pr);
    2851             :   GEN perm;
    2852       16345 :   if (f == nf_get_degree(nf)) return identity_perm(f);
    2853       14721 :   perm = cgetg(f+1, t_VECSMALL);
    2854       14721 :   perm[1] = 1;
    2855       14721 :   if (f > 1)
    2856             :   {
    2857        2541 :     GEN H = pr_hnf(nf,pr);
    2858             :     long i, k;
    2859        9779 :     for (i = k = 2; k <= f; i++)
    2860        7238 :       if (!equali1(gcoeff(H,i,i))) perm[k++] = i;
    2861             :   }
    2862       14721 :   return perm;
    2863             : }
    2864             : 
    2865             : /* \sum U[i]*y[i], U[i] ZM, y[i] ZC. We allow lg(y) > lg(U). */
    2866             : static GEN
    2867      358489 : ZMV_ZCV_mul(GEN U, GEN y)
    2868             : {
    2869      358489 :   long i, l = lg(U);
    2870      358489 :   GEN z = NULL;
    2871      358489 :   if (l == 1) return cgetg(1,t_COL);
    2872      997382 :   for (i = 1; i < l; i++)
    2873             :   {
    2874      638893 :     GEN u = ZM_ZC_mul(gel(U,i), gel(y,i));
    2875      638893 :     z = z? ZC_add(z, u): u;
    2876             :   }
    2877      358489 :   return z;
    2878             : }
    2879             : /* A * (U[1], ..., U[d] */
    2880             : static GEN
    2881         518 : ZM_ZMV_mul(GEN A, GEN U)
    2882             : {
    2883             :   long i, l;
    2884         518 :   GEN V = cgetg_copy(U,&l);
    2885        1057 :   for (i = 1; i < l; i++) gel(V,i) = ZM_mul(A,gel(U,i));
    2886         518 :   return V;
    2887             : }
    2888             : 
    2889             : /* a = 1 mod pr, sprk mod pr^e, e >= 1 */
    2890             : static GEN
    2891       49637 : sprk_log_prk1_2(GEN nf, GEN a, GEN sprk)
    2892             : {
    2893       49637 :   GEN U1, U2, y, L2 = sprk_get_L2(sprk);
    2894       49637 :   sprk_get_U2(sprk, &U1,&U2);
    2895       49637 :   y = ZM_ZC_mul(U2, log_prk1(nf, a, lg(U2)-1, L2, sprk_get_prk(sprk)));
    2896       49637 :   return vecmodii(y, sprk_get_cyc(sprk));
    2897             : }
    2898             : /* assume e >= 2 */
    2899             : static GEN
    2900       34370 : sprk_log_gen_pr2(GEN nf, GEN sprk, long e)
    2901             : {
    2902       34370 :   GEN M, G, pr = sprk_get_pr(sprk);
    2903             :   long i, l;
    2904       34370 :   if (e == 2)
    2905             :   {
    2906       18277 :     GEN L2 = sprk_get_L2(sprk), L = gel(L2,1);
    2907       18277 :     G = gel(L,2); l = lg(G);
    2908             :   }
    2909             :   else
    2910             :   {
    2911       16093 :     GEN perm = pr_basis_perm(nf,pr), PI = nfpow_u(nf, pr_get_gen(pr), e-1);
    2912       16093 :     l = lg(perm);
    2913       16093 :     G = cgetg(l, t_VEC);
    2914       16093 :     if (typ(PI) == t_INT)
    2915             :     { /* zk_ei_mul doesn't allow t_INT */
    2916        1617 :       long N = nf_get_degree(nf);
    2917        1617 :       gel(G,1) = addiu(PI,1);
    2918        2135 :       for (i = 2; i < l; i++)
    2919             :       {
    2920         518 :         GEN z = col_ei(N, 1);
    2921         518 :         gel(G,i) = z; gel(z, perm[i]) = PI;
    2922             :       }
    2923             :     }
    2924             :     else
    2925             :     {
    2926       14476 :       gel(G,1) = nfadd(nf, gen_1, PI);
    2927       20832 :       for (i = 2; i < l; i++)
    2928        6356 :         gel(G,i) = nfadd(nf, gen_1, zk_ei_mul(nf, PI, perm[i]));
    2929             :     }
    2930             :   }
    2931       34370 :   M = cgetg(l, t_MAT);
    2932       78099 :   for (i = 1; i < l; i++) gel(M,i) = sprk_log_prk1_2(nf, gel(G,i), sprk);
    2933       34370 :   return M;
    2934             : }
    2935             : /* Log on bid.gen of generators of P_{1,I pr^{e-1}} / P_{1,I pr^e} (I,pr) = 1,
    2936             :  * defined implicitly via CRT. 'ind' is the index of pr in modulus
    2937             :  * factorization */
    2938             : GEN
    2939       95445 : log_gen_pr(zlog_S *S, long ind, GEN nf, long e)
    2940             : {
    2941       95445 :   GEN Uind = gel(S->U, ind);
    2942       95445 :   if (e == 1) retmkmat( gel(Uind,1) );
    2943       32333 :   return ZM_mul(Uind, sprk_log_gen_pr2(nf, gel(S->sprk,ind), e));
    2944             : }
    2945             : GEN
    2946        2037 : sprk_log_gen_pr(GEN nf, GEN sprk, long e)
    2947             : {
    2948        2037 :   if (e == 1)
    2949             :   {
    2950           0 :     long n = lg(sprk_get_cyc(sprk))-1;
    2951           0 :     retmkmat(col_ei(n, 1));
    2952             :   }
    2953        2037 :   return sprk_log_gen_pr2(nf, sprk, e);
    2954             : }
    2955             : /* a = 1 mod pr */
    2956             : GEN
    2957        5908 : sprk_log_prk1(GEN nf, GEN a, GEN sprk)
    2958             : {
    2959        5908 :   if (lg(sprk) == 5) return mkcol(gen_0); /* mod pr */
    2960        5908 :   return sprk_log_prk1_2(nf, a, sprk);
    2961             : }
    2962             : /* Log on bid.gen of generator of P_{1,f} / P_{1,f v[index]}
    2963             :  * v = vector of r1 real places */
    2964             : GEN
    2965       16562 : log_gen_arch(zlog_S *S, long index)
    2966             : {
    2967       16562 :   GEN U = gel(S->U, lg(S->U)-1);
    2968       16562 :   return gel(U, index);
    2969             : }
    2970             : 
    2971             : /* compute bid.clgp: [h,cyc] or [h,cyc,gen] */
    2972             : static GEN
    2973       45787 : bid_grp(GEN nf, GEN U, GEN cyc, GEN g, GEN F, GEN sarch)
    2974             : {
    2975       45787 :   GEN G, h = ZV_prod(cyc);
    2976             :   long c;
    2977       45787 :   if (!U) return mkvec2(h,cyc);
    2978       45507 :   c = lg(U);
    2979       45507 :   G = cgetg(c,t_VEC);
    2980       45507 :   if (c > 1)
    2981             :   {
    2982       38332 :     GEN U0, Uoo, EX = cyc_get_expo(cyc); /* exponent of bid */
    2983       38332 :     long i, hU = nbrows(U), nba = lg(sarch_get_cyc(sarch))-1; /* #f_oo */
    2984       38332 :     if (!nba) { U0 = U; Uoo = NULL; }
    2985       20671 :     else if (nba == hU) { U0 = NULL; Uoo = U; }
    2986             :     else
    2987             :     {
    2988       16884 :       U0 = rowslice(U, 1, hU-nba);
    2989       16884 :       Uoo = rowslice(U, hU-nba+1, hU);
    2990             :     }
    2991      110558 :     for (i = 1; i < c; i++)
    2992             :     {
    2993       72226 :       GEN t = gen_1;
    2994       72226 :       if (U0) t = famat_to_nf_modideal_coprime(nf, g, gel(U0,i), F, EX);
    2995       72226 :       if (Uoo) t = set_sign_mod_divisor(nf, ZV_to_Flv(gel(Uoo,i),2), t, sarch);
    2996       72226 :       gel(G,i) = t;
    2997             :     }
    2998             :   }
    2999       45507 :   return mkvec3(h, cyc, G);
    3000             : }
    3001             : 
    3002             : /* remove prime ideals of norm 2 with exponent 1 from factorization */
    3003             : static GEN
    3004       45486 : famat_strip2(GEN fa)
    3005             : {
    3006       45486 :   GEN P = gel(fa,1), E = gel(fa,2), Q, F;
    3007       45486 :   long l = lg(P), i, j;
    3008       45486 :   Q = cgetg(l, t_COL);
    3009       45486 :   F = cgetg(l, t_COL);
    3010       98630 :   for (i = j = 1; i < l; i++)
    3011             :   {
    3012       53144 :     GEN pr = gel(P,i), e = gel(E,i);
    3013       53144 :     if (!absequaliu(pr_get_p(pr), 2) || itou(e) != 1 || pr_get_f(pr) != 1)
    3014             :     {
    3015       45101 :       gel(Q,j) = pr;
    3016       45101 :       gel(F,j) = e; j++;
    3017             :     }
    3018             :   }
    3019       45486 :   setlg(Q,j);
    3020       45486 :   setlg(F,j); return mkmat2(Q,F);
    3021             : }
    3022             : 
    3023             : /* Compute [[ideal,arch], [h,[cyc],[gen]], idealfact, [liste], U]
    3024             :    flag may include nf_GEN | nf_INIT */
    3025             : static GEN
    3026       45507 : Idealstarmod_i(GEN nf, GEN ideal, long flag, GEN MOD)
    3027             : {
    3028             :   long i, k, nbp, R1;
    3029       45507 :   GEN y, cyc, U, u1 = NULL, fa, fa2, sprk, x, arch, archp, E, P, sarch, gen;
    3030             : 
    3031       45507 :   nf = checknf(nf);
    3032       45507 :   R1 = nf_get_r1(nf);
    3033       45507 :   if (typ(ideal) == t_VEC && lg(ideal) == 3)
    3034             :   {
    3035       22372 :     arch = gel(ideal,2);
    3036       22372 :     ideal= gel(ideal,1);
    3037       44737 :     switch(typ(arch))
    3038             :     {
    3039       21861 :       case t_VEC:
    3040       21861 :         if (lg(arch) != R1+1)
    3041           0 :           pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
    3042       21861 :         archp = vec01_to_indices(arch);
    3043       21861 :         break;
    3044         511 :       case t_VECSMALL:
    3045         511 :         archp = arch;
    3046         511 :         k = lg(archp)-1;
    3047         511 :         if (k && archp[k] > R1)
    3048           7 :           pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
    3049         504 :         arch = indices_to_vec01(archp, R1);
    3050         504 :         break;
    3051           0 :       default:
    3052           0 :         pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
    3053             :         return NULL;/*LCOV_EXCL_LINE*/
    3054             :     }
    3055             :   }
    3056             :   else
    3057             :   {
    3058       23135 :     arch = zerovec(R1);
    3059       23135 :     archp = cgetg(1, t_VECSMALL);
    3060             :   }
    3061       45500 :   if (MOD)
    3062             :   {
    3063        2506 :     if (typ(MOD) != t_INT) pari_err_TYPE("bnrinit [incorrect cycmod]", MOD);
    3064        2506 :     if (mpodd(MOD) && lg(archp) != 1)
    3065         161 :       MOD = shifti(MOD, 1); /* ensure elements of G^MOD are >> 0 */
    3066             :   }
    3067       45500 :   if (is_nf_factor(ideal))
    3068             :   {
    3069        2919 :     fa = ideal;
    3070        2919 :     x = factorbackprime(nf, gel(fa,1), gel(fa,2));
    3071             :   }
    3072             :   else
    3073             :   {
    3074       42581 :     fa = idealfactor(nf, ideal);
    3075       42574 :     x = ideal;
    3076             :   }
    3077       45493 :   if (typ(x) != t_MAT) x = idealhnf_shallow(nf, x);
    3078       45493 :   if (lg(x) == 1) pari_err_DOMAIN("Idealstar", "ideal","=",gen_0,x);
    3079       45493 :   if (typ(gcoeff(x,1,1)) != t_INT)
    3080           7 :     pari_err_DOMAIN("Idealstar","denominator(ideal)", "!=",gen_1,x);
    3081       45486 :   sarch = nfarchstar(nf, x, archp);
    3082       45486 :   fa2 = famat_strip2(fa);
    3083       45486 :   P = gel(fa2,1);
    3084       45486 :   E = gel(fa2,2);
    3085       45486 :   nbp = lg(P)-1;
    3086       45486 :   sprk = cgetg(nbp+1,t_VEC);
    3087       45486 :   if (nbp)
    3088             :   {
    3089       34482 :     GEN t = (lg(gel(fa,1))==2)? NULL: x; /* beware fa != fa2 */
    3090       34482 :     cyc = cgetg(nbp+2,t_VEC);
    3091       34482 :     gen = cgetg(nbp+1,t_VEC);
    3092       79583 :     for (i = 1; i <= nbp; i++)
    3093             :     {
    3094       45101 :       GEN L = sprkinit(nf, gel(P,i), itou(gel(E,i)), t, MOD);
    3095       45101 :       gel(sprk,i) = L;
    3096       45101 :       gel(cyc,i) = sprk_get_cyc(L);
    3097             :       /* true gens are congruent to those mod x AND positive at archp */
    3098       45101 :       gel(gen,i) = sprk_get_gen(L);
    3099             :     }
    3100       34482 :     gel(cyc,i) = sarch_get_cyc(sarch);
    3101       34482 :     cyc = shallowconcat1(cyc);
    3102       34482 :     gen = shallowconcat1(gen);
    3103       34482 :     cyc = ZV_snf_group(cyc, &U, (flag & nf_GEN)? &u1: NULL);
    3104             :   }
    3105             :   else
    3106             :   {
    3107       11004 :     cyc = sarch_get_cyc(sarch);
    3108       11004 :     gen = cgetg(1,t_VEC);
    3109       11004 :     U = matid(lg(cyc)-1);
    3110       11004 :     if (flag & nf_GEN) u1 = U;
    3111             :   }
    3112       45486 :   y = bid_grp(nf, u1, cyc, gen, x, sarch);
    3113       45486 :   if (!(flag & nf_INIT)) return y;
    3114       44674 :   U = split_U(U, sprk);
    3115       44674 :   return mkvec5(mkvec2(x, arch), y, mkvec2(fa,fa2), mkvec2(sprk, sarch), U);
    3116             : }
    3117             : GEN
    3118       45234 : Idealstarmod(GEN nf, GEN ideal, long flag, GEN MOD)
    3119             : {
    3120       45234 :   pari_sp av = avma;
    3121       45234 :   if (!nf) nf = nfinit(pol_x(0), DEFAULTPREC);
    3122       45234 :   return gerepilecopy(av, Idealstarmod_i(nf, ideal, flag, MOD));
    3123             : }
    3124             : GEN
    3125         952 : Idealstar(GEN nf, GEN ideal, long flag) { return Idealstarmod(nf, ideal, flag, NULL); }
    3126             : GEN
    3127         273 : Idealstarprk(GEN nf, GEN pr, long k, long flag)
    3128             : {
    3129         273 :   pari_sp av = avma;
    3130         273 :   GEN z = Idealstarmod_i(nf, mkmat2(mkcol(pr),mkcols(k)), flag, NULL);
    3131         273 :   return gerepilecopy(av, z);
    3132             : }
    3133             : 
    3134             : /* FIXME: obsolete */
    3135             : GEN
    3136           0 : zidealstarinitgen(GEN nf, GEN ideal)
    3137           0 : { return Idealstar(nf,ideal, nf_INIT|nf_GEN); }
    3138             : GEN
    3139           0 : zidealstarinit(GEN nf, GEN ideal)
    3140           0 : { return Idealstar(nf,ideal, nf_INIT); }
    3141             : GEN
    3142           0 : zidealstar(GEN nf, GEN ideal)
    3143           0 : { return Idealstar(nf,ideal, nf_GEN); }
    3144             : 
    3145             : GEN
    3146          70 : idealstarmod(GEN nf, GEN ideal, long flag, GEN MOD)
    3147             : {
    3148          70 :   switch(flag)
    3149             :   {
    3150           0 :     case 0: return Idealstarmod(nf,ideal, nf_GEN, MOD);
    3151          56 :     case 1: return Idealstarmod(nf,ideal, nf_INIT, MOD);
    3152          14 :     case 2: return Idealstarmod(nf,ideal, nf_INIT|nf_GEN, MOD);
    3153           0 :     default: pari_err_FLAG("idealstar");
    3154             :   }
    3155             :   return NULL; /* LCOV_EXCL_LINE */
    3156             : }
    3157             : GEN
    3158           0 : idealstar0(GEN nf, GEN ideal,long flag) { return idealstarmod(nf, ideal, flag, NULL); }
    3159             : 
    3160             : void
    3161       30765 : check_nfelt(GEN x, GEN *den)
    3162             : {
    3163       30765 :   long l = lg(x), i;
    3164       30765 :   GEN t, d = NULL;
    3165       30765 :   if (typ(x) != t_COL) pari_err_TYPE("check_nfelt", x);
    3166      108311 :   for (i=1; i<l; i++)
    3167             :   {
    3168       77546 :     t = gel(x,i);
    3169       77546 :     switch (typ(t))
    3170             :     {
    3171       77343 :       case t_INT: break;
    3172         203 :       case t_FRAC:
    3173         203 :         if (!d) d = gel(t,2); else d = lcmii(d, gel(t,2));
    3174         203 :         break;
    3175           0 :       default: pari_err_TYPE("check_nfelt", x);
    3176             :     }
    3177             :   }
    3178       30765 :   *den = d;
    3179       30765 : }
    3180             : 
    3181             : GEN
    3182     1719296 : vecmodii(GEN x, GEN y)
    3183     4845216 : { pari_APPLY_same(modii(gel(x,i), gel(y,i))) }
    3184             : GEN
    3185      150032 : ZV_snf_gcd(GEN x, GEN mod)
    3186      441886 : { pari_APPLY_same(gcdii(gel(x,i), mod)); }
    3187             : 
    3188             : GEN
    3189       27230 : vecmoduu(GEN x, GEN y)
    3190       88158 : { pari_APPLY_ulong(uel(x,i) % uel(y,i)) }
    3191             : 
    3192             : /* assume a true bnf and bid */
    3193             : GEN
    3194       37282 : ideallog_units0(GEN bnf, GEN bid, GEN MOD)
    3195             : {
    3196       37282 :   GEN nf = bnf_get_nf(bnf), D, y, C, cyc;
    3197       37282 :   long j, lU = lg(bnf_get_logfu(bnf)); /* r1+r2 */
    3198             :   zlog_S S;
    3199       37282 :   init_zlog_mod(&S, bid, MOD);
    3200       37282 :   if (!S.hU) return zeromat(0,lU);
    3201       37282 :   cyc = bid_get_cyc(bid);
    3202       37282 :   if (MOD) cyc = ZV_snf_gcd(cyc, MOD);
    3203       37282 :   D = nfsign_fu(bnf, bid_get_archp(bid));
    3204       37282 :   y = cgetg(lU, t_MAT);
    3205       37282 :   if ((C = bnf_build_cheapfu(bnf)))
    3206       60396 :   { for (j = 1; j < lU; j++) gel(y,j) = zlog(nf, gel(C,j), gel(D,j), &S); }
    3207             :   else
    3208             :   {
    3209          21 :     long i, l = lg(S.U), l0 = lg(S.sprk);
    3210             :     GEN X, U;
    3211          21 :     if (!(C = bnf_compactfu_mat(bnf))) bnf_build_units(bnf); /* error */
    3212          21 :     X = gel(C,1); U = gel(C,2);
    3213          42 :     for (j = 1; j < lU; j++) gel(y,j) = cgetg(l, t_COL);
    3214          42 :     for (i = 1; i < l0; i++)
    3215             :     {
    3216          21 :       GEN sprk = gel(S.sprk, i);
    3217          21 :       GEN Xi = sunits_makecoprime(X, sprk_get_pr(sprk), sprk_get_prk(sprk));
    3218          42 :       for (j = 1; j < lU; j++)
    3219          21 :         gcoeff(y,i,j) = famat_zlog_pr_coprime(nf, Xi, gel(U,j), sprk, MOD);
    3220             :     }
    3221          21 :     if (l0 != l)
    3222          14 :       for (j = 1; j < lU; j++) gcoeff(y,l0,j) = Flc_to_ZC(gel(D,j));
    3223             :   }
    3224       37282 :   y = vec_prepend(y, zlog(nf, bnf_get_tuU(bnf), nfsign_tu(bnf, S.archp), &S));
    3225       97720 :   for (j = 1; j <= lU; j++)
    3226       60438 :     gel(y,j) = vecmodii(ZMV_ZCV_mul(S.U, gel(y,j)), cyc);
    3227       37282 :   return y;
    3228             : }
    3229             : GEN
    3230          84 : ideallog_units(GEN bnf, GEN bid)
    3231          84 : { return ideallog_units0(bnf, bid, NULL); }
    3232             : GEN
    3233         518 : log_prk_units(GEN nf, GEN D, GEN sprk)
    3234             : {
    3235         518 :   GEN L, Ltu = log_prk(nf, gel(D,1), sprk, NULL);
    3236         518 :   D = gel(D,2);
    3237         518 :   if (lg(D) != 3 || typ(gel(D,2)) != t_MAT) L = veclog_prk(nf, D, sprk);
    3238             :   else
    3239             :   {
    3240          21 :     GEN X = gel(D,1), U = gel(D,2);
    3241          21 :     long j, lU = lg(U);
    3242          21 :     X = sunits_makecoprime(X, sprk_get_pr(sprk), sprk_get_prk(sprk));
    3243          21 :     L = cgetg(lU, t_MAT);
    3244          63 :     for (j = 1; j < lU; j++)
    3245          42 :       gel(L,j) = famat_zlog_pr_coprime(nf, X, gel(U,j), sprk, NULL);
    3246             :   }
    3247         518 :   return vec_prepend(L, Ltu);
    3248             : }
    3249             : 
    3250             : static GEN
    3251      300088 : ideallog_i(GEN nf, GEN x, zlog_S *S)
    3252             : {
    3253      300088 :   pari_sp av = avma;
    3254             :   GEN y;
    3255      300088 :   if (!S->hU) return cgetg(1, t_COL);
    3256      298058 :   if (typ(x) == t_MAT)
    3257      292255 :     y = famat_zlog(nf, x, NULL, S);
    3258             :   else
    3259        5803 :     y = zlog(nf, x, NULL, S);
    3260      298051 :   y = ZMV_ZCV_mul(S->U, y);
    3261      298051 :   return gerepileupto(av, vecmodii(y, bid_get_cyc(S->bid)));
    3262             : }
    3263             : GEN
    3264      306766 : ideallogmod(GEN nf, GEN x, GEN bid, GEN mod)
    3265             : {
    3266             :   zlog_S S;
    3267      306766 :   if (!nf)
    3268             :   {
    3269        6671 :     if (mod) pari_err_IMPL("Zideallogmod");
    3270        6671 :     return Zideallog(bid, x);
    3271             :   }
    3272      300095 :   checkbid(bid); init_zlog_mod(&S, bid, mod);
    3273      300088 :   return ideallog_i(checknf(nf), x, &S);
    3274             : }
    3275             : GEN
    3276       13258 : ideallog(GEN nf, GEN x, GEN bid) { return ideallogmod(nf, x, bid, NULL); }
    3277             : 
    3278             : /*************************************************************************/
    3279             : /**                                                                     **/
    3280             : /**               JOIN BID STRUCTURES, IDEAL LISTS                      **/
    3281             : /**                                                                     **/
    3282             : /*************************************************************************/
    3283             : /* bid1, bid2: for coprime modules m1 and m2 (without arch. part).
    3284             :  * Output: bid for m1 m2 */
    3285             : static GEN
    3286         469 : join_bid(GEN nf, GEN bid1, GEN bid2)
    3287             : {
    3288         469 :   pari_sp av = avma;
    3289             :   long nbgen, l1,l2;
    3290             :   GEN I1,I2, G1,G2, sprk1,sprk2, cyc1,cyc2, sarch;
    3291         469 :   GEN sprk, fa,fa2, U, cyc, y, u1 = NULL, x, gen;
    3292             : 
    3293         469 :   I1 = bid_get_ideal(bid1);
    3294         469 :   I2 = bid_get_ideal(bid2);
    3295         469 :   if (gequal1(gcoeff(I1,1,1))) return bid2; /* frequent trivial case */
    3296         259 :   G1 = bid_get_grp(bid1);
    3297         259 :   G2 = bid_get_grp(bid2);
    3298         259 :   x = idealmul(nf, I1,I2);
    3299         259 :   fa = famat_mul_shallow(bid_get_fact(bid1), bid_get_fact(bid2));
    3300         259 :   fa2= famat_mul_shallow(bid_get_fact2(bid1), bid_get_fact2(bid2));
    3301         259 :   sprk1 = bid_get_sprk(bid1);
    3302         259 :   sprk2 = bid_get_sprk(bid2);
    3303         259 :   sprk = shallowconcat(sprk1, sprk2);
    3304             : 
    3305         259 :   cyc1 = abgrp_get_cyc(G1); l1 = lg(cyc1);
    3306         259 :   cyc2 = abgrp_get_cyc(G2); l2 = lg(cyc2);
    3307         259 :   gen = (lg(G1)>3 && lg(G2)>3)? gen_1: NULL;
    3308         259 :   nbgen = l1+l2-2;
    3309         259 :   cyc = ZV_snf_group(shallowconcat(cyc1,cyc2), &U, gen? &u1: NULL);
    3310         259 :   if (nbgen)
    3311             :   {
    3312         259 :     GEN U1 = bid_get_U(bid1), U2 = bid_get_U(bid2);
    3313           0 :     U1 = l1==1? const_vec(lg(sprk1), cgetg(1,t_MAT))
    3314         259 :               : ZM_ZMV_mul(vecslice(U, 1, l1-1),   U1);
    3315           0 :     U2 = l2==1? const_vec(lg(sprk2), cgetg(1,t_MAT))
    3316         259 :               : ZM_ZMV_mul(vecslice(U, l1, nbgen), U2);
    3317         259 :     U = shallowconcat(U1, U2);
    3318             :   }
    3319             :   else
    3320           0 :     U = const_vec(lg(sprk), cgetg(1,t_MAT));
    3321             : 
    3322         259 :   if (gen)
    3323             :   {
    3324         259 :     GEN uv = zkchinese1init2(nf, I2, I1, x);
    3325         518 :     gen = shallowconcat(zkVchinese1(gel(uv,1), abgrp_get_gen(G1)),
    3326         259 :                         zkVchinese1(gel(uv,2), abgrp_get_gen(G2)));
    3327             :   }
    3328         259 :   sarch = bid_get_sarch(bid1); /* trivial */
    3329         259 :   y = bid_grp(nf, u1, cyc, gen, x, sarch);
    3330         259 :   x = mkvec2(x, bid_get_arch(bid1));
    3331         259 :   y = mkvec5(x, y, mkvec2(fa, fa2), mkvec2(sprk, sarch), U);
    3332         259 :   return gerepilecopy(av,y);
    3333             : }
    3334             : 
    3335             : typedef struct _ideal_data {
    3336             :   GEN nf, emb, L, pr, prL, sgnU, archp;
    3337             : } ideal_data;
    3338             : 
    3339             : /* z <- ( z | f(v[i])_{i=1..#v} ) */
    3340             : static void
    3341      119252 : concat_join(GEN *pz, GEN v, GEN (*f)(ideal_data*,GEN), ideal_data *data)
    3342             : {
    3343      119252 :   long i, nz, lv = lg(v);
    3344             :   GEN z, Z;
    3345      119252 :   if (lv == 1) return;
    3346       52542 :   z = *pz; nz = lg(z)-1;
    3347       52542 :   *pz = Z = cgetg(lv + nz, typ(z));
    3348       86870 :   for (i = 1; i <=nz; i++) gel(Z,i) = gel(z,i);
    3349       52542 :   Z += nz;
    3350      113050 :   for (i = 1; i < lv; i++) gel(Z,i) = f(data, gel(v,i));
    3351             : }
    3352             : static GEN
    3353         469 : join_idealinit(ideal_data *D, GEN x)
    3354         469 : { return join_bid(D->nf, x, D->prL); }
    3355             : static GEN
    3356       60039 : join_ideal(ideal_data *D, GEN x)
    3357       60039 : { return idealmulpowprime(D->nf, x, D->pr, D->L); }
    3358             : static GEN
    3359         448 : join_unit(ideal_data *D, GEN x)
    3360             : {
    3361         448 :   GEN bid = join_idealinit(D, gel(x,1)), u = gel(x,2), v = mkvec(D->emb);
    3362         448 :   if (lg(u) != 1) v = shallowconcat(u, v);
    3363         448 :   return mkvec2(bid, v);
    3364             : }
    3365             : 
    3366             : GEN
    3367          49 : log_prk_units_init(GEN bnf)
    3368             : {
    3369          49 :   GEN U = bnf_has_fu(bnf);
    3370          49 :   if (U) U = matalgtobasis(bnf_get_nf(bnf), U);
    3371          21 :   else if (!(U = bnf_compactfu_mat(bnf))) (void)bnf_build_units(bnf);
    3372          49 :   return mkvec2(bnf_get_tuU(bnf), U);
    3373             : }
    3374             : /*  flag & nf_GEN : generators, otherwise no
    3375             :  *  flag &2 : units, otherwise no
    3376             :  *  flag &4 : ideals in HNF, otherwise bid
    3377             :  *  flag &8 : omit ideals which cannot be conductors (pr^1 with Npr=2) */
    3378             : static GEN
    3379        6041 : Ideallist(GEN bnf, ulong bound, long flag)
    3380             : {
    3381        6041 :   const long cond = flag & 8;
    3382        6041 :   const long do_units = flag & 2, big_id = !(flag & 4);
    3383        6041 :   const long istar_flag = (flag & nf_GEN) | nf_INIT;
    3384        6041 :   pari_sp av, av0 = avma;
    3385             :   long i, j;
    3386        6041 :   GEN nf, z, p, fa, id, BOUND, U, empty = cgetg(1,t_VEC);
    3387             :   forprime_t S;
    3388             :   ideal_data ID;
    3389        6041 :   GEN (*join_z)(ideal_data*, GEN) =
    3390             :     do_units? &join_unit
    3391        6041 :               : (big_id? &join_idealinit: &join_ideal);
    3392             : 
    3393        6041 :   if (do_units)
    3394             :   {
    3395          21 :     bnf = checkbnf(bnf);
    3396          21 :     nf = bnf_get_nf(bnf);
    3397             :   }
    3398             :   else
    3399        6020 :     nf = checknf(bnf);
    3400        6041 :   if ((long)bound <= 0) return empty;
    3401        6041 :   id = matid(nf_get_degree(nf));
    3402        6041 :   if (big_id) id = Idealstar(nf,id, istar_flag);
    3403             : 
    3404             :   /* z[i] will contain all "objects" of norm i. Depending on flag, this means
    3405             :    * an ideal, a bid, or a couple [bid, log(units)]. Such objects are stored
    3406             :    * in vectors, computed one primary component at a time; join_z
    3407             :    * reconstructs the global object */
    3408        6041 :   BOUND = utoipos(bound);
    3409        6041 :   z = const_vec(bound, empty);
    3410        6041 :   U = do_units? log_prk_units_init(bnf): NULL;
    3411        6041 :   gel(z,1) = mkvec(U? mkvec2(id, cgetg(1,t_VEC)): id);
    3412        6041 :   ID.nf = nf;
    3413             : 
    3414        6041 :   p = cgetipos(3);
    3415        6041 :   u_forprime_init(&S, 2, bound);
    3416        6041 :   av = avma;
    3417       27454 :   while ((p[2] = u_forprime_next(&S)))
    3418             :   {
    3419       21413 :     if (DEBUGLEVEL>1) err_printf("%ld ",p[2]);
    3420       21413 :     fa = idealprimedec_limit_norm(nf, p, BOUND);
    3421       43652 :     for (j = 1; j < lg(fa); j++)
    3422             :     {
    3423       22239 :       GEN pr = gel(fa,j), z2 = leafcopy(z);
    3424       22239 :       ulong Q, q = upr_norm(pr);
    3425             :       long l;
    3426       22239 :       ID.pr = ID.prL = pr;
    3427       22239 :       if (cond && q == 2) { l = 2; Q = 4; } else { l = 1; Q = q; }
    3428       53795 :       for (; Q <= bound; l++, Q *= q) /* add pr^l */
    3429             :       {
    3430             :         ulong iQ;
    3431       31556 :         ID.L = utoipos(l);
    3432       31556 :         if (big_id) {
    3433         210 :           ID.prL = Idealstarprk(nf, pr, l, istar_flag);
    3434         210 :           if (U)
    3435         189 :             ID.emb = Q == 2? cgetg(1,t_VEC)
    3436         189 :                            : log_prk_units(nf, U, gel(bid_get_sprk(ID.prL),1));
    3437             :         }
    3438      150808 :         for (iQ = Q,i = 1; iQ <= bound; iQ += Q,i++)
    3439      119252 :           concat_join(&gel(z,iQ), gel(z2,i), join_z, &ID);
    3440             :       }
    3441             :     }
    3442       21413 :     if (gc_needed(av,1))
    3443             :     {
    3444           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Ideallist");
    3445           0 :       z = gerepilecopy(av, z);
    3446             :     }
    3447             :   }
    3448        6041 :   return gerepilecopy(av0, z);
    3449             : }
    3450             : GEN
    3451         357 : ideallist0(GEN bnf,long bound, long flag) {
    3452         357 :   if (flag<0 || flag>15) pari_err_FLAG("ideallist");
    3453         357 :   return Ideallist(bnf,bound,flag);
    3454             : }
    3455             : GEN
    3456        5684 : ideallist(GEN bnf,long bound) { return Ideallist(bnf,bound,4); }
    3457             : 
    3458             : /* bid = for module m (without arch. part), arch = archimedean part.
    3459             :  * Output: bid for [m,arch] */
    3460             : static GEN
    3461          42 : join_bid_arch(GEN nf, GEN bid, GEN archp)
    3462             : {
    3463          42 :   pari_sp av = avma;
    3464             :   GEN G, U;
    3465          42 :   GEN sprk, cyc, y, u1 = NULL, x, sarch, gen;
    3466             : 
    3467          42 :   checkbid(bid);
    3468          42 :   G = bid_get_grp(bid);
    3469          42 :   x = bid_get_ideal(bid);
    3470          42 :   sarch = nfarchstar(nf, bid_get_ideal(bid), archp);
    3471          42 :   sprk = bid_get_sprk(bid);
    3472             : 
    3473          42 :   gen = (lg(G)>3)? gel(G,3): NULL;
    3474          42 :   cyc = diagonal_shallow(shallowconcat(gel(G,2), sarch_get_cyc(sarch)));
    3475          42 :   cyc = ZM_snf_group(cyc, &U, gen? &u1: NULL);
    3476          42 :   y = bid_grp(nf, u1, cyc, gen, x, sarch);
    3477          42 :   U = split_U(U, sprk);
    3478          42 :   y = mkvec5(mkvec2(x, archp), y, gel(bid,3), mkvec2(sprk, sarch), U);
    3479          42 :   return gerepilecopy(av,y);
    3480             : }
    3481             : static GEN
    3482          42 : join_arch(ideal_data *D, GEN x) {
    3483          42 :   return join_bid_arch(D->nf, x, D->archp);
    3484             : }
    3485             : static GEN
    3486          14 : join_archunit(ideal_data *D, GEN x) {
    3487          14 :   GEN bid = join_arch(D, gel(x,1)), u = gel(x,2), v = mkvec(D->emb);
    3488          14 :   if (lg(u) != 1) v = shallowconcat(u, v);
    3489          14 :   return mkvec2(bid, v);
    3490             : }
    3491             : 
    3492             : /* L from ideallist, add archimedean part */
    3493             : GEN
    3494          14 : ideallistarch(GEN bnf, GEN L, GEN arch)
    3495             : {
    3496             :   pari_sp av;
    3497          14 :   long i, j, l = lg(L), lz;
    3498             :   GEN v, z, V;
    3499             :   ideal_data ID;
    3500             :   GEN (*join_z)(ideal_data*, GEN);
    3501             : 
    3502          14 :   if (typ(L) != t_VEC) pari_err_TYPE("ideallistarch",L);
    3503          14 :   if (l == 1) return cgetg(1,t_VEC);
    3504          14 :   z = gel(L,1);
    3505          14 :   if (typ(z) != t_VEC) pari_err_TYPE("ideallistarch",z);
    3506          14 :   z = gel(z,1); /* either a bid or [bid,U] */
    3507          14 :   ID.nf = checknf(bnf);
    3508          14 :   ID.archp = vec01_to_indices(arch);
    3509          14 :   if (lg(z) == 3) { /* the latter: do units */
    3510           7 :     if (typ(z) != t_VEC) pari_err_TYPE("ideallistarch",z);
    3511           7 :     ID.emb = zm_to_ZM( rowpermute(nfsign_units(bnf,NULL,1), ID.archp) );
    3512           7 :     join_z = &join_archunit;
    3513             :   } else
    3514           7 :     join_z = &join_arch;
    3515          14 :   av = avma; V = cgetg(l, t_VEC);
    3516          63 :   for (i = 1; i < l; i++)
    3517             :   {
    3518          49 :     z = gel(L,i); lz = lg(z);
    3519          49 :     gel(V,i) = v = cgetg(lz,t_VEC);
    3520          91 :     for (j=1; j<lz; j++) gel(v,j) = join_z(&ID, gel(z,j));
    3521             :   }
    3522          14 :   return gerepilecopy(av,V);
    3523             : }

Generated by: LCOV version 1.13