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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - modules - stark.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16624-25b9976) Lines: 1789 1988 90.0 %
Date: 2014-06-24 Functions: 115 122 94.3 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 880 1195 73.6 %

           Branch data     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                 :            : /*        COMPUTATION OF STARK UNITS OF TOTALLY REAL FIELDS        */
      17                 :            : /*                                                                 */
      18                 :            : /*******************************************************************/
      19                 :            : #include "pari.h"
      20                 :            : #include "paripriv.h"
      21                 :            : 
      22                 :            : #define EXTRA_PREC (DEFAULTPREC-1)
      23                 :            : #define ADD_PREC   (DEFAULTPREC-2)*3
      24                 :            : 
      25                 :            : /* ComputeCoeff */
      26                 :            : typedef struct {
      27                 :            :   GEN L0, L1, L11, L2; /* VECSMALL of p */
      28                 :            :   GEN L1ray, L11ray; /* precomputed isprincipalray(pr), pr | p */
      29                 :            :   GEN rayZ; /* precomputed isprincipalray(i), i < condZ */
      30                 :            :   long condZ; /* generates cond(bnr) \cap Z, assumed small */
      31                 :            : } LISTray;
      32                 :            : 
      33                 :            : /* Char evaluation */
      34                 :            : typedef struct {
      35                 :            :   long ord;
      36                 :            :   GEN *val, chi;
      37                 :            : } CHI_t;
      38                 :            : 
      39                 :            : /* RecCoeff */
      40                 :            : typedef struct {
      41                 :            :   GEN M, beta, B, U, nB;
      42                 :            :   long v, G, N;
      43                 :            : } RC_data;
      44                 :            : 
      45                 :            : /********************************************************************/
      46                 :            : /*                    Miscellaneous functions                       */
      47                 :            : /********************************************************************/
      48                 :            : /* exp(2iPi/den), assume den a t_INT */
      49                 :            : static GEN
      50                 :       2260 : InitRU(GEN den, long prec)
      51                 :            : {
      52                 :            :   GEN c, s;
      53         [ +  + ]:       2260 :   if (equaliu(den, 2)) return gen_m1;
      54                 :       1565 :   gsincos(divri(Pi2n(1, prec), den), &s, &c, prec);
      55                 :       2260 :   return mkcomplex(c, s);
      56                 :            : }
      57                 :            : /* Compute the image of logelt by character chi, as a complex number */
      58                 :            : static GEN
      59                 :      10285 : ComputeImagebyChar(GEN chi, GEN logelt)
      60                 :            : {
      61                 :      10285 :   GEN gn = ZV_dotproduct(gel(chi,1), logelt), x = gel(chi,2);
      62                 :      10285 :   long d = itos(gel(chi,3)), n = smodis(gn, d);
      63                 :            :   /* x^d = 1 and, if d even, x^(d/2) = -1 */
      64         [ +  - ]:      10285 :   if ((d & 1) == 0)
      65                 :            :   {
      66                 :      10285 :     d /= 2;
      67         [ +  + ]:      10285 :     if (n >= d) return gneg(gpowgs(x, n-d));
      68                 :            :   }
      69                 :      10285 :   return gpowgs(x, n);
      70                 :            : }
      71                 :            : 
      72                 :            : /* return n such that C(elt) = z^n */
      73                 :            : static ulong
      74                 :     594266 : EvalChar_n(CHI_t *C, GEN logelt)
      75                 :            : {
      76                 :     594266 :   GEN n = ZV_dotproduct(C->chi, logelt);
      77                 :     594266 :   return umodiu(n, C->ord);
      78                 :            : }
      79                 :            : /* return C(elt) */
      80                 :            : static GEN
      81                 :     593146 : EvalChar(CHI_t *C, GEN logelt)
      82                 :            : {
      83                 :     593146 :   return C->val[EvalChar_n(C, logelt)];
      84                 :            : }
      85                 :            : 
      86                 :            : static void
      87                 :       2665 : init_CHI(CHI_t *c, GEN CHI, GEN z)
      88                 :            : {
      89                 :       2665 :   long i, d = itos(gel(CHI,3));
      90                 :       2665 :   GEN *v = (GEN*)new_chunk(d);
      91                 :       2665 :   v[0] = gen_1;
      92                 :       2665 :   v[1] = z;
      93         [ +  + ]:      22200 :   for (i=2; i<d; i++) v[i] = gmul(v[i-1], z);
      94                 :       2665 :   c->chi = gel(CHI,1);
      95                 :       2665 :   c->ord = d;
      96                 :       2665 :   c->val = v;
      97                 :       2665 : }
      98                 :            : /* as t_POLMOD */
      99                 :            : static void
     100                 :       1655 : init_CHI_alg(CHI_t *c, GEN CHI) {
     101                 :       1655 :   long d = itos(gel(CHI,3));
     102                 :            :   GEN z;
     103      [ -  +  + ]:       1655 :   switch(d)
     104                 :            :   {
     105                 :          0 :     case 1: z = gen_1; break;
     106                 :        640 :     case 2: z = gen_m1; break;
     107                 :       1015 :     default: z = mkpolmod(pol_x(0), polcyclo(d,0));
     108                 :            :   }
     109                 :       1655 :   init_CHI(c,CHI, z);
     110                 :       1655 : }
     111                 :            : /* as t_COMPLEX */
     112                 :            : static void
     113                 :       1010 : init_CHI_C(CHI_t *c, GEN CHI) {
     114                 :       1010 :   init_CHI(c,CHI, gel(CHI,2));
     115                 :       1010 : }
     116                 :            : 
     117                 :            : /* Compute the conjugate character [ZV] */
     118                 :            : static GEN
     119                 :        620 : ConjChar(GEN chi, GEN cyc)
     120                 :            : {
     121                 :        620 :   long i, l = lg(chi);
     122                 :        620 :   GEN z = cgetg(l, t_VEC);
     123         [ +  + ]:       1455 :   for (i = 1; i < l; i++)
     124         [ +  + ]:        835 :     gel(z,i) = signe(gel(chi,i))? subii(gel(cyc,i), gel(chi,i)): gen_0;
     125                 :        620 :   return z;
     126                 :            : }
     127                 :            : 
     128                 :            : typedef struct {
     129                 :            :   long r; /* rank = lg(gen) */
     130                 :            :   GEN j; /* current elt is gen[1]^j[1] ... gen[r]^j[r] */
     131                 :            :   GEN cyc; /* t_VECSMALL of elementary divisors */
     132                 :            : } GROUP_t;
     133                 :            : 
     134                 :            : static int
     135                 :      37840 : NextElt(GROUP_t *G)
     136                 :            : {
     137                 :      37840 :   long i = 1;
     138         [ +  + ]:      37840 :   if (G->r == 0) return 0; /* no more elt */
     139         [ +  + ]:      38090 :   while (++G->j[i] == G->cyc[i]) /* from 0 to cyc[i]-1 */
     140                 :            :   {
     141                 :        750 :     G->j[i] = 0;
     142         [ +  + ]:        750 :     if (++i > G->r) return 0; /* no more elt */
     143                 :            :   }
     144                 :      37840 :   return i; /* we have multiplied by gen[i] */
     145                 :            : }
     146                 :            : 
     147                 :            : /* Compute all the elements of a group given by its SNF */
     148                 :            : static GEN
     149                 :        710 : EltsOfGroup(long order, GEN cyc)
     150                 :            : {
     151                 :            :   long i;
     152                 :            :   GEN rep;
     153                 :            :   GROUP_t G;
     154                 :            : 
     155                 :        710 :   G.cyc = gtovecsmall(cyc);
     156                 :        710 :   G.r = lg(cyc)-1;
     157                 :        710 :   G.j = zero_zv(G.r);
     158                 :            : 
     159                 :        710 :   rep = cgetg(order + 1, t_VEC);
     160                 :        710 :   gel(rep,order) = vecsmall_to_col(G.j);
     161                 :            : 
     162         [ +  + ]:       5025 :   for  (i = 1; i < order; i++)
     163                 :            :   {
     164                 :       4315 :     (void)NextElt(&G);
     165                 :       4315 :     gel(rep,i) = vecsmall_to_col(G.j);
     166                 :            :   }
     167                 :        710 :   return rep;
     168                 :            : }
     169                 :            : 
     170                 :            : /* Let dataC as given by InitQuotient, compute a system of
     171                 :            :    representatives of the quotient */
     172                 :            : static GEN
     173                 :        455 : ComputeLift(GEN dataC)
     174                 :            : {
     175                 :            :   long order, i;
     176                 :        455 :   pari_sp av = avma;
     177                 :            :   GEN cyc, surj, eltq, elt;
     178                 :            : 
     179                 :        455 :   order = itos(gel(dataC,1));
     180                 :        455 :   cyc   = gel(dataC,2);
     181                 :        455 :   surj  = gel(dataC,3);
     182                 :            : 
     183                 :        455 :   eltq = EltsOfGroup(order, cyc);
     184                 :        455 :   elt = cgetg(order + 1, t_VEC);
     185         [ +  + ]:       2895 :   for (i = 1; i <= order; i++) gel(elt,i) = inverseimage(surj, gel(eltq,i));
     186                 :            : 
     187                 :        455 :   return gerepileupto(av, elt);
     188                 :            : }
     189                 :            : 
     190                 :            : /* Return c[1],  [c[1]/c[1] = 1,...,c[1]/c[n]] */
     191                 :            : static GEN
     192                 :        625 : init_get_chic(GEN c)
     193                 :            : {
     194                 :        625 :   long i, l = lg(c); /* > 1 */
     195                 :        625 :   GEN C, D = cgetg(l, t_VEC);
     196         [ -  + ]:        625 :   if (l == 1) C = gen_1;
     197                 :            :   else
     198                 :            :   {
     199                 :        625 :     C = gel(c,1); gel(D,1) = gen_1;
     200         [ +  + ]:        970 :     for (i = 2; i < l; i++) gel(D,i) = diviiexact(C, gel(c,i));
     201                 :            :   }
     202                 :        625 :   return mkvec2(C, D);
     203                 :            : }
     204                 :            : 
     205                 :            : static GEN
     206                 :       1710 : get_chic(GEN chi, GEN D)
     207                 :            : {
     208                 :       1710 :   long i, l = lg(chi);
     209                 :       1710 :   GEN chic = cgetg(l, t_VEC);
     210                 :       1710 :   gel(chic,1) = gel(chi,1);
     211         [ +  + ]:       2585 :   for (i = 2; i < l; i++) gel(chic,i) = mulii(gel(chi,i), gel(D,i));
     212                 :       1710 :   return chic;
     213                 :            : }
     214                 :            : 
     215                 :            : /* A character is given by a vector [(c_i), z, d] such that
     216                 :            :    chi(id) = z ^ sum(c_i * a_i) where
     217                 :            :      a_i= log(id) on the generators of bnr
     218                 :            :      z  = exp(2i * Pi / d) */
     219                 :            : /* U is NULL or a ZM */
     220                 :            : static GEN
     221                 :       1710 : get_Char(GEN chi, GEN initc, GEN U, long prec)
     222                 :            : {
     223                 :       1710 :   GEN d, chic = get_chic(chi, gel(initc,2));
     224         [ +  + ]:       1710 :   if (U) chic = ZV_ZM_mul(chic, U);
     225                 :       1710 :   d = ZV_content(chic);
     226         [ +  + ]:       1710 :   if (is_pm1(d)) d = gel(initc,1);
     227                 :            :   else
     228                 :            :   {
     229                 :       1136 :     GEN t = gred_frac2(gel(initc,1), d);
     230                 :       1136 :     chic = ZC_Z_divexact(chic, d);
     231         [ +  + ]:       1136 :     if (typ(t) == t_INT)
     232                 :        788 :       d = t;
     233                 :            :     else
     234                 :            :     {
     235                 :        348 :       d = gel(t,1);
     236                 :        348 :       chic = gmul(gel(t,2), chic);
     237                 :            :     }
     238                 :            :   }
     239                 :       1710 :   return mkvec3(chic, InitRU(d, prec), d);
     240                 :            : }
     241                 :            : 
     242                 :            : /* prime divisors of conductor */
     243                 :            : static GEN
     244                 :        600 : divcond(GEN bnr) { GEN bid = bnr_get_bid(bnr); return gmael(bid,3,1); }
     245                 :            : 
     246                 :            : /* vector of prime ideals dividing bnr but not bnrc */
     247                 :            : static GEN
     248                 :        110 : get_prdiff(GEN bnr, GEN condc)
     249                 :            : {
     250                 :        110 :   GEN prdiff, M = gel(condc,1), D = divcond(bnr), nf = bnr_get_nf(bnr);
     251                 :        110 :   long nd, i, l  = lg(D);
     252                 :        110 :   prdiff = cgetg(l, t_COL);
     253         [ +  + ]:        320 :   for (nd=1, i=1; i < l; i++)
     254         [ +  + ]:        210 :     if (!idealval(nf, M, gel(D,i))) gel(prdiff,nd++) = gel(D,i);
     255                 :        110 :   setlg(prdiff, nd); return prdiff;
     256                 :            : }
     257                 :            : 
     258                 :            : /* Let chi a character defined over bnr and primitive over bnrc, compute the
     259                 :            :  * corresponding primitive character. Returns NULL if bnr = bnrc */
     260                 :            : static GEN
     261                 :       1570 : GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
     262                 :            : {
     263                 :            :   long l;
     264                 :       1570 :   pari_sp av = avma;
     265                 :            :   GEN U, M, cond, condc, initc, Mrc;
     266                 :            : 
     267                 :       1570 :   cond  = bnr_get_mod(bnr);
     268         [ +  + ]:       1570 :   condc = bnr_get_mod(bnrc); if (gequal(cond, condc)) return NULL;
     269                 :            : 
     270                 :        145 :   initc = init_get_chic(bnr_get_cyc(bnr));
     271                 :        145 :   Mrc   = diagonal_shallow(bnr_get_cyc(bnrc));
     272                 :        145 :   M = bnrsurjection(bnr, bnrc);
     273                 :        145 :   (void)ZM_hnfall(shallowconcat(M, Mrc), &U, 1);
     274                 :        145 :   l = lg(M);
     275                 :        145 :   U = rowslice(vecslice(U, l, lg(U)-1), 1, l-1);
     276                 :       1570 :   return gerepilecopy(av, get_Char(chi, initc, U, prec));
     277                 :            : }
     278                 :            : 
     279                 :            : #define ch_chi(x)  gel(x,1)
     280                 :            : #define ch_C(x)    gel(x,2)
     281                 :            : #define ch_bnr(x)  gel(x,3)
     282                 :            : #define ch_4(x)    gel(x,4)
     283                 :            : #define ch_CHI(x)  gel(x,5)
     284                 :            : #define ch_diff(x) gel(x,6)
     285                 :            : #define ch_cond(x) gel(x,7)
     286                 :            : #define ch_CHI0(x) gel(x,8)
     287                 :            : #define ch_comp(x) gel(x,9)
     288                 :            : 
     289                 :            : static GEN
     290                 :        725 : GetDeg(GEN dataCR)
     291                 :            : {
     292                 :        725 :   long i, l = lg(dataCR);
     293                 :        725 :   GEN degs = cgetg(l, t_VECSMALL);
     294                 :            : 
     295         [ +  + ]:       3125 :   for (i = 1; i < l; i++) degs[i] = eulerphiu(itou(gel(ch_CHI(gel(dataCR, i)), 3)));
     296                 :        725 :   return degs;
     297                 :            : }
     298                 :            : 
     299                 :            : /********************************************************************/
     300                 :            : /*                    1rst part: find the field K                   */
     301                 :            : /********************************************************************/
     302                 :            : static GEN AllStark(GEN data, GEN nf, long flag, long prec);
     303                 :            : 
     304                 :            : /* Columns of C [HNF] give the generators of a subgroup of the finite abelian
     305                 :            :  * group A [ in terms of implicit generators ], compute data to work in A/C:
     306                 :            :  * 1) order
     307                 :            :  * 2) structure
     308                 :            :  * 3) the matrix A ->> A/C
     309                 :            :  * 4) the group C */
     310                 :            : static GEN
     311                 :       1815 : InitQuotient(GEN C)
     312                 :            : {
     313                 :            :   long junk;
     314                 :       1815 :   GEN U, D = ZM_snfall_i(C, &U, NULL, 1), h = detcyc(D, &junk);
     315                 :       1815 :   return mkvec4(h, D, U, C);
     316                 :            : }
     317                 :            : 
     318                 :            : /* Let s: A -> B given by P, and let cycA, cycB be the cyclic structure of
     319                 :            :  * A and B, compute the kernel of s. */
     320                 :            : static GEN
     321                 :        345 : ComputeKernel0(GEN P, GEN cycA, GEN cycB)
     322                 :            : {
     323                 :        345 :   pari_sp av = avma;
     324                 :        345 :   long nbA = lg(cycA)-1, rk;
     325                 :        345 :   GEN U, DB = diagonal_shallow(cycB);
     326                 :            : 
     327                 :        345 :   rk = nbA + lg(cycB) - lg(ZM_hnfall(shallowconcat(P, DB), &U, 1));
     328                 :        345 :   U = vecslice(U, 1,rk);
     329                 :        345 :   U = rowslice(U, 1,nbA);
     330                 :        345 :   return gerepileupto(av, ZM_hnfmodid(U, cycA));
     331                 :            : }
     332                 :            : 
     333                 :            : /* Let m and n be two moduli such that n|m and let C be a congruence
     334                 :            :    group modulo n, compute the corresponding congruence group modulo m
     335                 :            :    ie the kernel of the map Clk(m) ->> Clk(n)/C */
     336                 :            : static GEN
     337                 :        345 : ComputeKernel(GEN bnrm, GEN bnrn, GEN dtQ)
     338                 :            : {
     339                 :        345 :   pari_sp av = avma;
     340                 :        345 :   GEN P = ZM_mul(gel(dtQ,3), bnrsurjection(bnrm, bnrn));
     341                 :        345 :   return gerepileupto(av, ComputeKernel0(P, bnr_get_cyc(bnrm), gel(dtQ,2)));
     342                 :            : }
     343                 :            : 
     344                 :            : static GEN
     345                 :        779 : Order(GEN cyc, GEN x)
     346                 :            : {
     347                 :        779 :   pari_sp av = avma;
     348                 :        779 :   long i, l = lg(cyc);
     349                 :        779 :   GEN c,o,f = gen_1;
     350         [ +  + ]:       1933 :   for (i = 1; i < l; i++)
     351                 :            :   {
     352                 :       1154 :     o = gel(cyc,i);
     353                 :       1154 :     c = gcdii(o, gel(x,i));
     354         [ +  + ]:       1154 :     if (!is_pm1(c)) o = diviiexact(o,c);
     355                 :       1154 :     f = lcmii(f, o);
     356                 :            :   }
     357                 :        779 :   return gerepileuptoint(av, f);
     358                 :            : }
     359                 :            : 
     360                 :            : /* Let H be a subgroup of cl(bnr)/sugbroup, return 1 if
     361                 :            :    cl(bnr)/subgoup/H is cyclic and the signature of the
     362                 :            :    corresponding field is equal to sig and no finite prime
     363                 :            :    dividing cond(bnr) is totally split in this field. Return 0
     364                 :            :    otherwise. */
     365                 :            : static long
     366                 :        781 : IsGoodSubgroup(GEN H, GEN bnr, GEN map)
     367                 :            : {
     368                 :        781 :   pari_sp av = avma;
     369                 :            :   long j, f;
     370                 :            :   GEN bnf, mod, mod0, mod1, modH, modH0, p1, p2, p3, p4;
     371                 :            :   GEN bnrH, cycH, iH, qH;
     372                 :            : 
     373                 :        781 :   p1 = InitQuotient(H);
     374                 :        781 :   p2 = gel(p1, 2);
     375                 :            :   /* quotient is non cyclic */
     376 [ +  - ][ +  + ]:        781 :   if ((lg(p2) > 2) && !gequal1(gel(p2, 2))) { avma = av; return 0; }
     377                 :            : 
     378                 :        321 :   bnf  = bnr_get_bnf(bnr);
     379                 :        321 :   mod  = bnr_get_mod(bnr);
     380                 :        321 :   mod0 = gel(mod,1);
     381                 :        321 :   mod1 = gel(mod,2);
     382                 :            : 
     383                 :        321 :   p1 = concat(map, H);
     384                 :        321 :   p2 = ZM_hnfall(p1, &p3, 0);
     385                 :        321 :   setlg(p3, lg(H));
     386         [ +  + ]:       1108 :   for (j = 1; j < lg(p3); j++) setlg(p3[j], lg(H));
     387                 :        321 :   p1 = ZM_hnfmodid(p3, bnr_get_cyc(bnr)); /* H as a subgroup of bnr */
     388                 :        321 :   modH = bnrconductor(bnr, p1, 0);
     389                 :            : 
     390                 :            :   /* is the signature correct? */
     391         [ +  + ]:        321 :   if (!gequal(gel(modH, 2), mod1)) { avma = av; return 0; }
     392                 :            : 
     393                 :            :   /* if the finite part are the same, then it's good */
     394         [ +  + ]:        185 :   if (gequal(gel(modH, 1), mod0)) { avma = av; return 1; }
     395                 :            : 
     396                 :            :   /* we need to check the splitting of primes dividing mod0/p2 */
     397                 :         49 :   modH0 = gel(modH, 1);
     398                 :         49 :   p3 = idealdivexact(bnf, mod0, modH0);
     399                 :         49 :   p4 = gel(idealfactor(bnf, p3), 1);
     400                 :            : 
     401                 :         49 :   bnrH = Buchray(bnf, mkvec2(modH, mod1), nf_INIT|nf_GEN);
     402                 :         49 :   cycH = bnr_get_cyc(bnrH);
     403                 :         49 :   p2 = ZM_mul(bnrsurjection(bnr, bnrH), p1);
     404                 :            :   /* H as a subgroup of bnrH */
     405                 :         49 :   iH = ZM_hnfmodid(p2, cycH);
     406                 :         49 :   qH = InitQuotient(iH);
     407                 :            : 
     408         [ +  + ]:         93 :   for (j = 1; j < lg(p4); j++)
     409                 :            :   {
     410                 :         64 :     GEN pr = gel(p4, j);
     411                 :            :     /* if pr divides modH0, it is ramified, so it's good */
     412         [ +  + ]:         64 :     if (!idealval(bnf, modH0, pr))
     413                 :            :     {
     414                 :            :       /* we compute the inertia degree of pr in bnr(modH)/H*/
     415                 :         44 :       p1 = ZM_ZC_mul(gel(qH, 3), isprincipalray(bnrH, pr));
     416                 :         44 :       p2 = gel(qH, 2);
     417                 :         44 :       f  = itos(Order(p2, p1));
     418         [ +  + ]:         44 :       if (f == 1) { avma = av; return 0; }
     419                 :            :     }
     420                 :            :   }
     421                 :            : 
     422                 :         29 :   avma = av;
     423                 :        781 :   return 1;
     424                 :            : }
     425                 :            : 
     426                 :            : static GEN get_listCR(GEN bnr, GEN dtQ);
     427                 :            : static GEN InitChar(GEN bnr, GEN listCR, long prec);
     428                 :            : 
     429                 :            : /* Given a conductor and a subgroups, return the corresponding
     430                 :            :    complexity and precision required using quickpol. Fill data[5] with
     431                 :            :    listCR */
     432                 :            : static long
     433                 :        225 : CplxModulus(GEN data, long *newprec)
     434                 :            : {
     435                 :        225 :   long pr, ex, dprec = DEFAULTPREC;
     436                 :            :   pari_sp av;
     437                 :        225 :   GEN pol, listCR, cpl, bnr = gel(data,1), nf = checknf(bnr);
     438                 :            : 
     439                 :        225 :   listCR = get_listCR(bnr, gel(data,3));
     440                 :        225 :   for (av = avma;; avma = av)
     441                 :            :   {
     442                 :        225 :     gel(data,5) = InitChar(bnr, listCR, dprec);
     443                 :        225 :     pol = AllStark(data, nf, -1, dprec);
     444                 :        225 :     pr = nbits2extraprec( gexpo(pol) );
     445         [ -  + ]:        225 :     if (pr < 0) pr = 0;
     446                 :        225 :     dprec = maxss(dprec, pr) + EXTRA_PREC;
     447         [ +  - ]:        225 :     if (!gequal0(leading_term(pol)))
     448                 :            :     {
     449                 :        225 :       cpl = RgX_fpnorml2(pol, DEFAULTPREC);
     450         [ +  - ]:        225 :       if (!gequal0(cpl)) break;
     451                 :            :     }
     452         [ #  # ]:          0 :     if (DEBUGLEVEL>1) pari_warn(warnprec, "CplxModulus", dprec);
     453                 :          0 :   }
     454                 :        225 :   ex = gexpo(cpl); avma = av;
     455         [ -  + ]:        225 :   if (DEBUGLEVEL>1) err_printf("cpl = 2^%ld\n", ex);
     456                 :            : 
     457                 :        225 :   gel(data,5) = listCR;
     458                 :        225 :   *newprec = dprec; return ex;
     459                 :            : }
     460                 :            : 
     461                 :            :  /* Let f be a conductor without infinite part and let C be a
     462                 :            :    congruence group modulo f, compute (m,D) such that D is a
     463                 :            :    congruence group of conductor m where m is a multiple of f
     464                 :            :    divisible by all the infinite places but one, D is a subgroup of
     465                 :            :    index 2 of Im(C) in Clk(m), and m is such that the intersection
     466                 :            :    of the subgroups H of Clk(n)/Im(C) such that the quotient is
     467                 :            :    cyclic and no prime divding m, but the one infinite prime, is
     468                 :            :    totally split in the extension corresponding to H is trival.
     469                 :            :    Return bnr(m), D, the quotient Ck(m)/D and Clk(m)/C */
     470                 :            : static GEN
     471                 :        225 : FindModulus(GEN bnr, GEN dtQ, long *newprec)
     472                 :            : {
     473                 :        225 :   const long limnorm = 400;
     474                 :            :   long n, i, narch, maxnorm, minnorm, N, nbidnn, s, c, j, k, nbcand;
     475                 :        225 :   long first = 1, pr, rb, oldcpl = -1, iscyc = 0;
     476                 :        225 :   pari_sp av = avma, av1;
     477                 :        225 :   GEN bnf, nf, f, arch, m, listid, idnormn, bnrm, ImC, rep = NULL;
     478                 :            :   GEN candD, p2;
     479                 :            : 
     480                 :        225 :   bnf = bnr_get_bnf(bnr);
     481                 :        225 :   nf  = bnf_get_nf(bnf);
     482                 :        225 :   N   = nf_get_degree(nf);
     483                 :        225 :   f   = gel(bnr_get_mod(bnr), 1);
     484                 :            : 
     485                 :            :   /* if cpl < rb, it is not necessary to try another modulus */
     486                 :        225 :   rb = expi( powii(mulii(nf_get_disc(nf), ZM_det_triangular(f)), gmul2n(bnr_get_no(bnr), 3)) );
     487                 :            : 
     488                 :            :   /* Initialization of the possible infinite part */
     489                 :        225 :   arch = const_vec(N, gen_1);
     490                 :            : 
     491                 :            :   /* narch = (N == 2)? 1: N; -- if N=2, only one case is necessary */
     492                 :        225 :   narch = N;
     493                 :        225 :   m = mkvec2(NULL, arch);
     494                 :            : 
     495                 :            :   /* go from minnorm up to maxnorm. If necessary, increase these values.
     496                 :            :    * If we cannot find a suitable conductor of norm < limnorm, stop */
     497                 :        225 :   maxnorm = 50;
     498                 :        225 :   minnorm = 1;
     499                 :            : 
     500                 :            :   /* if the extension is cyclic then we _must_ find a suitable conductor */
     501 [ +  + ][ -  + ]:        225 :   if (lg(gel(dtQ,2)) == 2 || gequal1(gmael(dtQ,2,2))) iscyc = 1;
     502                 :            : 
     503         [ -  + ]:        225 :   if (DEBUGLEVEL>1)
     504                 :          0 :     err_printf("Looking for a modulus of norm: ");
     505                 :            : 
     506                 :            :   for(;;)
     507                 :            :   {
     508                 :        225 :     listid = ideallist(nf, maxnorm); /* all ideals of norm <= maxnorm */
     509                 :            : 
     510                 :        225 :     av1 = avma;
     511         [ +  - ]:       1030 :     for (n = minnorm; n <= maxnorm; n++)
     512                 :            :     {
     513         [ -  + ]:       1030 :       if (DEBUGLEVEL>1) err_printf(" %ld", n);
     514                 :       1030 :       avma = av1;
     515                 :            : 
     516                 :       1030 :       idnormn = gel(listid,n);
     517                 :       1030 :       nbidnn  = lg(idnormn) - 1;
     518         [ +  + ]:       1965 :       for (i = 1; i <= nbidnn; i++)
     519                 :            :       { /* finite part of the conductor */
     520                 :       1160 :         gel(m,1) = idealmul(nf, f, gel(idnormn,i));
     521                 :            : 
     522         [ +  + ]:       3365 :         for (s = 1; s <= narch; s++)
     523                 :            :         { /* infinite part */
     524                 :       2430 :           gel(arch,N+1-s) = gen_0;
     525                 :            : 
     526                 :            :           /* compute Clk(m), check if m is a conductor */
     527                 :       2430 :           bnrm = Buchray(bnf, m, nf_INIT|nf_GEN);
     528                 :       2430 :           c = bnrisconductor(bnrm, NULL);
     529                 :       2430 :           gel(arch,N+1-s) = gen_1;
     530         [ +  + ]:       2430 :           if (!c) continue;
     531                 :            : 
     532                 :            :           /* compute Im(C) in Clk(m)... */
     533                 :        345 :           ImC = ComputeKernel(bnrm, bnr, dtQ);
     534                 :            : 
     535                 :            :           /* ... and its subgroups of index 2 with conductor m */
     536                 :        345 :           candD = subgrouplist_cond_sub(bnrm, ImC, mkvec(gen_2));
     537                 :        345 :           nbcand = lg(candD) - 1;
     538         [ +  + ]:        400 :           for (c = 1; c <= nbcand; c++)
     539                 :            :           {
     540                 :        280 :             GEN D  = gel(candD,c);
     541                 :            :             long cpl;
     542                 :            : 
     543                 :            :             /* check if the conductor is suitable */
     544                 :            :             {
     545                 :        280 :               GEN p1 = InitQuotient(D), p2, ord = gel(p1, 1);
     546                 :        280 :               GEN map = gel(p1, 3), cyc = gel(p1, 2), H;
     547                 :        280 :               GEN lH = subgrouplist(cyc, NULL), IK = NULL;
     548                 :        280 :               long ok = 0;
     549                 :            : 
     550                 :            :               /* if the extension is cyclic, then it's suitable */
     551 [ +  + ][ +  + ]:        280 :               if ((lg(cyc) > 2) && !gequal1(gel(cyc, 2)))
     552                 :            :               {
     553         [ +  + ]:       1375 :                 for (j = 1; j < lg(lH); j++)
     554                 :            :                 {
     555                 :       1320 :                   H = gel(lH, j);
     556                 :            :                   /* if IK != NULL and H > IK, no need to test H */
     557         [ +  + ]:       1320 :                   if (IK)
     558                 :            :                   {
     559                 :        828 :                     p1 = RgM_mul(IK, RgM_inv_upper(H));
     560         [ +  + ]:        828 :                     if (RgM_is_ZM(p1)) continue;
     561                 :            :                   }
     562         [ +  + ]:        781 :                   if (IsGoodSubgroup(H, bnrm, map))
     563                 :            :                   {
     564         [ +  + ]:        165 :                     if (!IK)
     565                 :        105 :                       IK = H;
     566                 :            :                     else
     567                 :            :                     {
     568                 :            :                       /* compute the intersection of IK and H */
     569                 :         60 :                       p1 = shallowconcat(IK, H);
     570                 :         60 :                       p1 = ZM_hnfall(p1, &p2, 1);
     571                 :         60 :                       setlg(p2, lg(IK));
     572         [ +  + ]:        210 :                       for (k = 1; k < lg(p2); k++) setlg(p2[k], lg(p1));
     573                 :         60 :                       IK = ZM_mul(IK, p2);
     574                 :         60 :                       IK = ZM_hnf(shallowconcat(IK, diagonal(cyc)));
     575                 :            :                     }
     576         [ +  + ]:        165 :                     if (gequal(ord, ZM_det_triangular(IK)))
     577                 :            :                     {
     578                 :         50 :                       ok = 1; break;
     579                 :            :                     }
     580                 :            :                   }
     581                 :            :                 }
     582         [ +  + ]:        105 :                 if (!ok) continue;
     583                 :            :               }
     584                 :            :             }
     585                 :            : 
     586                 :        225 :             p2 = cgetg(6, t_VEC); /* p2[5] filled in CplxModulus */
     587                 :        225 :             gel(p2,1) = bnrm;
     588                 :        225 :             gel(p2,2) = D;
     589                 :        225 :             gel(p2,3) = InitQuotient(D);
     590                 :        225 :             gel(p2,4) = InitQuotient(ImC);
     591         [ -  + ]:        225 :             if (DEBUGLEVEL>1)
     592                 :          0 :               err_printf("\nTrying modulus = %Ps and subgroup = %Ps\n",
     593                 :            :                          bnr_get_mod(bnrm), D);
     594                 :        225 :             cpl = CplxModulus(p2, &pr);
     595 [ -  + ][ #  # ]:        225 :             if (oldcpl < 0 || cpl < oldcpl)
     596                 :            :             {
     597                 :        225 :               *newprec = pr;
     598         [ -  + ]:        225 :               if (rep) gunclone(rep);
     599                 :        225 :               rep    = gclone(p2);
     600                 :        225 :               oldcpl = cpl;
     601                 :            :             }
     602         [ +  - ]:        225 :             if (oldcpl < rb) goto END; /* OK */
     603                 :            : 
     604         [ #  # ]:          0 :             if (DEBUGLEVEL>1) err_printf("Trying to find another modulus...");
     605                 :          0 :             first = 0;
     606                 :            :           }
     607                 :            :         }
     608         [ -  + ]:        935 :         if (!first) goto END; /* OK */
     609                 :            :       }
     610                 :            :     }
     611                 :            :     /* if necessary compute more ideals */
     612                 :          0 :     minnorm = maxnorm;
     613                 :          0 :     maxnorm <<= 1;
     614 [ #  # ][ #  # ]:          0 :     if (!iscyc && maxnorm > limnorm) return NULL;
     615                 :            : 
     616                 :          0 :   }
     617                 :            : END:
     618         [ -  + ]:        225 :   if (DEBUGLEVEL>1)
     619                 :          0 :     err_printf("No, we're done!\nModulus = %Ps and subgroup = %Ps\n",
     620                 :          0 :                bnr_get_mod(gel(rep,1)), gel(rep,2));
     621                 :        225 :   gel(rep,5) = InitChar(gel(rep,1), gel(rep,5), *newprec);
     622                 :        225 :   return gerepilecopy(av, rep);
     623                 :            : }
     624                 :            : 
     625                 :            : /********************************************************************/
     626                 :            : /*                      2nd part: compute W(X)                      */
     627                 :            : /********************************************************************/
     628                 :            : 
     629                 :            : /* compute the list of W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*),
     630                 :            :  * for all chi in LCHI. All chi have the same conductor (= cond(bnr)).
     631                 :            :  * if check == 0 do not check the result */
     632                 :            : static GEN
     633                 :        585 : ArtinNumber(GEN bnr, GEN LCHI, long check, long prec)
     634                 :            : {
     635                 :        585 :   long ic, i, j, nz, nChar = lg(LCHI)-1;
     636                 :        585 :   pari_sp av = avma, av2, lim;
     637                 :            :   GEN sqrtnc, dc, cond, condZ, cond0, cond1, lambda, nf, T;
     638                 :            :   GEN cyc, vN, vB, diff, vt, idg, idh, zid, gen, z, nchi;
     639                 :            :   GEN indW, W, classe, s0, s, den, muslambda, sarch;
     640                 :            :   CHI_t **lC;
     641                 :            :   GROUP_t G;
     642                 :            : 
     643                 :        585 :   lC = (CHI_t**)new_chunk(nChar + 1);
     644                 :        585 :   indW = cgetg(nChar + 1, t_VECSMALL);
     645                 :        585 :   W = cgetg(nChar + 1, t_VEC);
     646         [ +  + ]:       2180 :   for (ic = 0, i = 1; i <= nChar; i++)
     647                 :            :   {
     648                 :       1595 :     GEN CHI = gel(LCHI,i);
     649         [ +  + ]:       1595 :     if (cmpui(2, gel(CHI,3)) >= 0) { gel(W,i) = gen_1; continue; } /* trivial case */
     650                 :       1010 :     ic++; indW[ic] = i;
     651                 :       1010 :     lC[ic] = (CHI_t*)new_chunk(sizeof(CHI_t));
     652                 :       1010 :     init_CHI_C(lC[ic], CHI);
     653                 :            :   }
     654         [ +  + ]:        585 :   if (!ic) return W;
     655                 :        500 :   nChar = ic;
     656                 :            : 
     657                 :        500 :   nf    = bnr_get_nf(bnr);
     658                 :        500 :   diff  = nf_get_diff(nf);
     659                 :        500 :   T     = nf_get_Tr(nf);
     660                 :        500 :   cond  = bnr_get_mod(bnr);
     661                 :        500 :   cond0 = gel(cond,1); condZ = gcoeff(cond0,1,1);
     662                 :        500 :   cond1 = vec01_to_indices(gel(cond,2));
     663                 :            : 
     664                 :        500 :   sqrtnc = gsqrt(idealnorm(nf, cond0), prec);
     665                 :        500 :   dc  = idealmul(nf, diff, cond0);
     666                 :            : 
     667                 :            :   /* compute a system of elements congruent to 1 mod cond0 and giving all
     668                 :            :      possible signatures for cond1 */
     669                 :        500 :   sarch = nfarchstar(nf, cond0, cond1);
     670                 :            : 
     671                 :            :   /* find lambda in diff.cond such that gcd(lambda.(diff.cond)^-1,cond0) = 1
     672                 :            :      and lambda >> 0 at cond1 */
     673                 :        500 :   lambda = idealappr(nf, dc);
     674                 :        500 :   lambda = set_sign_mod_divisor(nf, NULL, lambda, cond,sarch);
     675                 :        500 :   idg = idealdivexact(nf, lambda, dc);
     676                 :            : 
     677                 :            :   /* find mu in idg such that idh=(mu) / idg is coprime with cond0 and
     678                 :            :      mu >> 0 at cond1 */
     679         [ +  + ]:        500 :   if (!gequal1(gcoeff(idg, 1, 1))) {
     680                 :        490 :     GEN P = divcond(bnr);
     681                 :        490 :     GEN f = famat_mul_shallow(idealfactor(nf, idg),
     682                 :        490 :                           mkmat2(P, zerocol(lg(P)-1)));
     683                 :        490 :     GEN mu = set_sign_mod_divisor(nf, NULL, idealapprfact(nf, f), cond,sarch);
     684                 :        490 :     idh = idealdivexact(nf, mu, idg);
     685                 :        490 :     muslambda = nfdiv(nf, mu, lambda);
     686                 :            :   } else { /* mu = 1 */
     687                 :         10 :     idh = idg;
     688                 :         10 :     muslambda = nfinv(nf, lambda);
     689                 :            :   }
     690                 :        500 :   muslambda = Q_remove_denom(muslambda, &den);
     691                 :        500 :   z = InitRU(den, prec);
     692                 :            : 
     693                 :            :   /* compute a system of generators of (Ok/cond)^*, we'll make them
     694                 :            :    * cond1-positive in the main loop */
     695                 :        500 :   zid = Idealstar(nf, cond0, nf_GEN);
     696                 :        500 :   cyc = gel(zid, 2);
     697                 :        500 :   gen = gel(zid, 3);
     698                 :        500 :   nz = lg(gen) - 1;
     699                 :            : 
     700                 :        500 :   nchi = cgetg(nChar+1, t_VEC);
     701         [ +  + ]:       1510 :   for (ic = 1; ic <= nChar; ic++) gel(nchi,ic) = cgetg(nz + 1, t_VECSMALL);
     702                 :            : 
     703         [ +  + ]:       1000 :   for (i = 1; i <= nz; i++)
     704                 :            :   {
     705         [ -  + ]:        500 :     if (is_bigint(gel(cyc,i)))
     706                 :          0 :       pari_err_OVERFLOW("ArtinNumber [conductor too large]");
     707                 :        500 :     gel(gen,i) = set_sign_mod_divisor(nf, NULL, gel(gen,i), cond,sarch);
     708                 :        500 :     classe = isprincipalray(bnr, gel(gen,i));
     709         [ +  + ]:       1620 :     for (ic = 1; ic <= nChar; ic++) {
     710                 :       1120 :       GEN n = gel(nchi,ic);
     711                 :       1120 :       n[i] = EvalChar_n(lC[ic], classe);
     712                 :            :     }
     713                 :            :   }
     714                 :            : 
     715                 :            :   /* Sum chi(beta) * exp(2i * Pi * Tr(beta * mu / lambda) where beta
     716                 :            :      runs through the classes of (Ok/cond0)^* and beta cond1-positive */
     717                 :            : 
     718                 :        500 :   vt = gel(T,1); /* ( Tr(w_i) )_i */
     719                 :        500 :   vt = ZV_ZM_mul(vt, zk_multable(nf, muslambda)); /*den (Tr(w_i mu/lambda))_i */
     720                 :        500 :   G.cyc = gtovecsmall(cyc);
     721                 :        500 :   G.r = nz;
     722                 :        500 :   G.j = zero_zv(nz);
     723                 :            : 
     724                 :        500 :   vN = cgetg(nChar+1, t_VEC);
     725         [ +  + ]:       1510 :   for (ic = 1; ic <= nChar; ic++) gel(vN,ic) = zero_zv(nz);
     726                 :            : 
     727                 :        500 :   av2 = avma; lim = stack_lim(av2, 1);
     728                 :        500 :   vB = const_vec(nz, gen_1);
     729                 :        500 :   s0 = powgi(z, modii(gel(vt,1), den)); /* for beta = 1 */
     730                 :        500 :   s = const_vec(nChar, s0);
     731                 :            : 
     732         [ +  + ]:      33525 :   while ( (i = NextElt(&G)) )
     733                 :            :   {
     734                 :      33025 :     gel(vB,i) = FpC_red(nfmuli(nf, gel(vB,i), gel(gen,i)), condZ);
     735         [ +  + ]:      33150 :     for (j=1; j<i; j++) gel(vB,j) = gel(vB,i);
     736                 :            : 
     737         [ +  + ]:     143620 :     for (ic = 1; ic <= nChar; ic++)
     738                 :            :     {
     739                 :     110595 :       GEN v = gel(vN,ic), n = gel(nchi,ic);
     740                 :     110595 :       v[i] = Fl_add(v[i], n[i], lC[ic]->ord);
     741         [ +  + ]:     110910 :       for (j=1; j<i; j++) v[j] = v[i];
     742                 :            :     }
     743                 :            : 
     744                 :      33025 :     gel(vB,i) = set_sign_mod_divisor(nf, NULL, gel(vB,i), cond,sarch);
     745                 :      33025 :     s0 = powgi(z, FpV_dotproduct(vt, gel(vB,i), den));
     746         [ +  + ]:     143620 :     for (ic = 1; ic <= nChar; ic++)
     747                 :            :     {
     748                 :     110595 :       GEN n = gel(vN,ic), val = lC[ic]->val[ n[i] ];
     749                 :     110595 :       gel(s,ic) = gadd(gel(s,ic), gmul(val, s0));
     750                 :            :     }
     751                 :            : 
     752         [ -  + ]:      33025 :     if (low_stack(lim, stack_lim(av2, 1)))
     753                 :            :     {
     754         [ #  # ]:          0 :       if (DEBUGMEM > 1) pari_warn(warnmem,"ArtinNumber");
     755                 :          0 :       gerepileall(av2, 2, &s, &vB);
     756                 :            :     }
     757                 :            :   }
     758                 :            : 
     759                 :        500 :   classe = isprincipalray(bnr, idh);
     760                 :        500 :   z = powIs(- (lg(cond1)-1));
     761                 :            : 
     762         [ +  + ]:       1510 :   for (ic = 1; ic <= nChar; ic++)
     763                 :            :   {
     764                 :       1010 :     s0 = gmul(gel(s,ic), EvalChar(lC[ic], classe));
     765                 :       1010 :     s0 = gdiv(s0, sqrtnc);
     766 [ +  + ][ -  + ]:       1010 :     if (check && - expo(subrs(gnorm(s0), 1)) < prec2nbits(prec) >> 1)
     767                 :          0 :       pari_err_BUG("ArtinNumber");
     768                 :       1010 :     gel(W, indW[ic]) = gmul(s0, z);
     769                 :            :   }
     770                 :        585 :   return gerepilecopy(av, W);
     771                 :            : }
     772                 :            : 
     773                 :            : static GEN
     774                 :        485 : ComputeAllArtinNumbers(GEN dataCR, GEN vChar, int check, long prec)
     775                 :            : {
     776                 :        485 :   long j, k, cl = lg(dataCR) - 1, J = lg(vChar)-1;
     777                 :        485 :   GEN W = cgetg(cl+1,t_VEC), WbyCond, LCHI;
     778                 :            : 
     779         [ +  + ]:       1065 :   for (j = 1; j <= J; j++)
     780                 :            :   {
     781                 :        580 :     GEN LChar = gel(vChar,j), ldata = vecpermute(dataCR, LChar);
     782                 :        580 :     GEN dtcr = gel(ldata,1), bnr = ch_bnr(dtcr);
     783                 :        580 :     long l = lg(LChar);
     784                 :            : 
     785         [ -  + ]:        580 :     if (DEBUGLEVEL>1)
     786                 :          0 :       err_printf("* Root Number: cond. no %ld/%ld (%ld chars)\n", j, J, l-1);
     787                 :        580 :     LCHI = cgetg(l, t_VEC);
     788         [ +  + ]:       2170 :     for (k = 1; k < l; k++) gel(LCHI,k) = ch_CHI0(gel(ldata,k));
     789                 :        580 :     WbyCond = ArtinNumber(bnr, LCHI, check, prec);
     790         [ +  + ]:       2170 :     for (k = 1; k < l; k++) W[LChar[k]] = WbyCond[k];
     791                 :            :   }
     792                 :        485 :   return W;
     793                 :            : }
     794                 :            : static GEN
     795                 :          5 : SingleArtinNumber(GEN bnr, GEN chi, long prec)
     796                 :          5 : { return gel(ArtinNumber(bnr, mkvec(chi), 1, prec), 1); }
     797                 :            : 
     798                 :            : /* compute the constant W of the functional equation of
     799                 :            :    Lambda(chi). If flag = 1 then chi is assumed to be primitive */
     800                 :            : GEN
     801                 :          5 : bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
     802                 :            : {
     803                 :            :   long l;
     804                 :          5 :   pari_sp av = avma;
     805                 :            :   GEN cond, condc, bnrc, CHI, cyc;
     806                 :            : 
     807 [ +  - ][ -  + ]:          5 :   if (flag < 0 || flag > 1) pari_err_FLAG("bnrrootnumber");
     808                 :            : 
     809                 :          5 :   checkbnr(bnr);
     810                 :          5 :   cyc = bnr_get_cyc(bnr);
     811                 :          5 :   cond = bnr_get_mod(bnr);
     812                 :          5 :   l    = lg(cyc);
     813                 :            : 
     814 [ +  - ][ -  + ]:          5 :   if (typ(chi) != t_VEC || lg(chi) != l)
     815                 :          0 :     pari_err_TYPE("bnrrootnumber [character]", chi);
     816                 :            : 
     817         [ -  + ]:          5 :   if (flag) condc = NULL;
     818                 :            :   else
     819                 :            :   {
     820                 :          5 :     condc = bnrconductorofchar(bnr, chi);
     821         [ -  + ]:          5 :     if (gequal(cond, condc)) flag = 1;
     822                 :            :   }
     823                 :            : 
     824         [ -  + ]:          5 :   if (flag)
     825                 :            :   {
     826                 :          0 :     GEN initc = init_get_chic(cyc);
     827                 :          0 :     bnrc = bnr;
     828                 :          0 :     CHI = get_Char(chi, initc, NULL, prec);
     829                 :            :   }
     830                 :            :   else
     831                 :            :   {
     832                 :          5 :     bnrc = Buchray(bnr_get_bnf(bnr), condc, nf_INIT|nf_GEN);
     833                 :          5 :     CHI = GetPrimChar(chi, bnr, bnrc, prec);
     834                 :            :   }
     835                 :          5 :   return gerepilecopy(av, SingleArtinNumber(bnrc, CHI, prec));
     836                 :            : }
     837                 :            : 
     838                 :            : /********************************************************************/
     839                 :            : /*               3rd part: initialize the characters                */
     840                 :            : /********************************************************************/
     841                 :            : /* returns a ZV */
     842                 :            : static GEN
     843                 :       1485 : LiftChar(GEN cyc, GEN Mat, GEN chi, GEN D)
     844                 :            : {
     845                 :       1485 :   long lm = lg(cyc), l  = lg(chi), i, j;
     846                 :       1485 :   GEN lchi = cgetg(lm, t_VEC);
     847         [ +  + ]:       3705 :   for (i = 1; i < lm; i++)
     848                 :            :   {
     849                 :       2220 :     pari_sp av = avma;
     850                 :       2220 :     GEN t, s  = mulii(gel(chi,1), gcoeff(Mat, 1, i));
     851         [ +  + ]:       4050 :     for (j = 2; j < l; j++)
     852                 :            :     { /* rarely exercised: D[1]/D[j] could be precomputed */
     853                 :       1830 :       t = mulii(gel(chi,j), diviiexact(gel(D,1), gel(D,j)));
     854                 :       1830 :       s = addii(s, mulii(t, gcoeff(Mat, j, i)));
     855                 :            :     }
     856                 :       2220 :     t = diviiexact(mulii(s, gel(cyc,i)), gel(D,1));
     857                 :       2220 :     gel(lchi,i) = gerepileuptoint(av, modii(t, gel(cyc,i)));
     858                 :            :   }
     859                 :       1485 :   return lchi;
     860                 :            : }
     861                 :            : 
     862                 :            : /* Let chi be a character, A(chi) corresponding to the primes dividing diff
     863                 :            :    at s = flag. If s = 0, returns [r, A] where r is the order of vanishing
     864                 :            :    at s = 0 corresponding to diff. No GC */
     865                 :            : static GEN
     866                 :       1485 : ComputeAChi(GEN dtcr, long *r, long flag, long prec)
     867                 :            : {
     868                 :            :   long l, i;
     869                 :            :   GEN A, diff, chi, bnrc;
     870                 :            : 
     871                 :       1485 :   bnrc = ch_bnr(dtcr);
     872                 :       1485 :   diff = ch_diff(dtcr); l = lg(diff);
     873                 :       1485 :   chi  = ch_CHI0(dtcr);
     874                 :            : 
     875                 :       1485 :   A = gen_1;
     876                 :       1485 :   *r = 0;
     877         [ +  + ]:       1560 :   for (i = 1; i < l; i++)
     878                 :            :   {
     879                 :         75 :     GEN pr = gel(diff,i), B;
     880                 :         75 :     GEN z = ComputeImagebyChar(chi, isprincipalray(bnrc, pr));
     881                 :            : 
     882         [ -  + ]:         75 :     if (flag)
     883                 :          0 :       B = gsubsg(1, gdiv(z, pr_norm(pr)));
     884         [ +  + ]:         75 :     else if (gequal1(z))
     885                 :            :     {
     886                 :         10 :       B = glog(pr_norm(pr), prec);
     887                 :         10 :       (*r)++;
     888                 :            :     }
     889                 :            :     else
     890                 :         65 :       B = gsubsg(1, z);
     891                 :         75 :     A = gmul(A, B);
     892                 :            :   }
     893                 :       1485 :   return A;
     894                 :            : }
     895                 :            : /* simplified version of ComputeAchi: return 1 if L(0,chi) = 0 */
     896                 :            : static int
     897                 :       1495 : L_vanishes_at_0(GEN dtcr)
     898                 :            : {
     899                 :            :   long l, i;
     900                 :            :   GEN diff, chi, bnrc;
     901                 :            : 
     902                 :       1495 :   bnrc = ch_bnr(dtcr);
     903                 :       1495 :   diff = ch_diff(dtcr); l = lg(diff);
     904                 :       1495 :   chi  = ch_CHI0(dtcr);
     905         [ +  + ]:       1555 :   for (i = 1; i < l; i++)
     906                 :            :   {
     907                 :         80 :     GEN pr = gel(diff,i);
     908                 :         80 :     GEN z = ComputeImagebyChar(chi, isprincipalray(bnrc, pr));
     909         [ +  + ]:         80 :     if (gequal1(z)) return 1;
     910                 :            :   }
     911                 :       1495 :   return 0;
     912                 :            : }
     913                 :            : 
     914                 :            : static GEN
     915                 :        580 : _data4(GEN arch, long r1, long r2)
     916                 :            : {
     917                 :        580 :   GEN z = cgetg(5, t_VECSMALL);
     918                 :        580 :   long i, b, q = 0;
     919                 :            : 
     920 [ +  + ][ +  + ]:       1940 :   for (i=1; i<=r1; i++) if (signe(gel(arch,i))) q++;
     921                 :        580 :   z[1] = q; b = r1 - q;
     922                 :        580 :   z[2] = b;
     923                 :        580 :   z[3] = r2;
     924                 :        580 :   z[4] = maxss(b+r2+1, r2+q);
     925                 :        580 :   return z;
     926                 :            : }
     927                 :            : 
     928                 :            : /* Given a list [chi, F = cond(chi)] of characters over Cl(bnr), compute a
     929                 :            :    vector dataCR containing for each character:
     930                 :            :    2: the constant C(F) [t_REAL]
     931                 :            :    3: bnr(F)
     932                 :            :    4: [q, r1 - q, r2, rc] where
     933                 :            :         q = number of real places in F
     934                 :            :         rc = max{r1 + r2 - q + 1, r2 + q}
     935                 :            :    6: diff(chi) primes dividing m but not F
     936                 :            :    7: finite part of F
     937                 :            : 
     938                 :            :    1: chi
     939                 :            :    5: [(c_i), z, d] in bnr(m)
     940                 :            :    8: [(c_i), z, d] in bnr(F)
     941                 :            :    9: if NULL then does not compute (for AllStark) */
     942                 :            : static GEN
     943                 :        480 : InitChar(GEN bnr, GEN listCR, long prec)
     944                 :            : {
     945                 :        480 :   GEN bnf = checkbnf(bnr), nf = bnf_get_nf(bnf);
     946                 :            :   GEN modul, dk, C, dataCR, chi, cond, initc;
     947                 :            :   long N, r1, r2, prec2, i, j, l;
     948                 :        480 :   pari_sp av = avma;
     949                 :            : 
     950                 :        480 :   modul = bnr_get_mod(bnr);
     951                 :        480 :   dk    = nf_get_disc(nf);
     952                 :        480 :   N     = nf_get_degree(nf);
     953                 :        480 :   nf_get_sign(nf, &r1,&r2);
     954                 :        480 :   prec2 = precdbl(prec) + EXTRA_PREC;
     955                 :        480 :   C     = gmul2n(sqrtr_abs(divir(dk, powru(mppi(prec2),N))), -r2);
     956                 :        480 :   initc = init_get_chic( bnr_get_cyc(bnr) );
     957                 :            : 
     958                 :        480 :   dataCR = cgetg_copy(listCR, &l);
     959         [ +  + ]:       2045 :   for (i = 1; i < l; i++)
     960                 :            :   {
     961                 :       1565 :     GEN olddtcr, dtcr = cgetg(10, t_VEC);
     962                 :       1565 :     gel(dataCR,i) = dtcr;
     963                 :            : 
     964                 :       1565 :     chi  = gmael(listCR, i, 1);
     965                 :       1565 :     cond = gmael(listCR, i, 2);
     966                 :            : 
     967                 :            :     /* do we already know the invariants of chi? */
     968                 :       1565 :     olddtcr = NULL;
     969         [ +  + ]:       1930 :     for (j = 1; j < i; j++)
     970         [ +  + ]:       1350 :       if (gequal(cond, gmael(listCR,j,2))) { olddtcr = gel(dataCR,j); break; }
     971                 :            : 
     972         [ +  + ]:       1565 :     if (!olddtcr)
     973                 :            :     {
     974                 :        580 :       ch_C(dtcr) = gmul(C, gsqrt(ZM_det_triangular(gel(cond,1)), prec2));
     975                 :        580 :       ch_4(dtcr) = _data4(gel(cond,2),r1,r2);
     976                 :        580 :       ch_cond(dtcr) = gel(cond,1);
     977         [ +  + ]:        580 :       if (gequal(cond,modul))
     978                 :            :       {
     979                 :        470 :         ch_bnr(dtcr) = bnr;
     980                 :        470 :         ch_diff(dtcr) = cgetg(1, t_VEC);
     981                 :            :       }
     982                 :            :       else
     983                 :            :       {
     984                 :        110 :         ch_bnr(dtcr) = Buchray(bnf, cond, nf_INIT|nf_GEN);
     985                 :        110 :         ch_diff(dtcr) = get_prdiff(bnr, cond);
     986                 :            :       }
     987                 :            :     }
     988                 :            :     else
     989                 :            :     {
     990                 :        985 :       ch_C(dtcr) = ch_C(olddtcr);
     991                 :        985 :       ch_bnr(dtcr) = ch_bnr(olddtcr);
     992                 :        985 :       ch_4(dtcr) = ch_4(olddtcr);
     993                 :        985 :       ch_diff(dtcr) = ch_diff(olddtcr);
     994                 :        985 :       ch_cond(dtcr) = ch_cond(olddtcr);
     995                 :            :     }
     996                 :            : 
     997                 :       1565 :     ch_chi(dtcr) = chi; /* the character */
     998                 :       1565 :     ch_CHI(dtcr) = get_Char(chi,initc,NULL,prec); /* associated to bnr(m) */
     999                 :       1565 :     ch_comp(dtcr) = gen_1; /* compute this character (by default) */
    1000                 :       1565 :     chi = GetPrimChar(chi, bnr, ch_bnr(dtcr), prec2);
    1001         [ +  + ]:       1565 :     if (!chi) chi = ch_CHI(dtcr);
    1002                 :       1565 :     ch_CHI0(dtcr) = chi;
    1003                 :            :   }
    1004                 :            : 
    1005                 :        480 :   return gerepilecopy(av, dataCR);
    1006                 :            : }
    1007                 :            : 
    1008                 :            : /* compute the list of characters to consider for AllStark and
    1009                 :            :    initialize precision-independent data to compute with them */
    1010                 :            : static GEN
    1011                 :        225 : get_listCR(GEN bnr, GEN dtQ)
    1012                 :            : {
    1013                 :            :   GEN MrD, listCR, vecchi, lchi, Surj, cond, Mr, d, allCR;
    1014                 :            :   long hD, h, nc, i, j, tnc;
    1015                 :            : 
    1016                 :        225 :   Surj = gel(dtQ,3);
    1017                 :        225 :   MrD  = gel(dtQ,2);
    1018                 :        225 :   Mr   = bnr_get_cyc(bnr);
    1019                 :        225 :   hD   = itos(gel(dtQ,1));
    1020                 :        225 :   h    = hD >> 1;
    1021                 :            : 
    1022                 :        225 :   listCR = cgetg(h + 1, t_VEC); /* non-conjugate characters */
    1023                 :        225 :   nc  = 1;
    1024                 :        225 :   allCR  = cgetg(h + 1, t_VEC); /* all characters, including conjugates */
    1025                 :        225 :   tnc = 1;
    1026                 :            : 
    1027                 :        225 :   vecchi = EltsOfGroup(hD, MrD);
    1028                 :            : 
    1029         [ +  + ]:       1555 :   for (i = 1; tnc <= h; i++)
    1030                 :            :   {
    1031                 :            :     /* lift a character of D in Clk(m) */
    1032                 :       1330 :     lchi = LiftChar(Mr, Surj, gel(vecchi,i), MrD);
    1033                 :            : 
    1034         [ +  + ]:       4853 :     for (j = 1; j < tnc; j++)
    1035         [ +  + ]:       3540 :       if (ZV_equal(lchi, gel(allCR,j))) break;
    1036         [ +  + ]:       1330 :     if (j != tnc) continue;
    1037                 :            : 
    1038                 :       1313 :     cond = bnrconductorofchar(bnr, lchi);
    1039         [ +  + ]:       1313 :     if (gequal0(gel(cond,2))) continue;
    1040                 :            : 
    1041                 :            :     /* the infinite part of chi is non trivial */
    1042                 :        735 :     gel(listCR,nc++) = mkvec2(lchi, cond);
    1043                 :        735 :     gel(allCR,tnc++) = lchi;
    1044                 :            : 
    1045                 :            :     /* if chi is not real, add its conjugate character to allCR */
    1046                 :        735 :     d = Order(Mr, lchi);
    1047         [ +  + ]:        735 :     if (!equaliu(d, 2))
    1048                 :        465 :       gel(allCR,tnc++) = ConjChar(lchi, Mr);
    1049                 :            :   }
    1050                 :        225 :   setlg(listCR, nc); return listCR;
    1051                 :            : }
    1052                 :            : 
    1053                 :            : /* recompute dataCR with the new precision */
    1054                 :            : static GEN
    1055                 :          5 : CharNewPrec(GEN dataCR, GEN nf, long prec)
    1056                 :            : {
    1057                 :            :   GEN dk, C, p1;
    1058                 :            :   long N, l, j, prec2;
    1059                 :            : 
    1060                 :          5 :   dk    =  nf_get_disc(nf);
    1061                 :          5 :   N     =  nf_get_degree(nf);
    1062                 :          5 :   prec2 = precdbl(prec) + EXTRA_PREC;
    1063                 :            : 
    1064                 :          5 :   C = sqrtr(divir(absi(dk), powru(mppi(prec2), N)));
    1065                 :            : 
    1066                 :          5 :   l = lg(dataCR);
    1067         [ +  + ]:         30 :   for (j = 1; j < l; j++)
    1068                 :            :   {
    1069                 :         25 :     GEN dtcr = gel(dataCR,j);
    1070                 :         25 :     ch_C(dtcr) = gmul(C, gsqrt(ZM_det_triangular(ch_cond(dtcr)), prec2));
    1071                 :            : 
    1072                 :         25 :     gmael(ch_bnr(dtcr), 1, 7) = nf;
    1073                 :            : 
    1074                 :         25 :     p1 = ch_CHI( dtcr); gel(p1,2) = InitRU(gel(p1,3), prec2);
    1075                 :         25 :     p1 = ch_CHI0(dtcr); gel(p1,2) = InitRU(gel(p1,3), prec2);
    1076                 :            :   }
    1077                 :            : 
    1078                 :          5 :   return dataCR;
    1079                 :            : }
    1080                 :            : 
    1081                 :            : /********************************************************************/
    1082                 :            : /*             4th part: compute the coefficients an(chi)           */
    1083                 :            : /*                                                                  */
    1084                 :            : /* matan entries are arrays of ints containing the coefficients of  */
    1085                 :            : /* an(chi) as a polmod modulo polcyclo(order(chi))                     */
    1086                 :            : /********************************************************************/
    1087                 :            : 
    1088                 :            : static void
    1089                 :     977416 : _0toCoeff(int *rep, long deg)
    1090                 :            : {
    1091                 :            :   long i;
    1092         [ +  + ]:    3917810 :   for (i=0; i<deg; i++) rep[i] = 0;
    1093                 :     977416 : }
    1094                 :            : 
    1095                 :            : /* transform a polmod into Coeff */
    1096                 :            : static void
    1097                 :     351715 : Polmod2Coeff(int *rep, GEN polmod, long deg)
    1098                 :            : {
    1099                 :            :   long i;
    1100         [ +  + ]:     351715 :   if (typ(polmod) == t_POLMOD)
    1101                 :            :   {
    1102                 :     246522 :     GEN pol = gel(polmod,2);
    1103                 :     246522 :     long d = degpol(pol);
    1104                 :            : 
    1105                 :     246522 :     pol += 2;
    1106         [ +  + ]:    1008105 :     for (i=0; i<=d; i++) rep[i] = itos(gel(pol,i));
    1107         [ +  + ]:     537500 :     for (   ; i<deg; i++) rep[i] = 0;
    1108                 :            :   }
    1109                 :            :   else
    1110                 :            :   {
    1111                 :     105193 :     rep[0] = itos(polmod);
    1112         [ +  + ]:     112638 :     for (i=1; i<deg; i++) rep[i] = 0;
    1113                 :            :   }
    1114                 :     351715 : }
    1115                 :            : 
    1116                 :            : /* initialize a deg * n matrix of ints */
    1117                 :            : static int**
    1118                 :       2715 : InitMatAn(long n, long deg, long flag)
    1119                 :            : {
    1120                 :            :   long i, j;
    1121                 :       2715 :   int *a, **A = (int**)pari_malloc((n+1)*sizeof(int*));
    1122                 :       2715 :   A[0] = NULL;
    1123         [ +  + ]:    6036331 :   for (i = 1; i <= n; i++)
    1124                 :            :   {
    1125                 :    6033616 :     a = (int*)pari_malloc(deg*sizeof(int));
    1126 [ +  + ][ +  + ]:    6033616 :     A[i] = a; a[0] = (i == 1 || flag);
    1127         [ +  + ]:   17002900 :     for (j = 1; j < deg; j++) a[j] = 0;
    1128                 :            :   }
    1129                 :       2715 :   return A;
    1130                 :            : }
    1131                 :            : 
    1132                 :            : static void
    1133                 :       4295 : FreeMat(int **A, long n)
    1134                 :            : {
    1135                 :            :   long i;
    1136         [ +  + ]:    6046171 :   for (i = 0; i <= n; i++)
    1137         [ +  + ]:    6041876 :     if (A[i]) pari_free((void*)A[i]);
    1138                 :       4295 :   pari_free((void*)A);
    1139                 :       4295 : }
    1140                 :            : 
    1141                 :            : /* initialize Coeff reduction */
    1142                 :            : static int**
    1143                 :       1580 : InitReduction(GEN CHI, long deg)
    1144                 :            : {
    1145                 :            :   long j;
    1146                 :       1580 :   pari_sp av = avma;
    1147                 :            :   int **A;
    1148                 :            :   GEN d, polmod, pol;
    1149                 :            : 
    1150                 :       1580 :   A   = (int**)pari_malloc(deg*sizeof(int*));
    1151                 :       1580 :   d   = gel(CHI,3);
    1152                 :       1580 :   pol = polcyclo(itos(d), 0);
    1153         [ +  + ]:       7125 :   for (j = 0; j < deg; j++)
    1154                 :            :   {
    1155                 :       5545 :     A[j] = (int*)pari_malloc(deg*sizeof(int));
    1156                 :       5545 :     polmod = gmodulo(monomial(gen_1, deg+j, 0), pol);
    1157                 :       5545 :     Polmod2Coeff(A[j], polmod, deg);
    1158                 :            :   }
    1159                 :            : 
    1160                 :       1580 :   avma = av; return A;
    1161                 :            : }
    1162                 :            : 
    1163                 :            : #if 0
    1164                 :            : void
    1165                 :            : pan(int **an, long n, long deg)
    1166                 :            : {
    1167                 :            :   long i,j;
    1168                 :            :   for (i = 1; i <= n; i++)
    1169                 :            :   {
    1170                 :            :     err_printf("n = %ld: ",i);
    1171                 :            :     for (j = 0; j < deg; j++) err_printf("%d ",an[i][j]);
    1172                 :            :     err_printf("\n");
    1173                 :            :   }
    1174                 :            : }
    1175                 :            : #endif
    1176                 :            : 
    1177                 :            : /* returns 0 if c is zero, 1 otherwise. */
    1178                 :            : static int
    1179                 :    7962731 : IsZero(int* c, long deg)
    1180                 :            : {
    1181                 :            :   long i;
    1182         [ +  + ]:   24666909 :   for (i = 0; i < deg; i++)
    1183         [ +  + ]:   18627767 :     if (c[i]) return 0;
    1184                 :    7962731 :   return 1;
    1185                 :            : }
    1186                 :            : 
    1187                 :            : /* set c0 <-- c0 * c1 */
    1188                 :            : static void
    1189                 :    1449358 : MulCoeff(int *c0, int* c1, int** reduc, long deg)
    1190                 :            : {
    1191                 :            :   long i,j;
    1192                 :            :   int c, *T;
    1193                 :            : 
    1194         [ +  + ]:    2198500 :   if (IsZero(c0,deg)) return;
    1195                 :            : 
    1196                 :     749142 :   T = (int*)new_chunk(2*deg);
    1197         [ +  + ]:    7121234 :   for (i = 0; i < 2*deg; i++)
    1198                 :            :   {
    1199                 :    6372092 :     c = 0;
    1200         [ +  + ]:   64812990 :     for (j = 0; j <= i; j++)
    1201 [ +  + ][ +  + ]:   58440898 :       if (j < deg && j > i - deg) c += c0[j] * c1[i-j];
    1202                 :    6372092 :     T[i] = c;
    1203                 :            :   }
    1204         [ +  + ]:    3935188 :   for (i = 0; i < deg; i++)
    1205                 :            :   {
    1206                 :    3186046 :     c = T[i];
    1207         [ +  + ]:   30813472 :     for (j = 0; j < deg; j++) c += reduc[j][i] * T[deg+j];
    1208                 :    3186046 :     c0[i] = c;
    1209                 :            :   }
    1210                 :            : }
    1211                 :            : 
    1212                 :            : /* c0 <- c0 + c1 * c2 */
    1213                 :            : static void
    1214                 :    6513373 : AddMulCoeff(int *c0, int *c1, int* c2, int** reduc, long deg)
    1215                 :            : {
    1216                 :            :   long i, j;
    1217                 :            :   pari_sp av;
    1218                 :            :   int c, *t;
    1219                 :            : 
    1220         [ +  + ]:    6513373 :   if (IsZero(c2,deg)) return;
    1221         [ +  + ]:    1174447 :   if (!c1) /* c1 == 1 */
    1222                 :            :   {
    1223         [ +  + ]:     763966 :     for (i = 0; i < deg; i++) c0[i] += c2[i];
    1224                 :     320756 :     return;
    1225                 :            :   }
    1226                 :     853691 :   av = avma;
    1227                 :     853691 :   t = (int*)new_chunk(2*deg); /* = c1 * c2, not reduced */
    1228         [ +  + ]:    5923873 :   for (i = 0; i < 2*deg; i++)
    1229                 :            :   {
    1230                 :    5070182 :     c = 0;
    1231         [ +  + ]:   34172335 :     for (j = 0; j <= i; j++)
    1232 [ +  + ][ +  + ]:   29102153 :       if (j < deg && j > i - deg) c += c1[j] * c2[i-j];
    1233                 :    5070182 :     t[i] = c;
    1234                 :            :   }
    1235         [ +  + ]:    3388782 :   for (i = 0; i < deg; i++)
    1236                 :            :   {
    1237                 :    2535091 :     c = t[i];
    1238         [ +  + ]:   15818622 :     for (j = 0; j < deg; j++) c += reduc[j][i] * t[deg+j];
    1239                 :    2535091 :     c0[i] += c;
    1240                 :            :   }
    1241                 :    6513373 :   avma = av;
    1242                 :            : }
    1243                 :            : 
    1244                 :            : /* evaluate the Coeff. No Garbage collector */
    1245                 :            : static GEN
    1246                 :    3596616 : EvalCoeff(GEN z, int* c, long deg)
    1247                 :            : {
    1248                 :            :   long i,j;
    1249                 :            :   GEN e, r;
    1250                 :            : 
    1251         [ -  + ]:    3596616 :   if (!c) return gen_0;
    1252                 :            : #if 0
    1253                 :            :   /* standard Horner */
    1254                 :            :   e = stoi(c[deg - 1]);
    1255                 :            :   for (i = deg - 2; i >= 0; i--)
    1256                 :            :     e = gadd(stoi(c[i]), gmul(z, e));
    1257                 :            : #else
    1258                 :            :   /* specific attention to sparse polynomials */
    1259                 :    3596616 :   e = NULL;
    1260         [ +  + ]:    5165422 :   for (i = deg-1; i >=0; i=j-1)
    1261                 :            :   {
    1262         [ +  + ]:   10524033 :     for (j=i; c[j] == 0; j--)
    1263         [ +  + ]:    8955227 :       if (j==0)
    1264                 :            :       {
    1265         [ +  + ]:    2880972 :         if (!e) return NULL;
    1266         [ +  + ]:     278849 :         if (i!=j) z = gpowgs(z,i-j+1);
    1267                 :     278849 :         return gmul(e,z);
    1268                 :            :       }
    1269         [ +  + ]:    1568806 :     if (e)
    1270                 :            :     {
    1271         [ +  + ]:     574313 :       r = (i==j)? z: gpowgs(z,i-j+1);
    1272                 :     574313 :       e = gadd(gmul(e,r), stoi(c[j]));
    1273                 :            :     }
    1274                 :            :     else
    1275                 :     994493 :       e = stoi(c[j]);
    1276                 :            :   }
    1277                 :            : #endif
    1278                 :    3596616 :   return e;
    1279                 :            : }
    1280                 :            : 
    1281                 :            : /* copy the n * (m+1) array matan */
    1282                 :            : static void
    1283                 :     326356 : CopyCoeff(int** a, int** a2, long n, long m)
    1284                 :            : {
    1285                 :            :   long i,j;
    1286                 :            : 
    1287         [ +  + ]:    5474668 :   for (i = 1; i <= n; i++)
    1288                 :            :   {
    1289                 :    5148312 :     int *b = a[i], *b2 = a2[i];
    1290         [ +  + ]:   18282198 :     for (j = 0; j < m; j++) b2[j] = b[j];
    1291                 :            :   }
    1292                 :     326356 : }
    1293                 :            : 
    1294                 :            : /* return q*p if <= n. Beware overflow */
    1295                 :            : static long
    1296                 :     544725 : next_pow(long q, long p, long n)
    1297                 :            : {
    1298                 :     544725 :   const GEN x = muluu((ulong)q, (ulong)p);
    1299                 :     544725 :   const ulong qp = uel(x,2);
    1300 [ +  - ][ +  + ]:     544725 :   return (lgefint(x) > 3 || qp > (ulong)n)? 0: qp;
    1301                 :            : }
    1302                 :            : 
    1303                 :            : static void
    1304                 :     326356 : an_AddMul(int **an,int **an2, long np, long n, long deg, GEN chi, int **reduc)
    1305                 :            : {
    1306                 :     326356 :   GEN chi2 = chi;
    1307                 :            :   long q, qk, k;
    1308                 :     326356 :   int *c, *c2 = (int*)new_chunk(deg);
    1309                 :            : 
    1310                 :     326356 :   CopyCoeff(an, an2, n/np, deg);
    1311                 :     326356 :   for (q=np;;)
    1312                 :            :   {
    1313         [ +  + ]:     350718 :     if (gequal1(chi2)) c = NULL; else { Polmod2Coeff(c2, chi2, deg); c = c2; }
    1314         [ +  + ]:    6864091 :     for(k = 1, qk = q; qk <= n; k++, qk += q)
    1315                 :    6513373 :       AddMulCoeff(an[qk], c, an2[k], reduc, deg);
    1316         [ +  + ]:     350718 :     if (! (q = next_pow(q,np, n)) ) break;
    1317                 :            : 
    1318                 :      24362 :     chi2 = gmul(chi2, chi);
    1319                 :      24362 :   }
    1320                 :     326356 : }
    1321                 :            : 
    1322                 :            : /* correct the coefficients an(chi) according with diff(chi) in place */
    1323                 :            : static void
    1324                 :       1580 : CorrectCoeff(GEN dtcr, int** an, int** reduc, long n, long deg)
    1325                 :            : {
    1326                 :       1580 :   pari_sp av = avma;
    1327                 :            :   long lg, j, np;
    1328                 :            :   pari_sp av1;
    1329                 :            :   int **an2;
    1330                 :            :   GEN bnrc, diff, chi, pr;
    1331                 :            :   CHI_t C;
    1332                 :            : 
    1333                 :       1580 :   diff = ch_diff(dtcr); lg = lg(diff) - 1;
    1334         [ +  + ]:       3160 :   if (!lg) return;
    1335                 :            : 
    1336         [ -  + ]:         75 :   if (DEBUGLEVEL>2) err_printf("diff(CHI) = %Ps", diff);
    1337                 :         75 :   bnrc =  ch_bnr(dtcr);
    1338                 :         75 :   init_CHI_alg(&C, ch_CHI0(dtcr));
    1339                 :            : 
    1340                 :         75 :   an2 = InitMatAn(n, deg, 0);
    1341                 :         75 :   av1 = avma;
    1342         [ +  + ]:        155 :   for (j = 1; j <= lg; j++)
    1343                 :            :   {
    1344                 :         80 :     pr = gel(diff,j);
    1345                 :         80 :     np = itos( pr_norm(pr) );
    1346                 :            : 
    1347                 :         80 :     chi  = EvalChar(&C, isprincipalray(bnrc, pr));
    1348                 :            : 
    1349                 :         80 :     an_AddMul(an,an2,np,n,deg,chi,reduc);
    1350                 :         80 :     avma = av1;
    1351                 :            :   }
    1352                 :       1580 :   FreeMat(an2, n); avma = av;
    1353                 :            : }
    1354                 :            : 
    1355                 :            : /* compute the coefficients an in the general case */
    1356                 :            : static int**
    1357                 :       1060 : ComputeCoeff(GEN dtcr, LISTray *R, long n, long deg)
    1358                 :            : {
    1359                 :       1060 :   pari_sp av = avma, av2;
    1360                 :            :   long i, l, np;
    1361                 :            :   int **an, **reduc, **an2;
    1362                 :            :   GEN L, CHI, chi;
    1363                 :            :   CHI_t C;
    1364                 :            : 
    1365                 :       1060 :   CHI = ch_CHI(dtcr); init_CHI_alg(&C, CHI);
    1366                 :            : 
    1367                 :       1060 :   an  = InitMatAn(n, deg, 0);
    1368                 :       1060 :   an2 = InitMatAn(n, deg, 0);
    1369                 :       1060 :   reduc  = InitReduction(CHI, deg);
    1370                 :       1060 :   av2 = avma;
    1371                 :            : 
    1372                 :       1060 :   L = R->L1; l = lg(L);
    1373         [ +  + ]:     327336 :   for (i=1; i<l; i++, avma = av2)
    1374                 :            :   {
    1375                 :     326276 :     np = L[i];
    1376                 :     326276 :     chi  = EvalChar(&C, gel(R->L1ray,i));
    1377                 :     326276 :     an_AddMul(an,an2,np,n,deg,chi,reduc);
    1378                 :            :   }
    1379                 :       1060 :   FreeMat(an2, n);
    1380                 :            : 
    1381                 :       1060 :   CorrectCoeff(dtcr, an, reduc, n, deg);
    1382                 :       1060 :   FreeMat(reduc, deg-1);
    1383                 :       1060 :   avma = av; return an;
    1384                 :            : }
    1385                 :            : 
    1386                 :            : /********************************************************************/
    1387                 :            : /*              5th part: compute L-functions at s=1                */
    1388                 :            : /********************************************************************/
    1389                 :            : static void
    1390                 :        335 : deg11(LISTray *R, long p, GEN bnr, GEN pr) {
    1391                 :        335 :   GEN z = isprincipalray(bnr, pr);
    1392                 :        335 :   vecsmalltrunc_append(R->L1, p);
    1393                 :        335 :   vectrunc_append(R->L1ray, z);
    1394                 :        335 : }
    1395                 :            : static void
    1396                 :      20859 : deg12(LISTray *R, long p, GEN bnr, GEN Lpr) {
    1397                 :      20859 :   GEN z = isprincipalray(bnr, gel(Lpr,1));
    1398                 :      20859 :   vecsmalltrunc_append(R->L11, p);
    1399                 :      20859 :   vectrunc_append(R->L11ray, z);
    1400                 :      20859 : }
    1401                 :            : static void
    1402                 :         25 : deg0(LISTray *R, long p) { vecsmalltrunc_append(R->L0, p); }
    1403                 :            : static void
    1404                 :      21868 : deg2(LISTray *R, long p) { vecsmalltrunc_append(R->L2, p); }
    1405                 :            : 
    1406                 :            : /* pi(x) <= ?? */
    1407                 :            : static long
    1408                 :        485 : PiBound(ulong x)
    1409                 :            : {
    1410                 :        485 :   double lx = log((double)x);
    1411                 :        485 :   return 1 + (long) (x/lx * (1 + 3/(2*lx)));
    1412                 :            : }
    1413                 :            : 
    1414                 :            : static void
    1415                 :        140 : InitPrimesQuad(GEN bnr, ulong N0, LISTray *R)
    1416                 :            : {
    1417                 :        140 :   pari_sp av = avma;
    1418                 :        140 :   GEN bnf = bnr_get_bnf(bnr), cond = gel(bnr_get_mod(bnr), 1);
    1419                 :        140 :   long p,i,l, condZ = itos(gcoeff(cond,1,1)), contZ = itos(content(cond));
    1420                 :        140 :   GEN prime, Lpr, nf = bnf_get_nf(bnf), dk = nf_get_disc(nf);
    1421                 :            :   forprime_t T;
    1422                 :            : 
    1423                 :        140 :   l = 1 + PiBound(N0);
    1424                 :        140 :   R->L0 = vecsmalltrunc_init(l);
    1425                 :        140 :   R->L2 = vecsmalltrunc_init(l); R->condZ = condZ;
    1426                 :        140 :   R->L1 = vecsmalltrunc_init(l); R->L1ray = vectrunc_init(l);
    1427                 :        140 :   R->L11= vecsmalltrunc_init(l); R->L11ray= vectrunc_init(l);
    1428                 :        140 :   prime = utoipos(2);
    1429                 :        140 :   u_forprime_init(&T, 2, N0);
    1430         [ +  + ]:      43227 :   while ( (p = u_forprime_next(&T)) )
    1431                 :            :   {
    1432                 :      43087 :     prime[2] = p;
    1433      [ +  +  + ]:      43087 :     switch (kroiu(dk, p))
    1434                 :            :     {
    1435                 :            :     case -1: /* inert */
    1436         [ +  + ]:      21873 :       if (condZ % p == 0) deg0(R,p); else deg2(R,p);
    1437                 :      21873 :       break;
    1438                 :            :     case 1: /* split */
    1439                 :      20994 :       Lpr = idealprimedec(nf, prime);
    1440         [ +  + ]:      20994 :       if      (condZ % p != 0) deg12(R, p, bnr, Lpr);
    1441         [ +  + ]:        135 :       else if (contZ % p == 0) deg0(R,p);
    1442                 :            :       else
    1443                 :            :       {
    1444         [ +  - ]:        130 :         GEN pr = idealval(nf, cond, gel(Lpr,1))? gel(Lpr,2): gel(Lpr,1);
    1445                 :        130 :         deg11(R, p, bnr, pr);
    1446                 :            :       }
    1447                 :      20994 :       break;
    1448                 :            :     default: /* ramified */
    1449         [ +  + ]:        220 :       if (condZ % p == 0) deg0(R,p);
    1450                 :            :       else {
    1451                 :        205 :         GEN pr = gel(idealprimedec(nf,prime),1);
    1452                 :        205 :         deg11(R, p, bnr, pr);
    1453                 :            :       }
    1454                 :        220 :       break;
    1455                 :            :     }
    1456                 :            :   }
    1457                 :            :   /* precompute isprincipalray(x), x in Z */
    1458                 :        140 :   R->rayZ = cgetg(condZ, t_VEC);
    1459         [ +  + ]:       1395 :   for (i=1; i<condZ; i++)
    1460         [ +  + ]:       1255 :     gel(R->rayZ,i) = (ugcd(i,condZ) == 1)? isprincipalray(bnr, utoipos(i)): gen_0;
    1461                 :        140 :   gerepileall(av, 7, &(R->L0), &(R->L2), &(R->rayZ),
    1462                 :            :               &(R->L1), &(R->L1ray), &(R->L11), &(R->L11ray) );
    1463                 :        140 : }
    1464                 :            : 
    1465                 :            : static void
    1466                 :        345 : InitPrimes(GEN bnr, ulong N0, LISTray *R)
    1467                 :            : {
    1468                 :        345 :   GEN bnf = bnr_get_bnf(bnr), cond = gel(bnr_get_mod(bnr), 1);
    1469                 :        345 :   long p,j,k,l, condZ = itos(gcoeff(cond,1,1)), N = lg(cond)-1;
    1470                 :        345 :   GEN tmpray, tabpr, prime, nf = bnf_get_nf(bnf);
    1471                 :            :   forprime_t T;
    1472                 :            : 
    1473                 :        345 :   R->condZ = condZ; l = PiBound(N0) * N;
    1474                 :        345 :   tmpray = cgetg(N+1, t_VEC);
    1475                 :        345 :   R->L1 = vecsmalltrunc_init(l);
    1476                 :        345 :   R->L1ray = vectrunc_init(l);
    1477                 :        345 :   u_forprime_init(&T, 2, N0);
    1478                 :        345 :   prime = utoipos(2);
    1479         [ +  + ]:     123115 :   while ( (p = u_forprime_next(&T)) )
    1480                 :            :   {
    1481                 :     122770 :     pari_sp av = avma;
    1482                 :     122770 :     prime[2] = p;
    1483 [ -  + ][ #  # ]:     122770 :     if (DEBUGLEVEL>1 && (p & 2047) == 1) err_printf("%ld ", p);
    1484                 :     122770 :     tabpr = idealprimedec(nf, prime);
    1485         [ +  + ]:     244477 :     for (j = 1; j < lg(tabpr); j++)
    1486                 :            :     {
    1487                 :     215119 :       GEN pr  = gel(tabpr,j);
    1488                 :     215119 :       ulong np = upowuu(p, pr_get_f(pr));
    1489 [ +  + ][ +  + ]:     215119 :       if (!np || np > N0) break;
    1490 [ +  + ][ +  + ]:     121707 :       if (condZ % p == 0 && idealval(nf, cond, pr))
    1491                 :            :       {
    1492                 :        315 :         gel(tmpray,j) = NULL; continue;
    1493                 :            :       }
    1494                 :            : 
    1495                 :     121392 :       vecsmalltrunc_append(R->L1, np);
    1496                 :     121392 :       gel(tmpray,j) = gclone( isprincipalray(bnr, pr) );
    1497                 :            :     }
    1498                 :     122770 :     avma = av;
    1499         [ +  + ]:     244477 :     for (k = 1; k < j; k++)
    1500                 :            :     {
    1501         [ +  + ]:     121707 :       if (!tmpray[k]) continue;
    1502                 :     121392 :       vectrunc_append(R->L1ray, ZC_copy(gel(tmpray,k)));
    1503                 :     121392 :       gunclone(gel(tmpray,k));
    1504                 :            :     }
    1505                 :            :   }
    1506                 :        345 : }
    1507                 :            : 
    1508                 :            : static GEN /* cf polcoeff */
    1509                 :     206540 : _sercoeff(GEN x, long n)
    1510                 :            : {
    1511                 :     206540 :   long i = n - valp(x);
    1512         [ +  + ]:     206540 :   return (i < 0)? gen_0: gel(x,i+2);
    1513                 :            : }
    1514                 :            : 
    1515                 :            : static void
    1516                 :     206540 : affect_coeff(GEN q, long n, GEN y)
    1517                 :            : {
    1518                 :     206540 :   GEN x = _sercoeff(q,-n);
    1519         [ +  + ]:     206540 :   if (x == gen_0) gel(y,n) = NULL; else affgr(x, gel(y,n));
    1520                 :     206540 : }
    1521                 :            : 
    1522                 :            : typedef struct {
    1523                 :            :   GEN c1, aij, bij, cS, cT, powracpi;
    1524                 :            :   long i0, a,b,c, r, rc1, rc2;
    1525                 :            : } ST_t;
    1526                 :            : 
    1527                 :            : /* compute the principal part at the integers s = 0, -1, -2, ..., -i0
    1528                 :            :    of Gamma((s+1)/2)^a Gamma(s/2)^b Gamma(s)^c / (s - z) with z = 0 and 1 */
    1529                 :            : /* NOTE: surely not the best way to do this, but it's fast enough! */
    1530                 :            : static void
    1531                 :        160 : ppgamma(ST_t *T, long prec)
    1532                 :            : {
    1533                 :            :   GEN eul, gam,gamun,gamdm, an,bn,cn_evn,cn_odd, x,x2,X,Y, cf, sqpi;
    1534                 :            :   GEN p1, p2, aij, bij;
    1535                 :        160 :   long a = T->a;
    1536                 :        160 :   long b = T->b;
    1537                 :        160 :   long c = T->c, r = T->r, i0 = T->i0;
    1538                 :            :   long i,j, s,t;
    1539                 :            :   pari_sp av;
    1540                 :            : 
    1541                 :        160 :   aij = cgetg(i0+1, t_VEC);
    1542                 :        160 :   bij = cgetg(i0+1, t_VEC);
    1543         [ +  + ]:      51030 :   for (i = 1; i <= i0; i++)
    1544                 :            :   {
    1545                 :      50870 :     gel(aij,i) = p1 = cgetg(r+1, t_VEC);
    1546                 :      50870 :     gel(bij,i) = p2 = cgetg(r+1, t_VEC);
    1547         [ +  + ]:     154140 :     for (j=1; j<=r; j++) { gel(p1,j) = cgetr(prec); gel(p2,j) = cgetr(prec); }
    1548                 :            :   }
    1549                 :        160 :   av = avma;
    1550                 :            : 
    1551                 :        160 :   x   = pol_x(0);
    1552                 :        160 :   x2  = gmul2n(x, -1); /* x/2 */
    1553                 :        160 :   eul = mpeuler(prec);
    1554                 :        160 :   sqpi= sqrtr_abs(mppi(prec)); /* Gamma(1/2) */
    1555                 :            : 
    1556                 :            :   /* expansion of log(Gamma(u)) at u = 1 */
    1557                 :        160 :   gamun = cgetg(r+3, t_SER);
    1558                 :        160 :   gamun[1] = evalsigne(1) | _evalvalp(0) | evalvarn(0);
    1559                 :        160 :   gel(gamun,2) = gen_0;
    1560                 :        160 :   gel(gamun,3) = gneg(eul);
    1561         [ +  + ]:        330 :   for (i = 2; i <= r; i++)
    1562         [ +  + ]:        170 :     gel(gamun,i+2) = divrs(szeta(i,prec), odd(i)? -i: i);
    1563                 :        160 :   gamun = gexp(gamun, prec); /* Gamma(1 + x) */
    1564                 :        160 :   gam = gdiv(gamun,x); /* Gamma(x) */
    1565                 :            : 
    1566                 :            :   /* expansion of log(Gamma(u) / Gamma(1/2)) at u = 1/2 */
    1567                 :        160 :   gamdm = cgetg(r+3, t_SER);
    1568                 :        160 :   gamdm[1] = evalsigne(1) | _evalvalp(0) | evalvarn(0);
    1569                 :        160 :   gel(gamdm,2) = gen_0;
    1570                 :        160 :   gel(gamdm,3) = gneg(gadd(gmul2n(mplog2(prec), 1), eul));
    1571         [ +  + ]:        330 :   for (i = 2; i <= r; i++)
    1572                 :        170 :     gel(gamdm,i+2) = mulri(gel(gamun,i+2), subis(int2n(i), 1));
    1573                 :        160 :   gamdm = gmul(sqpi, gexp(gamdm, prec)); /* Gamma(1/2 + x) */
    1574                 :            : 
    1575                 :            :  /* We simplify to get one of the following two expressions
    1576                 :            :   * if (b > a) : sqrt{Pi}^a 2^{a-au} Gamma(u)^{a+c} Gamma(  u/2  )^{|b-a|}
    1577                 :            :   * if (b <= a): sqrt{Pi}^b 2^{b-bu} Gamma(u)^{b+c} Gamma((u+1)/2)^{|b-a|} */
    1578         [ +  + ]:        160 :   if (b > a)
    1579                 :            :   {
    1580                 :         20 :     t = a; s = b; X = x2; Y = gsub(x2,ghalf);
    1581                 :         20 :     p1 = ser_unscale(gam, ghalf);
    1582                 :         20 :     p2 = gdiv(ser_unscale(gamdm,ghalf), Y); /* Gamma((x-1)/2) */
    1583                 :            :   }
    1584                 :            :   else
    1585                 :            :   {
    1586                 :        140 :     t = b; s = a; X = gadd(x2,ghalf); Y = x2;
    1587                 :        140 :     p1 = ser_unscale(gamdm,ghalf);
    1588                 :        140 :     p2 = ser_unscale(gam,ghalf);
    1589                 :            :   }
    1590                 :        160 :   cf = powru(sqpi, t);
    1591                 :        160 :   an = gpowgs(gpow(gen_2, gsubsg(1,x), prec), t); /* 2^{t-tx} */
    1592                 :        160 :   bn = gpowgs(gam, t+c); /* Gamma(x)^{t+c} */
    1593                 :        160 :   cn_evn = gpowgs(p1, s-t); /* Gamma(X)^{s-t} */
    1594                 :        160 :   cn_odd = gpowgs(p2, s-t); /* Gamma(Y)^{s-t} */
    1595         [ +  + ]:      25595 :   for (i = 0; i < i0/2; i++)
    1596                 :            :   {
    1597                 :      25435 :     GEN C1,q1, A1 = gel(aij,2*i+1), B1 = gel(bij,2*i+1);
    1598                 :      25435 :     GEN C2,q2, A2 = gel(aij,2*i+2), B2 = gel(bij,2*i+2);
    1599                 :            : 
    1600                 :      25435 :     C1 = gmul(cf, gmul(bn, gmul(an, cn_evn)));
    1601                 :      25435 :     p1 = gdiv(C1, gsubgs(x, 2*i));
    1602                 :      25435 :     q1 = gdiv(C1, gsubgs(x, 2*i+1));
    1603                 :            : 
    1604                 :            :     /* an(x-u-1) = 2^t an(x-u) */
    1605                 :      25435 :     an = gmul2n(an, t);
    1606                 :            :     /* bn(x-u-1) = bn(x-u) / (x-u-1)^{t+c} */
    1607                 :      25435 :     bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+1), t+c));
    1608                 :            : 
    1609                 :      25435 :     C2 = gmul(cf, gmul(bn, gmul(an, cn_odd)));
    1610                 :      25435 :     p2 = gdiv(C2, gsubgs(x, 2*i+1));
    1611                 :      25435 :     q2 = gdiv(C2, gsubgs(x, 2*i+2));
    1612         [ +  + ]:      77070 :     for (j = 1; j <= r; j++)
    1613                 :            :     {
    1614                 :      51635 :       affect_coeff(p1, j, A1); affect_coeff(q1, j, B1);
    1615                 :      51635 :       affect_coeff(p2, j, A2); affect_coeff(q2, j, B2);
    1616                 :            :     }
    1617                 :            : 
    1618                 :      25435 :     an = gmul2n(an, t);
    1619                 :      25435 :     bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+2), t+c));
    1620                 :            :     /* cn_evn(x-2i-2) = cn_evn(x-2i)  / (X - (i+1))^{s-t} */
    1621                 :            :     /* cn_odd(x-2i-3) = cn_odd(x-2i-1)/ (Y - (i+1))^{s-t} */
    1622                 :      25435 :     cn_evn = gdiv(cn_evn, gpowgs(gsubgs(X,i+1), s-t));
    1623                 :      25435 :     cn_odd = gdiv(cn_odd, gpowgs(gsubgs(Y,i+1), s-t));
    1624                 :            :   }
    1625                 :        160 :   T->aij = aij;
    1626                 :        160 :   T->bij = bij; avma = av;
    1627                 :        160 : }
    1628                 :            : 
    1629                 :            : static GEN
    1630                 :       1565 : _cond(GEN dtcr) { return mkvec2(ch_cond(dtcr), ch_4(dtcr)); }
    1631                 :            : 
    1632                 :            : /* sort chars according to conductor */
    1633                 :            : static GEN
    1634                 :        480 : sortChars(GEN dataCR)
    1635                 :            : {
    1636                 :        480 :   const long cl = lg(dataCR) - 1;
    1637                 :        480 :   GEN vCond  = cgetg(cl+1, t_VEC);
    1638                 :        480 :   GEN CC     = cgetg(cl+1, t_VECSMALL);
    1639                 :        480 :   GEN nvCond = cgetg(cl+1, t_VECSMALL);
    1640                 :            :   long j,k, ncond;
    1641                 :            :   GEN vChar;
    1642                 :            : 
    1643         [ +  + ]:       2045 :   for (j = 1; j <= cl; j++) nvCond[j] = 0;
    1644                 :            : 
    1645                 :        480 :   ncond = 0;
    1646         [ +  + ]:       2045 :   for (j = 1; j <= cl; j++)
    1647                 :            :   {
    1648                 :       1565 :     GEN cond = _cond(gel(dataCR,j));
    1649         [ +  + ]:       1759 :     for (k = 1; k <= ncond; k++)
    1650         [ +  + ]:       1184 :       if (gequal(cond, gel(vCond,k))) break;
    1651         [ +  + ]:       1565 :     if (k > ncond) gel(vCond,++ncond) = cond;
    1652                 :       1565 :     nvCond[k]++; CC[j] = k; /* char j has conductor number k */
    1653                 :            :   }
    1654                 :        480 :   vChar = cgetg(ncond+1, t_VEC);
    1655         [ +  + ]:       1055 :   for (k = 1; k <= ncond; k++)
    1656                 :            :   {
    1657                 :        575 :     gel(vChar,k) = cgetg(nvCond[k]+1, t_VECSMALL);
    1658                 :        575 :     nvCond[k] = 0;
    1659                 :            :   }
    1660         [ +  + ]:       2045 :   for (j = 1; j <= cl; j++)
    1661                 :            :   {
    1662                 :       1565 :     k = CC[j]; nvCond[k]++;
    1663                 :       1565 :     mael(vChar,k,nvCond[k]) = j;
    1664                 :            :   }
    1665                 :        480 :   return vChar;
    1666                 :            : }
    1667                 :            : 
    1668                 :            : /* Given W(chi), S(chi) and T(chi), return L(1, chi) if fl & 1, else
    1669                 :            :    [r(chi), c(chi)] where L(s, chi) ~ c(chi) s^r(chi) at s = 0.
    1670                 :            :    If fl & 2, adjust the value to get L_S(s, chi). */
    1671                 :            : static GEN
    1672                 :        845 : GetValue(GEN dtcr, GEN W, GEN S, GEN T, long fl, long prec)
    1673                 :            : {
    1674                 :        845 :   pari_sp av = avma;
    1675                 :            :   GEN cf, z, p1;
    1676                 :            :   long q, b, c, r;
    1677                 :        845 :   int isreal = (itos(gel(ch_CHI0(dtcr), 3)) <= 2);
    1678                 :            : 
    1679                 :        845 :   p1 = ch_4(dtcr);
    1680                 :        845 :   q = p1[1];
    1681                 :        845 :   b = p1[2];
    1682                 :        845 :   c = p1[3];
    1683                 :            : 
    1684         [ -  + ]:        845 :   if (fl & 1)
    1685                 :            :   { /* S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */
    1686                 :          0 :     cf = gmul(ch_C(dtcr), powruhalf(mppi(prec), b));
    1687                 :            : 
    1688                 :          0 :     z = gadd(S, gmul(W, T));
    1689         [ #  # ]:          0 :     if (isreal) z = real_i(z);
    1690                 :          0 :     z = gdiv(z, cf);
    1691         [ #  # ]:          0 :     if (fl & 2) z = gmul(z, ComputeAChi(dtcr, &r, 1, prec));
    1692                 :            :   }
    1693                 :            :   else
    1694                 :            :   { /* (W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */
    1695                 :        845 :     cf = gmul2n(powruhalf(mppi(prec), q), b);
    1696                 :            : 
    1697                 :        845 :     z = gadd(gmul(W, gconj(S)), gconj(T));
    1698         [ +  + ]:        845 :     if (isreal) z = real_i(z);
    1699                 :        845 :     z = gdiv(z, cf); r = 0;
    1700         [ +  + ]:        845 :     if (fl & 2) z = gmul(z, ComputeAChi(dtcr, &r, 0, prec));
    1701                 :        845 :     z = mkvec2(utoi(b + c + r), z);
    1702                 :            :   }
    1703                 :        845 :   return gerepilecopy(av, z);
    1704                 :            : }
    1705                 :            : 
    1706                 :            : /* return the order and the first non-zero term of L(s, chi0)
    1707                 :            :    at s = 0. If flag != 0, adjust the value to get L_S(s, chi0). */
    1708                 :            : static GEN
    1709                 :         30 : GetValue1(GEN bnr, long flag, long prec)
    1710                 :            : {
    1711                 :         30 :   GEN bnf = checkbnf(bnr), nf = bnf_get_nf(bnf);
    1712                 :            :   GEN h, R, c, diff;
    1713                 :            :   long i, l, r, r1, r2;
    1714                 :         30 :   pari_sp av = avma;
    1715                 :            : 
    1716                 :         30 :   nf_get_sign(nf, &r1,&r2);
    1717                 :         30 :   h = bnf_get_no(bnf);
    1718                 :         30 :   R = bnf_get_reg(bnf);
    1719                 :            : 
    1720                 :         30 :   c = gneg_i(gdivgs(mpmul(h, R), bnf_get_tuN(bnf)));
    1721                 :         30 :   r = r1 + r2 - 1;
    1722                 :            : 
    1723         [ -  + ]:         30 :   if (flag)
    1724                 :            :   {
    1725                 :          0 :     diff = divcond(bnr);
    1726                 :          0 :     l = lg(diff) - 1; r += l;
    1727         [ #  # ]:          0 :     for (i = 1; i <= l; i++)
    1728                 :          0 :       c = gmul(c, glog(pr_norm(gel(diff,i)), prec));
    1729                 :            :   }
    1730                 :         30 :   return gerepilecopy(av, mkvec2(stoi(r), c));
    1731                 :            : }
    1732                 :            : 
    1733                 :            : /********************************************************************/
    1734                 :            : /*                6th part: recover the coefficients                */
    1735                 :            : /********************************************************************/
    1736                 :            : static long
    1737                 :       1776 : TestOne(GEN plg, RC_data *d)
    1738                 :            : {
    1739                 :       1776 :   long j, v = d->v;
    1740                 :       1776 :   GEN z = gsub(d->beta, gel(plg,v));
    1741         [ +  + ]:       1776 :   if (expo(z) >= d->G) return 0;
    1742         [ +  + ]:       4659 :   for (j = 1; j < lg(plg); j++)
    1743 [ +  + ][ +  + ]:       3246 :     if (j != v && mpcmp(d->B, mpabs(gel(plg,j))) < 0) return 0;
    1744                 :       1776 :   return 1;
    1745                 :            : }
    1746                 :            : 
    1747                 :            : static GEN
    1748                 :        303 : chk_reccoeff_init(FP_chk_fun *chk, GEN r, GEN mat)
    1749                 :            : {
    1750                 :        303 :   RC_data *d = (RC_data*)chk->data;
    1751                 :        303 :   (void)r; d->U = mat; return d->nB;
    1752                 :            : }
    1753                 :            : 
    1754                 :            : static GEN
    1755                 :        266 : chk_reccoeff(void *data, GEN x)
    1756                 :            : {
    1757                 :        266 :   RC_data *d = (RC_data*)data;
    1758                 :        266 :   GEN v = gmul(d->U, x), z = gel(v,1);
    1759                 :            : 
    1760         [ -  + ]:        266 :   if (!gequal1(z)) return NULL;
    1761                 :        266 :   *++v = evaltyp(t_COL) | evallg( lg(d->M) );
    1762         [ +  - ]:        266 :   if (TestOne(gmul(d->M, v), d)) return v;
    1763                 :        266 :   return NULL;
    1764                 :            : }
    1765                 :            : 
    1766                 :            : /* Using Cohen's method */
    1767                 :            : static GEN
    1768                 :        298 : RecCoeff3(GEN nf, RC_data *d, long prec)
    1769                 :            : {
    1770                 :            :   GEN A, M, nB, cand, p1, B2, C2, tB, beta2, nf2, Bd;
    1771                 :        298 :   GEN beta = d->beta, B = d->B;
    1772                 :        298 :   long N = d->N, v = d->v, e, BIG;
    1773                 :        298 :   long i, j, k, l, ct = 0, prec2;
    1774                 :        298 :   FP_chk_fun chk = { &chk_reccoeff, &chk_reccoeff_init, NULL, NULL, 0 };
    1775                 :        298 :   chk.data = (void*)d;
    1776                 :            : 
    1777                 :        298 :   d->G = minss(-10, -prec2nbits(prec) >> 4);
    1778                 :        298 :   BIG = maxss(32, -(d->G << 1));
    1779                 :        298 :   tB  = sqrtnr(real2n(BIG-N,DEFAULTPREC), N-1);
    1780                 :        298 :   Bd  = grndtoi(gmin(B, tB), &e);
    1781         [ -  + ]:        298 :   if (e > 0) return NULL; /* failure */
    1782                 :        298 :   Bd = addis(Bd, 1);
    1783                 :        298 :   prec2 = nbits2prec( expi(Bd) + 192 );
    1784                 :        298 :   prec2 = maxss(precdbl(prec), prec2);
    1785                 :        298 :   B2 = sqri(Bd);
    1786                 :        298 :   C2 = shifti(B2, BIG<<1);
    1787                 :            : 
    1788                 :        303 : LABrcf: ct++;
    1789                 :        303 :   beta2 = gprec_w(beta, prec2);
    1790                 :        303 :   nf2 = nfnewprec_shallow(nf, prec2);
    1791                 :        303 :   d->M = M = nf_get_M(nf2);
    1792                 :            : 
    1793                 :        303 :   A = cgetg(N+2, t_MAT);
    1794         [ +  + ]:       1252 :   for (i = 1; i <= N+1; i++) gel(A,i) = cgetg(N+2, t_COL);
    1795                 :            : 
    1796                 :        303 :   gcoeff(A, 1, 1) = gadd(gmul(C2, gsqr(beta2)), B2);
    1797         [ +  + ]:        949 :   for (j = 2; j <= N+1; j++)
    1798                 :            :   {
    1799                 :        646 :     p1 = gmul(C2, gmul(gneg_i(beta2), gcoeff(M, v, j-1)));
    1800                 :        646 :     gcoeff(A, 1, j) = gcoeff(A, j, 1) = p1;
    1801                 :            :   }
    1802         [ +  + ]:        949 :   for (i = 2; i <= N+1; i++)
    1803         [ +  + ]:       1675 :     for (j = i; j <= N+1; j++)
    1804                 :            :     {
    1805                 :       1029 :       p1 = gen_0;
    1806         [ +  + ]:       3327 :       for (k = 1; k <= N; k++)
    1807                 :            :       {
    1808                 :       2298 :         GEN p2 = gmul(gcoeff(M, k, j-1), gcoeff(M, k, i-1));
    1809         [ +  + ]:       2298 :         if (k == v) p2 = gmul(C2, p2);
    1810                 :       2298 :         p1 = gadd(p1,p2);
    1811                 :            :       }
    1812                 :       1029 :       gcoeff(A, i, j) = gcoeff(A, j, i) = p1;
    1813                 :            :     }
    1814                 :            : 
    1815                 :        303 :   nB = mului(N+1, B2);
    1816                 :        303 :   d->nB = nB;
    1817                 :        303 :   cand = fincke_pohst(A, nB, -1, prec2, &chk);
    1818                 :            : 
    1819         [ +  + ]:        303 :   if (!cand)
    1820                 :            :   {
    1821         [ -  + ]:          5 :     if (ct > 3) return NULL;
    1822                 :          5 :     prec2 = precdbl(prec2);
    1823         [ -  + ]:          5 :     if (DEBUGLEVEL>1) pari_warn(warnprec,"RecCoeff", prec2);
    1824                 :          5 :     goto LABrcf;
    1825                 :            :   }
    1826                 :            : 
    1827                 :        298 :   cand = gel(cand,1); l = lg(cand) - 1;
    1828         [ +  + ]:        298 :   if (l == 1)
    1829                 :        133 :     return coltoalg(nf, gel(cand,1));
    1830                 :            : 
    1831         [ -  + ]:        165 :   if (DEBUGLEVEL>1) err_printf("RecCoeff3: no solution found!\n");
    1832                 :        298 :   return NULL;
    1833                 :            : }
    1834                 :            : 
    1835                 :            : /* Using linear dependance relations */
    1836                 :            : static GEN
    1837                 :       1445 : RecCoeff2(GEN nf,  RC_data *d,  long prec)
    1838                 :            : {
    1839                 :            :   pari_sp av;
    1840                 :       1445 :   GEN vec, M = nf_get_M(nf), beta = d->beta;
    1841                 :       1445 :   long i, imin, imax, lM = lg(M);
    1842                 :            : 
    1843                 :       1445 :   d->G = minss(-20, -prec2nbits(prec) >> 4);
    1844                 :            : 
    1845                 :       1445 :   vec  = shallowconcat(mkvec(gneg(beta)), row(M, d->v));
    1846                 :       1445 :   imin = (long)prec2nbits_mul(prec, .225);
    1847                 :       1445 :   imax = (long)prec2nbits_mul(prec, .315);
    1848                 :            : 
    1849                 :       1445 :   av = avma;
    1850         [ +  + ]:       2006 :   for (i = imax; i >= imin; i-=16, avma = av)
    1851                 :            :   {
    1852                 :            :     long e;
    1853                 :       1708 :     GEN v = lindep2(vec, i), z = gel(v,1);
    1854         [ +  + ]:       1708 :     if (!signe(z)) continue;
    1855                 :       1510 :     *++v = evaltyp(t_COL) | evallg(lM);
    1856                 :       1510 :     v = grndtoi(gdiv(v, z), &e);
    1857         [ +  - ]:       1510 :     if (e > 0) break;
    1858         [ +  + ]:       1510 :     if (TestOne(gmul(M, v), d)) return coltoalg(nf, v);
    1859                 :            :   }
    1860                 :            :   /* failure */
    1861                 :       1445 :   return RecCoeff3(nf,d,prec);
    1862                 :            : }
    1863                 :            : 
    1864                 :            : /* Attempts to find a polynomial with coefficients in nf such that
    1865                 :            :    its coefficients are close to those of pol at the place v and
    1866                 :            :    less than B at all the other places */
    1867                 :            : static GEN
    1868                 :        390 : RecCoeff(GEN nf,  GEN pol,  long v, long prec)
    1869                 :            : {
    1870                 :        390 :   long j, md, cl = degpol(pol);
    1871                 :        390 :   pari_sp av = avma;
    1872                 :            :   RC_data d;
    1873                 :            : 
    1874                 :            :   /* if precision(pol) is too low, abort */
    1875         [ +  + ]:       2595 :   for (j = 2; j <= cl+1; j++)
    1876                 :            :   {
    1877                 :       2205 :     GEN t = gel(pol, j);
    1878         [ -  + ]:       2205 :     if (prec2nbits(gprecision(t)) - gexpo(t) < 34) return NULL;
    1879                 :            :   }
    1880                 :            : 
    1881                 :        390 :   md = cl/2;
    1882                 :        390 :   pol = leafcopy(pol);
    1883                 :            : 
    1884                 :        390 :   d.N = nf_get_degree(nf);
    1885                 :        390 :   d.v = v;
    1886                 :            : 
    1887         [ +  + ]:       1670 :   for (j = 1; j <= cl; j++)
    1888                 :            :   { /* start with the coefficients in the middle,
    1889                 :            :        since they are the harder to recognize! */
    1890         [ +  + ]:       1445 :     long cf = md + (j%2? j/2: -j/2);
    1891                 :       1445 :     GEN t, bound = shifti(binomial(utoipos(cl), cf), cl-cf);
    1892                 :            : 
    1893         [ -  + ]:       1445 :     if (DEBUGLEVEL>1) err_printf("RecCoeff (cf = %ld, B = %Ps)\n", cf, bound);
    1894                 :       1445 :     d.beta = real_i( gel(pol,cf+2) );
    1895                 :       1445 :     d.B    = bound;
    1896         [ +  + ]:       1445 :     if (! (t = RecCoeff2(nf, &d, prec)) ) return NULL;
    1897                 :       1280 :     gel(pol, cf+2) = t;
    1898                 :            :   }
    1899                 :        225 :   gel(pol,cl+2) = gen_1;
    1900                 :        390 :   return gerepilecopy(av, pol);
    1901                 :            : }
    1902                 :            : 
    1903                 :            : /* an[q * i] *= chi for all (i,p)=1 */
    1904                 :            : static void
    1905                 :     102494 : an_mul(int **an, long p, long q, long n, long deg, GEN chi, int **reduc)
    1906                 :            : {
    1907                 :            :   pari_sp av;
    1908                 :            :   long c,i;
    1909                 :            :   int *T;
    1910                 :            : 
    1911         [ +  + ]:     196505 :   if (gequal1(chi)) return;
    1912                 :      94011 :   av = avma;
    1913                 :      94011 :   T = (int*)new_chunk(deg); Polmod2Coeff(T,chi, deg);
    1914         [ +  + ]:    2033383 :   for (c = 1, i = q; i <= n; i += q, c++)
    1915         [ +  + ]:    1939372 :     if (c == p) c = 0; else MulCoeff(an[i], T, reduc, deg);
    1916                 :      94011 :   avma = av;
    1917                 :            : }
    1918                 :            : /* an[q * i] = 0 for all (i,p)=1 */
    1919                 :            : static void
    1920                 :      91513 : an_set0_coprime(int **an, long p, long q, long n, long deg)
    1921                 :            : {
    1922                 :            :   long c,i;
    1923         [ +  + ]:    1063531 :   for (c = 1, i = q; i <= n; i += q, c++)
    1924         [ +  + ]:     972018 :     if (c == p) c = 0; else _0toCoeff(an[i], deg);
    1925                 :      91513 : }
    1926                 :            : /* an[q * i] = 0 for all i */
    1927                 :            : static void
    1928                 :        105 : an_set0(int **an, long p, long n, long deg)
    1929                 :            : {
    1930                 :            :   long i;
    1931         [ +  + ]:     123190 :   for (i = p; i <= n; i += p) _0toCoeff(an[i], deg);
    1932                 :        105 : }
    1933                 :            : 
    1934                 :            : /* compute the coefficients an for the quadratic case */
    1935                 :            : static int**
    1936                 :        520 : computean(GEN dtcr, LISTray *R, long n, long deg)
    1937                 :            : {
    1938                 :        520 :   pari_sp av = avma, av2;
    1939                 :            :   long i, p, q, condZ, l;
    1940                 :            :   int **an, **reduc;
    1941                 :            :   GEN L, CHI, chi, chi1;
    1942                 :            :   CHI_t C;
    1943                 :            : 
    1944                 :        520 :   CHI = ch_CHI(dtcr); init_CHI_alg(&C, CHI);
    1945                 :        520 :   condZ= R->condZ;
    1946                 :            : 
    1947                 :        520 :   an = InitMatAn(n, deg, 1);
    1948                 :        520 :   reduc = InitReduction(CHI, deg);
    1949                 :        520 :   av2 = avma;
    1950                 :            : 
    1951                 :            :   /* all pr | p divide cond */
    1952                 :        520 :   L = R->L0; l = lg(L);
    1953         [ +  + ]:        625 :   for (i=1; i<l; i++) an_set0(an,L[i],n,deg);
    1954                 :            : 
    1955                 :            :   /* 1 prime of degree 2 */
    1956                 :        520 :   L = R->L2; l = lg(L);
    1957         [ +  + ]:      91011 :   for (i=1; i<l; i++, avma = av2)
    1958                 :            :   {
    1959                 :      90491 :     p = L[i];
    1960         [ +  + ]:      90491 :     if (condZ == 1) chi = C.val[0]; /* 1 */
    1961                 :      90381 :     else            chi = EvalChar(&C, gel(R->rayZ, p%condZ));
    1962                 :      90491 :     chi1 = chi;
    1963                 :      90491 :     for (q=p;;)
    1964                 :            :     {
    1965                 :      91513 :       an_set0_coprime(an, p,q,n,deg); /* v_p(q) odd */
    1966         [ +  + ]:      91513 :       if (! (q = next_pow(q,p, n)) ) break;
    1967                 :            : 
    1968                 :       3733 :       an_mul(an,p,q,n,deg,chi,reduc);
    1969         [ +  + ]:       3733 :       if (! (q = next_pow(q,p, n)) ) break;
    1970                 :       1022 :       chi = gmul(chi, chi1);
    1971                 :       1022 :     }
    1972                 :            :   }
    1973                 :            : 
    1974                 :            :   /* 1 prime of degree 1 */
    1975                 :        520 :   L = R->L1; l = lg(L);
    1976         [ +  + ]:       1835 :   for (i=1; i<l; i++, avma = av2)
    1977                 :            :   {
    1978                 :       1315 :     p = L[i];
    1979                 :       1315 :     chi = EvalChar(&C, gel(R->L1ray,i));
    1980                 :       1315 :     chi1 = chi;
    1981                 :       1315 :     for(q=p;;)
    1982                 :            :     {
    1983                 :       5765 :       an_mul(an,p,q,n,deg,chi,reduc);
    1984         [ +  + ]:       5765 :       if (! (q = next_pow(q,p, n)) ) break;
    1985                 :       4450 :       chi = gmul(chi, chi1);
    1986                 :       4450 :     }
    1987                 :            :   }
    1988                 :            : 
    1989                 :            :   /* 2 primes of degree 1 */
    1990                 :        520 :   L = R->L11; l = lg(L);
    1991         [ +  + ]:      87562 :   for (i=1; i<l; i++, avma = av2)
    1992                 :            :   {
    1993                 :            :     GEN ray1, ray2, chi11, chi12, chi2;
    1994                 :            : 
    1995                 :      87042 :     p = L[i]; ray1 = gel(R->L11ray,i); /* use pr1 pr2 = (p) */
    1996         [ +  + ]:      87042 :     if (condZ == 1)
    1997                 :         80 :       ray2 = ZC_neg(ray1);
    1998                 :            :     else
    1999                 :      86962 :       ray2 = ZC_sub(gel(R->rayZ, p%condZ),  ray1);
    2000                 :      87042 :     chi11 = EvalChar(&C, ray1);
    2001                 :      87042 :     chi12 = EvalChar(&C, ray2);
    2002                 :            : 
    2003                 :      87042 :     chi1 = gadd(chi11, chi12);
    2004                 :      87042 :     chi2 = chi12;
    2005                 :      87042 :     for(q=p;;)
    2006                 :            :     {
    2007                 :      92996 :       an_mul(an,p,q,n,deg,chi1,reduc);
    2008         [ +  + ]:      92996 :       if (! (q = next_pow(q,p, n)) ) break;
    2009                 :       5954 :       chi2 = gmul(chi2, chi12);
    2010                 :       5954 :       chi1 = gadd(chi2, gmul(chi1, chi11));
    2011                 :       5954 :     }
    2012                 :            :   }
    2013                 :            : 
    2014                 :        520 :   CorrectCoeff(dtcr, an, reduc, n, deg);
    2015                 :        520 :   FreeMat(reduc, deg-1);
    2016                 :        520 :   avma = av; return an;
    2017                 :            : }
    2018                 :            : 
    2019                 :            : /* return the vector of A^i/i for i = 1...n */
    2020                 :            : static GEN
    2021                 :        160 : mpvecpowdiv(GEN A, long n)
    2022                 :            : {
    2023                 :        160 :   pari_sp av = avma;
    2024                 :            :   long i;
    2025                 :        160 :   GEN v = powersr(A, n);
    2026                 :        160 :   GEN w = cgetg(n+1, t_VEC);
    2027                 :        160 :   gel(w,1) = rcopy(gel(v,2));
    2028         [ +  + ]:     314747 :   for (i=2; i<=n; i++) gel(w,i) = divru(gel(v,i+1), i);
    2029                 :        160 :   return gerepileupto(av, w);
    2030                 :            : }
    2031                 :            : 
    2032                 :            : static void GetST0(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec);
    2033                 :            : 
    2034                 :            : /* compute S and T for the quadratic case. The following cases (cs) are:
    2035                 :            :    1) bnr complex;
    2036                 :            :    2) bnr real and no infinite place divide cond_chi (TBD);
    2037                 :            :    3) bnr real and one infinite place divide cond_chi;
    2038                 :            :    4) bnr real and both infinite places divide cond_chi (TBD) */
    2039                 :            : static void
    2040                 :        155 : QuadGetST(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec)
    2041                 :            : {
    2042                 :        155 :   pari_sp av = avma, av1, av2;
    2043                 :            :   long ncond, n, j, k, n0;
    2044                 :        155 :   GEN N0, C, T = *pT, S = *pS, an, degs, cs;
    2045                 :            :   LISTray LIST;
    2046                 :            : 
    2047                 :            :   /* initializations */
    2048                 :        155 :   degs = GetDeg(dataCR);
    2049                 :        155 :   ncond = lg(vChar)-1;
    2050                 :        155 :   C    = cgetg(ncond+1, t_VEC);
    2051                 :        155 :   N0   = cgetg(ncond+1, t_VECSMALL);
    2052                 :        155 :   cs   = cgetg(ncond+1, t_VECSMALL);
    2053                 :        155 :   n0 = 0;
    2054         [ +  + ]:        325 :   for (j = 1; j <= ncond; j++)
    2055                 :            :   {
    2056                 :            :     /* FIXME: make sure that this value of c is correct for the general case */
    2057                 :            :     long r1, r2, q;
    2058                 :        185 :     GEN dtcr = gel(dataCR, mael(vChar,j,1)), p1 = ch_4(dtcr), c = ch_C(dtcr);
    2059                 :            : 
    2060                 :        185 :     gel(C,j) = c;
    2061                 :        185 :     q = p1[1];
    2062                 :            : 
    2063                 :        185 :     nf_get_sign(bnr_get_nf(gel(dtcr, 3)), &r1, &r2);
    2064         [ +  + ]:        185 :     if (r1 == 2) /* real quadratic */
    2065                 :            :     {
    2066                 :        175 :       cs[j] = 2 + q;
    2067                 :            :       /* FIXME:
    2068                 :            :          make sure that this value of N0 is correct for the general case */
    2069                 :        175 :       N0[j] = (long)prec2nbits_mul(prec, 0.35 * gtodouble(c));
    2070 [ +  + ][ -  + ]:        175 :       if (cs[j] == 2 || cs[j] == 4) /* NOT IMPLEMENTED YET */
    2071                 :            :       {
    2072                 :         15 :         GetST0(bnr, pS, pT, dataCR, vChar, prec);
    2073                 :        155 :         return;
    2074                 :            :       }
    2075                 :            :     }
    2076                 :            :     else /* complex quadratic */
    2077                 :            :     {
    2078                 :         10 :       cs[j] = 1;
    2079                 :         10 :       N0[j] = (long)prec2nbits_mul(prec, 0.7 * gtodouble(c));
    2080                 :            :     }
    2081         [ +  + ]:        170 :     if (n0 < N0[j]) n0 = N0[j];
    2082                 :            :   }
    2083         [ -  + ]:        140 :   if (DEBUGLEVEL>1) err_printf("N0 = %ld\n", n0);
    2084                 :        140 :   InitPrimesQuad(bnr, n0, &LIST);
    2085                 :            : 
    2086                 :        140 :   av1 = avma;
    2087                 :            :   /* loop over conductors */
    2088         [ +  + ]:        300 :   for (j = 1; j <= ncond; j++)
    2089                 :            :   {
    2090                 :        160 :     GEN c0 = gel(C,j), c1 = divur(1, c0), c2 = divur(2, c0);
    2091                 :        160 :     GEN ec1 = mpexp(c1), ec2 = mpexp(c2), LChar = gel(vChar,j);
    2092                 :            :     GEN vf0, vf1, cf0, cf1;
    2093                 :        160 :     const long nChar = lg(LChar)-1, NN = N0[j];
    2094                 :            : 
    2095         [ -  + ]:        160 :     if (DEBUGLEVEL>1)
    2096                 :          0 :       err_printf("* conductor no %ld/%ld (N = %ld)\n\tInit: ", j,ncond,NN);
    2097         [ +  - ]:        160 :     if (realprec(ec1) > prec) ec1 = rtor(ec1, prec);
    2098         [ +  - ]:        160 :     if (realprec(ec2) > prec) ec2 = rtor(ec2, prec);
    2099      [ +  +  - ]:        160 :     switch(cs[j])
    2100                 :            :     {
    2101                 :            :     case 1:
    2102                 :         10 :       cf0 = gen_1;
    2103                 :         10 :       cf1 = c0;
    2104                 :         10 :       vf0 = mpveceint1(rtor(c1, prec), ec1, NN);
    2105                 :         10 :       vf1 = mpvecpowdiv(invr(ec1), NN); break;
    2106                 :            : 
    2107                 :            :     case 3:
    2108                 :        150 :       cf0 = sqrtr(mppi(prec));
    2109                 :        150 :       cf1 = gmul2n(cf0, 1);
    2110                 :        150 :       cf0 = gmul(cf0, c0);
    2111                 :        150 :       vf0 = mpvecpowdiv(invr(ec2), NN);
    2112                 :        150 :       vf1 = mpveceint1(rtor(c2, prec), ec2, NN); break;
    2113                 :            : 
    2114                 :            :     default:
    2115                 :          0 :       cf0 = cf1 = NULL; /* FIXME: not implemented */
    2116                 :          0 :       vf0 = vf1 = NULL;
    2117                 :            :     }
    2118         [ +  + ]:        685 :     for (k = 1; k <= nChar; k++)
    2119                 :            :     {
    2120                 :        525 :       const long t = LChar[k], d = degs[t];
    2121                 :        525 :       const GEN dtcr = gel(dataCR, t), z = gel(ch_CHI(dtcr), 2);
    2122                 :        525 :       GEN p1 = gen_0, p2 = gen_0;
    2123                 :            :       int **matan;
    2124                 :        525 :       long c = 0;
    2125                 :            : 
    2126         [ -  + ]:        525 :       if (DEBUGLEVEL>1)
    2127                 :          0 :         err_printf("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar);
    2128         [ +  + ]:        525 :       if (isintzero( ch_comp(gel(dataCR, t)) ))
    2129                 :            :       {
    2130         [ -  + ]:          5 :         if (DEBUGLEVEL>1) err_printf("\t  no need to compute this character\n");
    2131                 :          5 :         continue;
    2132                 :            :       }
    2133                 :        520 :       av2 = avma;
    2134                 :        520 :       matan = computean(gel(dataCR,t), &LIST, NN, d);
    2135         [ +  + ]:    1221437 :       for (n = 1; n <= NN; n++)
    2136         [ +  + ]:    1220917 :         if ((an = EvalCoeff(z, matan[n], d)))
    2137                 :            :         {
    2138                 :     309367 :           p1 = gadd(p1, gmul(an, gel(vf0,n)));
    2139                 :     309367 :           p2 = gadd(p2, gmul(an, gel(vf1,n)));
    2140         [ +  + ]:     309367 :           if (++c == 256) { gerepileall(av2,2, &p1,&p2); c = 0; }
    2141                 :            :         }
    2142                 :        520 :       gaffect(gmul(cf0, p1), gel(S,t));
    2143                 :        520 :       gaffect(gmul(cf1,  gconj(p2)), gel(T,t));
    2144                 :        520 :       FreeMat(matan,NN); avma = av2;
    2145                 :            :     }
    2146         [ -  + ]:        160 :     if (DEBUGLEVEL>1) err_printf("\n");
    2147                 :        160 :     avma = av1;
    2148                 :            :   }
    2149                 :        140 :   avma = av;
    2150                 :            : }
    2151                 :            : 
    2152                 :            : /* s += t*u. All 3 of them t_REAL, except we allow s or u = NULL (for 0) */
    2153                 :            : static GEN
    2154                 :   85404734 : _addmulrr(GEN s, GEN t, GEN u)
    2155                 :            : {
    2156         [ +  + ]:   85404734 :   if (u)
    2157                 :            :   {
    2158                 :   85111871 :     GEN v = mulrr(t, u);
    2159         [ +  + ]:   85111871 :     return s? addrr(s, v): v;
    2160                 :            :   }
    2161                 :   85404734 :   return s;
    2162                 :            : }
    2163                 :            : /* s += t. Both real, except we allow s or t = NULL (for exact 0) */
    2164                 :            : static GEN
    2165                 :  170392664 : _addrr(GEN s, GEN t)
    2166 [ +  + ][ +  + ]:  170392664 : { return t? (s? addrr(s, t): t) : s; }
    2167                 :            : 
    2168                 :            : /* S & T for the general case. This is time-critical: optimize */
    2169                 :            : static void
    2170                 :     604741 : get_cS_cT(ST_t *T, long n)
    2171                 :            : {
    2172                 :            :   pari_sp av;
    2173                 :            :   GEN csurn, nsurc, lncsurn, A, B, s, t, Z, aij, bij;
    2174                 :            :   long i, j, r, i0;
    2175                 :            : 
    2176         [ +  + ]:     890009 :   if (T->cS[n]) return;
    2177                 :            : 
    2178                 :     285268 :   av = avma;
    2179                 :     285268 :   aij = T->aij; i0= T->i0;
    2180                 :     285268 :   bij = T->bij; r = T->r;
    2181                 :     285268 :   Z = cgetg(r+1, t_VEC);
    2182                 :     285268 :   gel(Z,1) = NULL; /* unused */
    2183                 :            : 
    2184                 :     285268 :   csurn = divru(T->c1, n);
    2185                 :     285268 :   nsurc = invr(csurn);
    2186                 :     285268 :   lncsurn = logr_abs(csurn);
    2187                 :            : 
    2188         [ +  + ]:     285268 :   if (r > 1)
    2189                 :            :   {
    2190                 :     285223 :     gel(Z,2) = lncsurn; /* r >= 2 */
    2191         [ +  + ]:     286678 :     for (i = 3; i <= r; i++)
    2192                 :       1455 :       gel(Z,i) = divru(mulrr(gel(Z,i-1), lncsurn), i-1);
    2193                 :            :     /* Z[i] = ln^(i-1)(c1/n) / (i-1)! */
    2194                 :            :   }
    2195                 :            : 
    2196                 :            :   /* i = i0 */
    2197                 :     285268 :     A = gel(aij,i0); t = _addrr(NULL, gel(A,1));
    2198                 :     285268 :     B = gel(bij,i0); s = _addrr(NULL, gel(B,1));
    2199         [ +  + ]:     571946 :     for (j = 2; j <= r; j++)
    2200                 :            :     {
    2201                 :     286678 :       s = _addmulrr(s, gel(Z,j),gel(B,j));
    2202                 :     286678 :       t = _addmulrr(t, gel(Z,j),gel(A,j));
    2203                 :            :     }
    2204         [ +  + ]:   84768430 :   for (i = i0 - 1; i > 1; i--)
    2205                 :            :   {
    2206         [ +  + ]:   84483162 :     A = gel(aij,i); if (t) t = mulrr(t, nsurc);
    2207         [ +  + ]:   84483162 :     B = gel(bij,i); if (s) s = mulrr(s, nsurc);
    2208 [ +  + ][ +  + ]:  126612173 :     for (j = odd(i)? T->rc2: T->rc1; j > 1; j--)
    2209                 :            :     {
    2210                 :   42129011 :       s = _addmulrr(s, gel(Z,j),gel(B,j));
    2211                 :   42129011 :       t = _addmulrr(t, gel(Z,j),gel(A,j));
    2212                 :            :     }
    2213                 :   84483162 :     s = _addrr(s, gel(B,1));
    2214                 :   84483162 :     t = _addrr(t, gel(A,1));
    2215                 :            :   }
    2216                 :            :   /* i = 1 */
    2217         [ +  - ]:     285268 :     A = gel(aij,1); if (t) t = mulrr(t, nsurc);
    2218         [ +  - ]:     285268 :     B = gel(bij,1); if (s) s = mulrr(s, nsurc);
    2219                 :     285268 :     s = _addrr(s, gel(B,1));
    2220                 :     285268 :     t = _addrr(t, gel(A,1));
    2221         [ +  + ]:     571946 :     for (j = 2; j <= r; j++)
    2222                 :            :     {
    2223                 :     286678 :       s = _addmulrr(s, gel(Z,j),gel(B,j));
    2224                 :     286678 :       t = _addmulrr(t, gel(Z,j),gel(A,j));
    2225                 :            :     }
    2226         [ +  + ]:     285268 :   s = _addrr(s, T->b? mulrr(csurn, gel(T->powracpi,T->b)): csurn);
    2227         [ -  + ]:     285268 :   if (!s) s = gen_0;
    2228         [ -  + ]:     285268 :   if (!t) t = gen_0;
    2229                 :     285268 :   gel(T->cS,n) = gclone(s);
    2230                 :     285268 :   gel(T->cT,n) = gclone(t); avma = av;
    2231                 :            : }
    2232                 :            : 
    2233                 :            : static void
    2234                 :        280 : clear_cScT(ST_t *T, long N)
    2235                 :            : {
    2236                 :        280 :   GEN cS = T->cS, cT = T->cT;
    2237                 :            :   long i;
    2238         [ +  + ]:    1813145 :   for (i=1; i<=N; i++)
    2239         [ +  + ]:    1812865 :     if (cS[i]) {
    2240                 :     285268 :       gunclone(gel(cS,i));
    2241                 :     285268 :       gunclone(gel(cT,i)); gel(cS,i) = gel(cT,i) = NULL;
    2242                 :            :     }
    2243                 :        280 : }
    2244                 :            : 
    2245                 :            : static void
    2246                 :        160 : init_cScT(ST_t *T, GEN dtcr, long N, long prec)
    2247                 :            : {
    2248                 :        160 :   GEN p1 = ch_4(dtcr);
    2249                 :        160 :   T->a = p1[1];
    2250                 :        160 :   T->b = p1[2];
    2251                 :        160 :   T->c = p1[3];
    2252                 :        160 :   T->rc1 = T->a + T->c;
    2253                 :        160 :   T->rc2 = T->b + T->c;
    2254                 :        160 :   T->r   = maxss(T->rc2+1, T->rc1); /* >= 2 */
    2255                 :        160 :   ppgamma(T, prec);
    2256                 :        160 :   clear_cScT(T, N);
    2257                 :        160 : }
    2258                 :            : 
    2259                 :            : static void
    2260                 :        120 : GetST0(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec)
    2261                 :            : {
    2262                 :        120 :   pari_sp av = avma, av1, av2;
    2263                 :            :   long ncond, n, j, k, jc, n0, prec2, i0, r1, r2;
    2264                 :        120 :   GEN nf = checknf(bnr), racpi, powracpi;
    2265                 :        120 :   GEN N0, C, T = *pT, S = *pS, an, degs, limx;
    2266                 :            :   LISTray LIST;
    2267                 :            :   ST_t cScT;
    2268                 :            : 
    2269                 :            :   /* initializations */
    2270                 :        120 :   degs = GetDeg(dataCR);
    2271                 :        120 :   ncond = lg(vChar)-1;
    2272                 :        120 :   nf_get_sign(nf,&r1,&r2);
    2273                 :            : 
    2274                 :        120 :   C  = cgetg(ncond+1, t_VEC);
    2275                 :        120 :   N0 = cgetg(ncond+1, t_VECSMALL);
    2276                 :        120 :   n0 = 0;
    2277                 :        120 :   limx = zeta_get_limx(r1, r2, prec2nbits(prec));
    2278         [ +  + ]:        280 :   for (j = 1; j <= ncond; j++)
    2279                 :            :   {
    2280                 :        160 :     GEN dtcr = gel(dataCR, mael(vChar,j,1)), c = ch_C(dtcr);
    2281                 :        160 :     gel(C,j) = c;
    2282                 :        160 :     N0[j] = zeta_get_N0(c, limx);
    2283         [ +  + ]:        160 :     if (n0 < N0[j]) n0  = N0[j];
    2284                 :            :   }
    2285                 :        120 :   i0 = zeta_get_i0(r1, r2, prec2nbits(prec), limx);
    2286                 :        120 :   InitPrimes(bnr, n0, &LIST);
    2287                 :            : 
    2288                 :        120 :   prec2 = precdbl(prec) + EXTRA_PREC;
    2289                 :        120 :   racpi = sqrtr(mppi(prec2));
    2290                 :        120 :   powracpi = cgetg(r1+2,t_VEC);
    2291                 :        120 :   gel(powracpi,1) = racpi;
    2292         [ +  + ]:        335 :   for (j=2; j<=r1; j++) gel(powracpi,j) = mulrr(gel(powracpi,j-1), racpi);
    2293                 :        120 :   cScT.powracpi = powracpi;
    2294                 :            : 
    2295                 :        120 :   cScT.cS = cgetg(n0+1, t_VEC);
    2296                 :        120 :   cScT.cT = cgetg(n0+1, t_VEC);
    2297         [ +  + ]:     887306 :   for (j=1; j<=n0; j++) gel(cScT.cS,j) = gel(cScT.cT,j) = NULL;
    2298                 :            : 
    2299                 :        120 :   cScT.i0 = i0;
    2300                 :            : 
    2301                 :        120 :   av1 = avma;
    2302         [ +  + ]:        280 :   for (jc = 1; jc <= ncond; jc++)
    2303                 :            :   {
    2304                 :        160 :     const GEN LChar = gel(vChar,jc);
    2305                 :        160 :     const long nChar = lg(LChar)-1, NN = N0[jc];
    2306                 :            : 
    2307         [ -  + ]:        160 :     if (DEBUGLEVEL>1)
    2308                 :          0 :       err_printf("* conductor no %ld/%ld (N = %ld)\n\tInit: ", jc,ncond,NN);
    2309                 :            : 
    2310                 :        160 :     cScT.c1 = gel(C,jc);
    2311                 :        160 :     init_cScT(&cScT, gel(dataCR, LChar[1]), NN, prec2);
    2312                 :        160 :     av2 = avma;
    2313         [ +  + ]:        490 :     for (k = 1; k <= nChar; k++)
    2314                 :            :     {
    2315                 :        330 :       const long t = LChar[k];
    2316         [ -  + ]:        330 :       if (DEBUGLEVEL>1)
    2317                 :          0 :         err_printf("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar);
    2318                 :            : 
    2319         [ +  + ]:        330 :       if (!isintzero( ch_comp(gel(dataCR, t)) ))
    2320                 :            :       {
    2321                 :        325 :         const long d = degs[t];
    2322                 :        325 :         const GEN dtcr = gel(dataCR, t), z = gel(ch_CHI(dtcr), 2);
    2323                 :        325 :         GEN p1 = gen_0, p2 = gen_0;
    2324                 :        325 :         long c = 0;
    2325                 :        325 :         int **matan = ComputeCoeff(gel(dataCR,t), &LIST, NN, d);
    2326         [ +  + ]:    2160274 :         for (n = 1; n <= NN; n++)
    2327         [ +  + ]:    2159949 :           if ((an = EvalCoeff(z, matan[n], d)))
    2328                 :            :           {
    2329                 :     604741 :            get_cS_cT(&cScT, n);
    2330                 :     604741 :            p1 = gadd(p1, gmul(an, gel(cScT.cS,n)));
    2331                 :     604741 :            p2 = gadd(p2, gmul(an, gel(cScT.cT,n)));
    2332         [ +  + ]:     604741 :            if (++c == 256) { gerepileall(av2,2, &p1,&p2); c = 0; }
    2333                 :            :           }
    2334                 :        325 :         gaffect(p1,        gel(S,t));
    2335                 :        325 :         gaffect(gconj(p2), gel(T,t));
    2336                 :        325 :         FreeMat(matan, NN); avma = av2;
    2337                 :            :       }
    2338         [ -  + ]:          5 :       else if (DEBUGLEVEL>1)
    2339                 :          0 :         err_printf("\t  no need to compute this character\n");
    2340                 :            :     }
    2341         [ -  + ]:        160 :     if (DEBUGLEVEL>1) err_printf("\n");
    2342                 :        160 :     avma = av1;
    2343                 :            :   }
    2344                 :        120 :   clear_cScT(&cScT, n0);
    2345                 :        120 :   avma = av;
    2346                 :        120 : }
    2347                 :            : 
    2348                 :            : static void
    2349                 :        260 : GetST(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec)
    2350                 :            : {
    2351                 :        260 :   const long cl = lg(dataCR) - 1;
    2352                 :        260 :   GEN S, T, nf  = checknf(bnr);
    2353                 :            :   long j;
    2354                 :            : 
    2355                 :            :   /* allocate memory for answer */
    2356                 :        260 :   *pS = S = cgetg(cl+1, t_VEC);
    2357                 :        260 :   *pT = T = cgetg(cl+1, t_VEC);
    2358         [ +  + ]:       1115 :   for (j = 1; j <= cl; j++)
    2359                 :            :   {
    2360                 :        855 :     gel(S,j) = cgetc(prec);
    2361                 :        855 :     gel(T,j) = cgetc(prec);
    2362                 :            :   }
    2363         [ +  + ]:        260 :   if (nf_get_degree(nf) == 2)
    2364                 :            :   {
    2365                 :        155 :     QuadGetST(bnr, pS, pT, dataCR, vChar, prec);
    2366                 :        260 :     return;
    2367                 :            :   }
    2368                 :        105 :   GetST0(bnr, pS, pT, dataCR, vChar, prec);
    2369                 :            : }
    2370                 :            : 
    2371                 :            : /*******************************************************************/
    2372                 :            : /*                                                                 */
    2373                 :            : /*     Class fields of real quadratic fields using Stark units     */
    2374                 :            : /*                                                                 */
    2375                 :            : /*******************************************************************/
    2376                 :            : /* compute the Hilbert class field using genus class field theory when
    2377                 :            :    the exponent of the class group is exactly 2 (trivial group not covered) */
    2378                 :            : /* Cf Herz, Construction of class fields, LNM 21, Theorem 1 (VII-6) */
    2379                 :            : static GEN
    2380                 :         10 : GenusFieldQuadReal(GEN disc)
    2381                 :            : {
    2382                 :         10 :   long i, i0 = 0, l;
    2383                 :         10 :   pari_sp av = avma;
    2384                 :         10 :   GEN T = NULL, p0 = NULL, P;
    2385                 :            : 
    2386                 :         10 :   P = gel(Z_factor(disc), 1);
    2387                 :         10 :   l = lg(P);
    2388         [ +  + ]:         30 :   for (i = 1; i < l; i++)
    2389                 :            :   {
    2390                 :         25 :     GEN p = gel(P,i);
    2391         [ +  + ]:         25 :     if (mod4(p) == 3) { p0 = p; i0 = i; break; }
    2392                 :            :   }
    2393                 :         10 :   l--; /* remove last prime */
    2394         [ -  + ]:         10 :   if (i0 == l) l--; /* ... remove p0 and last prime */
    2395         [ +  + ]:         35 :   for (i = 1; i < l; i++)
    2396                 :            :   {
    2397                 :         25 :     GEN p = gel(P,i), d, t;
    2398         [ +  + ]:         25 :     if (i == i0) continue;
    2399         [ +  + ]:         20 :     if (equaliu(p, 2))
    2400      [ +  -  - ]:         10 :       switch (mod32(disc))
    2401                 :            :       {
    2402                 :         10 :         case  8: d = gen_2; break;
    2403                 :          0 :         case 24: d = shifti(p0, 1); break;
    2404                 :         10 :         default: d = p0; break;
    2405                 :            :       }
    2406                 :            :     else
    2407         [ -  + ]:         10 :       d = (mod4(p) == 1)? p: mulii(p0, p);
    2408                 :         20 :     t = mkpoln(3, gen_1, gen_0, negi(d)); /* x^2 - d */
    2409         [ +  + ]:         20 :     T = T? ZX_compositum_disjoint(T, t): t;
    2410                 :            :   }
    2411                 :         10 :   return gerepileupto(av, polredbest(T, 0));
    2412                 :            : }
    2413                 :            : static GEN
    2414                 :        280 : GenusFieldQuadImag(GEN disc)
    2415                 :            : {
    2416                 :            :   long i, l;
    2417                 :        280 :   pari_sp av = avma;
    2418                 :        280 :   GEN T = NULL, P;
    2419                 :            : 
    2420                 :        280 :   P = gel(absi_factor(disc), 1);
    2421                 :        280 :   l = lg(P);
    2422                 :        280 :   l--; /* remove last prime */
    2423         [ +  + ]:        825 :   for (i = 1; i < l; i++)
    2424                 :            :   {
    2425                 :        545 :     GEN p = gel(P,i), d, t;
    2426         [ +  + ]:        545 :     if (equaliu(p, 2))
    2427      [ +  +  + ]:        160 :       switch (mod32(disc))
    2428                 :            :       {
    2429                 :         40 :         case 24: d = gen_2; break;  /* disc = 8 mod 32 */
    2430                 :         30 :         case  8: d = gen_m2; break; /* disc =-8 mod 32 */
    2431                 :        160 :         default: d = gen_m1; break;
    2432                 :            :       }
    2433                 :            :     else
    2434         [ +  + ]:        385 :       d = (mod4(p) == 1)? p: negi(p);
    2435                 :        545 :     t = mkpoln(3, gen_1, gen_0, negi(d)); /* x^2 - d */
    2436         [ +  + ]:        545 :     T = T? ZX_compositum_disjoint(T, t): t;
    2437                 :            :   }
    2438                 :        280 :   return gerepileupto(av, polredbest(T, 0));
    2439                 :            : }
    2440                 :            : 
    2441                 :            : /* if flag != 0, computes a fast and crude approximation of the result */
    2442                 :            : static GEN
    2443                 :        450 : AllStark(GEN data,  GEN nf,  long flag,  long newprec)
    2444                 :            : {
    2445                 :        450 :   const long BND = 300;
    2446                 :        450 :   long cl, i, j, cpt = 0, N, h, v, n, r1, r2, den;
    2447                 :            :   pari_sp av, av2;
    2448                 :            :   int **matan;
    2449                 :        450 :   GEN bnr = gel(data,1), p1, p2, S, T, polrelnum, polrel, Lp, W, veczeta;
    2450                 :            :   GEN vChar, degs, C, dataCR, cond1, L1, an;
    2451                 :            :   LISTray LIST;
    2452                 :            :   pari_timer ti;
    2453                 :            : 
    2454                 :        450 :   nf_get_sign(nf, &r1,&r2);
    2455                 :        450 :   N     = nf_get_degree(nf);
    2456                 :        450 :   cond1 = gel(bnr_get_mod(bnr), 2);
    2457                 :        450 :   dataCR = gel(data,5);
    2458                 :        450 :   vChar = sortChars(dataCR);
    2459                 :            : 
    2460                 :        450 :   v = 1;
    2461         [ +  + ]:        930 :   while (gequal1(gel(cond1,v))) v++;
    2462                 :            : 
    2463                 :        450 :   cl = lg(dataCR)-1;
    2464                 :        450 :   degs = GetDeg(dataCR);
    2465                 :        450 :   h  = itos(ZM_det_triangular(gel(data,2))) >> 1;
    2466                 :            : 
    2467                 :            : LABDOUB:
    2468         [ -  + ]:        455 :   if (DEBUGLEVEL) timer_start(&ti);
    2469                 :        455 :   av = avma;
    2470                 :            : 
    2471                 :            :   /* characters with rank > 1 should not be computed */
    2472         [ +  + ]:       1950 :   for (i = 1; i <= cl; i++)
    2473                 :            :   {
    2474                 :       1495 :     GEN chi = gel(dataCR, i);
    2475         [ +  + ]:       1495 :     if (L_vanishes_at_0(chi)) ch_comp(chi) = gen_0;
    2476                 :            :   }
    2477                 :            : 
    2478                 :        455 :   W = ComputeAllArtinNumbers(dataCR, vChar, (flag >= 0), newprec);
    2479         [ -  + ]:        455 :   if (DEBUGLEVEL) timer_printf(&ti,"Compute W");
    2480                 :        455 :   Lp = cgetg(cl + 1, t_VEC);
    2481         [ +  + ]:        455 :   if (!flag)
    2482                 :            :   {
    2483                 :        230 :     GetST(bnr, &S, &T, dataCR, vChar, newprec);
    2484         [ -  + ]:        230 :     if (DEBUGLEVEL) timer_printf(&ti, "S&T");
    2485         [ +  + ]:        990 :     for (i = 1; i <= cl; i++)
    2486                 :            :     {
    2487                 :        760 :       GEN chi = gel(dataCR, i), v = gen_0;
    2488         [ +  + ]:        760 :       if (!isintzero( ch_comp(chi) ))
    2489                 :        750 :         v = gel(GetValue(chi, gel(W,i), gel(S,i), gel(T,i), 2, newprec), 2);
    2490                 :        760 :       gel(Lp, i) = v;
    2491                 :            :     }
    2492                 :            :   }
    2493                 :            :   else
    2494                 :            :   { /* compute a crude approximation of the result */
    2495                 :        225 :     C = cgetg(cl + 1, t_VEC);
    2496         [ +  + ]:        960 :     for (i = 1; i <= cl; i++) gel(C,i) = ch_C(gel(dataCR, i));
    2497                 :        225 :     n = zeta_get_N0(vecmax(C), zeta_get_limx(r1, r2, prec2nbits(newprec)));
    2498         [ +  + ]:        225 :     if (n > BND) n = BND;
    2499         [ -  + ]:        225 :     if (DEBUGLEVEL) err_printf("N0 in QuickPol: %ld \n", n);
    2500                 :        225 :     InitPrimes(bnr, n, &LIST);
    2501                 :            : 
    2502                 :        225 :     L1 = cgetg(cl+1, t_VEC);
    2503                 :            :     /* use L(1) = sum (an / n) */
    2504         [ +  + ]:        960 :     for (i = 1; i <= cl; i++)
    2505                 :            :     {
    2506                 :        735 :       GEN dtcr = gel(dataCR,i);
    2507                 :        735 :       matan = ComputeCoeff(dtcr, &LIST, n, degs[i]);
    2508                 :        735 :       av2 = avma;
    2509                 :        735 :       p1 = real_0(newprec); p2 = gel(ch_CHI(dtcr), 2);
    2510         [ +  + ]:     216485 :       for (j = 1; j <= n; j++)
    2511         [ +  + ]:     215750 :         if ( (an = EvalCoeff(p2, matan[j], degs[i])) )
    2512                 :      80385 :           p1 = gadd(p1, gdivgs(an, j));
    2513                 :        735 :       gel(L1,i) = gerepileupto(av2, p1);
    2514                 :        735 :       FreeMat(matan, n);
    2515                 :            :     }
    2516                 :        225 :     p1 = gmul2n(powruhalf(mppi(newprec), N-2), 1);
    2517                 :            : 
    2518         [ +  + ]:        960 :     for (i = 1; i <= cl; i++)
    2519                 :            :     {
    2520                 :            :       long r;
    2521                 :        735 :       GEN WW, A = ComputeAChi(gel(dataCR,i), &r, 0, newprec);
    2522                 :        735 :       WW = gmul(gel(C,i), gmul(A, gel(W,i)));
    2523                 :        735 :       gel(Lp,i) = gdiv(gmul(WW, gconj(gel(L1,i))), p1);
    2524                 :            :     }
    2525                 :            :   }
    2526                 :            : 
    2527                 :        455 :   p1 = ComputeLift(gel(data,4));
    2528                 :            : 
    2529         [ +  + ]:        455 :   den = flag ? h: 2*h;
    2530                 :        455 :   veczeta = cgetg(h + 1, t_VEC);
    2531         [ +  + ]:       2895 :   for (i = 1; i <= h; i++)
    2532                 :            :   {
    2533                 :       2440 :     GEN z = gen_0, sig = gel(p1,i);
    2534         [ +  + ]:      12570 :     for (j = 1; j <= cl; j++)
    2535                 :            :     {
    2536                 :      10130 :       GEN dtcr = gel(dataCR,j), CHI = ch_CHI(dtcr);
    2537                 :      10130 :       GEN t = mulreal(gel(Lp,j), ComputeImagebyChar(CHI, sig));
    2538         [ +  + ]:      10130 :       if (itos(gel(CHI,3)) != 2) t = gmul2n(t, 1); /* character not real */
    2539                 :      10130 :       z = gadd(z, t);
    2540                 :            :     }
    2541                 :       2440 :     gel(veczeta,i) = gdivgs(z, den);
    2542                 :            :   }
    2543         [ +  + ]:       2895 :   for (j = 1; j <= h; j++)
    2544                 :       2440 :     gel(veczeta,j) = gmul2n(gcosh(gel(veczeta,j), newprec), 1);
    2545                 :        455 :   polrelnum = roots_to_pol(veczeta, 0);
    2546         [ -  + ]:        455 :   if (DEBUGLEVEL)
    2547                 :            :   {
    2548         [ #  # ]:          0 :     if (DEBUGLEVEL>1) {
    2549                 :          0 :       err_printf("polrelnum = %Ps\n", polrelnum);
    2550                 :          0 :       err_printf("zetavalues = %Ps\n", veczeta);
    2551         [ #  # ]:          0 :       if (!flag)
    2552                 :          0 :         err_printf("Checking the square-root of the Stark unit...\n");
    2553                 :            :     }
    2554         [ #  # ]:          0 :     timer_printf(&ti, "Compute %s", flag? "quickpol": "polrelnum");
    2555                 :            :   }
    2556                 :            : 
    2557         [ +  + ]:        455 :   if (flag)
    2558                 :        225 :     return gerepilecopy(av, polrelnum);
    2559                 :            : 
    2560                 :            :   /* try to recognize this polynomial */
    2561                 :        230 :   polrel = RecCoeff(nf, polrelnum, v, newprec);
    2562         [ +  + ]:        230 :   if (!polrel)
    2563                 :            :   {
    2564         [ +  + ]:       1125 :     for (j = 1; j <= h; j++)
    2565                 :        965 :       gel(veczeta,j) = gsubgs(gsqr(gel(veczeta,j)), 2);
    2566                 :        160 :     polrelnum = roots_to_pol(veczeta, 0);
    2567         [ -  + ]:        160 :     if (DEBUGLEVEL)
    2568                 :            :     {
    2569         [ #  # ]:          0 :       if (DEBUGLEVEL>1) {
    2570                 :          0 :         err_printf("It's not a square...\n");
    2571                 :          0 :         err_printf("polrelnum = %Ps\n", polrelnum);
    2572                 :            :       }
    2573                 :          0 :       timer_printf(&ti, "Compute polrelnum");
    2574                 :            :     }
    2575                 :        160 :     polrel = RecCoeff(nf, polrelnum, v, newprec);
    2576                 :            :   }
    2577         [ +  + ]:        230 :   if (!polrel) /* FAILED */
    2578                 :            :   {
    2579                 :            :     long incr_pr;
    2580         [ -  + ]:          5 :     if (++cpt >= 3) pari_err_PREC( "stark (computation impossible)");
    2581                 :            : 
    2582                 :            :     /* compute the precision, we need
    2583                 :            :           a) get at least EXTRA_PREC fractional digits if there is none;
    2584                 :            :        or b) double the fractional digits.
    2585                 :            :     */
    2586                 :          5 :     incr_pr = nbits2extraprec( prec2nbits(gprecision(polrelnum))- gexpo(polrelnum) );
    2587         [ -  + ]:          5 :     if (incr_pr < 0) incr_pr = -incr_pr + EXTRA_PREC;
    2588                 :          5 :     newprec = newprec + maxss(ADD_PREC, cpt*incr_pr);
    2589         [ -  + ]:          5 :     if (DEBUGLEVEL) pari_warn(warnprec, "AllStark", newprec);
    2590                 :            : 
    2591                 :          5 :     nf = nfnewprec_shallow(nf, newprec);
    2592                 :          5 :     dataCR = CharNewPrec(dataCR, nf, newprec);
    2593                 :            : 
    2594                 :          5 :     gerepileall(av, 2, &nf, &dataCR);
    2595                 :          5 :     goto LABDOUB;
    2596                 :            :   }
    2597                 :            : 
    2598         [ -  + ]:        225 :   if (DEBUGLEVEL) {
    2599         [ #  # ]:          0 :     if (DEBUGLEVEL>1) err_printf("polrel = %Ps\n", polrel);
    2600                 :          0 :     timer_printf(&ti, "Recpolnum");
    2601                 :            :   }
    2602                 :        450 :   return gerepilecopy(av, polrel);
    2603                 :            : }
    2604                 :            : 
    2605                 :            : /********************************************************************/
    2606                 :            : /*                        Main functions                            */
    2607                 :            : /********************************************************************/
    2608                 :            : 
    2609                 :            : static GEN
    2610                 :        175 : get_subgroup(GEN H, GEN cyc, const char *s)
    2611                 :            : {
    2612 [ +  + ][ +  + ]:        175 :   if (!H || gequal0(H)) return diagonal_shallow(cyc);
    2613         [ +  + ]:         10 :   if (typ(H) != t_MAT) pari_err_TYPE(stack_strcat(s," [subgroup]"), H);
    2614                 :          5 :   RgM_check_ZM(H, s);
    2615                 :        170 :   return ZM_hnfmodid(H, cyc);
    2616                 :            : }
    2617                 :            : 
    2618                 :            : GEN
    2619                 :        145 : bnrstark(GEN bnr, GEN subgrp, long prec)
    2620                 :            : {
    2621                 :            :   long N, newprec;
    2622                 :        145 :   pari_sp av = avma;
    2623                 :            :   GEN bnf, p1, cycbnr, nf, data, dtQ;
    2624                 :            : 
    2625                 :            :   /* check the bnr */
    2626                 :        145 :   checkbnr(bnr);
    2627                 :        145 :   bnf = checkbnf(bnr);
    2628                 :        145 :   nf  = bnf_get_nf(bnf);
    2629                 :        145 :   N   = nf_get_degree(nf);
    2630         [ -  + ]:        145 :   if (N == 1) return galoissubcyclo(bnr, subgrp, 0, 0);
    2631                 :            : 
    2632                 :            :   /* check the bnf */
    2633         [ -  + ]:        145 :   if (!nf_get_varn(nf))
    2634                 :          0 :     pari_err_PRIORITY("bnrstark", nf_get_pol(nf), "=", 0);
    2635         [ +  + ]:        145 :   if (nf_get_r2(nf)) pari_err_DOMAIN("bnrstark", "r2", "!=", gen_0, nf);
    2636                 :        140 :   subgrp = get_subgroup(subgrp,bnr_get_cyc(bnr),"bnrstark");
    2637                 :            : 
    2638                 :            :   /* compute bnr(conductor) */
    2639                 :        140 :   p1     = bnrconductor(bnr, subgrp, 2);
    2640                 :        140 :   bnr    = gel(p1,2); cycbnr = bnr_get_cyc(bnr);
    2641                 :        140 :   subgrp = gel(p1,3);
    2642         [ -  + ]:        140 :   if (gequal1( ZM_det_triangular(subgrp) )) { avma = av; return pol_x(0); }
    2643                 :            : 
    2644                 :            :   /* check the class field */
    2645         [ +  + ]:        140 :   if (!gequal0(gel(bnr_get_mod(bnr), 2)))
    2646                 :          5 :     pari_err_DOMAIN("bnrstark", "r2(class field)", "!=", gen_0, bnr);
    2647                 :            : 
    2648                 :            :   /* find a suitable extension N */
    2649                 :        135 :   dtQ = InitQuotient(subgrp);
    2650                 :        135 :   data  = FindModulus(bnr, dtQ, &newprec);
    2651         [ -  + ]:        135 :   if (!data)
    2652                 :            :   {
    2653                 :          0 :     GEN vec, H, cyc = gel(dtQ,2), U = gel(dtQ,3), M = RgM_inv(U);
    2654                 :          0 :     long i, j = 1, l = lg(M);
    2655                 :            : 
    2656                 :            :     /* M = indep. generators of Cl_f/subgp, restrict to cyclic components */
    2657                 :          0 :     vec = cgetg(l, t_VEC);
    2658         [ #  # ]:          0 :     for (i = 1; i < l; i++)
    2659                 :            :     {
    2660         [ #  # ]:          0 :       if (is_pm1(gel(cyc,i))) continue;
    2661                 :          0 :       H = ZM_hnfmodid(vecsplice(M,i), cycbnr);
    2662                 :          0 :       gel(vec,j++) = bnrstark(bnr, H, prec);
    2663                 :            :     }
    2664                 :          0 :     setlg(vec, j); return gerepilecopy(av, vec);
    2665                 :            :   }
    2666                 :            : 
    2667         [ +  - ]:        135 :   if (newprec > prec)
    2668                 :            :   {
    2669         [ -  + ]:        135 :     if (DEBUGLEVEL>1) err_printf("new precision: %ld\n", newprec);
    2670                 :        135 :     nf = nfnewprec_shallow(nf, newprec);
    2671                 :            :   }
    2672                 :        135 :   return gerepileupto(av, AllStark(data, nf, 0, newprec));
    2673                 :            : }
    2674                 :            : 
    2675                 :            : /* For each character of Cl(bnr)/subgp, compute L(1, chi) (or equivalently
    2676                 :            :  * the first non-zero term c(chi) of the expansion at s = 0).
    2677                 :            :  * If flag & 1: compute the value at s = 1 (for non-trivial characters),
    2678                 :            :  * else compute the term c(chi) and return [r(chi), c(chi)] where r(chi) is
    2679                 :            :  *   the order of L(s, chi) at s = 0.
    2680                 :            :  * If flag & 2: compute the value of the L-function L_S(s, chi) where S is the
    2681                 :            :  *   set of places dividing the modulus of bnr (and the infinite places),
    2682                 :            :  * else
    2683                 :            :  *   compute the value of the primitive L-function associated to chi,
    2684                 :            :  * If flag & 4: return also the character */
    2685                 :            : GEN
    2686                 :         35 : bnrL1(GEN bnr, GEN subgp, long flag, long prec)
    2687                 :            : {
    2688                 :            :   GEN cyc, L1, allCR, listCR;
    2689                 :            :   GEN indCR, invCR, Qt;
    2690                 :            :   long cl, i, nc;
    2691                 :         35 :   pari_sp av = avma;
    2692                 :            : 
    2693                 :         35 :   checkbnr(bnr);
    2694 [ +  - ][ -  + ]:         35 :   if (flag < 0 || flag > 8) pari_err_FLAG("bnrL1");
    2695                 :            : 
    2696                 :            :   /* compute bnr(conductor) */
    2697         [ +  - ]:         35 :   if (!(flag & 2)) bnr = gel(bnrconductor(bnr, NULL, 2),2);
    2698                 :         35 :   cyc  = bnr_get_cyc(bnr);
    2699                 :         35 :   subgp = get_subgroup(subgp, cyc, "bnrL1");
    2700                 :            : 
    2701                 :         30 :   cl = itou( ZM_det_triangular(subgp) );
    2702                 :         30 :   Qt = InitQuotient(subgp);
    2703                 :            : 
    2704                 :            :   /* compute all characters */
    2705                 :         30 :   allCR = EltsOfGroup(cl, gel(Qt,2));
    2706                 :            : 
    2707                 :            :   /* make a list of all non-trivial characters modulo conjugation */
    2708                 :         30 :   listCR = cgetg(cl, t_VEC);
    2709                 :         30 :   indCR = new_chunk(cl);
    2710                 :         30 :   invCR = new_chunk(cl); nc = 0;
    2711         [ +  + ]:        185 :   for (i = 1; i < cl; i++)
    2712                 :            :   {
    2713                 :            :     /* lift to a character on Cl(bnr) */
    2714                 :        155 :     GEN lchi = LiftChar(cyc, gel(Qt,3), gel(allCR,i), gel(Qt,2));
    2715                 :        155 :     GEN clchi = ConjChar(lchi, cyc);
    2716                 :        155 :     long j, a = i;
    2717         [ +  + ]:        495 :     for (j = 1; j <= nc; j++)
    2718         [ +  + ]:        400 :       if (ZV_equal(gmael(listCR, j, 1), clchi)) { a = -j; break; }
    2719                 :            : 
    2720         [ +  + ]:        155 :     if (a > 0)
    2721                 :            :     {
    2722                 :         95 :       nc++;
    2723                 :         95 :       gel(listCR,nc) = mkvec2(lchi, bnrconductorofchar(bnr, lchi));
    2724                 :         95 :       indCR[i]  = nc;
    2725                 :         95 :       invCR[nc] = i;
    2726                 :            :     }
    2727                 :            :     else
    2728                 :         60 :       indCR[i] = -invCR[-a];
    2729                 :            : 
    2730                 :        155 :     gel(allCR,i) = lchi;
    2731                 :            :   }
    2732                 :         30 :   settyp(allCR[cl], t_VEC); /* set correct type for trivial character */
    2733                 :            : 
    2734                 :         30 :   setlg(listCR, nc + 1);
    2735         [ +  - ]:         30 :   L1 = cgetg((flag&1)? cl: cl+1, t_VEC);
    2736         [ +  - ]:         30 :   if (nc)
    2737                 :            :   {
    2738                 :         30 :     GEN dataCR = InitChar(bnr, listCR, prec);
    2739                 :         30 :     GEN W, S, T, vChar = sortChars(dataCR);
    2740                 :         30 :     GetST(bnr, &S, &T, dataCR, vChar, prec);
    2741                 :         30 :     W = ComputeAllArtinNumbers(dataCR, vChar, 1, prec);
    2742         [ +  + ]:        185 :     for (i = 1; i < cl; i++)
    2743                 :            :     {
    2744                 :        155 :       long a = indCR[i];
    2745         [ +  + ]:        155 :       if (a > 0)
    2746                 :         95 :         gel(L1,i) = GetValue(gel(dataCR,a), gel(W,a), gel(S,a), gel(T,a),
    2747                 :            :                              flag, prec);
    2748                 :            :       else
    2749                 :         60 :         gel(L1,i) = gconj(gel(L1,-a));
    2750                 :            :     }
    2751                 :            :   }
    2752         [ +  - ]:         30 :   if (!(flag & 1))
    2753                 :         30 :     gel(L1,cl) = GetValue1(bnr, flag & 2, prec);
    2754                 :            :   else
    2755                 :          0 :     cl--;
    2756                 :            : 
    2757         [ -  + ]:         30 :   if (flag & 4) {
    2758         [ #  # ]:          0 :     for (i = 1; i <= cl; i++) gel(L1,i) = mkvec2(gel(allCR,i), gel(L1,i));
    2759                 :            :   }
    2760                 :         30 :   return gerepilecopy(av, L1);
    2761                 :            : }
    2762                 :            : 
    2763                 :            : /*******************************************************************/
    2764                 :            : /*                                                                 */
    2765                 :            : /*       Hilbert and Ray Class field using Stark                   */
    2766                 :            : /*                                                                 */
    2767                 :            : /*******************************************************************/
    2768                 :            : /* P in A[x,y], deg_y P < 2, return P0 and P1 in A[x] such that P = P0 + P1 y */
    2769                 :            : static void
    2770                 :         90 : split_pol_quad(GEN P, GEN *gP0, GEN *gP1)
    2771                 :            : {
    2772                 :         90 :   long i, l = lg(P);
    2773                 :         90 :   GEN P0 = cgetg(l, t_POL), P1 = cgetg(l, t_POL);
    2774                 :         90 :   P0[1] = P1[1] = P[1];
    2775         [ +  + ]:        835 :   for (i = 2; i < l; i++)
    2776                 :            :   {
    2777                 :        745 :     GEN c = gel(P,i), c0 = c, c1 = gen_0;
    2778         [ +  + ]:        745 :     if (typ(c) == t_POL) /* write c = c1 y + c0 */
    2779      [ -  +  - ]:        655 :       switch(degpol(c))
    2780                 :            :       {
    2781                 :          0 :         case -1: c0 = gen_0; break;
    2782                 :        655 :         default: c1 = gel(c,3); /* fall through */
    2783                 :        655 :         case  0: c0 = gel(c,2); break;
    2784                 :            :       }
    2785                 :        745 :     gel(P0,i) = c0; gel(P1,i) = c1;
    2786                 :            :   }
    2787                 :         90 :   *gP0 = normalizepol_lg(P0, l);
    2788                 :         90 :   *gP1 = normalizepol_lg(P1, l);
    2789                 :         90 : }
    2790                 :            : 
    2791                 :            : /* k = nf quadratic field, P relative equation of H_k (Hilbert class field)
    2792                 :            :  * return T in Z[X], such that H_k / Q is the compositum of Q[X]/(T) and k */
    2793                 :            : static GEN
    2794                 :         90 : makescind(GEN nf, GEN P)
    2795                 :            : {
    2796                 :         90 :   GEN Pp, p, pol, G, L, a, roo, P0,P1, Ny,Try, nfpol = nf_get_pol(nf);
    2797                 :            :   long i, is_P;
    2798                 :            : 
    2799                 :         90 :   split_pol_quad(lift_intern(P), &P0, &P1);
    2800                 :            :   /* P = P0 + y P1, Norm_{k/Q}(P) = P0^2 + Tr y P0P1 + Ny P1^2, irreducible/Q */
    2801                 :         90 :   Ny = gel(nfpol, 2);
    2802                 :         90 :   Try = negi(gel(nfpol, 3));
    2803                 :         90 :   pol = RgX_add(RgX_sqr(P0), RgX_Rg_mul(RgX_sqr(P1), Ny));
    2804         [ +  + ]:         90 :   if (signe(Try)) pol = RgX_add(pol, RgX_Rg_mul(RgX_mul(P0,P1), Try));
    2805                 :            :   /* pol = rnfequation(nf, P); */
    2806                 :         90 :   G = galoisinit(pol, NULL);
    2807                 :         90 :   L = gal_get_group(G);
    2808                 :         90 :   p = gal_get_p(G);
    2809                 :         90 :   a = FpX_quad_root(nfpol, p, 0);
    2810                 :            :   /* P mod a prime \wp above p (which splits) */
    2811                 :         90 :   Pp = FpXY_evalx(P, a, p);
    2812                 :         90 :   roo = gal_get_roots(G);
    2813                 :         90 :   is_P = gequal0( FpX_eval(Pp, remii(gel(roo,1),p), p) );
    2814                 :            :   /* each roo[i] mod p is a root of P or (exclusive) tau(P) mod \wp */
    2815                 :            :   /* record whether roo[1] is a root of P or tau(P) */
    2816                 :            : 
    2817         [ +  - ]:        715 :   for (i = 1; i < lg(L); i++)
    2818                 :            :   {
    2819                 :        715 :     GEN perm = gel(L,i);
    2820         [ +  + ]:        715 :     long k = perm[1]; if (k == 1) continue;
    2821                 :        625 :     k = gequal0( FpX_eval(Pp, remii(gel(roo,k),p), p) );
    2822                 :            :     /* roo[k] is a root of the other polynomial */
    2823         [ +  + ]:        625 :     if (k != is_P)
    2824                 :            :     {
    2825                 :         90 :       long o = perm_order(perm);
    2826         [ -  + ]:         90 :       if (o != 2) perm = perm_pow(perm, o >> 1);
    2827                 :            :       /* perm has order two and doesn't belong to Gal(H_k/k) */
    2828                 :         90 :       return galoisfixedfield(G, perm, 1, varn(P));
    2829                 :            :     }
    2830                 :            :   }
    2831                 :          0 :   pari_err_BUG("makescind");
    2832                 :         90 :   return NULL; /*not reached*/
    2833                 :            : }
    2834                 :            : 
    2835                 :            : /* pbnf = NULL if no bnf is needed, f = NULL may be passed for a trivial
    2836                 :            :  * conductor */
    2837                 :            : static void
    2838                 :        555 : quadray_init(GEN *pD, GEN f, GEN *pbnf, long prec)
    2839                 :            : {
    2840                 :        555 :   GEN D = *pD, nf, bnf = NULL;
    2841         [ +  + ]:        555 :   if (typ(D) == t_INT)
    2842                 :            :   {
    2843                 :            :     int isfund;
    2844         [ +  + ]:        540 :     if (pbnf) {
    2845         [ +  + ]:        145 :       long v = f? gvar(f): NO_VARIABLE;
    2846         [ +  - ]:        145 :       if (v == NO_VARIABLE) v = fetch_user_var("y");
    2847                 :        145 :       bnf = Buchall(quadpoly0(D, v), nf_FORCE, prec);
    2848                 :        145 :       nf = bnf_get_nf(bnf);
    2849                 :        145 :       isfund = equalii(D, nf_get_disc(nf));
    2850                 :            :     }
    2851                 :            :     else
    2852                 :        395 :       isfund = Z_isfundamental(D);
    2853         [ +  + ]:        540 :     if (!isfund) pari_err_DOMAIN("quadray", "isfundamental(D)", "=",gen_0, D);
    2854                 :            :   }
    2855                 :            :   else
    2856                 :            :   {
    2857                 :         15 :     bnf = checkbnf(D);
    2858                 :         15 :     nf = bnf_get_nf(bnf);
    2859         [ +  + ]:         15 :     if (nf_get_degree(nf) != 2)
    2860                 :          5 :       pari_err_DOMAIN("quadray", "degree", "!=", gen_2, nf_get_pol(nf));
    2861                 :         10 :     D = nf_get_disc(nf);
    2862                 :            :   }
    2863         [ +  + ]:        545 :   if (pbnf) *pbnf = bnf;
    2864                 :        545 :   *pD = D;
    2865                 :        545 : }
    2866                 :            : 
    2867                 :            : /* compute the polynomial over Q of the Hilbert class field of
    2868                 :            :    Q(sqrt(D)) where D is a positive fundamental discriminant */
    2869                 :            : static GEN
    2870                 :        100 : quadhilbertreal(GEN D, long prec)
    2871                 :            : {
    2872                 :        100 :   pari_sp av = avma;
    2873                 :            :   long newprec;
    2874                 :            :   GEN bnf;
    2875                 :            :   VOLATILE GEN bnr, dtQ, data, nf, cyc, M;
    2876                 :            :   pari_timer ti;
    2877         [ -  + ]:        100 :   if (DEBUGLEVEL) timer_start(&ti);
    2878                 :            : 
    2879                 :            :   (void)&prec; /* prevent longjmp clobbering it */
    2880                 :            :   (void)&bnf;  /* prevent longjmp clobbering it, avoid warning due to
    2881                 :            :                 * quadray_init call : discards qualifiers from pointer type */
    2882                 :        100 :   quadray_init(&D, NULL, &bnf, prec);
    2883                 :        100 :   cyc = bnf_get_cyc(bnf);
    2884         [ -  + ]:        100 :   if (lg(cyc) == 1) { avma = av; return pol_x(0); }
    2885                 :            :   /* if the exponent of the class group is 2, use Genus Theory */
    2886         [ +  + ]:        100 :   if (equaliu(gel(cyc,1), 2)) return gerepileupto(av, GenusFieldQuadReal(D));
    2887                 :            : 
    2888                 :         90 :   bnr  = Buchray(bnf, gen_1, nf_INIT|nf_GEN);
    2889                 :         90 :   M = diagonal_shallow(bnr_get_cyc(bnr));
    2890                 :         90 :   dtQ = InitQuotient(M);
    2891                 :         90 :   nf  = bnf_get_nf(bnf);
    2892                 :            : 
    2893                 :            :   for(;;) {
    2894                 :         90 :     VOLATILE GEN pol = NULL;
    2895 [ -  + ][ #  # ]:         90 :     pari_CATCH(e_PREC) {
                 [ #  # ]
    2896                 :          0 :       prec += EXTRA_PREC;
    2897         [ #  # ]:          0 :       if (DEBUGLEVEL) pari_warn(warnprec, "quadhilbertreal", prec);
    2898                 :          0 :       bnr = bnrnewprec_shallow(bnr, prec);
    2899                 :          0 :       bnf = bnr_get_bnf(bnr);
    2900                 :          0 :       nf  = bnf_get_nf(bnf);
    2901                 :            :     } pari_TRY {
    2902                 :            :       /* find the modulus defining N */
    2903                 :            :       pari_timer T;
    2904         [ -  + ]:         90 :       if (DEBUGLEVEL) timer_start(&T);
    2905                 :         90 :       data = FindModulus(bnr, dtQ, &newprec);
    2906         [ -  + ]:         90 :       if (DEBUGLEVEL) timer_printf(&T,"FindModulus");
    2907         [ -  + ]:         90 :       if (!data)
    2908                 :            :       {
    2909                 :          0 :         long i, l = lg(M);
    2910                 :          0 :         GEN vec = cgetg(l, t_VEC);
    2911         [ #  # ]:          0 :         for (i = 1; i < l; i++)
    2912                 :            :         {
    2913                 :          0 :           GEN t = gcoeff(M,i,i);
    2914                 :          0 :           gcoeff(M,i,i) = gen_1;
    2915                 :          0 :           gel(vec,i) = bnrstark(bnr, M, prec);
    2916                 :          0 :           gcoeff(M,i,i) = t;
    2917                 :            :         }
    2918                 :          0 :         return gerepileupto(av, vec);
    2919                 :            :       }
    2920                 :            : 
    2921         [ +  - ]:         90 :       if (newprec > prec)
    2922                 :            :       {
    2923         [ -  + ]:         90 :         if (DEBUGLEVEL>1) err_printf("new precision: %ld\n", newprec);
    2924                 :         90 :         nf = nfnewprec_shallow(nf, newprec);
    2925                 :            :       }
    2926                 :         90 :       pol = AllStark(data, nf, 0, newprec);
    2927                 :         90 :     } pari_ENDCATCH;
    2928         [ +  - ]:         90 :     if (pol) {
    2929                 :         90 :       pol = makescind(nf, pol);
    2930                 :         90 :       return gerepileupto(av, polredbest(pol, 0));
    2931                 :            :     }
    2932                 :        100 :   }
    2933                 :            : }
    2934                 :            : 
    2935                 :            : /*******************************************************************/
    2936                 :            : /*                                                                 */
    2937                 :            : /*       Hilbert and Ray Class field using CM (Schertz)            */
    2938                 :            : /*                                                                 */
    2939                 :            : /*******************************************************************/
    2940                 :            : /* form^2 = 1 ? */
    2941                 :            : static int
    2942                 :        186 : hasexp2(GEN form)
    2943                 :            : {
    2944                 :        186 :   GEN a = gel(form,1), b = gel(form,2), c = gel(form,3);
    2945 [ +  - ][ +  + ]:        186 :   return !signe(b) || absi_equal(a,b) || equalii(a,c);
                 [ -  + ]
    2946                 :            : }
    2947                 :            : static int
    2948                 :        910 : uhasexp2(GEN form)
    2949                 :            : {
    2950                 :        910 :   long a = form[1], b = form[2], c = form[3];
    2951 [ +  + ][ +  + ]:        910 :   return !b || a == labs(b) || a == c;
                 [ +  - ]
    2952                 :            : }
    2953                 :            : 
    2954                 :            : GEN
    2955                 :        295 : qfbforms(GEN D)
    2956                 :            : {
    2957                 :        295 :   ulong d = itou(D), dover3 = d/3, t, b2, a, b, c, h;
    2958                 :        295 :   GEN L = cgetg((long)(sqrt(d) * log2(d)), t_VEC);
    2959                 :        295 :   b2 = b = (d&1); h = 0;
    2960         [ +  + ]:        295 :   if (!b) /* b = 0 treated separately to avoid special cases */
    2961                 :            :   {
    2962                 :        165 :     t = d >> 2; /* (b^2 - D) / 4*/
    2963         [ +  + ]:       2035 :     for (a=1; a*a<=t; a++)
    2964         [ +  + ]:       1870 :       if (c = t/a, t == c*a) gel(L,++h) = mkvecsmall3(a,0,c);
    2965                 :        165 :     b = 2; b2 = 4;
    2966                 :            :   }
    2967                 :            :   /* now b > 0, b = D (mod 2) */
    2968         [ +  + ]:       5675 :   for ( ; b2 <= dover3; b += 2, b2 = b*b)
    2969                 :            :   {
    2970                 :       5380 :     t = (b2 + d) >> 2; /* (b^2 - D) / 4*/
    2971                 :            :     /* b = a */
    2972         [ +  + ]:       5380 :     if (c = t/b, t == c*b) gel(L,++h) = mkvecsmall3(b,b,c);
    2973                 :            :     /* b < a < c */
    2974         [ +  + ]:    1365545 :     for (a = b+1; a*a < t; a++)
    2975         [ +  + ]:    1360165 :       if (c = t/a, t == c*a)
    2976                 :            :       {
    2977                 :        700 :         gel(L,++h) = mkvecsmall3(a, b,c);
    2978                 :        700 :         gel(L,++h) = mkvecsmall3(a,-b,c);
    2979                 :            :       }
    2980                 :            :     /* a = c */
    2981         [ +  + ]:       5380 :     if (a * a == t) gel(L,++h) = mkvecsmall3(a,b,a);
    2982                 :            :   }
    2983                 :        295 :   setlg(L,h+1); return L;
    2984                 :            : }
    2985                 :            : 
    2986                 :            : /* gcd(n, 24) */
    2987                 :            : static long
    2988                 :        186 : GCD24(long n)
    2989                 :            : {
    2990   [ +  +  +  -  :        186 :   switch(n % 24)
          +  -  +  -  -  
          -  +  -  +  -  
          -  -  +  -  +  
          -  -  -  +  -  
                      - ]
    2991                 :            :   {
    2992                 :          5 :     case 0: return 24;
    2993                 :          5 :     case 1: return 1;
    2994                 :          5 :     case 2: return 2;
    2995                 :          0 :     case 3: return 3;
    2996                 :         25 :     case 4: return 4;
    2997                 :          0 :     case 5: return 1;
    2998                 :         15 :     case 6: return 6;
    2999                 :          0 :     case 7: return 1;
    3000                 :          0 :     case 8: return 8;
    3001                 :          0 :     case 9: return 3;
    3002                 :         35 :     case 10: return 2;
    3003                 :          0 :     case 11: return 1;
    3004                 :         25 :     case 12: return 12;
    3005                 :          0 :     case 13: return 1;
    3006                 :          0 :     case 14: return 2;
    3007                 :          0 :     case 15: return 3;
    3008                 :         20 :     case 16: return 8;
    3009                 :          0 :     case 17: return 1;
    3010                 :         31 :     case 18: return 6;
    3011                 :          0 :     case 19: return 1;
    3012                 :          0 :     case 20: return 4;
    3013                 :          0 :     case 21: return 3;
    3014                 :         20 :     case 22: return 2;
    3015                 :          0 :     case 23: return 1;
    3016                 :        186 :     default: return 0;
    3017                 :            :   }
    3018                 :            : }
    3019                 :            : 
    3020                 :            : struct gpq_data {
    3021                 :            :   long p, q;
    3022                 :            :   GEN sqd; /* sqrt(D), t_REAL */
    3023                 :            :   GEN u, D;
    3024                 :            :   GEN pq, pq2; /* p*q, 2*p*q */
    3025                 :            :   GEN qfpq ; /* class of \P * \Q */
    3026                 :            : };
    3027                 :            : 
    3028                 :            : /* find P and Q two non principal prime ideals (above p <= q) such that
    3029                 :            :  *   cl(P) = cl(Q) if P,Q have order 2 in Cl(K).
    3030                 :            :  *   Ensure that e = 24 / gcd(24, (p-1)(q-1)) = 1 */
    3031                 :            : /* D t_INT, discriminant */
    3032                 :            : static void
    3033                 :         15 : init_pq(GEN D, struct gpq_data *T)
    3034                 :            : {
    3035                 :         15 :   const long Np = 6547; /* N.B. primepi(50000) = 5133 */
    3036                 :         15 :   const ulong maxq = 50000;
    3037                 :         15 :   GEN listp = cgetg(Np + 1, t_VECSMALL); /* primes p */
    3038                 :         15 :   GEN listP = cgetg(Np + 1, t_VEC); /* primeform(p) if of order 2, else NULL */
    3039                 :         15 :   GEN gcd24 = cgetg(Np + 1, t_VECSMALL); /* gcd(p-1, 24) */
    3040                 :            :   forprime_t S;
    3041                 :         15 :   long l = 1;
    3042                 :         15 :   double best = 0.;
    3043                 :            :   ulong q;
    3044                 :            : 
    3045                 :         15 :   u_forprime_init(&S, 2, ULONG_MAX);
    3046                 :         15 :   T->D = D;
    3047                 :         15 :   T->p = T->q = 0;
    3048                 :            :   for(;;)
    3049                 :            :   {
    3050                 :            :     GEN Q;
    3051                 :            :     long i, gcdq, mod;
    3052                 :            :     int order2, store;
    3053                 :            :     double t;
    3054                 :            : 
    3055                 :        371 :     q = u_forprime_next(&S);
    3056 [ +  + ][ -  + ]:        371 :     if (best > 0 && q >= maxq)
    3057                 :            :     {
    3058         [ #  # ]:          0 :       if (DEBUGLEVEL)
    3059                 :          0 :         pari_warn(warner,"possibly suboptimal (p,q) for D = %Ps", D);
    3060                 :          0 :       break;
    3061                 :            :     }
    3062         [ +  + ]:        371 :     if (kroiu(D, q) < 0) continue; /* inert */
    3063                 :        186 :     Q = redimag(primeform_u(D, q));
    3064         [ -  + ]:        186 :     if (is_pm1(gel(Q,1))) continue; /* Q | q is principal */
    3065                 :            : 
    3066                 :        186 :     store = 1;
    3067                 :        186 :     order2 = hasexp2(Q);
    3068                 :        186 :     gcd24[l] = gcdq = GCD24(q-1);
    3069                 :        186 :     mod = 24 / gcdq; /* mod must divide p-1 otherwise e > 1 */
    3070                 :        186 :     listp[l] = q;
    3071         [ +  + ]:        186 :     gel(listP,l) = order2 ? Q : NULL;
    3072                 :        186 :     t = (q+1)/(double)(q-1);
    3073         [ +  + ]:        441 :     for (i = 1; i < l; i++) /* try all (p, q), p < q in listp */
    3074                 :            :     {
    3075                 :        326 :       long p = listp[i], gcdp = gcd24[i];
    3076                 :            :       double b;
    3077                 :            :       /* P,Q order 2 => cl(Q) = cl(P) */
    3078 [ +  + ][ +  - ]:        326 :       if (order2 && gel(listP,i) && !gequal(gel(listP,i), Q)) continue;
                 [ +  - ]
    3079         [ +  + ]:        321 :       if (gcdp % gcdq == 0) store = 0; /* already a better one in the list */
    3080         [ +  + ]:        321 :       if ((p-1) % mod) continue;
    3081                 :            : 
    3082                 :         71 :       b = (t*(p+1)) / (p-1); /* (p+1)(q+1) / (p-1)(q-1) */
    3083         [ +  + ]:         71 :       if (b > best) {
    3084                 :         25 :         store = 0; /* (p,q) always better than (q,r) for r >= q */
    3085                 :         25 :         best = b; T->q = q; T->p = p;
    3086         [ -  + ]:         25 :         if (DEBUGLEVEL>2) err_printf("p,q = %ld,%ld\n", p, q);
    3087                 :            :       }
    3088                 :            :       /* won't improve with this q as largest member */
    3089         [ +  - ]:         71 :       if (best > 0) break;
    3090                 :            :     }
    3091                 :            :     /* if !store or (q,r) won't improve on current best pair, forget that q */
    3092 [ +  + ][ +  + ]:        186 :     if (store && t*t > best)
    3093         [ -  + ]:         30 :       if (++l >= Np) pari_err_BUG("quadhilbert (not enough primes)");
    3094         [ +  + ]:        186 :     if (!best) /* (p,q) with p < q always better than (q,q) */
    3095                 :            :     { /* try (q,q) */
    3096 [ +  + ][ +  - ]:         35 :       if (gcdq >= 12 && umodiu(D, q)) /* e = 1 and unramified */
    3097                 :            :       {
    3098                 :          5 :         double b = (t*q) / (q-1); /* q(q+1) / (q-1)^2 */
    3099         [ +  - ]:          5 :         if (b > best) {
    3100                 :          5 :           best = b; T->q = T->p = q;
    3101         [ -  + ]:          5 :           if (DEBUGLEVEL>2) err_printf("p,q = %ld,%ld\n", q, q);
    3102                 :            :         }
    3103                 :            :       }
    3104                 :            :     }
    3105                 :            :     /* If (p1+1)(q+1) / (p1-1)(q-1) <= best, we can no longer improve
    3106                 :            :      * even with best p : stop */
    3107         [ +  + ]:        186 :     if ((listp[1]+1)*t <= (listp[1]-1)*best) break;
    3108                 :        356 :   }
    3109         [ -  + ]:         15 :   if (DEBUGLEVEL>1)
    3110                 :          0 :     err_printf("(p, q) = %ld, %ld; gain = %f\n", T->p, T->q, 12*best);
    3111                 :         15 : }
    3112                 :            : 
    3113                 :            : static GEN
    3114                 :       2780 : gpq(GEN form, struct gpq_data *T)
    3115                 :            : {
    3116                 :       2780 :   pari_sp av = avma;
    3117                 :       2780 :   long a = form[1], b = form[2], c = form[3];
    3118                 :       2780 :   long p = T->p, q = T->q;
    3119                 :            :   GEN form2, w, z;
    3120                 :       2780 :   int fl, real = 0;
    3121                 :            : 
    3122                 :       2780 :   form2 = qficomp(T->qfpq, mkvec3s(a, -b, c));
    3123                 :            :   /* form2 and form yield complex conjugate roots : only compute for the
    3124                 :            :    * lexicographically smallest of the 2 */
    3125                 :       2780 :   fl = cmpis(gel(form2,1), a);
    3126         [ +  + ]:       2780 :   if (fl <= 0)
    3127                 :            :   {
    3128         [ +  + ]:       1405 :     if (fl < 0) return NULL;
    3129                 :         30 :     fl = cmpis(gel(form2,2), b);
    3130         [ +  - ]:         30 :     if (fl <= 0)
    3131                 :            :     {
    3132         [ -  + ]:         30 :       if (fl < 0) return NULL;
    3133                 :            :       /* form == form2 : real root */
    3134                 :         30 :       real = 1;
    3135                 :            :     }
    3136                 :            :   }
    3137                 :            : 
    3138         [ +  + ]:       1405 :   if (p == 2) { /* (a,b,c) = (1,1,0) mod 2 ? */
    3139 [ -  + ][ #  # ]:         55 :     if (a % q == 0 && (a & b & 1) && !(c & 1))
                 [ #  # ]
    3140                 :            :     { /* apply S : make sure that (a,b,c) represents odd values */
    3141                 :          0 :       lswap(a,c); b = -b;
    3142                 :            :     }
    3143                 :            :   }
    3144 [ +  + ][ +  + ]:       1405 :   if (a % p == 0 || a % q == 0)
    3145                 :            :   { /* apply T^k, look for c' = a k^2 + b k + c coprime to N */
    3146 [ +  + ][ +  + ]:        355 :     while (c % p == 0 || c % q == 0)
    3147                 :            :     {
    3148                 :         50 :       c += a + b;
    3149                 :         50 :       b += a << 1;
    3150                 :            :     }
    3151                 :        305 :     lswap(a, c); b = -b; /* apply S */
    3152                 :            :   }
    3153                 :            :   /* now (a,b,c) ~ form and (a,pq) = 1 */
    3154                 :            : 
    3155                 :            :   /* gcd(2a, u) = 2,  w = u mod 2pq, -b mod 2a */
    3156                 :       1405 :   w = Z_chinese(T->u, stoi(-b), T->pq2, utoipos(a << 1));
    3157                 :       1405 :   z = double_eta_quotient(utoipos(a), w, T->D, T->p, T->q, T->pq, T->sqd);
    3158 [ +  + ][ +  - ]:       1405 :   if (real && typ(z) == t_COMPLEX) z = gcopy(gel(z, 1));
    3159                 :       2780 :   return gerepileupto(av, z);
    3160                 :            : }
    3161                 :            : 
    3162                 :            : /* returns an equation for the Hilbert class field of Q(sqrt(D)), D < 0
    3163                 :            :  * fundamental discriminant */
    3164                 :            : static GEN
    3165                 :        295 : quadhilbertimag(GEN D)
    3166                 :            : {
    3167                 :            :   GEN L, P, Pi, Pr, qfp, u;
    3168                 :        295 :   pari_sp av = avma;
    3169                 :            :   long h, i, prec;
    3170                 :            :   struct gpq_data T;
    3171                 :            :   pari_timer ti;
    3172                 :            : 
    3173         [ -  + ]:        295 :   if (DEBUGLEVEL>1) timer_start(&ti);
    3174         [ +  - ]:        295 :   if (lgefint(D) == 3)
    3175         [ -  + ]:        295 :     switch (D[2]) { /* = |D|; special cases where e > 1 */
    3176                 :            :       case 3:
    3177                 :            :       case 4:
    3178                 :            :       case 7:
    3179                 :            :       case 8:
    3180                 :            :       case 11:
    3181                 :            :       case 19:
    3182                 :            :       case 43:
    3183                 :            :       case 67:
    3184                 :          0 :       case 163: return pol_x(0);
    3185                 :            :     }
    3186                 :        295 :   L = qfbforms(D);
    3187                 :        295 :   h = lg(L)-1;
    3188         [ +  + ]:        295 :   if ((1L << vals(h)) == h) /* power of 2 */
    3189                 :            :   { /* check whether > |Cl|/2 elements have order <= 2 ==> 2-elementary */
    3190                 :        280 :     long lim = (h>>1) + 1;
    3191         [ +  + ]:       1190 :     for (i=1; i <= lim; i++)
    3192         [ -  + ]:        910 :       if (!uhasexp2(gel(L,i))) break;
    3193         [ +  - ]:        280 :     if (i > lim) return GenusFieldQuadImag(D);
    3194                 :            :   }
    3195         [ -  + ]:         15 :   if (DEBUGLEVEL>1) timer_printf(&ti,"class number = %ld",h);
    3196                 :         15 :   init_pq(D, &T);
    3197                 :         15 :   qfp = primeform_u(D, T.p);
    3198                 :         15 :   T.pq =  muluu(T.p, T.q);
    3199                 :         15 :   T.pq2 = shifti(T.pq,1);
    3200         [ -  + ]:         15 :   if (T.p == T.q)
    3201                 :            :   {
    3202                 :          0 :     GEN qfbp2 = qficompraw(qfp, qfp);
    3203                 :          0 :     u = gel(qfbp2,2);
    3204                 :          0 :     T.u = modii(u, T.pq2);
    3205                 :          0 :     T.qfpq = redimag(qfbp2);
    3206                 :            :   }
    3207                 :            :   else
    3208                 :            :   {
    3209                 :         15 :     GEN qfq = primeform_u(D, T.q), bp = gel(qfp,2), bq = gel(qfq,2);
    3210                 :         15 :     T.u = Z_chinese(bp, bq, utoipos(T.p << 1), utoipos(T.q << 1));
    3211                 :            :     /* T.u = bp (mod 2p), T.u = bq (mod 2q) */
    3212                 :         15 :     T.qfpq = qficomp(qfp, qfq);
    3213                 :            :   }
    3214                 :            :   /* u modulo 2pq */
    3215                 :         15 :   prec = LOWDEFAULTPREC;
    3216                 :         15 :   Pr = cgetg(h+1,t_VEC);
    3217                 :         15 :   Pi = cgetg(h+1,t_VEC);
    3218                 :            :   for(;;)
    3219                 :            :   {
    3220                 :         25 :     long ex, exmax = 0, r1 = 0, r2 = 0;
    3221                 :         25 :     pari_sp av0 = avma;
    3222                 :         25 :     T.sqd = sqrtr_abs(itor(D, prec));
    3223         [ +  + ]:       2805 :     for (i=1; i<=h; i++)
    3224                 :            :     {
    3225                 :       2780 :       GEN s = gpq(gel(L,i), &T);
    3226         [ -  + ]:       2780 :       if (DEBUGLEVEL>3) err_printf("%ld ", i);
    3227         [ +  + ]:       2780 :       if (!s) continue;
    3228         [ +  + ]:       1405 :       if (typ(s) != t_COMPLEX) gel(Pr, ++r1) = s; /* real root */
    3229                 :       1375 :       else                     gel(Pi, ++r2) = s;
    3230         [ +  + ]:       1405 :       ex = gexpo(s); if (ex > 0) exmax += ex;
    3231                 :            :     }
    3232         [ -  + ]:         25 :     if (DEBUGLEVEL>1) timer_printf(&ti,"roots");
    3233                 :         25 :     setlg(Pr, r1+1);
    3234                 :         25 :     setlg(Pi, r2+1);
    3235                 :         25 :     P = roots_to_pol_r1(shallowconcat(Pr,Pi), 0, r1);
    3236                 :         25 :     P = grndtoi(P,&exmax);
    3237         [ -  + ]:         25 :     if (DEBUGLEVEL>1) timer_printf(&ti,"product, error bits = %ld",exmax);
    3238         [ +  + ]:         25 :     if (exmax <= -10) break;
    3239                 :         10 :     avma = av0; prec += (DEFAULTPREC-2) + nbits2extraprec(exmax);
    3240         [ -  + ]:         10 :     if (DEBUGLEVEL) pari_warn(warnprec,"quadhilbertimag",prec);
    3241                 :         10 :   }
    3242                 :        295 :   return gerepileupto(av,P);
    3243                 :            : }
    3244                 :            : 
    3245                 :            : GEN
    3246                 :        405 : quadhilbert(GEN D, long prec)
    3247                 :            : {
    3248                 :        405 :   GEN d = D;
    3249                 :        405 :   quadray_init(&d, NULL, NULL, 0);
    3250                 :        395 :   return (signe(d)>0)? quadhilbertreal(D,prec)
    3251         [ +  + ]:        395 :                      : quadhilbertimag(d);
    3252                 :            : }
    3253                 :            : 
    3254                 :            : /* return a vector of all roots of 1 in bnf [not necessarily quadratic] */
    3255                 :            : static GEN
    3256                 :         45 : getallrootsof1(GEN bnf)
    3257                 :            : {
    3258                 :         45 :   GEN T, u, nf = bnf_get_nf(bnf), tu;
    3259                 :         45 :   long i, n = bnf_get_tuN(bnf);
    3260                 :            : 
    3261         [ +  + ]:         45 :   if (n == 2) {
    3262                 :         35 :     long N = nf_get_degree(nf);
    3263                 :         35 :     return mkvec2(scalarcol_shallow(gen_m1, N),
    3264                 :            :                   scalarcol_shallow(gen_1, N));
    3265                 :            :   }
    3266                 :         10 :   tu = poltobasis(nf, bnf_get_tuU(bnf));
    3267                 :         10 :   T = zk_multable(nf, tu);
    3268                 :         10 :   u = cgetg(n+1, t_VEC); gel(u,1) = tu;
    3269         [ +  + ]:         40 :   for (i=2; i <= n; i++) gel(u,i) = ZM_ZC_mul(T, gel(u,i-1));
    3270                 :         45 :   return u;
    3271                 :            : }
    3272                 :            : /* assume bnr has the right conductor */
    3273                 :            : static GEN
    3274                 :         45 : get_lambda(GEN bnr)
    3275                 :            : {
    3276                 :         45 :   GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf), pol = nf_get_pol(nf);
    3277                 :         45 :   GEN f = gel(bnr_get_mod(bnr), 1), labas, lamodf, u;
    3278                 :         45 :   long a, b, f2, i, lu, v = varn(pol);
    3279                 :            : 
    3280                 :         45 :   f2 = 2 * itos(gcoeff(f,1,1));
    3281                 :         45 :   u = getallrootsof1(bnf); lu = lg(u);
    3282         [ +  + ]:        155 :   for (i=1; i<lu; i++)
    3283                 :        110 :     gel(u,i) = ZC_hnfrem(gel(u,i), f); /* roots of 1, mod f */
    3284         [ -  + ]:         45 :   if (DEBUGLEVEL>1)
    3285                 :          0 :     err_printf("quadray: looking for [a,b] != unit mod 2f\n[a,b] = ");
    3286         [ +  - ]:        105 :   for (a=0; a<f2; a++)
    3287         [ +  + ]:       1660 :     for (b=0; b<f2; b++)
    3288                 :            :     {
    3289                 :       1600 :       GEN la = deg1pol_shallow(stoi(a), stoi(b), v); /* ax + b */
    3290         [ +  + ]:       1600 :       if (umodiu(gnorm(mkpolmod(la, pol)), f2) != 1) continue;
    3291         [ -  + ]:        145 :       if (DEBUGLEVEL>1) err_printf("[%ld,%ld] ",a,b);
    3292                 :            : 
    3293                 :        145 :       labas = poltobasis(nf, la);
    3294                 :        145 :       lamodf = ZC_hnfrem(labas, f);
    3295         [ +  + ]:        305 :       for (i=1; i<lu; i++)
    3296         [ +  + ]:        260 :         if (ZV_equal(lamodf, gel(u,i))) break;
    3297         [ +  + ]:        145 :       if (i < lu) continue; /* la = unit mod f */
    3298         [ -  + ]:         45 :       if (DEBUGLEVEL)
    3299                 :            :       {
    3300         [ #  # ]:          0 :         if (DEBUGLEVEL>1) err_printf("\n");
    3301                 :          0 :         err_printf("lambda = %Ps\n",la);
    3302                 :            :       }
    3303                 :         45 :       return labas;
    3304                 :            :     }
    3305                 :          0 :   pari_err_BUG("get_lambda");
    3306                 :         45 :   return NULL;
    3307                 :            : }
    3308                 :            : 
    3309                 :            : static GEN
    3310                 :       6055 : to_approx(GEN nf, GEN a)
    3311                 :            : {
    3312                 :       6055 :   GEN M = nf_get_M(nf);
    3313                 :       6055 :   return gadd(gel(a,1), gmul(gcoeff(M,1,2),gel(a,2)));
    3314                 :            : }
    3315                 :            : /* Z-basis for a (over C) */
    3316                 :            : static GEN
    3317                 :       3005 : get_om(GEN nf, GEN a) {
    3318                 :       3005 :   return mkvec2(to_approx(nf,gel(a,2)),
    3319                 :       3005 :                 to_approx(nf,gel(a,1)));
    3320                 :            : }
    3321                 :            : 
    3322                 :            : /* Compute all elts in class group G = [|G|,c,g], c=cyclic factors, g=gens.
    3323                 :            :  * Set list[j + 1] = g1^e1...gk^ek where j is the integer
    3324                 :            :  *   ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */
    3325                 :            : static GEN
    3326                 :         45 : getallelts(GEN bnr)
    3327                 :            : {
    3328                 :            :   GEN nf, C, c, g, list, pows, gk;
    3329                 :            :   long lc, i, j, no;
    3330                 :            : 
    3331                 :         45 :   nf = bnr_get_nf(bnr);
    3332                 :         45 :   no = itos( bnr_get_no(bnr) );
    3333                 :         45 :   c = bnr_get_cyc(bnr);
    3334                 :         45 :   g = bnr_get_gen_nocheck(bnr); lc = lg(c)-1;
    3335                 :         45 :   list = cgetg(no+1,t_VEC);
    3336                 :         45 :   gel(list,1) = matid(nf_get_degree(nf)); /* (1) */
    3337         [ -  + ]:         45 :   if (!no) return list;
    3338                 :            : 
    3339                 :         45 :   pows = cgetg(lc+1,t_VEC);
    3340                 :         45 :   c = leafcopy(c); settyp(c, t_VECSMALL);
    3341         [ +  + ]:         90 :   for (i=1; i<=lc; i++)
    3342                 :            :   {
    3343                 :         45 :     long k = itos(gel(c,i));
    3344                 :         45 :     c[i] = k;
    3345                 :         45 :     gk = cgetg(k, t_VEC); gel(gk,1) = gel(g,i);
    3346         [ +  + ]:       2960 :     for (j=2; j<k; j++)
    3347                 :       2915 :       gel(gk,j) = idealmoddivisor(bnr, idealmul(nf, gel(gk,j-1), gel(gk,1)));
    3348                 :         45 :     gel(pows,i) = gk; /* powers of g[i] */
    3349                 :            :   }
    3350                 :            : 
    3351                 :         45 :   C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc];
    3352         [ -  + ]:         45 :   for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1];
    3353                 :            :   /* C[i] = c(k-i+1) * ... * ck */
    3354                 :            :   /* j < C[i+1] <==> j only involves g(k-i)...gk */
    3355                 :         45 :   i = 1;
    3356         [ +  + ]:       3005 :   for (j=1; j < C[1]; j++)
    3357                 :       2960 :     gel(list, j+1) = gmael(pows,lc,j);
    3358         [ -  + ]:         45 :   while(j<no)
    3359                 :            :   {
    3360                 :            :     long k;
    3361                 :            :     GEN a;
    3362         [ #  # ]:          0 :     if (j == C[i+1]) i++;
    3363                 :          0 :     a = gmael(pows,lc-i,j/C[i]);
    3364                 :          0 :     k = j%C[i] + 1;
    3365         [ #  # ]:          0 :     if (k > 1) a = idealmoddivisor(bnr, idealmul(nf, a, gel(list,k)));
    3366                 :          0 :     gel(list, ++j) = a;
    3367                 :            :   }
    3368                 :         45 :   return list;
    3369                 :            : }
    3370                 :            : 
    3371                 :            : /* x quadratic integer (approximate), recognize it. If error return NULL */
    3372                 :            : static GEN
    3373                 :       3050 : findbezk(GEN nf, GEN x)
    3374                 :            : {
    3375                 :       3050 :   GEN a,b, M = nf_get_M(nf), u = gcoeff(M,1,2);
    3376                 :            :   long ea, eb;
    3377                 :            : 
    3378                 :            :   /* u t_COMPLEX generator of nf.zk, write x ~ a + b u, a,b in Z */
    3379                 :       3050 :   b = grndtoi(mpdiv(imag_i(x), gel(u,2)), &eb);
    3380         [ -  + ]:       3050 :   if (eb > -20) return NULL;
    3381                 :       3050 :   a = grndtoi(mpsub(real_i(x), mpmul(b,gel(u,1))), &ea);
    3382         [ -  + ]:       3050 :   if (ea > -20) return NULL;
    3383         [ +  + ]:       3050 :   return signe(b)? coltoalg(nf, mkcol2(a,b)): a;
    3384                 :            : }
    3385                 :            : 
    3386                 :            : static GEN
    3387                 :         45 : findbezk_pol(GEN nf, GEN x)
    3388                 :            : {
    3389                 :         45 :   long i, lx = lg(x);
    3390                 :         45 :   GEN y = cgetg(lx,t_POL);
    3391         [ +  + ]:       3095 :   for (i=2; i<lx; i++)
    3392         [ -  + ]:       3050 :     if (! (gel(y,i) = findbezk(nf,gel(x,i))) ) return NULL;
    3393                 :         45 :   y[1] = x[1]; return y;
    3394                 :            : }
    3395                 :            : 
    3396                 :            : /* allow t_QF[IR], and t_VEC/t_COL with 3 components */
    3397                 :            : GEN
    3398                 :          0 : form_to_ideal(GEN x)
    3399                 :            : {
    3400                 :          0 :   long tx = typ(x);
    3401                 :            :   GEN b;
    3402 [ #  # ][ #  # ]:          0 :   if ((is_vec_t(tx) || lg(x) != 4)
    3403 [ #  # ][ #  # ]:          0 :        && tx != t_QFR && tx != t_QFI) pari_err_TYPE("form_to_ideal",x);
    3404         [ #  # ]:          0 :   b = negi(gel(x,2)); if (mpodd(b)) b = addis(b,1);
    3405                 :          0 :   return mkmat2( mkcol2(gel(x,1), gen_0),
    3406                 :            :                  mkcol2(shifti(b,-1), gen_1) );
    3407                 :            : }
    3408                 :            : 
    3409                 :            : /* P approximation computed at initial precision prec. Compute needed prec
    3410                 :            :  * to know P with 1 word worth of trailing decimals */
    3411                 :            : static long
    3412                 :          0 : get_prec(GEN P, long prec)
    3413                 :            : {
    3414                 :          0 :   long k = gprecision(P);
    3415         [ #  # ]:          0 :   if (k == 3) return precdbl(prec); /* approximation not trustworthy */
    3416                 :          0 :   k = prec - k; /* lost precision when computing P */
    3417         [ #  # ]:          0 :   if (k < 0) k = 0;
    3418                 :          0 :   k += nbits2prec(gexpo(P) + 128);
    3419         [ #  # ]:          0 :   if (k <= prec) k = precdbl(prec); /* dubious: old prec should have worked */
    3420                 :          0 :   return k;
    3421                 :            : }
    3422                 :            : 
    3423                 :            : /* Compute data for ellphist */
    3424                 :            : static GEN
    3425                 :       3005 : ellphistinit(GEN om, long prec)
    3426                 :            : {
    3427                 :       3005 :   GEN res,om1b,om2b, om1 = gel(om,1), om2 = gel(om,2);
    3428                 :            : 
    3429         [ -  + ]:       3005 :   if (gsigne(imag_i(gdiv(om1,om2))) < 0) { swap(om1,om2); om = mkvec2(om1,om2); }
    3430                 :       3005 :   om1b = gconj(om1);
    3431                 :       3005 :   om2b = gconj(om2); res = cgetg(4,t_VEC);
    3432                 :       3005 :   gel(res,1) = gdivgs(elleisnum(om,2,0,prec),12);
    3433                 :       3005 :   gel(res,2) = gdiv(PiI2(prec), gmul(om2, imag_i(gmul(om1b,om2))));
    3434                 :       3005 :   gel(res,3) = om2b; return res;
    3435                 :            : }
    3436                 :            : 
    3437                 :            : /* Computes log(phi^*(z,om)), using res computed by ellphistinit */
    3438                 :            : static GEN
    3439                 :       6010 : ellphist(GEN om, GEN res, GEN z, long prec)
    3440                 :            : {
    3441                 :       6010 :   GEN u = imag_i(gmul(z, gel(res,3)));
    3442                 :       6010 :   GEN zst = gsub(gmul(u, gel(res,2)), gmul(z,gel(res,1)));
    3443                 :       6010 :   return gsub(ellsigma(om,z,1,prec),gmul2n(gmul(z,zst),-1));
    3444                 :            : }
    3445                 :            : 
    3446                 :            : /* Computes phi^*(la,om)/phi^*(1,om) where (1,om) is an oriented basis of the
    3447                 :            :    ideal gf*gc^{-1} */
    3448                 :            : static GEN
    3449                 :       3005 : computeth2(GEN om, GEN la, long prec)
    3450                 :            : {
    3451                 :       3005 :   GEN p1,p2,res = ellphistinit(om,prec);
    3452                 :            : 
    3453                 :       3005 :   p1 = gsub(ellphist(om,res,la,prec), ellphist(om,res,gen_1,prec));
    3454                 :       3005 :   p2 = imag_i(p1);
    3455 [ +  - ][ -  + ]:       3005 :   if (gexpo(real_i(p1))>20 || gexpo(p2)> prec2nbits(minss(prec,realprec(p2)))-10)
    3456                 :          0 :     return NULL;
    3457                 :       3005 :   return gexp(p1,prec);
    3458                 :            : }
    3459                 :            : 
    3460                 :            : /* Computes P_2(X)=polynomial in Z_K[X] closest to prod_gc(X-th2(gc)) where
    3461                 :            :    the product is over the ray class group bnr.*/
    3462                 :            : static GEN
    3463                 :         45 : computeP2(GEN bnr, long prec)
    3464                 :            : {
    3465                 :         45 :   long clrayno, i, first = 1;
    3466                 :         45 :   pari_sp av=avma, av2;
    3467                 :         45 :   GEN listray, P0, P, lanum, la = get_lambda(bnr);
    3468                 :         45 :   GEN nf = bnr_get_nf(bnr), f = gel(bnr_get_mod(bnr), 1);
    3469                 :         45 :   listray = getallelts(bnr);
    3470                 :         45 :   clrayno = lg(listray)-1; av2 = avma;
    3471                 :            : PRECPB:
    3472         [ -  + ]:         45 :   if (!first)
    3473                 :            :   {
    3474         [ #  # ]:          0 :     if (DEBUGLEVEL) pari_warn(warnprec,"computeP2",prec);
    3475                 :          0 :     nf = gerepilecopy(av2, nfnewprec_shallow(checknf(bnr),prec));
    3476                 :            :   }
    3477                 :         45 :   first = 0; lanum = to_approx(nf,la);
    3478                 :         45 :   P = cgetg(clrayno+1,t_VEC);
    3479         [ +  + ]:       3050 :   for (i=1; i<=clrayno; i++)
    3480                 :            :   {
    3481                 :       3005 :     GEN om = get_om(nf, idealdiv(nf,f,gel(listray,i)));
    3482                 :       3005 :     GEN s = computeth2(om,lanum,prec);
    3483         [ -  + ]:       3005 :     if (!s) { prec = precdbl(prec); goto PRECPB; }
    3484                 :       3005 :     gel(P,i) = s;
    3485                 :            :   }
    3486                 :         45 :   P0 = roots_to_pol(P, 0);
    3487                 :         45 :   P = findbezk_pol(nf, P0);
    3488         [ -  + ]:         45 :   if (!P) { prec = get_prec(P0, prec); goto PRECPB; }
    3489                 :         45 :   return gerepilecopy(av, P);
    3490                 :            : }
    3491                 :            : 
    3492                 :            : #define nexta(a) (a>0 ? -a : 1-a)
    3493                 :            : static GEN
    3494                 :          0 : do_compo(GEN x, GEN y)
    3495                 :            : {
    3496                 :          0 :   long a, i, l = lg(y);
    3497                 :            :   GEN z;
    3498                 :          0 :   y = leafcopy(y); /* y := t^deg(y) y(#/t) */
    3499         [ #  # ]:          0 :   for (i = 2; i < l; i++)
    3500         [ #  # ]:          0 :     if (signe(gel(y,i))) gel(y,i) = monomial(gel(y,i), l-i-1, MAXVARN);
    3501         [ #  # ]:          0 :   for  (a = 0;; a = nexta(a))
    3502                 :            :   {
    3503         [ #  # ]:          0 :     if (a) x = gsubst(x, 0, gaddsg(a, pol_x(0)));
    3504                 :          0 :     z = gsubst(resultant(x,y), MAXVARN, pol_x(0));
    3505         [ #  # ]:          0 :     if (issquarefree(z)) return z;
    3506                 :          0 :   }
    3507                 :            : }
    3508                 :            : #undef nexta
    3509                 :            : 
    3510                 :            : static GEN
    3511                 :          0 : galoisapplypol(GEN nf, GEN s, GEN x)
    3512                 :            : {
    3513                 :          0 :   long i, lx = lg(x);
    3514                 :          0 :   GEN y = cgetg(lx,t_POL);
    3515                 :            : 
    3516         [ #  # ]:          0 :   for (i=2; i<lx; i++) gel(y,i) = galoisapply(nf,s,gel(x,i));
    3517                 :          0 :   y[1] = x[1]; return y;
    3518                 :            : }
    3519                 :            : /* x quadratic, write it as ua + v, u,v rational */
    3520                 :            : static GEN
    3521                 :          0 : findquad(GEN a, GEN x, GEN p)
    3522                 :            : {
    3523                 :            :   long tu, tv;
    3524                 :          0 :   pari_sp av = avma;
    3525                 :            :   GEN u,v;
    3526         [ #  # ]:          0 :   if (typ(x) == t_POLMOD) x = gel(x,2);
    3527         [ #  # ]:          0 :   if (typ(a) == t_POLMOD) a = gel(a,2);
    3528                 :          0 :   u = poldivrem(x, a, &v);
    3529                 :          0 :   u = simplify_shallow(u); tu = typ(u);
    3530                 :          0 :   v = simplify_shallow(v); tv = typ(v);
    3531         [ #  # ]:          0 :   if (!is_scalar_t(tu)) pari_err_TYPE("findquad", u);
    3532         [ #  # ]:          0 :   if (!is_scalar_t(tv)) pari_err_TYPE("findquad", v);
    3533                 :          0 :   x = deg1pol(v, u, varn(a));
    3534         [ #  # ]:          0 :   if (typ(x) == t_POL) x = gmodulo(x,p);
    3535                 :          0 :   return gerepileupto(av, x);
    3536                 :            : }
    3537                 :            : static GEN
    3538                 :          0 : findquad_pol(GEN p, GEN a, GEN x)
    3539                 :            : {
    3540                 :          0 :   long i, lx = lg(x);
    3541                 :          0 :   GEN y = cgetg(lx,t_POL);
    3542         [ #  # ]:          0 :   for (i=2; i<lx; i++) gel(y,i) = findquad(a, gel(x,i), p);
    3543                 :          0 :   y[1] = x[1]; return y;
    3544                 :            : }
    3545                 :            : static GEN
    3546                 :          0 : compocyclo(GEN nf, long m, long d)
    3547                 :            : {
    3548                 :          0 :   GEN sb,a,b,s,p1,p2,p3,res,polL,polLK,nfL, D = nf_get_disc(nf);
    3549                 :            :   long ell,vx;
    3550                 :            : 
    3551                 :          0 :   p1 = quadhilbertimag(D);
    3552                 :          0 :   p2 = polcyclo(m,0);
    3553         [ #  # ]:          0 :   if (d==1) return do_compo(p1,p2);
    3554                 :            : 
    3555         [ #  # ]:          0 :   ell = m&1 ? m : (m>>2);
    3556         [ #  # ]:          0 :   if (equalui(ell,D)) /* ell = |D| */
    3557                 :            :   {
    3558                 :          0 :     p2 = gcoeff(nffactor(nf,p2),1,1);
    3559                 :          0 :     return do_compo(p1,p2);
    3560                 :            :   }
    3561         [ #  # ]:          0 :   if (ell%4 == 3) ell = -ell;
    3562                 :            :   /* nf = K = Q(a), L = K(b) quadratic extension = Q(t) */
    3563                 :          0 :   polLK = quadpoly(stoi(ell)); /* relative polynomial */
    3564                 :          0 :   res = rnfequation2(nf, polLK);
    3565                 :          0 :   vx = nf_get_varn(nf);
    3566                 :          0 :   polL = gsubst(gel(res,1),0,pol_x(vx)); /* = charpoly(t) */
    3567                 :          0 :   a = gsubst(lift(gel(res,2)), 0,pol_x(vx));
    3568                 :          0 :   b = gsub(pol_x(vx), gmul(gel(res,3), a));
    3569                 :          0 :   nfL = nfinit(polL, DEFAULTPREC);
    3570                 :          0 :   p1 = gcoeff(nffactor(nfL,p1),1,1);
    3571                 :          0 :   p2 = gcoeff(nffactor(nfL,p2),1,1);
    3572                 :          0 :   p3 = do_compo(p1,p2); /* relative equation over L */
    3573                 :            :   /* compute non trivial s in Gal(L / K) */
    3574                 :          0 :   sb= gneg(gadd(b, RgX_coeff(polLK,1))); /* s(b) = Tr(b) - b */
    3575                 :          0 :   s = gadd(pol_x(vx), gsub(sb, b)); /* s(t) = t + s(b) - b */
    3576                 :          0 :   p3 = gmul(p3, galoisapplypol(nfL, s, p3));
    3577                 :          0 :   return findquad_pol(nf_get_pol(nf), a, p3);
    3578                 :            : }
    3579                 :            : 
    3580                 :            : /* I integral ideal in HNF. (x) = I, x small in Z ? */
    3581                 :            : static long
    3582                 :         45 : isZ(GEN I)
    3583                 :            : {
    3584                 :         45 :   GEN x = gcoeff(I,1,1);
    3585 [ +  - ][ -  + ]:         45 :   if (signe(gcoeff(I,1,2)) || !equalii(x, gcoeff(I,2,2))) return 0;
    3586         [ +  - ]:         45 :   return is_bigint(x)? -1: itos(x);
    3587                 :            : }
    3588                 :            : 
    3589                 :            : /* Treat special cases directly. return NULL if not special case */
    3590                 :            : static GEN
    3591                 :         45 : treatspecialsigma(GEN bnr)
    3592                 :            : {
    3593                 :         45 :   GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf);
    3594                 :         45 :   GEN f = gel(bnr_get_mod(bnr), 1),  D = nf_get_disc(nf);
    3595                 :            :   GEN p1, p2;
    3596                 :         45 :   long Ds, fl, tryf, i = isZ(f);
    3597                 :            : 
    3598         [ -  + ]:         45 :   if (i == 1) return quadhilbertimag(D); /* f = 1 */
    3599                 :            : 
    3600         [ -  + ]:         45 :   if (equaliu(D,3)) /* Q(j) */
    3601                 :            :   {
    3602 [ #  # ][ #  # ]:          0 :     if (i == 4 || i == 5 || i == 7) return polcyclo(i,0);
                 [ #  # ]
    3603 [ #  # ][ #  # ]:          0 :     if (!equaliu(gcoeff(f,1,1),9) || !equaliu(Q_content(f),3)) return NULL;
    3604                 :            :     /* f = P_3^3 */
    3605                 :          0 :     p1 = mkpolmod(bnf_get_tuU(bnf), nf_get_pol(nf));
    3606                 :          0 :     return gadd(monomial(gen_1,3,0), p1); /* x^3+j */
    3607                 :            :   }
    3608         [ +  + ]:         45 :   if (equaliu(D,4)) /* Q(i) */
    3609                 :            :   {
    3610 [ +  - ][ -  + ]:         10 :     if (i == 3 || i == 5) return polcyclo(i,0);
    3611         [ +  - ]:         10 :     if (i != 4) return NULL;
    3612                 :          0 :     p1 = mkpolmod(bnf_get_tuU(bnf), nf_get_pol(nf));
    3613                 :          0 :     return gadd(monomial(gen_1,2,0), p1); /* x^2+i */
    3614                 :            :   }
    3615                 :         35 :   Ds = smodis(D,48);
    3616         [ +  - ]:         35 :   if (i)
    3617                 :            :   {
    3618 [ +  + ][ -  + ]:         35 :     if (i==2 && Ds%16== 8) return compocyclo(nf, 4,1);
    3619 [ +  + ][ -  + ]:         35 :     if (i==3 && Ds% 3== 1) return compocyclo(nf, 3,1);
    3620 [ -  + ][ #  # ]:         35 :     if (i==4 && Ds% 8== 1) return compocyclo(nf, 4,1);
    3621 [ -  + ][ #  # ]:         35 :     if (i==6 && Ds   ==40) return compocyclo(nf,12,1);
    3622                 :         35 :     return NULL;
    3623                 :            :   }
    3624                 :            : 
    3625                 :          0 :   p1 = gcoeff(f,1,1); /* integer > 0 */
    3626         [ #  # ]:          0 :   if (is_bigint(p1)) return NULL;
    3627                 :          0 :   tryf = p1[2];
    3628                 :          0 :   p2 = gcoeff(f,2,2); /* integer > 0 */
    3629         [ #  # ]:          0 :   if (is_pm1(p2)) fl = 0;
    3630                 :            :   else {
    3631 [ #  # ][ #  # ]:          0 :     if (Ds % 16 != 8 || !equaliu(Q_content(f),2)) return NULL;
    3632                 :          0 :     fl = 1; tryf >>= 1;
    3633                 :            :   }
    3634 [ #  # ][ #  # ]:          0 :   if (tryf <= 3 || umodiu(D, tryf) || !uisprime(tryf)) return NULL;
                 [ #  # ]
    3635         [ #  # ]:          0 :   if (fl) tryf <<= 2;
    3636                 :         45 :   return compocyclo(nf,tryf,2);
    3637                 :            : }
    3638                 :            : 
    3639                 :            : GEN
    3640                 :         75 : quadray(GEN D, GEN f, long prec)
    3641                 :            : {
    3642                 :            :   GEN bnr, y, bnf;
    3643                 :         75 :   pari_sp av = avma;
    3644                 :            : 
    3645         [ +  + ]:         75 :   if (isint1(f)) return quadhilbert(D, prec);
    3646                 :         50 :   quadray_init(&D, f, &bnf, prec);
    3647                 :         50 :   bnr = Buchray(bnf, f, nf_INIT|nf_GEN);
    3648         [ -  + ]:         50 :   if (is_pm1(bnr_get_no(bnr))) { avma = av; return pol_x(0); }
    3649         [ +  + ]:         50 :   if (signe(D) > 0)
    3650                 :          5 :     y = bnrstark(bnr,NULL,prec);
    3651                 :            :   else
    3652                 :            :   {
    3653                 :         45 :     bnr = gel(bnrconductor(bnr,NULL,2), 2);
    3654                 :         45 :     y = treatspecialsigma(bnr);
    3655         [ +  - ]:         45 :     if (!y) y = computeP2(bnr, prec);
    3656                 :            :   }
    3657                 :         65 :   return gerepileupto(av, y);
    3658                 :            : }

Generated by: LCOV version 1.9