Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

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

LCOV - code coverage report
Current view: top level - basemath - quad.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 631 684 92.3 %
Date: 2022-07-03 07:33:15 Functions: 52 54 96.3 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_arith
      19             : 
      20             : /*********************************************************************/
      21             : /**                                                                 **/
      22             : /**                    FUNDAMENTAL DISCRIMINANTS                    **/
      23             : /**                                                                 **/
      24             : /*********************************************************************/
      25             : static long
      26        1407 : fa_isfundamental(GEN F)
      27             : {
      28        1407 :   GEN P = gel(F,1), E = gel(F,2);
      29        1407 :   long i, s, l = lg(P);
      30             : 
      31        1407 :   if (l == 1) return 1;
      32        1400 :   s = signe(gel(P,1)); /* = signe(x) */
      33        1400 :   if (!s) return 0;
      34        1393 :   if (s < 0) { l--; P = vecslice(P,2,l); E = vecslice(E,2,l); }
      35        1393 :   if (l == 1) return 0;
      36        1386 :   if (!absequaliu(gel(P,1), 2))
      37         686 :     i = 1; /* need x = 1 mod 4 */
      38             :   else
      39             :   {
      40         700 :     i = 2;
      41         700 :     switch(itou(gel(E,1)))
      42             :     {
      43         182 :       case 2: s = -s; break; /* need x/4 = 3 mod 4 */
      44          84 :       case 3: s = 0; break; /* no condition mod 4 */
      45         434 :       default: return 0;
      46             :     }
      47             :   }
      48        1974 :   for(; i < l; i++)
      49             :   {
      50        1190 :     if (!equali1(gel(E,i))) return 0;
      51        1022 :     if (s && Mod4(gel(P,i)) == 3) s = -s;
      52             :   }
      53         784 :   return s >= 0;
      54             : }
      55             : long
      56       20482 : isfundamental(GEN x)
      57             : {
      58       20482 :   if (typ(x) != t_INT)
      59             :   {
      60        1407 :     pari_sp av = avma;
      61        1407 :     long v = fa_isfundamental(check_arith_all(x,"isfundamental"));
      62        1407 :     return gc_long(av,v);
      63             :   }
      64       19075 :   return Z_isfundamental(x);
      65             : }
      66             : 
      67             : /* x fundamental ? */
      68             : long
      69       16652 : uposisfundamental(ulong x)
      70             : {
      71       16652 :   ulong r = x & 15; /* x mod 16 */
      72       16652 :   if (!r) return 0;
      73       15840 :   switch(r & 3)
      74             :   { /* x mod 4 */
      75        3431 :     case 0: return (r == 4)? 0: uissquarefree(x >> 2);
      76        5958 :     case 1: return uissquarefree(x);
      77        6451 :     default: return 0;
      78             :   }
      79             : }
      80             : /* -x fundamental ? */
      81             : long
      82       29337 : unegisfundamental(ulong x)
      83             : {
      84       29337 :   ulong r = x & 15; /* x mod 16 */
      85       29337 :   if (!r) return 0;
      86       27692 :   switch(r & 3)
      87             :   { /* x mod 4 */
      88        6938 :     case 0: return (r == 12)? 0: uissquarefree(x >> 2);
      89       10425 :     case 3: return uissquarefree(x);
      90       10329 :     default: return 0;
      91             :   }
      92             : }
      93             : long
      94       25053 : sisfundamental(long x)
      95       25053 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
      96             : 
      97             : long
      98       19642 : Z_isfundamental(GEN x)
      99             : {
     100             :   long r;
     101       19642 :   switch(lgefint(x))
     102             :   {
     103           7 :     case 2: return 0;
     104        9233 :     case 3: return signe(x) < 0? unegisfundamental(x[2])
     105       26858 :                                : uposisfundamental(x[2]);
     106             :   }
     107        2010 :   r = mod16(x);
     108        2010 :   if (!r) return 0;
     109        1884 :   if ((r & 3) == 0)
     110             :   {
     111             :     pari_sp av;
     112         376 :     r >>= 2; /* |x|/4 mod 4 */
     113         376 :     if (signe(x) < 0) r = 4-r;
     114         376 :     if (r == 1) return 0;
     115         250 :     av = avma;
     116         250 :     r = Z_issquarefree( shifti(x,-2) );
     117         250 :     return gc_long(av, r);
     118             :   }
     119        1508 :   r &= 3; /* |x| mod 4 */
     120        1508 :   if (signe(x) < 0) r = 4-r;
     121        1508 :   return (r==1) ? Z_issquarefree(x) : 0;
     122             : }
     123             : 
     124             : static GEN
     125        2058 : fa_quaddisc(GEN f)
     126             : {
     127        2058 :   GEN P = gel(f,1), E = gel(f,2), s = gen_1;
     128        2058 :   long i, l = lg(P);
     129        6783 :   for (i = 1; i < l; i++) /* possibly including -1 */
     130        4725 :     if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
     131        2058 :   if (Mod4(s) > 1) s = shifti(s,2);
     132        2058 :   return s;
     133             : }
     134             : 
     135             : GEN
     136        2912 : quaddisc(GEN x)
     137             : {
     138        2912 :   const pari_sp av = avma;
     139        2912 :   long tx = typ(x);
     140             :   GEN F;
     141        2912 :   if (is_rational_t(tx)) F = factor(x);
     142             :   else
     143             :   {
     144        1407 :     F = check_arith_all(x,"quaddisc");
     145        1407 :     if (tx == t_VEC && typ(gel(x,1)) == t_INT
     146        1407 :                     && Z_issquarefree_fact(gel(x,2)))
     147             :     {
     148         854 :       x = gel(x,1);
     149         854 :       return (Mod4(x) > 1)? shifti(x, 2): icopy(x);
     150             :     }
     151             :   }
     152        2058 :   return gerepileuptoint(av, fa_quaddisc(F));
     153             : }
     154             : 
     155             : 
     156             : /***********************************************************************/
     157             : /**                                                                   **/
     158             : /**         FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS)         **/
     159             : /**                                                                   **/
     160             : /***********************************************************************/
     161             : /* replace f by f * [u,1; 1,0] */
     162             : static void
     163    22156268 : update_f(GEN f, GEN u)
     164             : {
     165    22156268 :   GEN a = gcoeff(f,1,1), b = gcoeff(f,1,2);
     166    22156268 :   GEN c = gcoeff(f,2,1), d = gcoeff(f,2,2);
     167    22156268 :   gcoeff(f,1,1) = addmulii(b, u,a); gcoeff(f,1,2) = a;
     168    22156268 :   gcoeff(f,2,1) = addmulii(d, u,c); gcoeff(f,2,2) = c;
     169    22156268 : }
     170             : 
     171             : /* f is a vector of matrices and i an index whose bits give the non-zero
     172             :  * entries; the product of the non zero entries is the actual result.
     173             :  * if i odd, f[1] is implicitely [f[1],1;1,0] */
     174             : 
     175             : static long
     176    22317344 : update_fm(GEN f, GEN a, long i)
     177             : {
     178             : #ifdef LONG_IS_64BIT
     179    19129152 :   const long LIM = 10;
     180             : #else
     181     3188192 :   const long LIM = 18;
     182             : #endif
     183    22317344 :   pari_sp av = avma;
     184             :   long k, v;
     185             :   GEN u;
     186    22317344 :   if (!odd(i)) { gel(f,1) = a; return i+1; }
     187    22236806 :   u = gel(f, 1);
     188    22236806 :   if (typ(u) == t_INT) /* [u,1;1,0] * [a,1;1,0] */
     189       80538 :   { gel(f,1) = mkmat22(addiu(mulii(a, u), 1), u, a, gen_1); return i; }
     190    22156268 :   update_f(u, a); if (lgefint(gcoeff(u,1,1)) < LIM) return i;
     191       80538 :   v = vals(i+1); gel(f,1) = gen_0;
     192      161003 :   for (k = 1; k < v; k++) { u = ZM2_mul(gel(f,k+1), u); gel(f,k+1) = gen_0; }
     193       80538 :   if (v != 1) u = gerepileupto(av, u);
     194       80538 :   gel(f,v+1) = u; return i+1;
     195             : }
     196             : /* \prod f[j]; if first only return column 1 */
     197             : static GEN
     198       69041 : prod_fm(GEN f, long i, long first)
     199             : {
     200       69041 :   long k, v = vals(i) + 1;
     201       69041 :   GEN u = gel(f, v);
     202             :   /* i a power of 2: f[1] can't be a t_INT */
     203       69041 :   if (!(i >>= v)) return first? gel(u,1): u;
     204         138 :   for (k = v+1; i; i >>= 1, k++)
     205         118 :     if (odd(i))
     206             :     {
     207          73 :       GEN w = gel(f,k);
     208          73 :       switch(typ(u))
     209             :       {
     210           0 :         case t_INT: update_f(w, u);
     211           0 :           u = first? gel(w,1): w; break;
     212           6 :         case t_COL: /* implies 'first' */
     213           6 :           u = ZM_ZC_mul(w, u); break;
     214          67 :         default: /* t_MAT */
     215          67 :           u = first? ZM_ZC_mul(w, gel(u,1)): ZM2_mul(w, u); break;
     216             :       }
     217          45 :     }
     218          20 :   return u;
     219             : }
     220             : 
     221             : GEN
     222       69048 : quadunit0(GEN x, long v)
     223             : {
     224       69048 :   GEN y = quadunit(x);
     225       69041 :   if (v==-1) v = fetch_user_var("w");
     226       69041 :   setvarn(gel(y,1), v); return y;
     227             : }
     228             : 
     229             : struct uimod { GEN N, T; };
     230             : static GEN
     231         973 : ui_pow(void *E, GEN x, GEN n)
     232         973 : { struct uimod *S = (struct uimod*)E; return FpXQ_pow(x, n, S->T, S->N); }
     233             : static int
     234         994 : ui_equal1(GEN x) { return degpol(x) < 1; }
     235             : static const struct bb_group
     236             : ui_group={ NULL,ui_pow,NULL,NULL,NULL,ui_equal1,NULL};
     237             : 
     238             : static void
     239          77 : quadunit_mod(GEN D, GEN d, GEN N, GEN *pu, GEN *pv)
     240             : {
     241             :   GEN u1, u2, v1, v2, p, q, q1, u, v;
     242          77 :   int m = mpodd(D);
     243          77 :   pari_sp av = avma;
     244          77 :   p = (mpodd(d) == m)? d: subiu(d, 1);
     245          77 :   u1 = negi(p); u2 = gen_2;
     246          77 :   v1 = gen_1; v2 = gen_0; q = gen_2;
     247          77 :   q1 = shifti(subii(D, sqri(p)), -1);
     248             :   for(;;)
     249         105 :   {
     250         182 :     GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
     251         182 :     p = subii(d, r);
     252         182 :     if (equalii(p1, p) && signe(v2))
     253             :     { /* even period */
     254           7 :       u = addmulii(sqri(u2), D, sqri(v2));
     255           7 :       v = shifti(mulii(u2,v2), 1);
     256           7 :       break;
     257             :     }
     258         175 :     t = Fp_addmul(u1, A, u2, N); u1 = u2; u2 = t;
     259         175 :     t = Fp_addmul(v1, A, v2, N); v1 = v2; v2 = t;
     260         175 :     t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
     261         175 :     if (equalii(q, t))
     262             :     { /* odd period */
     263          70 :       u = addmulii(mulii(u1,u2), D, mulii(v1,v2));
     264          70 :       v = addmulii(mulii(u1,v2), u2, v1);
     265          70 :       break;
     266             :     }
     267         105 :     if (gc_needed(av, 2))
     268             :     {
     269           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadunit_mod");
     270           0 :       gerepileall(av, 7, &p, &u1,&u2,&v1,&v2, &q,&q1);
     271             :     }
     272             :   }
     273          77 :   *pu = modii(u, N);
     274          77 :   *pv = modii(v, N); if (m) *pu = Fp_sub(*pu, *pv, N);
     275          77 : }
     276             : GEN
     277           0 : quadunit_basecase(GEN D)
     278             : {
     279           0 :   pari_sp av0 = avma;
     280             :   GEN d, u1, u2, v1, v2, p, q, q1, u, v, a, b, c;
     281             :   int m;
     282           0 :   long first = 1;
     283           0 :   check_quaddisc_real(D, NULL, "quadunit");
     284           0 :   m = mpodd(D); d = sqrtremi(D, &a);
     285           0 :   p = d; q1 = shifti(a, -1); q = gen_2;
     286           0 :   if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
     287           0 :   u1 = gen_2; u2 = p;
     288           0 :   v1 = gen_0; v2 = gen_1;
     289             :   for(;;)
     290           0 :   {
     291           0 :     GEN t = q;
     292           0 :     if (first) { first = 0; q = q1; }
     293             :     else
     294             :     {
     295           0 :       GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
     296           0 :       p = subii(d, r);
     297           0 :       if (equalii(p1, p)) /* even period */
     298           0 :       { a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break; }
     299           0 :       r = addmulii(u1, A, u2); u1 = u2; u2 = r;
     300           0 :       r = addmulii(v1, A, v2); v1 = v2; v2 = r;
     301           0 :       q = submulii(q1, A, subii(p, p1));
     302             :     }
     303           0 :     q1 = t;
     304           0 :     if (equalii(q, t))
     305             :     { /* odd period */
     306           0 :       a = mulii(u1, u2); b = mulii(v1, v2);
     307           0 :       c = mulii(addii(u1, v1), addii(u2, v2)); break;
     308             :     }
     309             :   }
     310           0 :   u = diviiexact(addmulii(a, D, b), q);
     311           0 :   v = diviiexact(subii(c, addii(a, b)), q);
     312           0 :   if (m == 1) u = subii(u, v);
     313           0 :   u = shifti(u, -1);
     314           0 :   return gerepilecopy(av0, mkquad(quadpoly_i(D), u, v));
     315             : }
     316             : 
     317             : GEN
     318       69048 : quadunit(GEN D)
     319             : {
     320       69048 :   pari_sp av0 = avma, av;
     321             :   GEN f, d, p, q, q1, u, v, a, b, c;
     322             :   int m;
     323       69048 :   long i = 0;
     324       69048 :   check_quaddisc_real(D, NULL, "quadunit");
     325       69041 :   m = mpodd(D); d = sqrtremi(D, &a); av = avma;
     326       69041 :   p = d; q1 = shifti(a, -1); q = gen_2;
     327       69041 :   if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
     328       69041 :   f = zerovec(2 + (expi(D)>>1));
     329       69041 :   gel(f,1) = mkmat22(p, gen_2, gen_1, gen_0);
     330             :   for(;;)
     331    22369158 :   {
     332    22438199 :     GEN t = q, u1,u2, v1,v2;
     333    22438199 :     if (!i) { i = 1; q = q1; }
     334             :     else
     335             :     {
     336    22369158 :       GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
     337    22369158 :       p = subii(d, r);
     338    22369158 :       if (equalii(p1, p))
     339             :       { /* even period */
     340       51814 :         f = prod_fm(f, i, 1); u2 = gel(f,1); v2 = gel(f,2);
     341       51814 :         a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break;
     342             :       }
     343    22317344 :       i = update_fm(f, A, i);
     344    22317344 :       q = submulii(q1, A, subii(p, p1));
     345             :     }
     346    22386385 :     q1 = t;
     347    22386385 :     if (equalii(q, t))
     348             :     { /* odd period */
     349       17227 :       f = prod_fm(f, i, 0);
     350       17227 :       u2 = gcoeff(f,1,1); u1 = gcoeff(f,1,2); a = mulii(u1, u2);
     351       17227 :       v2 = gcoeff(f,2,1); v1 = gcoeff(f,2,2); b = mulii(v1, v2);
     352       17227 :       c = mulii(addii(u1, v1), addii(u2, v2)); break;
     353             :     }
     354    22369158 :     if (gc_needed(av, 2))
     355             :     {
     356          83 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadunit (%ld)", i);
     357          83 :       gerepileall(av, 4, &p, &f, &q,&q1);
     358             :     }
     359             :   }
     360       69041 :   u = diviiexact(addmulii(a, D, b), q);
     361       69041 :   v = diviiexact(subii(c, addii(a, b)), q);
     362       69041 :   if (m == 1) u = subii(u, v);
     363       69041 :   u = shifti(u, -1);
     364       69041 :   return gerepilecopy(av0, mkquad(quadpoly_i(D), u, v));
     365             : }
     366             : /* D > 0, d = sqrti(D) */
     367             : static GEN
     368       69097 : quadunit_q(GEN D, GEN d, long *pN)
     369             : {
     370       69097 :   pari_sp av = avma;
     371             :   GEN p, q, q1;
     372       69097 :   long first = 1;
     373       69097 :   p = (Mod2(d) == Mod2(D))? d: subiu(d, 1);
     374       69097 :   q = gen_2;
     375       69097 :   q1 = shifti(subii(D, sqri(p)), -1);
     376             :   for(;;)
     377     1000930 :   {
     378     1070027 :     GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
     379     1070027 :     p = subii(d, r);
     380     1087317 :     if (!first && equalii(p1, p)) { *pN = 1; return q; } /* even period */
     381     1018220 :     first = 0;
     382     1018220 :     t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
     383     1018220 :     if (equalii(q, t)) { *pN = -1; return q; } /* odd period */
     384     1000930 :     if (gc_needed(av, 2))
     385             :     {
     386           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadunitnorm");
     387           0 :       gerepileall(av, 3, &p, &q, &q1);
     388             :     }
     389             :   }
     390             : }
     391             : long
     392       69027 : quadunitnorm(GEN D)
     393             : {
     394       69027 :   pari_sp av = avma;
     395             :   long s, r;
     396       69027 :   check_quaddisc(D, &s, &r, "quadunitnorm");
     397       69020 :   if (s < 0) return 1;
     398       69020 :   (void)quadunit_q(D, sqrti(D), &s); return gc_long(av, s);
     399             : }
     400             : GEN
     401         119 : quadunitindex(GEN D, GEN N)
     402             : {
     403         119 :   pari_sp av = avma, av2;
     404             :   GEN y, u, v, q, d, P, E, F, a, faH, H;
     405             :   struct uimod S;
     406             :   long r, s;
     407             : 
     408         119 :   check_quaddisc(D, &s, &r, "quadunitindex");
     409         112 :   if ((F = check_arith_pos(N,"quadunitindex")))
     410          35 :     N = typ(N) == t_VEC? gel(N,1): factorback(F);
     411         105 :   if (equali1(N)) return gen_1;
     412          98 :   if (s < 0) switch(itos_or_0(D)) {
     413           7 :     case -3: return utoipos(3);
     414           7 :     case -4: return utoipos(2);
     415           7 :     default: return gen_1;
     416             :   }
     417          77 :   if (!F) F = Z_factor(N);
     418          77 :   d = sqrti(D); av2 = avma;
     419          77 :   q = gerepileuptoint(av2, quadunit_q(D, d, &s));
     420          77 :   if (mpodd(N) && equali1(gcdii(q, N)))
     421             :   {
     422           7 :     quadunit_mod(D, d, N, &u, &v);
     423           7 :     q = Fp_inv(shifti(q, 1), N);
     424           7 :     u = Fp_mul(u, q, N);
     425           7 :     v = Fp_mul(v, q, N); v = modii(shifti(v, 1), N);
     426             :   }
     427             :   else
     428             :   {
     429          70 :     GEN M = shifti(mulii(q, N), 1);
     430          70 :     quadunit_mod(D, d, M, &u, &v);
     431          70 :     u = diviiexact(u, q);
     432          70 :     v = diviiexact(v, q); u = shifti(u,-1);
     433             :   }
     434             :   /* fundamental unit = y mod N */
     435          77 :   S.N = N; S.T = quadpoly_i(D); y = deg1pol_shallow(v, u, 0);
     436          77 :   H = quadclassnoF_fact(D, gel(F,1), gel(F,2));
     437          77 :   a = Z_smoothen(H, gel(F,1), &P, &E);
     438          77 :   faH = mkmat2(P, E);
     439          77 :   if (a) faH = merge_factor(Z_factor(a), faH,(void*)&cmpii,cmp_nodata);
     440          77 :   y = gen_order(y, mkvec2(H, faH), (void*)&S, &ui_group);
     441          77 :   return gerepileupto(av, y);
     442             : }
     443             : 
     444             : GEN
     445          42 : quadregulator(GEN x, long prec)
     446             : {
     447          42 :   pari_sp av = avma, av2;
     448             :   GEN R, rsqd, u, v, sqd;
     449             :   long r, e;
     450             : 
     451          42 :   check_quaddisc_real(x, &r, "quadregulator");
     452          42 :   sqd = sqrti(x);
     453          42 :   rsqd = gsqrt(x,prec); av2 = avma;
     454          42 :   e = 0; R = real2n(1, prec); u = utoi(r); v = gen_2;
     455             :   for(;;)
     456          49 :   {
     457          91 :     GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
     458          91 :     GEN v1 = divii(subii(x,sqri(u1)),v);
     459          91 :     if (equalii(v,v1)) { R = mulrr(sqrr(R), divri(addir(u1,rsqd),v)); break; }
     460          63 :     if (equalii(u,u1)) { R = sqrr(R); break; }
     461          49 :     R = mulrr(R, divri(addir(u1,rsqd),v));
     462          49 :     e += expo(R); setexpo(R,0);
     463          49 :     u = u1; v = v1;
     464          49 :     if (e & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
     465          49 :     if (gc_needed(av2,2))
     466             :     {
     467           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
     468           0 :       gerepileall(av2,3, &R,&u,&v);
     469             :     }
     470             :   }
     471          42 :   R = divri(R, v); e = 2*e - 1;
     472             :   /* avoid loss of accuracy */
     473          42 :   if (!((e + expo(R)) & ~EXPOBITS)) { setexpo(R, e + expo(R)); e = 0; }
     474          42 :   R = logr_abs(R);
     475          42 :   if (e) R = addrr(R, mulsr(e, mplog2(prec)));
     476          42 :   return gerepileuptoleaf(av, R);
     477             : }
     478             : 
     479             : /*************************************************************************/
     480             : /**                                                                     **/
     481             : /**                            CLASS NUMBER                             **/
     482             : /**                                                                     **/
     483             : /*************************************************************************/
     484             : 
     485             : int
     486     8243821 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
     487             : 
     488     9198757 : static GEN qfi_pow(void *E, GEN f, GEN n)
     489     9198757 : { return E? nupow(f,n,(GEN)E): qfbpow_i(f,n); }
     490     6350443 : static GEN qfi_comp(void *E, GEN f, GEN g)
     491     6350443 : { return E? nucomp(f,g,(GEN)E): qfbcomp_i(f,g); }
     492             : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
     493             :                                          gidentical,qfb_equal1,NULL};
     494             : 
     495             : GEN
     496     3022086 : qfi_order(GEN q, GEN o)
     497     3022086 : { return gen_order(q, o, NULL, &qfi_group); }
     498             : 
     499             : GEN
     500           0 : qfi_log(GEN a, GEN g, GEN o)
     501           0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
     502             : 
     503             : GEN
     504      646322 : qfi_Shanks(GEN a, GEN g, long n)
     505             : {
     506      646322 :   pari_sp av = avma;
     507             :   GEN T, X;
     508             :   long rt_n, c;
     509             : 
     510      646322 :   a = qfbred_i(a);
     511      646322 :   g = qfbred_i(g);
     512             : 
     513      646322 :   rt_n = sqrt((double)n);
     514      646322 :   c = n / rt_n;
     515      646322 :   c = (c * rt_n < n + 1) ? c + 1 : c;
     516             : 
     517      646322 :   T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
     518      646322 :   X = gen_Shanks(T, a, c, NULL, &qfi_group);
     519      646322 :   return X? gerepileuptoint(av, X): gc_NULL(av);
     520             : }
     521             : 
     522             : GEN
     523         455 : qfbclassno0(GEN x,long flag)
     524             : {
     525         455 :   switch(flag)
     526             :   {
     527         441 :     case 0: return map_proto_G(classno,x);
     528          14 :     case 1: return map_proto_G(classno2,x);
     529           0 :     default: pari_err_FLAG("qfbclassno");
     530             :   }
     531             :   return NULL; /* LCOV_EXCL_LINE */
     532             : }
     533             : 
     534             : /* f^h = 1, return order(f). Set *pfao to its factorization */
     535             : static GEN
     536        5159 : find_order(void *E, GEN f, GEN h, GEN *pfao)
     537             : {
     538        5159 :   GEN v = gen_factored_order(f, h, E, &qfi_group);
     539        5159 :   *pfao = gel(v,2); return gel(v,1);
     540             : }
     541             : 
     542             : static int
     543         301 : ok_q(GEN q, GEN h, GEN d2, long r2)
     544             : {
     545         301 :   if (d2)
     546             :   {
     547         245 :     if (r2 <= 2 && !mpodd(q)) return 0;
     548         245 :     return is_pm1(Z_ppo(q,d2));
     549             :   }
     550             :   else
     551             :   {
     552          56 :     if (r2 <= 1 && !mpodd(q)) return 0;
     553          56 :     return is_pm1(Z_ppo(q,h));
     554             :   }
     555             : }
     556             : 
     557             : /* a,b given by their factorizations. Return factorization of lcm(a,b).
     558             :  * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
     559             : static GEN
     560         182 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
     561             : {
     562         182 :   GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
     563         182 :   GEN A = gen_1, B = gen_1;
     564         182 :   long i, l = lg(P);
     565         182 :   GEN E = cgetg(l, t_COL);
     566         707 :   for (i=1; i<l; i++)
     567             :   {
     568         525 :     GEN p = gel(P,i);
     569         525 :     long va = Z_pval(a,p);
     570         525 :     long vb = Z_pval(b,p);
     571         525 :     if (va < vb)
     572             :     {
     573         203 :       B = mulii(B,powiu(p,vb));
     574         203 :       gel(E,i) = utoi(vb);
     575             :     }
     576             :     else
     577             :     {
     578         322 :       A = mulii(A,powiu(p,va));
     579         322 :       gel(E,i) = utoi(va);
     580             :     }
     581             :   }
     582         182 :   *pA = A;
     583         182 :   *pB = B; return mkmat2(P,E);
     584             : }
     585             : 
     586             : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
     587             : static void
     588         182 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
     589             : {
     590         182 :   GEN A,B, g1 = *pg1, d1 = *pd1;
     591         182 :   *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
     592         182 :   *pg1 = gmul(qfbpow_i(g1, diviiexact(d1,A)),  qfbpow_i(f, diviiexact(o,B)));
     593         182 :   *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
     594         182 : }
     595             : 
     596             : /* *pD = coredisc(x), *pR = regulator (x > 0) or NULL */
     597             : GEN
     598         455 : quadclassnoF(GEN x, GEN *pD)
     599             : {
     600             :   GEN D, H, P, E;
     601         455 :   if (lgefint(x) == 3)
     602             :   {
     603             :     ulong d, h;
     604         407 :     if (signe(x) < 0)
     605             :     {
     606         380 :       h = unegquadclassnoF(x[2], &d);
     607         380 :       if (pD) *pD = utoineg(d);
     608             :     }
     609             :     else
     610             :     {
     611          27 :       h = uposquadclassnoF(x[2], &d);
     612          27 :       if (pD) *pD = utoipos(d);
     613             :     }
     614         407 :     return utoi(h);
     615             :   }
     616          48 :   D = coredisc2_fact(absZ_factor(x), signe(x), &P, &E);
     617          48 :   H = quadclassnoF_fact(D, P, E);
     618             :   /* divide by [ O_K^* : O^* ] */
     619          48 :   if (signe(x) < 0) switch(itou_or_0(D))
     620             :   { /* |x| >= 2^BIL, hence x != D */
     621           0 :     case 4: H = shifti(H,-1); break;
     622           0 :     case 3: H = divis(H,3); break;
     623             :   }
     624          48 :   else if (!equalii(x,D))
     625           8 :     H = diviiexact(H, quadunitindex(D, mkmat2(P, zc_to_ZC(E))));
     626          48 :   if (pD) *pD = D; return H;
     627             : }
     628             : 
     629             : /* f \prod_{p|f}  [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ];
     630             :  * s = 1 or -1; D = s * d; assume Df^2 fits in an ulong */
     631             : ulong
     632        4358 : uquadclassnoF_fact(ulong d, long s, GEN P, GEN E)
     633             : {
     634        4358 :   long i, l = lg(P);
     635        4358 :   ulong H = 1;
     636        4567 :   for (i = 1; i < l; i++)
     637             :   {
     638         209 :     ulong p = P[i], e = E[i];
     639         209 :     long D = (long)(p == 2? d & 7: d % p), a;
     640         209 :     if (s < 0) D = -D;
     641         209 :     a = kross(D,p);
     642         209 :     if (!a)
     643          56 :       H *= upowuu(p, e);
     644             :     else
     645             :     {
     646         153 :       H *= p - a;
     647         153 :       if (e >= 2) H *= upowuu(p, e-1);
     648             :     }
     649             :   }
     650        4358 :   return H;
     651             : }
     652             : /* f \prod_{p|f}  [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
     653             : GEN
     654         125 : quadclassnoF_fact(GEN D, GEN P, GEN E)
     655             : {
     656         125 :   long i, l = lg(P);
     657             :   GEN H;
     658         125 :   if (typ(E) != t_VECSMALL) E = vec_to_vecsmall(E);
     659         245 :   for (i = 1, H = gen_1; i < l; i++)
     660             :   {
     661         120 :     GEN p = gel(P,i);
     662         120 :     long e = E[i], s = kronecker(D,p);
     663         120 :     if (!s)
     664          42 :       H = mulii(H, e == 1? p: powiu(p, e));
     665             :     else
     666             :     {
     667          78 :       H = mulii(H, subis(p, s));
     668          78 :       if (e >= 2) H = mulii(H, e == 2? p: powiu(p,e-1));
     669             :     }
     670             :   }
     671         125 :   return H;
     672             : }
     673             : ulong
     674         380 : unegquadclassnoF(ulong x, ulong *pD)
     675             : {
     676         380 :   pari_sp av = avma;
     677             :   GEN E, P;
     678         380 :   ulong D = coredisc2u_fact(factoru(x), -1, &P, &E);
     679         380 :   long H = uquadclassnoF_fact(D, -1, P, E);
     680         380 :   switch(D)
     681             :   { /* divide by [ O_K^* : O^* ] */
     682           0 :     case 4: H >>= 1; break;
     683           0 :     case 3: H /= 3; break;
     684             :   }
     685         380 :   *pD = D; return gc_ulong(av, H);
     686             : }
     687             : ulong
     688          27 : uposquadclassnoF(ulong x, ulong *pD)
     689             : {
     690             :   GEN P, E;
     691          27 :   ulong D = coredisc2u_fact(factoru(x), 1, &P, &E);
     692          27 :   ulong H = uquadclassnoF_fact(D, 1, P, E);
     693          27 :   if (x != D)
     694             :   {
     695          13 :     GEN F = mkvec2(utoipos(usqrt(x / D)), mkmat2(zc_to_ZC(P), zc_to_ZC(E)));
     696          13 :     H /= itou(quadunitindex(utoipos(D), F));
     697             :   }
     698          27 :   *pD = D; return H;
     699             : }
     700             : 
     701             : static long
     702         413 : two_rank(GEN x)
     703             : {
     704         413 :   GEN p = gel(absZ_factor(x),1);
     705         413 :   long l = lg(p)-1;
     706             : #if 0 /* positive disc not needed */
     707             :   if (signe(x) > 0)
     708             :   {
     709             :     long i;
     710             :     for (i=1; i<=l; i++)
     711             :       if (mod4(gel(p,i)) == 3) { l--; break; }
     712             :   }
     713             : #endif
     714         413 :   return l-1;
     715             : }
     716             : 
     717             : static GEN
     718        7847 : sqr_primeform(GEN x, ulong p) { return qfbsqr_i(primeform_u(x, p)); }
     719             : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
     720             : static GEN
     721         413 : get_forms(GEN D, GEN *pL)
     722             : {
     723         413 :   const long MAXFORM = 20;
     724         413 :   GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
     725         413 :   GEN forms = vectrunc_init(MAXFORM+1);
     726         413 :   long s, nforms = 0;
     727             :   ulong p;
     728             :   forprime_t S;
     729         413 :   L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
     730         413 :   s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
     731         413 :   if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
     732         413 :   if      (s < 10)   s = 200;
     733         413 :   else if (s < 20)   s = 1000;
     734         413 :   else if (s < 5000) s = 5000;
     735         413 :   u_forprime_init(&S, 2, s);
     736    13458732 :   while ( (p = u_forprime_next(&S)) )
     737             :   {
     738    13458319 :     long d, k = kroiu(D,p);
     739             :     pari_sp av2;
     740    13458319 :     if (!k) continue;
     741    13457605 :     if (k > 0)
     742             :     {
     743     6730451 :       if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
     744     6730451 :       d = p-1;
     745             :     }
     746             :     else
     747     6727154 :       d = p+1;
     748    13457605 :     av2 = avma; affrr(divru(mulur(p,L),d), L); set_avma(av2);
     749             :   }
     750         413 :   *pL = L; return forms;
     751             : }
     752             : 
     753             : /* h ~ #G, return o = order of f, set fao = its factorization */
     754             : static  GEN
     755         511 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
     756             : {
     757         511 :   long s = minss(itos(sqrti(h)), 10000);
     758         511 :   GEN T = gen_Shanks_init(f, s, E, &qfi_group);
     759         511 :   GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
     760         511 :   return find_order(E, f, addiu(v,1), pfao);
     761             : }
     762             : 
     763             : /* if g = 1 in  G/<f> ? */
     764             : static int
     765        5684 : equal1(void *E, GEN T, ulong N, GEN g)
     766        5684 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
     767             : 
     768             : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
     769             :  * FIXME: should be gen_order, but equal1 has the wrong prototype */
     770             : static GEN
     771         350 : relative_order(void *E, GEN a, GEN o, ulong N,  GEN T)
     772             : {
     773         350 :   pari_sp av = avma;
     774             :   long i, l;
     775             :   GEN m;
     776             : 
     777         350 :   m = get_arith_ZZM(o);
     778         350 :   if (!m) pari_err_TYPE("gen_order [missing order]",a);
     779         350 :   o = gel(m,1);
     780         350 :   m = gel(m,2); l = lgcols(m);
     781        1148 :   for (i = l-1; i; i--)
     782             :   {
     783         798 :     GEN t, y, p = gcoeff(m,i,1);
     784         798 :     long j, e = itos(gcoeff(m,i,2));
     785         798 :     if (l == 2) {
     786          49 :       t = gen_1;
     787          49 :       y = a;
     788             :     } else {
     789         749 :       t = diviiexact(o, powiu(p,e));
     790         749 :       y = powgi(a, t);
     791             :     }
     792         798 :     if (equal1(E, T,N,y)) o = t;
     793             :     else {
     794         364 :       for (j = 1; j < e; j++)
     795             :       {
     796          84 :         y = powgi(y, p);
     797          84 :         if (equal1(E, T,N,y)) break;
     798             :       }
     799         357 :       if (j < e) {
     800          77 :         if (j > 1) p = powiu(p, j);
     801          77 :         o = mulii(t, p);
     802             :       }
     803             :     }
     804             :   }
     805         350 :   return gerepilecopy(av, o);
     806             : }
     807             : 
     808             : /* h(x) for x<0 using Baby Step/Giant Step.
     809             :  * Assumes G is not too far from being cyclic.
     810             :  *
     811             :  * Compute G^2 instead of G so as to kill most of the noncyclicity */
     812             : GEN
     813         441 : classno(GEN x)
     814             : {
     815         441 :   pari_sp av = avma;
     816             :   long r2, k, s, i, l;
     817             :   GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
     818             :   void *E;
     819             : 
     820         441 :   if (signe(x) >= 0) return classno2(x);
     821             : 
     822         413 :   check_quaddisc(x, &s, &k, "classno");
     823         413 :   if (abscmpiu(x,12) <= 0) return gen_1;
     824             : 
     825         413 :   Hf = quadclassnoF(x, &D);
     826         413 :   if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
     827         413 :   forms =  get_forms(D, &L);
     828         413 :   r2 = two_rank(D);
     829         413 :   hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
     830             : 
     831         413 :   l = lg(forms);
     832         413 :   order_bound = const_vec(l-1, NULL);
     833         413 :   E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
     834         413 :   g1 = gel(forms,1);
     835         413 :   gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
     836         413 :   q = diviiround(hin, d1); /* approximate order of G/<g1> */
     837         413 :   d2 = NULL; /* not computed yet */
     838         413 :   if (is_pm1(q)) goto END;
     839        7315 :   for (i=2; i < l; i++)
     840             :   {
     841        6951 :     GEN o, fao, a, F, fd, f = gel(forms,i);
     842        6951 :     fd = qfbpow_i(f, d1); if (is_pm1(gel(fd,1))) continue;
     843         182 :     F = qfbpow_i(fd, q);
     844         182 :     a = gel(F,1);
     845         182 :     o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
     846             :     /* f^(d1 q) = 1 */
     847         182 :     fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
     848         182 :     o = find_order(E, f, fao, &fao);
     849         182 :     gel(order_bound,i) = o;
     850             :     /* o = order of f, fao = factor(o) */
     851         182 :     update_g1(&g1,&d1,&fad1, f,o,fao);
     852         182 :     q = diviiround(hin, d1);
     853         182 :     if (is_pm1(q)) goto END;
     854             :   }
     855             :   /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
     856         364 :   if (expi(q) > 3)
     857             :   { /* q large: compute d2, 2nd elt divisor */
     858         308 :     ulong N, n = 2*itou(sqrti(d1));
     859         308 :     GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
     860         308 :     d2 = gen_1;
     861         308 :     N = itou( gceil(gdivgu(d1,n)) ); /* order(g1) <= n*N */
     862        5047 :     for (i = 1; i < l; i++)
     863             :     {
     864        4802 :       GEN d, f = gel(forms,i), B = gel(order_bound,i);
     865        4802 :       if (!B) B = find_order(E, f, fad1, /*junk*/&d);
     866        4802 :       f = qfbpow_i(f,d2);
     867        4802 :       if (equal1(E, T,N,f)) continue;
     868         350 :       B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
     869             :       /* f^B = 1 */
     870         350 :       d = relative_order(E, f, B, N,T);
     871         350 :       d2= mulii(d,d2);
     872         350 :       D = mulii(d1,d2);
     873         350 :       q = diviiround(hin,D);
     874         350 :       if (is_pm1(q)) { d1 = D; goto END; }
     875             :     }
     876             :     /* very probably, d2 is the 2nd elementary divisor */
     877         245 :     d1 = D; /* product of first two elt divisors */
     878             :   }
     879             :   /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
     880             :    * 2-rank */
     881         301 :   if (!ok_q(q,d1,d2,r2))
     882             :   {
     883           0 :     GEN q0 = q;
     884             :     long d;
     885           0 :     if (cmpii(mulii(q,d1), hin) < 0)
     886             :     { /* try q = q0+1,-1,+2,-2 */
     887           0 :       d = 1;
     888           0 :       do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
     889             :     }
     890             :     else
     891             :     { /* q0-1,+1,-2,+2  */
     892           0 :       d = -1;
     893           0 :       do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
     894             :     }
     895             :   }
     896         301 :   d1 = mulii(d1,q);
     897             : 
     898         413 : END:
     899         413 :   return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
     900             : }
     901             : 
     902             : /* use Euler products */
     903             : GEN
     904          42 : classno2(GEN x)
     905             : {
     906          42 :   pari_sp av = avma;
     907          42 :   const long prec = DEFAULTPREC;
     908             :   long n, i, s;
     909          42 :   GEN p1, p2, S, p4, p5, p7, Hf, Pi, logd, sqrtd, D, half, reg = NULL;
     910             : 
     911          42 :   check_quaddisc(x, &s, NULL, "classno2");
     912          42 :   if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
     913             : 
     914          42 :   Hf = quadclassnoF(x, &D);
     915          42 :   if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
     916             : 
     917          42 :   Pi = mppi(prec);
     918          42 :   sqrtd = sqrtr_abs(itor(D, prec));
     919          42 :   logd = logr_abs(sqrtd); shiftr_inplace(logd,1);
     920          42 :   p1 = sqrtr_abs(divrr(mulir(D,logd), gmul2n(Pi,1)));
     921          42 :   if (s > 0)
     922             :   {
     923          35 :     GEN invlogd = invr(logd);
     924          35 :     reg = quadregulator(D, prec);
     925          35 :     p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
     926          35 :     if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
     927             :   }
     928          42 :   n = itos_or_0( mptrunc(p1) );
     929          42 :   if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
     930             : 
     931          42 :   p4 = divri(Pi, D); setsigne(p4, 1);
     932          42 :   p7 = invr(sqrtr_abs(Pi));
     933          42 :   half = real2n(-1, prec);
     934          42 :   if (s > 0)
     935             :   { /* i = 1, shortcut */
     936          35 :     p1 = sqrtd;
     937          35 :     p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
     938          35 :     S = addrr(mulrr(p1,p5), eint1(p4,prec));
     939         588 :     for (i=2; i<=n; i++)
     940             :     {
     941         553 :       long k = kroiu(D,i); if (!k) continue;
     942         455 :       p2 = mulir(sqru(i), p4);
     943         455 :       p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
     944         455 :       p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
     945         455 :       S = (k>0)? addrr(S,p5): subrr(S,p5);
     946             :     }
     947          35 :     S = shiftr(divrr(S,reg),-1);
     948             :   }
     949             :   else
     950             :   { /* i = 1, shortcut */
     951           7 :     p1 = gdiv(sqrtd, Pi);
     952           7 :     p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
     953           7 :     S = addrr(p5, divrr(p1, mpexp(p4)));
     954         952 :     for (i=2; i<=n; i++)
     955             :     {
     956         945 :       long k = kroiu(D,i); if (!k) continue;
     957         945 :       p2 = mulir(sqru(i), p4);
     958         945 :       p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
     959         945 :       p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
     960         945 :       S = (k>0)? addrr(S,p5): subrr(S,p5);
     961             :     }
     962             :   }
     963          42 :   return gerepileuptoint(av, mulii(Hf, roundr(S)));
     964             : }
     965             : 
     966             : /* 1 + q + ... + q^v, v > 0 */
     967             : static GEN
     968         973 : geomsumu(ulong q, long v)
     969             : {
     970         973 :   GEN u = utoipos(1+q);
     971        1267 :   for (; v > 1; v--) u = addui(1, mului(q, u));
     972         973 :   return u;
     973             : }
     974             : static GEN
     975         973 : geomsum(GEN q, long v)
     976             : {
     977             :   GEN u;
     978         973 :   if (lgefint(q) == 3) return geomsumu(q[2], v);
     979           0 :   u = addiu(q,1);
     980           0 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
     981           0 :   return u;
     982             : }
     983             : 
     984             : /* 1+p+...+p^(e-1), e >= 1; assuming result fits in an ulong */
     985             : static ulong
     986        2433 : usumpow(ulong p, long e)
     987             : {
     988             :   ulong q;
     989             :   long i;
     990        2433 :   if (p == 2) return (1UL << e) - 1; /* also OK if e = BITS_IN_LONG */
     991        1523 :   e--; for (i = 1, q = p + 1; i < e; i++) q = p*q + 1;
     992        1488 :   return q;
     993             : }
     994             : long
     995       44025 : uhclassnoF_fact(GEN faF, long D)
     996             : {
     997       44025 :   GEN P = gel(faF,1), E = gel(faF,2);
     998       44025 :   long i, t, l = lg(P);
     999      101954 :   for (i = t = 1; i < l; i++)
    1000             :   {
    1001       57929 :     long p = P[i], e = E[i], s = kross(D,p);
    1002       57929 :     if (e == 1) { t *= 1 + p - s; continue; }
    1003       15137 :     if (s == 1) { t *= upowuu(p,e); continue; }
    1004        2433 :     t *= 1 + usumpow(p, e) * (p - s);
    1005             :   }
    1006       44025 :   return t;
    1007             : }
    1008             : /* Hurwitz(D F^2)/ Hurwitz(D)
    1009             :  * = \sum_{f|F}  f \prod_{p|f} (1-kro(D/p)/p)
    1010             :  * = \prod_{p^e || F} (1 + (p^e-1) / (p-1) * (p-kro(D/p))) */
    1011             : GEN
    1012       60723 : hclassnoF_fact(GEN P, GEN E, GEN D)
    1013             : {
    1014             :   GEN H;
    1015       60723 :   long i, l = lg(P);
    1016       60723 :   if (l == 1) return gen_1;
    1017       56169 :   for (i = 1, H = NULL; i < l; i++)
    1018             :   {
    1019       30242 :     GEN t, p = gel(P,i);
    1020       30242 :     long e = E[i], s = kronecker(D,p);
    1021       30243 :     if (e == 1) t = addiu(p, 1-s);
    1022        1555 :     else if (s == 1) t = powiu(p, e);
    1023         973 :     else t = addui(1, mulii(subis(p, s), geomsum(p, e-1)));
    1024       30244 :     H = H? mulii(H,t): t;
    1025             :   }
    1026       25927 :   return H;
    1027             : }
    1028             : static GEN
    1029       60726 : hclassno6_large(GEN x)
    1030             : {
    1031       60726 :   GEN H = NULL, P, E, D = coredisc2_fact(absZ_factor(x), -1, &P, &E);
    1032       60728 :   long l = lg(P);
    1033             : 
    1034       60728 :   if (l > 1 && lgefint(x) == 3)
    1035             :   { /* F != 1, second chance */
    1036       25928 :     ulong h = hclassno6u_no_cache(x[2]);
    1037       25928 :     if (h) H = utoipos(h);
    1038             :   }
    1039       60728 :   if (!H)
    1040             :   {
    1041       60728 :     H = quadclassno(D);
    1042       60727 :     switch(itou_or_0(D))
    1043             :     {
    1044          49 :       case 3: H = shifti(H,1);break;
    1045           7 :       case 4: H = muliu(H,3); break;
    1046       60671 :       default:H = muliu(H,6); break;
    1047             :     }
    1048           0 :   }
    1049       60723 :   return mulii(H, hclassnoF_fact(P, E, D));
    1050             : }
    1051             : 
    1052             : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
    1053             : GEN
    1054      132477 : hclassno6(GEN x)
    1055             : {
    1056      132477 :   ulong d = itou_or_0(x);
    1057      132478 :   if (d)
    1058             :   { /* create cache if d small */
    1059      132477 :     ulong h = d < 500000 ? hclassno6u(d): hclassno6u_no_cache(d);
    1060      132479 :     if (h) return utoipos(h);
    1061             :   }
    1062       60726 :   return hclassno6_large(x);
    1063             : }
    1064             : 
    1065             : GEN
    1066          49 : hclassno(GEN x)
    1067             : {
    1068             :   long a, s;
    1069          49 :   if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
    1070          49 :   s = signe(x);
    1071          49 :   if (s < 0) return gen_0;
    1072          49 :   if (!s) return gdivgs(gen_1, -12);
    1073          49 :   a = mod4(x); if (a == 1 || a == 2) return gen_0;
    1074          49 :   return gdivgu(hclassno6(x), 6);
    1075             : }
    1076             : 
    1077             : /* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
    1078             : GEN
    1079     2774818 : Zn_quad_roots(GEN N, GEN B, GEN C)
    1080             : {
    1081     2774818 :   pari_sp av = avma;
    1082     2774818 :   GEN fa = NULL, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
    1083             :   long l, i, j, ct;
    1084             : 
    1085     2774818 :   if ((fa = check_arith_non0(N,"Zn_quad_roots")))
    1086             :   {
    1087       35994 :     N = typ(N) == t_VEC? gel(N,1): factorback(N);
    1088       35994 :     fa = clean_Z_factor(fa);
    1089             :   }
    1090     2774827 :   N = absi_shallow(N);
    1091     2774826 :   N4 = shifti(N,2);
    1092     2774794 :   D = modii(subii(sqri(B), shifti(C,2)), N4);
    1093     2774776 :   if (!signe(D))
    1094             :   { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
    1095         812 :     if (!fa) fa = Z_factor(N);
    1096         812 :     P = gel(fa,1);
    1097         812 :     E = ZV_to_zv(gel(fa,2));
    1098         812 :     l = lg(P);
    1099        1757 :     for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
    1100         812 :     Np = factorback2(P, E); /* x = -B mod N' */
    1101         812 :     B = shifti(B,-1);
    1102         812 :     return gerepilecopy(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
    1103             :   }
    1104     2773964 :   if (!fa)
    1105     2738204 :     fa = Z_factor(N4);
    1106             :   else  /* convert to factorization of N4 = 4*N */
    1107       35760 :     fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
    1108     2774005 :   P = gel(fa,1); l = lg(P);
    1109     2774005 :   E = ZV_to_zv(gel(fa,2));
    1110     2774001 :   F = cgetg(l, t_VEC);
    1111     2773998 :   mF= cgetg(l, t_VEC); F0 = gen_0;
    1112     2773998 :   Q = cgetg(l, t_VEC); Q0 = gen_1;
    1113     6677279 :   for (i = j = 1, ct = 0; i < l; i++)
    1114             :   {
    1115     6039191 :     GEN p = gel(P,i), q, f, mf, D0;
    1116     6039191 :     long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
    1117     6039207 :     if (d <= 0)
    1118             :     {
    1119     1377859 :       q = powiu(p, (s+1)>>1);
    1120     2385187 :       Q0 = mulii(Q0, q); continue;
    1121             :     }
    1122             :     /* d > 0 */
    1123     6611737 :     if (odd(t)) return NULL;
    1124     4475903 :     t2 = t >> 1;
    1125     4475903 :     if (i > 1)
    1126             :     { /* p > 2 */
    1127     2823211 :       if (kronecker(D0, p) == -1) return NULL;
    1128     1375067 :       q = powiu(p, s - t2);
    1129     1375069 :       f = Zp_sqrt(D0, p, d);
    1130     1375079 :       if (!f) return NULL; /* p was not actually prime... */
    1131     1375065 :       if (t2) f = mulii(powiu(p,t2), f);
    1132     1375065 :       mf = Fp_neg(f, q);
    1133             :     }
    1134             :     else
    1135             :     { /* p = 2 */
    1136     1652692 :       if (d <= 3)
    1137             :       {
    1138     1282770 :         if (d == 3 && Mod8(D0) != 1) return NULL;
    1139     1058231 :         if (d == 2 && Mod4(D0) != 1) return NULL;
    1140     1007327 :         Q0 = int2n(1+t2); F0 = NULL; continue;
    1141             :       }
    1142      369922 :       if (Mod8(D0) != 1) return NULL;
    1143      143143 :       q = int2n(d - 1 + t2);
    1144      143143 :       f = Z2_sqrt(D0, d);
    1145      143143 :       if (t2) f = shifti(f, t2);
    1146      143143 :       mf = Fp_neg(f, q);
    1147             :     }
    1148     1518167 :     gel(Q,j) = q;
    1149     1518167 :     gel(F,j) = f;
    1150     1518167 :     gel(mF,j)= mf; j++;
    1151             :   }
    1152      638088 :   setlg(Q,j);
    1153      638097 :   setlg(F,j);
    1154      638095 :   setlg(mF,j);
    1155      638096 :   if (is_pm1(Q0)) A = leafcopy(F);
    1156             :   else
    1157             :   { /* append the fixed congruence (F0 mod Q0) */
    1158      603181 :     if (!F0) F0 = shifti(Q0,-1);
    1159      603181 :     A = shallowconcat(F, F0);
    1160      603200 :     Q = shallowconcat(Q, Q0);
    1161             :   }
    1162      638124 :   ct = 1 << (j-1);
    1163      638124 :   T = ZV_producttree(Q);
    1164      638083 :   R = ZV_chinesetree(Q,T);
    1165      638093 :   Np = gmael(T, lg(T)-1, 1);
    1166      638093 :   B = modii(B, Np);
    1167      638104 :   if (!signe(B)) B = NULL;
    1168      638104 :   Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
    1169      638091 :   w = cgetg(3, t_VEC);
    1170      638110 :   gel(w,1) = icopy(Np);
    1171      638131 :   gel(w,2) = v = cgetg(ct+1, t_VEC);
    1172      638133 :   l = lg(F);
    1173     2737040 :   for (j = 1; j <= ct; j++)
    1174             :   {
    1175     2098926 :     pari_sp av2 = avma;
    1176     2098926 :     long m = j - 1;
    1177             :     GEN u;
    1178     6101649 :     for (i = 1; i < l; i++)
    1179             :     {
    1180     4002723 :       gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
    1181     4002723 :       m >>= 1;
    1182             :     }
    1183     2098926 :     u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
    1184     2098886 :     if (B) u = subii(u,B);
    1185     2098886 :     gel(v,j) = gerepileuptoint(av2, modii(shifti(u,-1), Np));
    1186             :   }
    1187      638114 :   return gerepileupto(av, w);
    1188             : }

Generated by: LCOV version 1.13