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.18.0 lcov report (development 29580-d8b649a849) Lines: 687 715 96.1 %
Date: 2024-10-04 09:07:15 Functions: 57 60 95.0 %
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       22862 : isfundamental(GEN x)
      57             : {
      58       22862 :   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       21455 :   return Z_isfundamental(x);
      65             : }
      66             : 
      67             : /* x fundamental ? */
      68             : long
      69       16953 : uposisfundamental(ulong x)
      70             : {
      71       16953 :   ulong r = x & 15; /* x mod 16 */
      72       16953 :   if (!r) return 0;
      73       16141 :   switch(r & 3)
      74             :   { /* x mod 4 */
      75        3487 :     case 0: return (r == 4)? 0: uissquarefree(x >> 2);
      76        6203 :     case 1: return uissquarefree(x);
      77        6451 :     default: return 0;
      78             :   }
      79             : }
      80             : /* -x fundamental ? */
      81             : long
      82       31472 : unegisfundamental(ulong x)
      83             : {
      84       31472 :   ulong r = x & 15; /* x mod 16 */
      85       31472 :   if (!r) return 0;
      86       29827 :   switch(r & 3)
      87             :   { /* x mod 4 */
      88        7827 :     case 0: return (r == 12)? 0: uissquarefree(x >> 2);
      89       11671 :     case 3: return uissquarefree(x);
      90       10329 :     default: return 0;
      91             :   }
      92             : }
      93             : long
      94       25109 : sisfundamental(long x)
      95       25109 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
      96             : 
      97             : long
      98       22022 : Z_isfundamental(GEN x)
      99             : {
     100             :   long r;
     101       22022 :   switch(lgefint(x))
     102             :   {
     103           7 :     case 2: return 0;
     104       11368 :     case 3: return signe(x) < 0? unegisfundamental(x[2])
     105       31373 :                                : 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        2079 : fa_quaddisc(GEN f)
     126             : {
     127        2079 :   GEN P = gel(f,1), E = gel(f,2), s = gen_1;
     128        2079 :   long i, l = lg(P);
     129        6881 :   for (i = 1; i < l; i++) /* possibly including -1 */
     130        4802 :     if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
     131        2079 :   if (Mod4(s) > 1) s = shifti(s,2);
     132        2079 :   return s;
     133             : }
     134             : 
     135             : GEN
     136        2933 : quaddisc(GEN x)
     137             : {
     138        2933 :   const pari_sp av = avma;
     139        2933 :   long tx = typ(x);
     140             :   GEN F;
     141        2933 :   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        2079 :   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    21199714 : update_f(GEN f, GEN u)
     164             : {
     165    21199714 :   GEN a = gcoeff(f,1,1), b = gcoeff(f,1,2);
     166    21199714 :   GEN c = gcoeff(f,2,1), d = gcoeff(f,2,2);
     167    21199714 :   gcoeff(f,1,1) = addmulii(b, u,a); gcoeff(f,1,2) = a;
     168    21199714 :   gcoeff(f,2,1) = addmulii(d, u,c); gcoeff(f,2,2) = c;
     169    21199714 : }
     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] may be an int: implicitely represent [f[1],1;1,0] */
     174             : static long
     175    21360738 : update_fm(GEN f, GEN a, long i)
     176             : {
     177             : #ifdef LONG_IS_64BIT
     178    18309204 :   const long LIM = 10;
     179             : #else
     180     3051534 :   const long LIM = 18;
     181             : #endif
     182    21360738 :   pari_sp av = avma;
     183             :   long k, v;
     184             :   GEN u;
     185    21360738 :   if (!odd(i)) { gel(f,1) = a; return i+1; }
     186    21280226 :   u = gel(f, 1);
     187    21280226 :   if (typ(u) == t_INT) /* [u,1;1,0] * [a,1;1,0] */
     188       80512 :   { gel(f,1) = mkmat22(addiu(mulii(a, u), 1), u, a, gen_1); return i; }
     189    21199714 :   update_f(u, a); if (lgefint(gcoeff(u,1,1)) < LIM) return i;
     190       80512 :   v = vals(i+1); gel(f,1) = gen_0;
     191      160970 :   for (k = 1; k < v; k++) { u = ZM2_mul(gel(f,k+1), u); gel(f,k+1) = gen_0; }
     192       80512 :   if (v != 1) u = gerepileupto(av, u);
     193       80512 :   gel(f,v+1) = u; return i+1;
     194             : }
     195             : /* \prod f[j]; if first only return column 1 */
     196             : static GEN
     197           7 : prod_fm(GEN f, long i, long first)
     198             : {
     199           7 :   long k, v = vals(i) + 1;
     200           7 :   GEN u = gel(f, v);
     201             :   /* i a power of 2: f[1] can't be a t_INT */
     202           7 :   if (!(i >>= v)) return first? gel(u,1): u;
     203         105 :   for (k = v+1; i; i >>= 1, k++)
     204          98 :     if (odd(i))
     205             :     {
     206          54 :       GEN w = gel(f,k);
     207          54 :       switch(typ(u))
     208             :       {
     209           0 :         case t_INT: update_f(w, u);
     210           0 :           u = first? gel(w,1): w; break;
     211           0 :         case t_COL: /* implies 'first' */
     212           0 :           u = ZM_ZC_mul(w, u); break;
     213          54 :         default: /* t_MAT */
     214          54 :           u = first? ZM_ZC_mul(w, gel(u,1)): ZM2_mul(w, u); break;
     215             :       }
     216             :     }
     217           7 :   return u;
     218             : }
     219             : 
     220             : GEN
     221       69048 : quadunit0(GEN x, long v)
     222             : {
     223       69048 :   GEN y = quadunit(x);
     224       69041 :   if (v==-1) v = fetch_user_var("w");
     225       69041 :   setvarn(gel(y,1), v); return y;
     226             : }
     227             : 
     228             : struct uimod { GEN N, T; };
     229             : static GEN
     230       21378 : ui_pow(void *E, GEN x, GEN n)
     231       21378 : { struct uimod *S = (struct uimod*)E; return FpXQ_pow(x, n, S->T, S->N); }
     232             : static int
     233       38612 : ui_equal1(GEN x) { return degpol(x) < 1; }
     234             : static const struct bb_group
     235             : ui_group={ NULL,ui_pow,NULL,NULL,NULL,ui_equal1,NULL};
     236             : 
     237             : static void
     238          98 : quadunit_uvmod(GEN D, GEN d, GEN N, GEN *pu, GEN *pv)
     239             : {
     240             :   GEN u1, u2, v1, v2, p, q, q1, u, v;
     241          98 :   int m = mpodd(D), first = 1;
     242          98 :   pari_sp av = avma;
     243          98 :   p = (mpodd(d) == m)? d: subiu(d, 1);
     244          98 :   u1 = negi(p); u2 = gen_2;
     245          98 :   v1 = gen_1; v2 = gen_0; q = gen_2;
     246          98 :   q1 = shifti(subii(D, sqri(p)), -1);
     247             :   for(;;)
     248         308 :   {
     249         406 :     GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
     250         406 :     p = subii(d, r);
     251         406 :     if (equalii(p1, p) && !first)
     252             :     { /* even period */
     253          14 :       u = addmulii(sqri(u2), D, sqri(v2));
     254          14 :       v = shifti(mulii(u2,v2), 1);
     255          14 :       break;
     256             :     }
     257         392 :     first = 0;
     258         392 :     t = Fp_addmul(u1, A, u2, N); u1 = u2; u2 = t;
     259         392 :     t = Fp_addmul(v1, A, v2, N); v1 = v2; v2 = t;
     260         392 :     t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
     261         392 :     if (equalii(q, t))
     262             :     { /* odd period */
     263          84 :       u = addmulii(mulii(u1,u2), D, mulii(v1,v2));
     264          84 :       v = addmulii(mulii(u1,v2), u2, v1);
     265          84 :       break;
     266             :     }
     267         308 :     if (gc_needed(av, 2))
     268             :     {
     269           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadunit_uvmod");
     270           0 :       gerepileall(av, 7, &p, &u1,&u2,&v1,&v2, &q,&q1);
     271             :     }
     272             :   }
     273          98 :   *pu = modii(u, N);
     274          98 :   *pv = modii(v, N); if (m) *pu = Fp_sub(*pu, *pv, N);
     275          98 : }
     276             : /* fundamental unit is u + vx mod quadpoly(D); always called with D
     277             :  * fundamental and relatively small but would work in all cases. Should be
     278             :  * called whenever the fundamental unit is so "small" that asymptotically
     279             :  * fast multiplication is not used in the continued fraction loop */
     280             : static void
     281       69034 : quadunit_uv_basecase(GEN D, GEN *pu, GEN *pv)
     282             : {
     283       69034 :   GEN u1, u2, v1, v2, p, q, q1, u, v, a, b, c, d = sqrtremi(D, &a);
     284       69034 :   int m = mpodd(D);
     285       69034 :   long first = 1;
     286             : 
     287       69034 :   p = d; q1 = shifti(a, -1); q = gen_2;
     288       69034 :   if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
     289       69034 :   u1 = gen_2; u2 = p;
     290       69034 :   v1 = gen_0; v2 = gen_1;
     291             :   for(;;)
     292      871052 :   {
     293      940086 :     GEN t = q;
     294      940086 :     if (first) { first = 0; q = q1; }
     295             :     else
     296             :     {
     297      871052 :       GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
     298      871052 :       p = subii(d, r);
     299      871052 :       if (equalii(p1, p)) /* even period */
     300       45738 :       { a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break; }
     301      825314 :       r = addmulii(u1, A, u2); u1 = u2; u2 = r;
     302      825314 :       r = addmulii(v1, A, v2); v1 = v2; v2 = r;
     303      825314 :       q = submulii(q1, A, subii(p, p1));
     304             :     }
     305      894348 :     q1 = t;
     306      894348 :     if (equalii(q, t))
     307             :     { /* odd period */
     308       23296 :       a = mulii(u1, u2); b = mulii(v1, v2);
     309       23296 :       c = mulii(addii(u1, v1), addii(u2, v2)); break;
     310             :     }
     311             :   }
     312       69034 :   u = diviiexact(addmulii(a, D, b), q);
     313       69034 :   v = diviiexact(subii(c, addii(a, b)), q);
     314       69034 :   if (m == 1) u = subii(u, v);
     315       69034 :   *pu = shifti(u, -1); *pv = v;
     316       69034 : }
     317             : 
     318             : /* D > 0, d = sqrti(D) */
     319             : static GEN
     320       69118 : quadunit_q(GEN D, GEN d, long *pN)
     321             : {
     322       69118 :   pari_sp av = avma;
     323             :   GEN p, q, q1;
     324       69118 :   long first = 1;
     325       69118 :   p = (Mod2(d) == Mod2(D))? d: subiu(d, 1);
     326       69118 :   q = gen_2;
     327       69118 :   q1 = shifti(subii(D, sqri(p)), -1);
     328             :   for(;;)
     329     1001133 :   {
     330     1070251 :     GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
     331     1070251 :     p = subii(d, r);
     332     1087555 :     if (!first && equalii(p1, p)) { *pN = 1; return q; } /* even period */
     333     1018437 :     first = 0;
     334     1018437 :     t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
     335     1018437 :     if (equalii(q, t)) { *pN = -1; return q; } /* odd period */
     336     1001133 :     if (gc_needed(av, 2))
     337             :     {
     338           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadunitnorm");
     339           0 :       gerepileall(av, 3, &p, &q, &q1);
     340             :     }
     341             :   }
     342             : }
     343             : /* fundamental unit mod N */
     344             : static GEN
     345          98 : quadunit_mod(GEN D, GEN N)
     346             : {
     347          98 :   GEN q, u, v, d = sqrti(D);
     348          98 :   pari_sp av = avma;
     349             :   long s;
     350          98 :   q = gerepileuptoint(av, quadunit_q(D, d, &s));
     351          98 :   if (mpodd(N) && equali1(gcdii(q, N)))
     352             :   {
     353          14 :     quadunit_uvmod(D, d, N, &u, &v);
     354          14 :     q = Fp_inv(shifti(q, 1), N);
     355          14 :     u = Fp_mul(u, q, N);
     356          14 :     v = Fp_mul(v, q, N); v = modii(shifti(v, 1), N);
     357             :   }
     358             :   else
     359             :   {
     360          84 :     GEN M = shifti(mulii(q, N), 1);
     361          84 :     quadunit_uvmod(D, d, M, &u, &v);
     362          84 :     u = diviiexact(u, q);
     363          84 :     v = modii(diviiexact(v, q), N); u = shifti(u,-1);
     364             :   }
     365          98 :   return deg1pol_shallow(v, u, 0);
     366             : }
     367             : 
     368             : /* f \prod_{p|f}  [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
     369             : static GEN
     370       26592 : quadclassnoEuler_fact(GEN D, GEN P, GEN E)
     371             : {
     372       26592 :   long i, l = lg(P);
     373             :   GEN H;
     374       26592 :   if (typ(E) != t_VECSMALL) E = vec_to_vecsmall(E);
     375       56525 :   for (i = 1, H = gen_1; i < l; i++)
     376             :   {
     377       29933 :     GEN p = gel(P,i);
     378       29933 :     long e = E[i], s = kronecker(D,p);
     379       29933 :     if (!s)
     380        8526 :       H = mulii(H, e == 1? p: powiu(p, e));
     381             :     else
     382             :     {
     383       21407 :       H = mulii(H, subis(p, s));
     384       21407 :       if (e >= 2) H = mulii(H, e == 2? p: powiu(p,e-1));
     385             :     }
     386             :   }
     387       26592 :   return H;
     388             : }
     389             : 
     390             : /* D > 0; y mod (N,T) congruent to fundamental unit of maximal order and
     391             :  * disc D. Return unit index of order of conductor N */
     392             : static GEN
     393       26544 : quadunitindex_ii(GEN D, GEN N, GEN F, GEN y, GEN T)
     394             : {
     395       26544 :   GEN H = quadclassnoEuler_fact(D, gel(F,1), gel(F,2));
     396       26544 :   GEN P, E, a = Z_smoothen(H, gel(F,1), &P, &E), faH = mkmat2(P, E);
     397             :   struct uimod S;
     398             : 
     399       26544 :   if (a) faH = ZM_merge_factor(Z_factor(a), faH);
     400             :   /* multiple of unit index, in [H, factor(H)] format */
     401       26544 :   S.N = N; S.T = FpX_red(T, N);
     402       26544 :   return gen_order(y, mkvec2(H,faH), (void*)&S, &ui_group);
     403             : }
     404             : static GEN
     405          98 : quadunitindex_i(GEN D, GEN N, GEN F)
     406          98 : { return quadunitindex_ii(D, N, F, quadunit_mod(D, N), quadpoly_i(D)); }
     407             : GEN
     408         112 : quadunitindex(GEN D, GEN N)
     409             : {
     410         112 :   pari_sp av = avma;
     411             :   long r, s;
     412             :   GEN F;
     413         112 :   check_quaddisc(D, &s, &r, "quadunitindex");
     414         105 :   if ((F = check_arith_pos(N,"quadunitindex")))
     415          14 :     N = typ(N) == t_VEC? gel(N,1): factorback(F);
     416          98 :   if (equali1(N)) return gen_1;
     417          91 :   if (s < 0) switch(itos_or_0(D)) {
     418           7 :     case -3: return utoipos(3);
     419           7 :     case -4: return utoipos(2);
     420           7 :     default: return gen_1;
     421             :   }
     422          70 :   return gerepileuptoint(av, quadunitindex_i(D, N, F? F: Z_factor(N)));
     423             : }
     424             : 
     425             : /* fundamental unit is u + vx mod quadpoly(D); always called with D
     426             :  * fundamental but would work in all cases. Same algorithm as basecase,
     427             :  * except we compute the product of elementary matrices with a product tree */
     428             : static void
     429           7 : quadunit_uv(GEN D, GEN *pu, GEN *pv)
     430             : {
     431           7 :   GEN a, b, c, u, v, p, q, q1, f, d = sqrtremi(D, &a);
     432           7 :   pari_sp av = avma;
     433           7 :   long i = 0;
     434           7 :   int m = mpodd(D);
     435             : 
     436           7 :   p = d; q1 = shifti(a, -1); q = gen_2;
     437           7 :   if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
     438           7 :   f = zerovec(2 + (expi(D)>>1));
     439           7 :   gel(f,1) = mkmat22(p, gen_2, gen_1, gen_0);
     440             :   for(;;)
     441    21360738 :   {
     442    21360745 :     GEN t = q, u1,u2, v1,v2;
     443    21360745 :     if (!i) { i = 1; q = q1; }
     444             :     else
     445             :     {
     446    21360738 :       GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
     447    21360738 :       p = subii(d, r);
     448    21360738 :       if (equalii(p1, p))
     449             :       { /* even period */
     450           0 :         f = prod_fm(f, i, 1); u2 = gel(f,1); v2 = gel(f,2);
     451           0 :         a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break;
     452             :       }
     453    21360738 :       i = update_fm(f, A, i);
     454    21360738 :       q = submulii(q1, A, subii(p, p1));
     455             :     }
     456    21360745 :     q1 = t;
     457    21360745 :     if (equalii(q, t))
     458             :     { /* odd period */
     459           7 :       f = prod_fm(f, i, 0);
     460           7 :       u2 = gcoeff(f,1,1); u1 = gcoeff(f,1,2); a = mulii(u1, u2);
     461           7 :       v2 = gcoeff(f,2,1); v1 = gcoeff(f,2,2); b = mulii(v1, v2);
     462           7 :       c = mulii(addii(u1, v1), addii(u2, v2)); break;
     463             :     }
     464    21360738 :     if (gc_needed(av, 2))
     465             :     {
     466          96 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadunit (%ld)", i);
     467          96 :       gerepileall(av, 4, &p, &f, &q,&q1);
     468             :     }
     469             :   }
     470           7 :   u = diviiexact(addmulii(a, D, b), q);
     471           7 :   v = diviiexact(subii(c, addii(a, b)), q);
     472           7 :   if (m == 1) u = subii(u, v);
     473           7 :   *pu = shifti(u, -1); *pv = v;
     474           7 : }
     475             : GEN
     476       69048 : quadunit(GEN D0)
     477             : {
     478       69048 :   pari_sp av = avma;
     479             :   GEN P, E, D, u, v;
     480       69048 :   long s = signe(D0);
     481             :   /* check_quaddisc_real omitting test for squares */
     482       69048 :   if (typ(D0) != t_INT) pari_err_TYPE("quadunit", D0);
     483       69041 :   if (s <= 0) pari_err_DOMAIN("quadunit", "disc","<=",gen_0,D0);
     484       69041 :   if (mod4(D0) > 1) pari_err_DOMAIN("quadunit","disc % 4",">", gen_1,D0);
     485       69041 :   D = coredisc2_fact(Z_factor(D0), s, &P, &E);
     486             :   /* test for squares done here for free */
     487       69041 :   if (equali1(D)) pari_err_DOMAIN("quadunit","issquare(disc)","=", gen_1,D0);
     488       69041 :   if (cmpiu(D, 2000000) < 0)
     489       69034 :     quadunit_uv_basecase(D, &u, &v);
     490             :   else
     491           7 :     quadunit_uv(D, &u, &v);
     492       69041 :   if (lg(P) != 1)
     493             :   { /* non-trivial conductor N > 1 */
     494       26446 :     GEN N = factorback2(P,E), qD = quadpoly_i(D);
     495       26446 :     GEN n, y = deg1pol_shallow(v, u, 0); /* maximal order fund unit */
     496       26446 :     n = quadunitindex_ii(D, N, mkvec2(P,E), FpX_red(y,N), qD); /* unit index */
     497       26446 :     y = ZXQ_powu(y, itou(n), qD); /* fund unit of order of conductor N */
     498       26446 :     v = gel(y,3); u = gel(y,2); /* u + v w_D */
     499       26446 :     if (mpodd(D))
     500             :     { /* w_D = (1+sqrt(D))/2 */
     501       17353 :       if (mpodd(D0))
     502             :       { /* w_D0 = (1 + N sqrt(D)) / 2 */
     503        6167 :         GEN v0 = v;
     504        6167 :         v = diviiexact(v, N);
     505        6167 :         u = addii(u, shifti(subii(v0, v), -1));
     506             :       }
     507             :       else
     508             :       { /* w_D0 = N sqrt(D)/2, N is even */
     509       11186 :         v = shifti(v, -1);
     510       11186 :         u = addii(u, v);
     511       11186 :         v = diviiexact(v, shifti(N,-1));
     512             :       }
     513             :     }
     514             :     else /* w_D = sqrt(D), w_D0 = N sqrt(D) */
     515        9093 :       v = diviiexact(v, N);
     516             :   }
     517       69041 :   return gerepilecopy(av, mkquad(quadpoly_i(D0), u, v));
     518             : }
     519             : long
     520       69027 : quadunitnorm(GEN D)
     521             : {
     522       69027 :   pari_sp av = avma;
     523             :   long s, r;
     524       69027 :   check_quaddisc(D, &s, &r, "quadunitnorm");
     525       69020 :   if (s < 0) return 1;
     526       69020 :   (void)quadunit_q(D, sqrti(D), &s); return gc_long(av, s);
     527             : }
     528             : 
     529             : GEN
     530          49 : quadregulator(GEN x, long prec)
     531             : {
     532          49 :   pari_sp av = avma, av2;
     533             :   GEN R, rsqd, u, v, sqd;
     534             :   long r, e;
     535             : 
     536          49 :   check_quaddisc_real(x, &r, "quadregulator");
     537          49 :   sqd = sqrti(x);
     538          49 :   rsqd = gsqrt(x,prec); av2 = avma;
     539          49 :   e = 0; R = real2n(1, prec); u = utoi(r); v = gen_2;
     540             :   for(;;)
     541         140 :   {
     542         189 :     GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
     543         189 :     GEN v1 = divii(subii(x,sqri(u1)),v);
     544         189 :     if (equalii(v,v1)) { R = mulrr(sqrr(R), divri(addir(u1,rsqd),v)); break; }
     545         154 :     if (equalii(u,u1)) { R = sqrr(R); break; }
     546         140 :     R = mulrr(R, divri(addir(u1,rsqd),v));
     547         140 :     e += expo(R); setexpo(R,0);
     548         140 :     u = u1; v = v1;
     549         140 :     if (e & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
     550         140 :     if (gc_needed(av2,2))
     551             :     {
     552           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
     553           0 :       gerepileall(av2,3, &R,&u,&v);
     554             :     }
     555             :   }
     556          49 :   R = divri(R, v); e = 2*e - 1;
     557             :   /* avoid loss of accuracy */
     558          49 :   if (!((e + expo(R)) & ~EXPOBITS)) { setexpo(R, e + expo(R)); e = 0; }
     559          49 :   R = logr_abs(R);
     560          49 :   if (e) R = addrr(R, mulsr(e, mplog2(prec)));
     561          49 :   return gerepileuptoleaf(av, R);
     562             : }
     563             : 
     564             : /*************************************************************************/
     565             : /**                                                                     **/
     566             : /**                            CLASS NUMBER                             **/
     567             : /**                                                                     **/
     568             : /*************************************************************************/
     569             : 
     570             : int
     571     9204105 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
     572             : 
     573    10288260 : static GEN qfi_pow(void *E, GEN f, GEN n)
     574    10288260 : { return E? nupow(f,n,(GEN)E): qfbpow_i(f,n); }
     575     7015563 : static GEN qfi_comp(void *E, GEN f, GEN g)
     576     7015563 : { return E? nucomp(f,g,(GEN)E): qfbcomp_i(f,g); }
     577             : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
     578             :                                          gidentical,qfb_equal1,NULL};
     579             : 
     580             : GEN
     581     3370157 : qfi_order(GEN q, GEN o)
     582     3370157 : { return gen_order(q, o, NULL, &qfi_group); }
     583             : 
     584             : GEN
     585           0 : qfi_log(GEN a, GEN g, GEN o)
     586           0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
     587             : 
     588             : GEN
     589      722080 : qfi_Shanks(GEN a, GEN g, long n)
     590             : {
     591      722080 :   pari_sp av = avma;
     592             :   GEN T, X;
     593             :   long rt_n, c;
     594             : 
     595      722080 :   a = qfbred_i(a);
     596      722080 :   g = qfbred_i(g);
     597             : 
     598      722080 :   rt_n = sqrt((double)n);
     599      722080 :   c = n / rt_n;
     600      722080 :   c = (c * rt_n < n + 1) ? c + 1 : c;
     601             : 
     602      722080 :   T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
     603      722080 :   X = gen_Shanks(T, a, c, NULL, &qfi_group);
     604      722080 :   return X? gerepileuptoint(av, X): gc_NULL(av);
     605             : }
     606             : 
     607             : GEN
     608         462 : qfbclassno0(GEN x,long flag)
     609             : {
     610         462 :   switch(flag)
     611             :   {
     612         448 :     case 0: return map_proto_G(classno,x);
     613          14 :     case 1: return map_proto_G(classno2,x);
     614           0 :     default: pari_err_FLAG("qfbclassno");
     615             :   }
     616             :   return NULL; /* LCOV_EXCL_LINE */
     617             : }
     618             : 
     619             : /* f^h = 1, return order(f). Set *pfao to its factorization */
     620             : static GEN
     621        5397 : find_order(void *E, GEN f, GEN h, GEN *pfao)
     622             : {
     623        5397 :   GEN v = gen_factored_order(f, h, E, &qfi_group);
     624        5397 :   *pfao = gel(v,2); return gel(v,1);
     625             : }
     626             : 
     627             : static int
     628         301 : ok_q(GEN q, GEN h, GEN d2, long r2)
     629             : {
     630         301 :   if (d2)
     631             :   {
     632         245 :     if (r2 <= 2 && !mpodd(q)) return 0;
     633         245 :     return is_pm1(Z_ppo(q,d2));
     634             :   }
     635             :   else
     636             :   {
     637          56 :     if (r2 <= 1 && !mpodd(q)) return 0;
     638          56 :     return is_pm1(Z_ppo(q,h));
     639             :   }
     640             : }
     641             : 
     642             : /* a,b given by their factorizations. Return factorization of lcm(a,b).
     643             :  * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
     644             : static GEN
     645         182 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
     646             : {
     647         182 :   GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
     648         182 :   GEN A = gen_1, B = gen_1;
     649         182 :   long i, l = lg(P);
     650         182 :   GEN E = cgetg(l, t_COL);
     651         707 :   for (i=1; i<l; i++)
     652             :   {
     653         525 :     GEN p = gel(P,i);
     654         525 :     long va = Z_pval(a,p);
     655         525 :     long vb = Z_pval(b,p);
     656         525 :     if (va < vb)
     657             :     {
     658         203 :       B = mulii(B,powiu(p,vb));
     659         203 :       gel(E,i) = utoi(vb);
     660             :     }
     661             :     else
     662             :     {
     663         322 :       A = mulii(A,powiu(p,va));
     664         322 :       gel(E,i) = utoi(va);
     665             :     }
     666             :   }
     667         182 :   *pA = A;
     668         182 :   *pB = B; return mkmat2(P,E);
     669             : }
     670             : 
     671             : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
     672             : static void
     673         182 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
     674             : {
     675         182 :   GEN A,B, g1 = *pg1, d1 = *pd1;
     676         182 :   *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
     677         182 :   *pg1 = gmul(qfbpow_i(g1, diviiexact(d1,A)),  qfbpow_i(f, diviiexact(o,B)));
     678         182 :   *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
     679         182 : }
     680             : 
     681             : /* Let s = 1 or -1; D = s * d; assume Df^2 fits in an ulong
     682             :  * Return  f / [O_{Df^2}^*:O_D^*] * \prod_{p|f}  [ 1 - (D/p) p^-1 ]
     683             :  * The Euler product is \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
     684             : ulong
     685        7081 : uquadclassnoF_fact(ulong d, long s, GEN P, GEN E)
     686             : {
     687        7081 :   long i, l = lg(P);
     688        7081 :   ulong H = 1;
     689        9054 :   for (i = 1; i < l; i++)
     690             :   {
     691        1973 :     ulong p = P[i], e = E[i];
     692        1973 :     long D = (long)(p == 2? d & 7: d % p), a;
     693        1973 :     if (s < 0) D = -D;
     694        1973 :     a = kross(D,p);
     695        1973 :     if (!a)
     696         504 :       H *= upowuu(p, e);
     697             :     else
     698             :     {
     699        1469 :       H *= p - a;
     700        1469 :       if (e >= 2) H *= upowuu(p, e-1);
     701             :     }
     702             :   }
     703        7081 :   if (l == 1) return H;
     704        1490 :   if (s < 0)
     705             :   {
     706        1470 :     switch(d)
     707             :     { /* divide by [ O_K^* : O^* ] */
     708         105 :       case 4: H >>= 1; break;
     709         287 :       case 3: H /= 3; break;
     710             :     }
     711             :   }
     712             :   else
     713             :   {
     714          20 :     GEN fa = mkmat2(zc_to_ZC(P), zc_to_ZC(E));
     715          20 :     H /= itou(quadunitindex_i(utoipos(d), factorback(fa), fa));
     716             :   }
     717        1490 :   return H;
     718             : }
     719             : GEN
     720          48 : quadclassnoF_fact(GEN D, GEN P, GEN E)
     721             : {
     722          48 :   GEN H = quadclassnoEuler_fact(D, P, E);
     723          48 :   if (lg(P) == 1) return H;
     724          15 :   if (signe(D) < 0)
     725             :   {
     726           7 :     switch(itou_or_0(D))
     727             :     { /* divide by [ O_K^* : O^* ] */
     728           0 :       case 4: H = shifti(H,-1); break;
     729           0 :       case 3: H = diviuexact(H,3); break;
     730             :     }
     731             :   }
     732             :   else
     733             :   {
     734           8 :     GEN fa = mkvec2(P, E);
     735           8 :     H = diviiexact(H, quadunitindex_i(D, factorback2(P, E), fa));
     736             :   }
     737          15 :   return H;
     738             : }
     739             : 
     740             : static ulong
     741         652 : quadclassnoF_u(ulong x, long s, ulong *pD)
     742             : {
     743         652 :   pari_sp av = avma;
     744             :   GEN P, E;
     745         652 :   ulong D = coredisc2u_fact(factoru(x), s, &P, &E);
     746         652 :   long H = uquadclassnoF_fact(D, s, P, E);
     747         652 :   *pD = D; return gc_ulong(av, H);
     748             : }
     749             : ulong
     750           0 : unegquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, -1, pD); }
     751             : ulong
     752           0 : uposquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, 1, pD); }
     753             : 
     754             : /* *pD = coredisc(x), *pR = regulator (x > 0) or NULL */
     755             : GEN
     756         700 : quadclassnoF(GEN x, GEN *pD)
     757             : {
     758             :   GEN D, P, E;
     759         700 :   if (lgefint(x) == 3)
     760             :   {
     761         652 :     long s = signe(x);
     762         652 :     ulong d, h = quadclassnoF_u(x[2], s, &d);
     763         652 :     if (pD) *pD = s > 0? utoipos(d): utoineg(d);
     764         652 :     return utoipos(h);
     765             :   }
     766          48 :   D = coredisc2_fact(absZ_factor(x), signe(x), &P, &E);
     767          48 :   if (pD) *pD = D;
     768          48 :   return quadclassnoF_fact(D, P, E);
     769             : }
     770             : 
     771             : static long
     772         651 : two_rank(GEN x)
     773             : {
     774         651 :   GEN p = gel(absZ_factor(x),1);
     775         651 :   long l = lg(p)-1;
     776             : #if 0 /* positive disc not needed */
     777             :   if (signe(x) > 0)
     778             :   {
     779             :     long i;
     780             :     for (i=1; i<=l; i++)
     781             :       if (mod4(gel(p,i)) == 3) { l--; break; }
     782             :   }
     783             : #endif
     784         651 :   return l-1;
     785             : }
     786             : 
     787             : static GEN
     788       12334 : sqr_primeform(GEN x, ulong p) { return qfbsqr_i(primeform_u(x, p)); }
     789             : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
     790             : static GEN
     791         651 : get_forms(GEN D, GEN *pL)
     792             : {
     793         651 :   const long MAXFORM = 20;
     794         651 :   GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
     795         651 :   GEN forms = vectrunc_init(MAXFORM+1);
     796         651 :   long s, nforms = 0;
     797             :   ulong p;
     798             :   forprime_t S;
     799         651 :   L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
     800         651 :   s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
     801         651 :   if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
     802         651 :   if      (s < 10)   s = 200;
     803         427 :   else if (s < 20)   s = 1000;
     804         413 :   else if (s < 5000) s = 5000;
     805         651 :   u_forprime_init(&S, 2, s);
     806    13471626 :   while ( (p = u_forprime_next(&S)) )
     807             :   {
     808    13470975 :     long d, k = kroiu(D,p);
     809             :     pari_sp av2;
     810    13470975 :     if (!k) continue;
     811    13469841 :     if (k > 0)
     812             :     {
     813     6736219 :       if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
     814     6736219 :       d = p-1;
     815             :     }
     816             :     else
     817     6733622 :       d = p+1;
     818    13469841 :     av2 = avma; affrr(divru(mulur(p,L),d), L); set_avma(av2);
     819             :   }
     820         651 :   *pL = L; return forms;
     821             : }
     822             : 
     823             : /* h ~ #G, return o = order of f, set fao = its factorization */
     824             : static  GEN
     825         749 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
     826             : {
     827         749 :   long s = minss(itos(sqrti(h)), 10000);
     828         749 :   GEN T = gen_Shanks_init(f, s, E, &qfi_group);
     829         749 :   GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
     830         749 :   return find_order(E, f, addiu(v,1), pfao);
     831             : }
     832             : 
     833             : /* if g = 1 in  G/<f> ? */
     834             : static int
     835        5684 : equal1(void *E, GEN T, ulong N, GEN g)
     836        5684 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
     837             : 
     838             : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
     839             :  * FIXME: should be gen_order, but equal1 has the wrong prototype */
     840             : static GEN
     841         350 : relative_order(void *E, GEN a, GEN o, ulong N,  GEN T)
     842             : {
     843         350 :   pari_sp av = avma;
     844             :   long i, l;
     845             :   GEN m;
     846             : 
     847         350 :   m = get_arith_ZZM(o);
     848         350 :   if (!m) pari_err_TYPE("gen_order [missing order]",a);
     849         350 :   o = gel(m,1);
     850         350 :   m = gel(m,2); l = lgcols(m);
     851        1148 :   for (i = l-1; i; i--)
     852             :   {
     853         798 :     GEN t, y, p = gcoeff(m,i,1);
     854         798 :     long j, e = itos(gcoeff(m,i,2));
     855         798 :     if (l == 2) {
     856          49 :       t = gen_1;
     857          49 :       y = a;
     858             :     } else {
     859         749 :       t = diviiexact(o, powiu(p,e));
     860         749 :       y = powgi(a, t);
     861             :     }
     862         798 :     if (equal1(E, T,N,y)) o = t;
     863             :     else {
     864         364 :       for (j = 1; j < e; j++)
     865             :       {
     866          84 :         y = powgi(y, p);
     867          84 :         if (equal1(E, T,N,y)) break;
     868             :       }
     869         357 :       if (j < e) {
     870          77 :         if (j > 1) p = powiu(p, j);
     871          77 :         o = mulii(t, p);
     872             :       }
     873             :     }
     874             :   }
     875         350 :   return gerepilecopy(av, o);
     876             : }
     877             : 
     878             : /* h(x) for x<0 using Baby Step/Giant Step.
     879             :  * Assumes G is not too far from being cyclic.
     880             :  *
     881             :  * Compute G^2 instead of G so as to kill most of the noncyclicity */
     882             : GEN
     883         798 : classno(GEN x)
     884             : {
     885         798 :   pari_sp av = avma;
     886             :   long r2, k, s, i, l;
     887             :   GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
     888             :   void *E;
     889             : 
     890         798 :   if (signe(x) >= 0) return classno2(x);
     891             : 
     892         763 :   check_quaddisc(x, &s, &k, "classno");
     893         763 :   if (abscmpiu(x,12) <= 0) return gen_1;
     894             : 
     895         651 :   Hf = quadclassnoF(x, &D);
     896         651 :   if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
     897         651 :   forms =  get_forms(D, &L);
     898         651 :   r2 = two_rank(D);
     899         651 :   hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
     900             : 
     901         651 :   l = lg(forms);
     902         651 :   order_bound = const_vec(l-1, NULL);
     903         651 :   E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
     904         651 :   g1 = gel(forms,1);
     905         651 :   gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
     906         651 :   q = diviiround(hin, d1); /* approximate order of G/<g1> */
     907         651 :   d2 = NULL; /* not computed yet */
     908         651 :   if (is_pm1(q)) goto END;
     909        7315 :   for (i=2; i < l; i++)
     910             :   {
     911        6951 :     GEN o, fao, a, F, fd, f = gel(forms,i);
     912        6951 :     fd = qfbpow_i(f, d1); if (is_pm1(gel(fd,1))) continue;
     913         182 :     F = qfbpow_i(fd, q);
     914         182 :     a = gel(F,1);
     915         182 :     o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
     916             :     /* f^(d1 q) = 1 */
     917         182 :     fao = ZM_merge_factor(fad1,fao);
     918         182 :     o = find_order(E, f, fao, &fao);
     919         182 :     gel(order_bound,i) = o;
     920             :     /* o = order of f, fao = factor(o) */
     921         182 :     update_g1(&g1,&d1,&fad1, f,o,fao);
     922         182 :     q = diviiround(hin, d1);
     923         182 :     if (is_pm1(q)) goto END;
     924             :   }
     925             :   /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
     926         364 :   if (expi(q) > 3)
     927             :   { /* q large: compute d2, 2nd elt divisor */
     928         308 :     ulong N, n = 2*itou(sqrti(d1));
     929         308 :     GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
     930         308 :     d2 = gen_1;
     931         308 :     N = itou( gceil(gdivgu(d1,n)) ); /* order(g1) <= n*N */
     932        5047 :     for (i = 1; i < l; i++)
     933             :     {
     934        4802 :       GEN d, f = gel(forms,i), B = gel(order_bound,i);
     935        4802 :       if (!B) B = find_order(E, f, fad1, /*junk*/&d);
     936        4802 :       f = qfbpow_i(f,d2);
     937        4802 :       if (equal1(E, T,N,f)) continue;
     938         350 :       B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
     939             :       /* f^B = 1 */
     940         350 :       d = relative_order(E, f, B, N,T);
     941         350 :       d2= mulii(d,d2);
     942         350 :       D = mulii(d1,d2);
     943         350 :       q = diviiround(hin,D);
     944         350 :       if (is_pm1(q)) { d1 = D; goto END; }
     945             :     }
     946             :     /* very probably, d2 is the 2nd elementary divisor */
     947         245 :     d1 = D; /* product of first two elt divisors */
     948             :   }
     949             :   /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
     950             :    * 2-rank */
     951         301 :   if (!ok_q(q,d1,d2,r2))
     952             :   {
     953           0 :     GEN q0 = q;
     954             :     long d;
     955           0 :     if (cmpii(mulii(q,d1), hin) < 0)
     956             :     { /* try q = q0+1,-1,+2,-2 */
     957           0 :       d = 1;
     958           0 :       do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
     959             :     }
     960             :     else
     961             :     { /* q0-1,+1,-2,+2  */
     962           0 :       d = -1;
     963           0 :       do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
     964             :     }
     965             :   }
     966         301 :   d1 = mulii(d1,q);
     967             : 
     968         651 : END:
     969         651 :   return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
     970             : }
     971             : 
     972             : /* use Euler products */
     973             : GEN
     974          49 : classno2(GEN x)
     975             : {
     976          49 :   pari_sp av = avma;
     977          49 :   const long prec = DEFAULTPREC;
     978             :   long n, i, s;
     979          49 :   GEN p1, p2, S, p4, p5, p7, Hf, Pi, logd, sqrtd, D, half, reg = NULL;
     980             : 
     981          49 :   check_quaddisc(x, &s, NULL, "classno2");
     982          49 :   if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
     983             : 
     984          49 :   Hf = quadclassnoF(x, &D);
     985          49 :   if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
     986             : 
     987          49 :   Pi = mppi(prec);
     988          49 :   sqrtd = sqrtr_abs(itor(D, prec));
     989          49 :   logd = logr_abs(sqrtd); shiftr_inplace(logd,1);
     990          49 :   p1 = sqrtr_abs(divrr(mulir(D,logd), gmul2n(Pi,1)));
     991          49 :   if (s > 0)
     992             :   {
     993          42 :     GEN invlogd = invr(logd);
     994          42 :     reg = quadregulator(D, prec);
     995          42 :     p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
     996          42 :     if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
     997             :   }
     998          49 :   n = itos_or_0( mptrunc(p1) );
     999          49 :   if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
    1000             : 
    1001          49 :   p4 = divri(Pi, D); setsigne(p4, 1);
    1002          49 :   p7 = invr(sqrtr_abs(Pi));
    1003          49 :   half = real2n(-1, prec);
    1004          49 :   if (s > 0)
    1005             :   { /* i = 1, shortcut */
    1006          42 :     p1 = sqrtd;
    1007          42 :     p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
    1008          42 :     S = addrr(mulrr(p1,p5), eint1(p4,prec));
    1009        1358 :     for (i=2; i<=n; i++)
    1010             :     {
    1011        1316 :       long k = kroiu(D,i); if (!k) continue;
    1012        1218 :       p2 = mulir(sqru(i), p4);
    1013        1218 :       p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
    1014        1218 :       p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
    1015        1218 :       S = (k>0)? addrr(S,p5): subrr(S,p5);
    1016             :     }
    1017          42 :     S = shiftr(divrr(S,reg),-1);
    1018             :   }
    1019             :   else
    1020             :   { /* i = 1, shortcut */
    1021           7 :     p1 = gdiv(sqrtd, Pi);
    1022           7 :     p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
    1023           7 :     S = addrr(p5, divrr(p1, mpexp(p4)));
    1024         952 :     for (i=2; i<=n; i++)
    1025             :     {
    1026         945 :       long k = kroiu(D,i); if (!k) continue;
    1027         945 :       p2 = mulir(sqru(i), p4);
    1028         945 :       p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
    1029         945 :       p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
    1030         945 :       S = (k>0)? addrr(S,p5): subrr(S,p5);
    1031             :     }
    1032             :   }
    1033          49 :   return gerepileuptoint(av, mulii(Hf, roundr(S)));
    1034             : }
    1035             : 
    1036             : /* 1 + q + ... + q^v, v > 0 */
    1037             : static GEN
    1038         973 : geomsumu(ulong q, long v)
    1039             : {
    1040         973 :   GEN u = utoipos(1+q);
    1041        1267 :   for (; v > 1; v--) u = addui(1, mului(q, u));
    1042         973 :   return u;
    1043             : }
    1044             : static GEN
    1045         973 : geomsum(GEN q, long v)
    1046             : {
    1047             :   GEN u;
    1048         973 :   if (lgefint(q) == 3) return geomsumu(q[2], v);
    1049           0 :   u = addiu(q,1);
    1050           0 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
    1051           0 :   return u;
    1052             : }
    1053             : 
    1054             : /* 1+p+...+p^(e-1), e >= 1; assuming result fits in an ulong */
    1055             : static ulong
    1056        9426 : usumpow(ulong p, long e)
    1057             : {
    1058             :   ulong q;
    1059             :   long i;
    1060        9426 :   if (p == 2) return (1UL << e) - 1; /* also OK if e = BITS_IN_LONG */
    1061        6643 :   e--; for (i = 1, q = p + 1; i < e; i++) q = p*q + 1;
    1062        6587 :   return q;
    1063             : }
    1064             : long
    1065      163463 : uhclassnoF_fact(GEN faF, long D)
    1066             : {
    1067      163463 :   GEN P = gel(faF,1), E = gel(faF,2);
    1068      163463 :   long i, t, l = lg(P);
    1069      357280 :   for (i = t = 1; i < l; i++)
    1070             :   {
    1071      193817 :     long p = P[i], e = E[i], s = kross(D,p);
    1072      193817 :     if (e == 1) { t *= 1 + p - s; continue; }
    1073       51718 :     if (s == 1) { t *= upowuu(p,e); continue; }
    1074        9426 :     t *= 1 + usumpow(p, e) * (p - s);
    1075             :   }
    1076      163463 :   return t;
    1077             : }
    1078             : /* Hurwitz(D F^2)/ Hurwitz(D)
    1079             :  * = \sum_{f|F}  f \prod_{p|f} (1-kro(D/p)/p)
    1080             :  * = \prod_{p^e || F} (1 + (p^e-1) / (p-1) * (p-kro(D/p))) */
    1081             : GEN
    1082       60723 : hclassnoF_fact(GEN P, GEN E, GEN D)
    1083             : {
    1084             :   GEN H;
    1085       60723 :   long i, l = lg(P);
    1086       60723 :   if (l == 1) return gen_1;
    1087       56167 :   for (i = 1, H = NULL; i < l; i++)
    1088             :   {
    1089       30243 :     GEN t, p = gel(P,i);
    1090       30243 :     long e = E[i], s = kronecker(D,p);
    1091       30244 :     if (e == 1) t = addiu(p, 1-s);
    1092        1555 :     else if (s == 1) t = powiu(p, e);
    1093         973 :     else t = addui(1, mulii(subis(p, s), geomsum(p, e-1)));
    1094       30243 :     H = H? mulii(H,t): t;
    1095             :   }
    1096       25924 :   return H;
    1097             : }
    1098             : static GEN
    1099       60724 : hclassno6_large(GEN x)
    1100             : {
    1101       60724 :   GEN H = NULL, P, E, D = coredisc2_fact(absZ_factor(x), -1, &P, &E);
    1102       60720 :   long l = lg(P);
    1103             : 
    1104       60720 :   if (l > 1 && lgefint(x) == 3)
    1105             :   { /* F != 1, second chance */
    1106       25926 :     ulong h = hclassno6u_no_cache(x[2]);
    1107       25927 :     if (h) H = utoipos(h);
    1108             :   }
    1109       60721 :   if (!H)
    1110             :   {
    1111       60721 :     H = quadclassno(D);
    1112       60728 :     switch(itou_or_0(D))
    1113             :     {
    1114          49 :       case 3: H = shifti(H,1);break;
    1115           7 :       case 4: H = muliu(H,3); break;
    1116       60671 :       default:H = muliu(H,6); break;
    1117             :     }
    1118             :   }
    1119       60723 :   return mulii(H, hclassnoF_fact(P, E, D));
    1120             : }
    1121             : 
    1122             : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
    1123             : GEN
    1124      132478 : hclassno6(GEN x)
    1125             : {
    1126      132478 :   ulong d = itou_or_0(x);
    1127      132478 :   if (d)
    1128             :   { /* create cache if d small */
    1129      132477 :     ulong h = d < 500000 ? hclassno6u(d): hclassno6u_no_cache(d);
    1130      132477 :     if (h) return utoipos(h);
    1131             :   }
    1132       60724 :   return hclassno6_large(x);
    1133             : }
    1134             : 
    1135             : GEN
    1136          49 : hclassno(GEN x)
    1137             : {
    1138             :   long a, s;
    1139          49 :   if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
    1140          49 :   s = signe(x);
    1141          49 :   if (s < 0) return gen_0;
    1142          49 :   if (!s) return gdivgs(gen_1, -12);
    1143          49 :   a = mod4(x); if (a == 1 || a == 2) return gen_0;
    1144          49 :   return gdivgu(hclassno6(x), 6);
    1145             : }
    1146             : 
    1147             : /* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
    1148             : GEN
    1149     2551766 : Zn_quad_roots(GEN N, GEN B, GEN C)
    1150             : {
    1151     2551766 :   pari_sp av = avma;
    1152     2551766 :   GEN fa = NULL, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
    1153             :   long l, i, j, ct;
    1154             : 
    1155     2551766 :   if ((fa = check_arith_non0(N,"Zn_quad_roots")))
    1156             :   {
    1157        7665 :     N = typ(N) == t_VEC? gel(N,1): factorback(N);
    1158        7665 :     fa = clean_Z_factor(fa);
    1159             :   }
    1160     2551766 :   N = absi_shallow(N);
    1161     2551766 :   N4 = shifti(N,2);
    1162     2551765 :   D = modii(subii(sqri(B), shifti(C,2)), N4);
    1163     2551763 :   if (!signe(D))
    1164             :   { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
    1165         693 :     if (!fa) fa = Z_factor(N);
    1166         693 :     P = gel(fa,1);
    1167         693 :     E = ZV_to_zv(gel(fa,2));
    1168         693 :     l = lg(P);
    1169        1519 :     for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
    1170         693 :     Np = factorback2(P, E); /* x = -B mod N' */
    1171         693 :     B = shifti(B,-1);
    1172         693 :     return gerepilecopy(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
    1173             :   }
    1174     2551070 :   if (!fa)
    1175     2543576 :     fa = Z_factor(N4);
    1176             :   else  /* convert to factorization of N4 = 4*N */
    1177        7494 :     fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
    1178     2551073 :   P = gel(fa,1); l = lg(P);
    1179     2551073 :   E = ZV_to_zv(gel(fa,2));
    1180     2551073 :   F = cgetg(l, t_VEC);
    1181     2551073 :   mF= cgetg(l, t_VEC); F0 = gen_0;
    1182     2551073 :   Q = cgetg(l, t_VEC); Q0 = gen_1;
    1183     6107395 :   for (i = j = 1, ct = 0; i < l; i++)
    1184             :   {
    1185     5596122 :     GEN p = gel(P,i), q, f, mf, D0;
    1186     5596122 :     long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
    1187     5596120 :     if (d <= 0)
    1188             :     {
    1189     1252292 :       q = powiu(p, (s+1)>>1);
    1190     2162306 :       Q0 = mulii(Q0, q); continue;
    1191             :     }
    1192             :     /* d > 0 */
    1193     6198128 :     if (odd(t)) return NULL;
    1194     4158331 :     t2 = t >> 1;
    1195     4158331 :     if (i > 1)
    1196             :     { /* p > 2 */
    1197     2603157 :       if (kronecker(D0, p) == -1) return NULL;
    1198     1250877 :       q = powiu(p, s - t2);
    1199     1250878 :       f = Zp_sqrt(D0, p, d);
    1200     1250879 :       if (!f) return NULL; /* p was not actually prime... */
    1201     1250879 :       if (t2) f = mulii(powiu(p,t2), f);
    1202     1250879 :       mf = Fp_neg(f, q);
    1203             :     }
    1204             :     else
    1205             :     { /* p = 2 */
    1206     1555174 :       if (d <= 3)
    1207             :       {
    1208     1185247 :         if (d == 3 && Mod8(D0) != 1) return NULL;
    1209      960918 :         if (d == 2 && Mod4(D0) != 1) return NULL;
    1210      910014 :         Q0 = int2n(1+t2); F0 = NULL; continue;
    1211             :       }
    1212      369927 :       if (Mod8(D0) != 1) return NULL;
    1213      143143 :       q = int2n(d - 1 + t2);
    1214      143143 :       f = Z2_sqrt(D0, d);
    1215      143143 :       if (t2) f = shifti(f, t2);
    1216      143143 :       mf = Fp_neg(f, q);
    1217             :     }
    1218     1394020 :     gel(Q,j) = q;
    1219     1394020 :     gel(F,j) = f;
    1220     1394020 :     gel(mF,j)= mf; j++;
    1221             :   }
    1222      511273 :   setlg(Q,j);
    1223      511273 :   setlg(F,j);
    1224      511273 :   setlg(mF,j);
    1225      511273 :   if (is_pm1(Q0)) A = leafcopy(F);
    1226             :   else
    1227             :   { /* append the fixed congruence (F0 mod Q0) */
    1228      476350 :     if (!F0) F0 = shifti(Q0,-1);
    1229      476350 :     A = shallowconcat(F, F0);
    1230      476350 :     Q = shallowconcat(Q, Q0);
    1231             :   }
    1232      511273 :   ct = 1 << (j-1);
    1233      511273 :   T = ZV_producttree(Q);
    1234      511272 :   R = ZV_chinesetree(Q,T);
    1235      511271 :   Np = gmael(T, lg(T)-1, 1);
    1236      511271 :   B = modii(B, Np);
    1237      511271 :   if (!signe(B)) B = NULL;
    1238      511271 :   Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
    1239      511271 :   w = cgetg(3, t_VEC);
    1240      511271 :   gel(w,1) = icopy(Np);
    1241      511272 :   gel(w,2) = v = cgetg(ct+1, t_VEC);
    1242      511272 :   l = lg(F);
    1243     2359177 :   for (j = 1; j <= ct; j++)
    1244             :   {
    1245     1847904 :     pari_sp av2 = avma;
    1246     1847904 :     long m = j - 1;
    1247             :     GEN u;
    1248     5602263 :     for (i = 1; i < l; i++)
    1249             :     {
    1250     3754359 :       gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
    1251     3754359 :       m >>= 1;
    1252             :     }
    1253     1847904 :     u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
    1254     1847902 :     if (B) u = subii(u,B);
    1255     1847902 :     gel(v,j) = gerepileuptoint(av2, modii(shifti(u,-1), Np));
    1256             :   }
    1257      511273 :   return gerepileupto(av, w);
    1258             : }

Generated by: LCOV version 1.16