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 - Qfb.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30316-0578a48613) Lines: 1095 1216 90.0 %
Date: 2025-05-31 09:20:01 Functions: 132 142 93.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000-2005  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             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : /*******************************************************************/
      17             : /*                                                                 */
      18             : /*         QUADRATIC POLYNOMIAL ASSOCIATED TO A DISCRIMINANT       */
      19             : /*                                                                 */
      20             : /*******************************************************************/
      21             : 
      22             : void
      23      151584 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
      24             : {
      25             :   long r;
      26      151584 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
      27      151570 :   *s = signe(x);
      28      151570 :   if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
      29      151570 :   r = mod4(x); if (*s < 0 && r) r = 4 - r;
      30      151570 :   if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
      31      151556 :   if (pr) *pr = r;
      32      151556 : }
      33             : void
      34        6916 : check_quaddisc_real(GEN x, long *r, const char *f)
      35             : {
      36        6916 :   long sx; check_quaddisc(x, &sx, r, f);
      37        6916 :   if (sx < 0) pari_err_DOMAIN(f, "disc","<",gen_0,x);
      38        6916 : }
      39             : void
      40        2177 : check_quaddisc_imag(GEN x, long *r, const char *f)
      41             : {
      42        2177 :   long sx; check_quaddisc(x, &sx, r, f);
      43        2170 :   if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
      44        2170 : }
      45             : 
      46             : /* X^2 + b X + c is the canonical quadratic t_POL of discriminant D.
      47             :  * Dodd is nonzero iff D is odd */
      48             : static void
      49      978317 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
      50             : {
      51      978317 :   if (Dodd)
      52             :   {
      53      879708 :     pari_sp av = avma;
      54      879708 :     *b = gen_m1;
      55      879708 :     *c = gc_INT(av, shifti(subui(1,D), -2));
      56             :   }
      57             :   else
      58             :   {
      59       98609 :     *b = gen_0;
      60       98609 :     *c = shifti(D,-2); togglesign(*c);
      61             :   }
      62      978317 : }
      63             : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
      64             : static GEN
      65      243341 : quadpoly_ii(GEN D, long Dmod4)
      66             : {
      67      243341 :   GEN b, c, y = cgetg(5,t_POL);
      68      243341 :   y[1] = evalsigne(1) | evalvarn(0);
      69      243341 :   quadpoly_bc(D, Dmod4, &b,&c);
      70      243341 :   gel(y,2) = c;
      71      243341 :   gel(y,3) = b;
      72      243341 :   gel(y,4) = gen_1; return y;
      73             : }
      74             : GEN
      75        2044 : quadpoly(GEN D)
      76             : {
      77             :   long s, Dmod4;
      78        2044 :   check_quaddisc(D, &s, &Dmod4, "quadpoly");
      79        2037 :   return quadpoly_ii(D, Dmod4);
      80             : }
      81             : GEN /* no checks */
      82      241304 : quadpoly_i(GEN D) { return quadpoly_ii(D, Mod4(D)); }
      83             : 
      84             : GEN
      85        1029 : quadpoly0(GEN x, long v)
      86             : {
      87        1029 :   GEN T = quadpoly(x);
      88        1022 :   if (v > 0) setvarn(T, v);
      89        1022 :   return T;
      90             : }
      91             : 
      92             : GEN
      93           0 : quadgen(GEN x)
      94           0 : { retmkquad(quadpoly(x), gen_0, gen_1); }
      95             : 
      96             : GEN
      97         574 : quadgen0(GEN x, long v)
      98             : {
      99         574 :   if (v==-1) v = fetch_user_var("w");
     100         574 :   retmkquad(quadpoly0(x, v), gen_0, gen_1);
     101             : }
     102             : 
     103             : /***********************************************************************/
     104             : /**                                                                   **/
     105             : /**                      BINARY QUADRATIC FORMS                       **/
     106             : /**                                                                   **/
     107             : /***********************************************************************/
     108             : static int
     109      814212 : is_qfi(GEN q) { return typ(q)==t_QFB && qfb_is_qfi(q); }
     110             : 
     111             : static GEN
     112     1846188 : check_qfbext(const char *fun, GEN x)
     113             : {
     114     1846188 :   long t = typ(x);
     115     1846188 :   if (t == t_QFB) return x;
     116         196 :   if (t == t_VEC && lg(x)==3)
     117             :   {
     118         196 :     GEN q = gel(x,1);
     119         196 :     if (!is_qfi(q) && typ(gel(x,2))==t_REAL) return q;
     120             :   }
     121           0 :   pari_err_TYPE(fun, x);
     122             :   return NULL;/* LCOV_EXCL_LINE */
     123             : }
     124             : 
     125             : static GEN
     126      139081 : qfb3(GEN x, GEN y, GEN z)
     127      139081 : { retmkqfb(icopy(x), icopy(y), icopy(z), qfb_disc3(x,y,z)); }
     128             : 
     129             : static int
     130    23783633 : qfb_equal(GEN x, GEN y)
     131             : {
     132    23783633 :   return equalii(gel(x,1),gel(y,1))
     133     1592911 :       && equalii(gel(x,2),gel(y,2))
     134    25376542 :       && equalii(gel(x,3),gel(y,3));
     135             : }
     136             : 
     137             : /* valid for t_QFB, qfr3, qfr5; shallow */
     138             : static GEN
     139      892304 : qfb_inv(GEN x) {
     140      892304 :   GEN z = shallowcopy(x);
     141      892303 :   gel(z,2) = negi(gel(z,2));
     142      892304 :   return z;
     143             : }
     144             : /* valid for t_QFB, GC clean */
     145             : static GEN
     146           7 : qfbinv(GEN x)
     147           7 : { retmkqfb(icopy(gel(x,1)),negi(gel(x,2)),icopy(gel(x,3)), icopy(gel(x,4))); }
     148             : 
     149             : GEN
     150       77231 : Qfb0(GEN a, GEN b, GEN c)
     151             : {
     152             :   GEN q, D;
     153       77231 :   if (!b)
     154             :   {
     155          49 :     if (c) pari_err_TYPE("Qfb",c);
     156          42 :     if (typ(a) == t_VEC && lg(a) == 4)
     157          21 :     { b = gel(a,2); c = gel(a,3); a = gel(a,1); }
     158          21 :     else if (typ(a) == t_POL && degpol(a) == 2)
     159           7 :     { b = gel(a,3); c = gel(a,2); a = gel(a,4); }
     160          14 :     else if (typ(a) == t_MAT && lg(a)==3 && lgcols(a)==3)
     161             :     {
     162           7 :       b = gadd(gcoeff(a,2,1), gcoeff(a,1,2));
     163           7 :       c = gcoeff(a,2,2); a = gcoeff(a,1,1);
     164             :     }
     165             :     else
     166           7 :       pari_err_TYPE("Qfb",a);
     167             :   }
     168       77182 :   else if (!c)
     169           7 :     pari_err_TYPE("Qfb",b);
     170       77210 :   if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
     171       77203 :   if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
     172       77203 :   if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
     173       77203 :   q = qfb3(a, b, c); D = qfb_disc(q);
     174       77203 :   if (signe(D) < 0)
     175       42392 :   { if (signe(a) < 0) pari_err_IMPL("negative definite t_QFB"); }
     176       34811 :   else if (Z_issquare(D)) pari_err_DOMAIN("Qfb","issquare(disc)","=", gen_1,q);
     177       77196 :   return q;
     178             : }
     179             : 
     180             : /***********************************************************************/
     181             : /**                                                                   **/
     182             : /**                         Reduction                                 **/
     183             : /**                                                                   **/
     184             : /***********************************************************************/
     185             : 
     186             : /* assume a > 0. Write b = q*2a + r, with -a < r <= a */
     187             : static GEN
     188    16738610 : dvmdii_round(GEN b, GEN a, GEN *r)
     189             : {
     190    16738610 :   GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
     191    16738207 :   if (signe(b) >= 0) {
     192     9197443 :     if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
     193             :   } else { /* r <= 0 */
     194     7540764 :     if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
     195             :   }
     196    16737786 :   return q;
     197             : }
     198             : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
     199             : static long
     200   111670109 : dvmdsu_round(long b, ulong a, long *r)
     201             : {
     202   111670109 :   ulong a2 = a << 1, q, ub, ur;
     203   111670109 :   if (b >= 0) {
     204    71207677 :     ub = b;
     205    71207677 :     q = ub / a2;
     206    71207677 :     ur = ub % a2;
     207    71207677 :     if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
     208    26050729 :     else *r = (long)ur;
     209    71207677 :     return (long)q;
     210             :   } else { /* r <= 0 */
     211    40462432 :     ub = (ulong)-b; /* |b| */
     212    40462432 :     q = ub / a2;
     213    40462432 :     ur = ub % a2;
     214    40462432 :     if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
     215    21571040 :     else *r = -(long)ur;
     216    40462432 :     return -(long)q;
     217             :   }
     218             : }
     219             : /* reduce b mod 2*a. Update b,c */
     220             : static void
     221     2783988 : REDB(GEN a, GEN *b, GEN *c)
     222             : {
     223     2783988 :   GEN r, q = dvmdii_round(*b, a, &r);
     224     2783988 :   if (!signe(q)) return;
     225     2714529 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     226     2714529 :   *b = r;
     227             : }
     228             : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
     229             : static void
     230   111659600 : sREDB(ulong a, long *b, ulong *c)
     231             : {
     232             :   long r, q;
     233             :   ulong uz;
     234   121213473 :   if (a > LONG_MAX) return; /* b already reduced */
     235   111659600 :   q = dvmdsu_round(*b, a, &r);
     236   111690088 :   if (q == 0) return;
     237             :   /* Final (a,r,c2) satisfies |r| <= |b| hence c2 <= c, c2 = c - q*z,
     238             :    * where z = (b+r) / 2, representable as long, has the same sign as q. */
     239   102136215 :   if (*b < 0)
     240             :   { /* uz = -z >= 0, q < 0 */
     241    35635890 :     if (r >= 0) /* different signs=>no overflow, exact division */
     242    18947614 :       uz = (ulong)-((*b + r)>>1);
     243             :     else
     244             :     {
     245    16688276 :       ulong ub = (ulong)-*b, ur = (ulong)-r;
     246    16688276 :       uz = (ub + ur) >> 1;
     247             :     }
     248    35635890 :     *c -= (-q) * uz; /* c -= qz */
     249             :   }
     250             :   else
     251             :   { /* uz = z >= 0, q > 0 */
     252    66500325 :     if (r <= 0)
     253    45218664 :       uz = (*b + r)>>1;
     254             :     else
     255             :     {
     256    21281661 :       ulong ub = (ulong)*b, ur = (ulong)r;
     257    21281661 :       uz = ((ub + ur) >> 1);
     258             :     }
     259    66500325 :     *c -= q * uz; /* c -= qz */
     260             :   }
     261   102136215 :   *b = r;
     262             : }
     263             : static void
     264    13954627 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
     265             : { /* REDB(a,b,c) */
     266    13954627 :   GEN r, q = dvmdii_round(*b, a, &r);
     267    13953782 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     268    13953638 :   *b = r;
     269    13953638 :   *u2 = subii(*u2, mulii(q, u1));
     270    13953754 : }
     271             : 
     272             : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
     273             : static GEN
     274     6633700 : qfi_redsl2_basecase(GEN q, GEN *U)
     275             : {
     276     6633700 :   pari_sp av = avma;
     277             :   GEN z, u1,u2,v1,v2,Q;
     278     6633700 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     279             :   long cmp;
     280     6633700 :   u1 = gen_1; u2 = gen_0;
     281     6633700 :   cmp = abscmpii(a, b);
     282     6633707 :   if (cmp < 0)
     283     2110359 :     REDBU(a,&b,&c, u1,&u2);
     284     4523348 :   else if (cmp == 0 && signe(b) < 0)
     285             :   { /* b = -a */
     286       11767 :     b = negi(b);
     287       11767 :     u2 = gen_1;
     288             :   }
     289             :   for(;;)
     290             :   {
     291    18477271 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     292    11843174 :     swap(a,c); b = negi(b);
     293    11844236 :     z = u1; u1 = u2; u2 = negi(z);
     294    11844263 :     REDBU(a,&b,&c, u1,&u2);
     295    11843579 :     if (gc_needed(av, 1)) {
     296           7 :       if (DEBUGMEM>1) pari_warn(warnmem, "qfbredsl2");
     297           7 :       (void)gc_all(av, 5, &a,&b,&c, &u1,&u2);
     298             :     }
     299             :   }
     300     6633661 :   if (cmp == 0 && signe(b) < 0)
     301             :   {
     302       17677 :     b = negi(b);
     303       17677 :     z = u1; u1 = u2; u2 = negi(z);
     304             :   }
     305             :   /* Let q = (A,B,C). q o [u1,u2; v1,v2] = Q implies
     306             :    * [v1,v2] = (1/C) [(b-B)/2 u1 - a u2, c u1 - (b+B)/2 u2] */
     307     6633661 :   z = shifti(subii(b, gel(q,2)), -1);
     308     6633649 :   v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
     309     6633644 :   z = subii(z, b);
     310     6633658 :   v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
     311     6633645 :   *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
     312     6633714 :   Q = mkqfb(a,b,c,gel(q,4));
     313     6633711 :   return gc_all(av, 2, &Q, U);
     314             : }
     315             : 
     316             : static GEN
     317     1063492 : setq_b0(ulong a, ulong c, GEN D)
     318     1063492 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
     319             : /* assume |sb| = 1 */
     320             : static GEN
     321    95089342 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
     322    95089342 : { retmkqfb(utoipos(a), sb==1? utoipos(b): utoineg(b), utoipos(c), icopy(D)); }
     323             : /* 0 < a, c < 2^BIL, b = 0 */
     324             : static GEN
     325      949110 : qfi_red_1_b0(ulong a, ulong c, GEN D)
     326      949110 : { return (a <= c)? setq_b0(a, c, D): setq_b0(c, a, D); }
     327             : 
     328             : /* 0 < a, c < 2^BIL: single word affair */
     329             : static GEN
     330    96248537 : qfi_red_1(pari_sp av, GEN a, GEN b, GEN c, GEN D)
     331             : {
     332             :   ulong ua, ub, uc;
     333             :   long sb;
     334             :   for(;;)
     335      141369 :   { /* at most twice */
     336    96248537 :     long lb = lgefint(b); /* <= 3 after first loop */
     337    96248537 :     if (lb == 2) return qfi_red_1_b0(a[2],c[2], D);
     338    95299427 :     if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
     339      139482 :     REDB(a,&b,&c);
     340      140504 :     if (uel(a,2) <= uel(c,2))
     341             :     { /* lg(b) <= 3 but may be too large for itos */
     342           0 :       long s = signe(b);
     343           0 :       set_avma(av);
     344           0 :       if (!s) return qfi_red_1_b0(a[2], c[2], D);
     345           0 :       if (a[2] == c[2]) s = 1;
     346           0 :       return setq(a[2], b[2], c[2], s, D);
     347             :     }
     348      140504 :     swap(a,c); b = negi(b);
     349             :   }
     350             :   /* b != 0 */
     351    95161019 :   set_avma(av);
     352    95160671 :   ua = a[2];
     353    95160671 :   ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
     354    95160671 :   uc = c[2];
     355    95160671 :   if (ua < ub)
     356    35930955 :     sREDB(ua, &sb, &uc);
     357    59229716 :   else if (ua == ub && sb < 0) sb = (long)ub;
     358   170954251 :   while(ua > uc)
     359             :   {
     360    75752686 :     lswap(ua,uc); sb = -sb;
     361    75752686 :     sREDB(ua, &sb, &uc);
     362             :   }
     363    95201565 :   if (!sb) return setq_b0(ua, uc, D);
     364             :   else
     365             :   {
     366    95087182 :     long s = 1;
     367    95087182 :     if (sb < 0)
     368             :     {
     369    35860482 :       sb = -sb;
     370    35860482 :       if (ua != uc) s = -1;
     371             :     }
     372    95087182 :     return setq(ua, sb, uc, s, D);
     373             :   }
     374             : }
     375             : 
     376             : static GEN
     377           7 : qfi_rho(GEN x)
     378             : {
     379           7 :   pari_sp av = avma;
     380           7 :   GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
     381           7 :   int fl = abscmpii(a, c);
     382           7 :   if (fl <= 0)
     383             :   {
     384           7 :     int fg = abscmpii(a, b);
     385           7 :     if (fg >= 0)
     386             :     {
     387           7 :       x = gcopy(x);
     388           7 :       if ((!fl || !fg) && signe(gel(x,2)) < 0) setsigne(gel(x,2), 1);
     389           7 :       return x;
     390             :     }
     391             :   }
     392           0 :   swap(a,c); b = negi(b);
     393           0 :   REDB(a, &b, &c);
     394           0 :   return gc_GEN(av, mkqfb(a,b,c, qfb_disc(x)));
     395             : }
     396             : 
     397             : /* qfr3 / qfr5 */
     398             : 
     399             : /* t_QFB are unusable: D, sqrtD, isqrtD are recomputed all the time and the
     400             :  * logarithmic Shanks's distance is costly and hard to control.
     401             :  * qfr3 / qfr5 routines take a container of t_INTs (e.g a t_VEC) as argument,
     402             :  * at least 3 (resp. 5) components [it is a feature that they do not check the
     403             :  * precise type or length of the input]. They return a vector of length 3
     404             :  * (resp. 5). A qfr3 [a,b,c] contains the form coeffs, in a qfr5 [a,b,c, e,d]
     405             :  * the t_INT e is a binary exponent, d a t_REAL, coding the distance in
     406             :  * multiplicative form: the true distance is obtained from qfr5_dist.
     407             :  * All other qfr routines are obsolete (inefficient) wrappers */
     408             : 
     409             : /* static functions are not stack-clean. Unless mentionned otherwise, public
     410             :  * functions are. */
     411             : 
     412             : #define EMAX 22
     413             : static void
     414    10217928 : fix_expo(GEN x)
     415             : {
     416    10217928 :   if (expo(gel(x,5)) >= (1L << EMAX)) {
     417           0 :     gel(x,4) = addiu(gel(x,4), 1);
     418           0 :     shiftr_inplace(gel(x,5), - (1L << EMAX));
     419             :   }
     420    10217928 : }
     421             : 
     422             : /* (1/2) log (|d| * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
     423             : GEN
     424      184688 : qfr5_dist(GEN e, GEN d, long prec)
     425             : {
     426      184688 :   GEN t = logr_abs(d);
     427      184688 :   if (signe(e)) {
     428           0 :     GEN u = mulir(e, mplog2(prec));
     429           0 :     shiftr_inplace(u, EMAX); t = addrr(t, u);
     430             :   }
     431      184688 :   shiftr_inplace(t, -1); return t;
     432             : }
     433             : 
     434             : static void
     435    14152531 : rho_get_BC(GEN *B, GEN *C, GEN a, GEN b, GEN c, struct qfr_data *S)
     436             : {
     437             :   GEN t, u, q;
     438    14152531 :   t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: absi_shallow(c);
     439    14152531 :   q = truedvmdii(addii(t, b), shifti(c,1), &u);
     440    14152531 :   *B = subii(t, u); /* t - ((t+b) % 2c) */
     441    14152531 :   *C = subii(a, mulii(q, subii(b, mulii(q,c))));
     442    14152531 : }
     443             : /* Not stack-clean */
     444             : GEN
     445     1139446 : qfr3_rho(GEN x, struct qfr_data *S)
     446             : {
     447     1139446 :   GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3);
     448     1139446 :   rho_get_BC(&B, &C, a, b, c, S);
     449     1139446 :   return mkvec3(c, B, C);
     450             : }
     451             : 
     452             : /* Not stack-clean */
     453             : GEN
     454     8486681 : qfr5_rho(GEN x, struct qfr_data *S)
     455             : {
     456     8486681 :   GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3), y;
     457     8486681 :   long sb = signe(b);
     458     8486681 :   rho_get_BC(&B, &C, a, b, c, S);
     459     8486681 :   y = mkvec5(c, B, C, gel(x,4), gel(x,5));
     460     8486681 :   if (sb) {
     461     8482698 :     GEN t = subii(sqri(b), S->D);
     462     8482698 :     if (sb < 0)
     463     2509500 :       t = divir(t, sqrr(subir(b,S->sqrtD)));
     464             :     else
     465     5973198 :       t = divri(sqrr(addir(b,S->sqrtD)), t);
     466             :     /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
     467     8482698 :     gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
     468        3983 :   } else gel(y,5) = negr(gel(y,5));
     469     8486681 :   return y;
     470             : }
     471             : 
     472             : /* Not stack-clean */
     473             : GEN
     474      217728 : qfr_to_qfr5(GEN x, long prec)
     475      217728 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
     476             : 
     477             : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
     478             : GEN
     479         483 : qfr5_to_qfr(GEN x, GEN D, GEN d0)
     480             : {
     481         483 :   if (d0)
     482             :   {
     483         140 :     GEN n = gel(x,4), d = absr(gel(x,5));
     484         140 :     if (signe(n))
     485             :     {
     486           0 :       n = addis(shifti(n, EMAX), expo(d));
     487           0 :       setexpo(d, 0); d = logr_abs(d);
     488           0 :       if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
     489           0 :       shiftr_inplace(d, -1);
     490           0 :       d0 = addrr(d0, d);
     491             :     }
     492         140 :     else if (!gequal1(d)) /* avoid loss of precision */
     493             :     {
     494          91 :       d = logr_abs(d);
     495          91 :       shiftr_inplace(d, -1);
     496          91 :       d0 = addrr(d0, d);
     497             :     }
     498             :   }
     499         483 :   x = qfr3_to_qfr(x, D);
     500         483 :   return d0 ? mkvec2(x,d0): x;
     501             : }
     502             : 
     503             : /* Not stack-clean */
     504             : GEN
     505       31920 : qfr3_to_qfr(GEN x, GEN d) { retmkqfb(gel(x,1), gel(x,2), gel(x,3), d); }
     506             : 
     507             : static int
     508    17712973 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
     509             : {
     510             :   GEN t;
     511    17712973 :   if (signe(b) <= 0 || abscmpii(b, isqrtD) > 0) return 0;
     512     5271643 :   t = addii_sign(isqrtD,1, shifti(a,1),-1); /* floor(sqrt(D)) - |2a| */
     513     1099055 :   return signe(t) < 0 ? abscmpii(b, t) >= 0
     514     6370698 :                       : abscmpii(b, t) > 0;
     515             : }
     516             : 
     517             : /* Not stack-clean */
     518             : GEN
     519     1952895 : qfr5_red(GEN x, struct qfr_data *S)
     520             : {
     521     1952895 :   pari_sp av = avma;
     522     8465338 :   while (!ab_isreduced(gel(x,1), gel(x,2), S->isqrtD))
     523             :   {
     524     6512443 :     x = qfr5_rho(x, S);
     525     6512443 :     if (gc_needed(av,2))
     526             :     {
     527           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
     528           0 :       x = gc_GEN(av, x);
     529             :     }
     530             :   }
     531     1952895 :   return x;
     532             : }
     533             : /* Not stack-clean */
     534             : GEN
     535     1172847 : qfr3_red(GEN x, struct qfr_data *S)
     536             : {
     537     1172847 :   pari_sp av = avma;
     538     1172847 :   GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
     539     5699251 :   while (!ab_isreduced(a, b, S->isqrtD))
     540             :   {
     541             :     GEN B, C;
     542     4526404 :     rho_get_BC(&B, &C, a, b, c, S);
     543     4526404 :     a = c; b = B; c = C;
     544     4526404 :     if (gc_needed(av,2))
     545             :     {
     546           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
     547           0 :       (void)gc_all(av, 3, &a, &b, &c);
     548             :     }
     549             :   }
     550     1172847 :   return mkvec3(a, b, c);
     551             : }
     552             : 
     553             : void
     554        2170 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
     555             : {
     556        2170 :   S->D = D;
     557        2170 :   S->sqrtD = sqrtr(itor(S->D,prec));
     558        2170 :   S->isqrtD = truncr(S->sqrtD);
     559        2170 : }
     560             : 
     561             : static GEN
     562         140 : qfr5_init(GEN x, GEN d, struct qfr_data *S)
     563             : {
     564         140 :   long prec = realprec(d), l = -expo(d);
     565         140 :   if (l < BITS_IN_LONG) l = BITS_IN_LONG;
     566         140 :   prec = maxss(prec, nbits2prec(l));
     567         140 :   S->D = qfb_disc(x);
     568         140 :   x = qfr_to_qfr5(x,prec);
     569         140 :   if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
     570           0 :   else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
     571             : 
     572         140 :   if (!S->isqrtD)
     573             :   {
     574         126 :     pari_sp av=avma;
     575             :     long e;
     576         126 :     S->isqrtD = gcvtoi(S->sqrtD,&e);
     577         126 :     if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
     578             :   }
     579          14 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
     580         140 :   return x;
     581             : }
     582             : static GEN
     583         371 : qfr3_init(GEN x, struct qfr_data *S)
     584             : {
     585         371 :   S->D = qfb_disc(x);
     586         371 :   if (!S->isqrtD) S->isqrtD = sqrti(S->D);
     587         280 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
     588         371 :   return x;
     589             : }
     590             : 
     591             : #define qf_NOD  2
     592             : #define qf_STEP 1
     593             : 
     594             : static GEN
     595         427 : qfr_red_basecase_i(GEN x, long flag, GEN isqrtD, GEN sqrtD)
     596             : {
     597             :   struct qfr_data S;
     598         427 :   GEN d = NULL, y;
     599         427 :   if (typ(x)==t_VEC) { d = gel(x,2); x = gel(x,1); } else flag |= qf_NOD;
     600         427 :   S.sqrtD = sqrtD;
     601         427 :   S.isqrtD = isqrtD;
     602         427 :   y = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, d, &S);
     603         427 :   switch(flag) {
     604          63 :     case 0:              y = qfr5_red(y,&S); break;
     605         336 :     case qf_NOD:         y = qfr3_red(y,&S); break;
     606          21 :     case qf_STEP:        y = qfr5_rho(y,&S); break;
     607           7 :     case qf_STEP|qf_NOD: y = qfr3_rho(y,&S); break;
     608           0 :     default: pari_err_FLAG("qfbred");
     609             :   }
     610         427 :   return qfr5_to_qfr(y, qfb_disc(x), d);
     611             : }
     612             : 
     613             : static void
     614    13379357 : qfr_rhosl2_i(GEN *pa, GEN *pb, GEN *pc, GEN *pu1, GEN *pu2, GEN *pv1,
     615             :              GEN *pv2, GEN rd)
     616             : {
     617    13379357 :   GEN C = mpabs_shallow(*pc), t = addii(*pb, gmax_shallow(rd,C));
     618    13379357 :   GEN r, q = truedvmdii(t, shifti(C,1), &r);
     619    13379357 :   GEN a = *pa, b = *pb, c = *pc;
     620    13379357 :   if (signe(c) < 0) togglesign(q);
     621    13379357 :   *pa = *pc;
     622    13379357 :   *pb = subii(t, addii(r, *pb));
     623    13379357 :   *pc = subii(a, mulii(q, subii(b, mulii(q,c))));
     624    13379357 :   r = *pu1; *pu1 = *pv1; *pv1 = subii(mulii(q, *pv1), r);
     625    13379357 :   r = *pu2; *pu2 = *pv2; *pv2 = subii(mulii(q, *pv2), r);
     626    13379357 : }
     627             : 
     628             : static GEN
     629    10810674 : qfr_rhosl2(GEN A, GEN rd)
     630             : {
     631    10810674 :   GEN V = gel(A,1), M = gel(A,2);
     632    10810674 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
     633    10810674 :   GEN u1 = gcoeff(M,1,1), v1 = gcoeff(M,1,2);
     634    10810674 :   GEN u2 = gcoeff(M,2,1), v2 = gcoeff(M,2,2);
     635    10810674 :   qfr_rhosl2_i(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
     636    10810674 :   return mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2));
     637             : }
     638             : 
     639             : static GEN
     640      979701 : qfr_redsl2_basecase(GEN V, GEN rd)
     641             : {
     642      979701 :   pari_sp av = avma;
     643      979701 :   GEN u1 = gen_1, u2 = gen_0, v1 = gen_0, v2 = gen_1;
     644      979701 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
     645     3548384 :   while (!ab_isreduced(a,b,rd))
     646             :   {
     647     2568683 :     qfr_rhosl2_i(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
     648     2568683 :     if (gc_needed(av, 1))
     649             :     {
     650           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfbredsl2");
     651           0 :       (void)gc_all(av, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
     652             :     }
     653             :   }
     654      979701 :   return gc_GEN(av, mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2)));
     655             : }
     656             : 
     657             : /* fast reduction of qfb with positive coefficients, based on
     658             : Arnold Schoenhage, Fast reduction and composition of binary quadratic forms,
     659             : Proc. of Intern. Symp. on Symbolic and Algebraic Computation (Bonn) (S. M.
     660             : Watt, ed.), ACM Press, 1991, pp. 128-133.
     661             : <https://dl.acm.org/doi/pdf/10.1145/120694.120711>
     662             : With thanks to Keegan Ryan
     663             : BA20230927
     664             : */
     665             : 
     666             : /* pqfb: qf with positive coefficients */
     667             : 
     668             : static int
     669     4940497 : lti2n(GEN a, long m) { return signe(a) < 0 || expi(a) < m;}
     670             : 
     671             : static GEN
     672     1921185 : pqfbred_1(GEN Q, long m, GEN U)
     673             : {
     674     1921185 :   GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     675     1921185 :   if (abscmpii(a, c) < 0)
     676             :   {
     677             :     GEN t, at, r;
     678      960317 :     GEN r2 = addii(shifti(a, m + 2), d);
     679      960317 :     long e2 = expi(r2);
     680      960317 :     r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
     681      960317 :     t = truedivii(subii(b, r), shifti(a,1));
     682      960317 :     if (signe(t)==0) pari_err_BUG("pqfbred_1");
     683      960317 :     at = mulii(a,t);
     684      960317 :     c = addii(subii(c, mulii(b, t)), mulii(at, t));
     685      960317 :     b = subii(b, shifti(at,1));
     686      960317 :     gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,1,1), t));
     687      960317 :     gcoeff(U,2,2) = subii( gcoeff(U,2,2), mulii(gcoeff(U,2,1), t));
     688             :   } else
     689             :   {
     690             :     GEN t, ct, r;
     691      960868 :     GEN r2 = addii(shifti(c, m + 2), d);
     692      960868 :     long e2 = expi(r2);
     693      960868 :     r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
     694      960868 :     t = truedivii(subii(b, r), shifti(c,1));
     695      960868 :     if (signe(t)==0) pari_err_BUG("pqfbred_1");
     696      960868 :     ct = mulii(c, t);
     697      960868 :     a = addii(subii(a, mulii(b, t)), mulii(ct, t));
     698      960868 :     b = subii(b, shifti(ct, 1));
     699      960868 :     gcoeff(U,1,1) = subii(gcoeff(U,1,1), mulii(gcoeff(U,1,2), t));
     700      960868 :     gcoeff(U,2,1) = subii(gcoeff(U,2,1), mulii(gcoeff(U,2,2), t));
     701             :   }
     702     1921185 :   return mkqfb(a,b,c,d);
     703             : }
     704             : 
     705             : static int
     706     2046352 : is_minimal(GEN Q, long m)
     707             : {
     708     2046352 :   pari_sp av = avma;
     709     2046352 :   GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3);
     710     4940497 :   return gc_bool(av, lti2n(addii(subii(a,b), c), m)
     711     1927231 :                  || (lti2n(subii(b, shifti(a,1)), m+1)
     712      966914 :                      && lti2n(subii(b, shifti(c,1)), m+1)));
     713             : }
     714             : 
     715             : static GEN
     716      124020 : pqfbred_iter_1(GEN Q, ulong m, GEN U)
     717             : {
     718      124020 :   pari_sp av = avma;
     719     1923484 :   while (!is_minimal(Q,m))
     720             :   {
     721     1799464 :     Q = pqfbred_1(Q, m, U);
     722     1799464 :     if (gc_needed(av, 1))
     723             :     {
     724           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"pqfbred_iter_1, lc = %ld", expi(gel(Q,3)));
     725           0 :       (void)gc_all(av, 3, &Q, &gel(U,1), &gel(U,2));
     726             :     }
     727             :   }
     728      124020 :   return Q;
     729             : }
     730             : 
     731             : static GEN
     732       61492 : pqfbred_basecase(GEN Q, ulong m, GEN *pt_U)
     733             : {
     734       61492 :   pari_sp av = avma;
     735       61492 :   GEN  U = matid(2);
     736       61492 :   Q = pqfbred_iter_1(Q, m, U);
     737       61492 :   *pt_U = U;
     738       61492 :   return gc_all(av, 2, &Q, pt_U);
     739             : }
     740             : 
     741             : static long
     742   101644731 : qfb_maxexpi(GEN Q)
     743   101644731 : { return 1+maxss(expi(gel(Q,1)), maxss(expi(gel(Q,2)), expi(gel(Q,3)))); }
     744             : 
     745             : /* use asymptotically fast reduction ? */
     746             : static int
     747   101334640 : qfi_red_fast(GEN Q)
     748             : {
     749   101334640 :   const long QFBRED_LIMIT = 9000;
     750   101334640 :   return 2*qfb_maxexpi(Q) - expi(gel(Q,4)) > QFBRED_LIMIT;
     751             : }
     752             : 
     753             : static long
     754      125056 : qfb_minexpi(GEN Q)
     755             : {
     756      125056 :   long m =  minss(expi(gel(Q,1)), minss(expi(gel(Q,2)), expi(gel(Q,3))));
     757      125056 :   return m < 0 ? 0: m;
     758             : }
     759             : 
     760             : static GEN
     761      124020 : pqfbred_rec(GEN Q, long m, GEN *pt_U)
     762             : {
     763      124020 :   pari_sp av = avma;
     764      124020 :   GEN U, Q0, Q1, QR, d = qfb_disc(Q);
     765      124020 :   long h, n = qfb_maxexpi(Q) - m;
     766      124020 :   int going_to_r8 = 0;
     767             : 
     768      124020 :   if (n < 170) return pqfbred_basecase(Q, m, pt_U);
     769       62528 :   if (qfb_minexpi(Q) <= m + 2) { U = matid(2); QR = Q; }
     770             :   else
     771             :   {
     772             :     long p, mm;
     773       62528 :     if (m <= n) { mm = m; p = 0; Q1 = Q; }
     774             :     else
     775             :     {
     776       61878 :       mm = n; p = m + 1 - n;
     777       61878 :       Q0 = mkvec3(remi2n(gel(Q,1),p), remi2n(gel(Q,2),p), remi2n(gel(Q,3),p));
     778       61878 :       Q1 = qfb3(shifti(gel(Q,1),-p), shifti(gel(Q,2),-p), shifti(gel(Q,3),-p));
     779             :     }
     780       62528 :     h = mm + (n>>1);
     781       62528 :     if (qfb_minexpi(Q1) <= h) { U = matid(2); QR = Q1; }
     782             :     else
     783       62321 :       QR = pqfbred_rec(Q1, h, &U);
     784      184249 :     while (qfb_maxexpi(QR) > h)
     785             :     {
     786      122868 :       if (is_minimal(QR, mm)) { going_to_r8 = 1; break; }
     787      121721 :       QR = pqfbred_1(QR, mm, U);
     788             :     }
     789       62528 :     if (!going_to_r8)
     790             :     {
     791             :       GEN V;
     792       61381 :       QR = pqfbred_rec(QR, mm, &V);
     793       61381 :       U = ZM2_mul(U,V);
     794             :     }
     795       62528 :     if (p > 0)
     796             :     {
     797       61878 :       GEN Q0U = qfb_ZM_apply(Q0,U);
     798      123756 :       QR = mkqfb(addii(shifti(gel(QR,1), p), gel(Q0U,1)),
     799       61878 :                  addii(shifti(gel(QR,2), p), gel(Q0U,2)),
     800       61878 :                  addii(shifti(gel(QR,3), p), gel(Q0U,3)), d);
     801             :     }
     802             :   }
     803       62528 :   QR = pqfbred_iter_1(QR, m, U);
     804       62528 :   *pt_U = U; return gc_all(av, 2, &QR, pt_U);
     805             : }
     806             : 
     807             : static GEN
     808      209575 : qfr_redsl2(GEN Q, GEN isqrtD)
     809             : {
     810      209575 :   pari_sp av = avma;
     811      209575 :   if (!qfi_red_fast(Q))
     812      209575 :     return qfr_redsl2_basecase(Q, isqrtD);
     813             :   else
     814             :   {
     815           0 :     GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     816           0 :     GEN Qf, Qr, W, U, t = NULL;
     817           0 :     long sa = signe(a), sb;
     818           0 :     if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
     819           0 :     if (signe(c) < 0)
     820             :     {
     821             :       GEN at;
     822           0 :       t  = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
     823           0 :       at = mulii(a,t);
     824           0 :       c = addii(subii(c, mulii(b, t)), mulii(at, t));
     825           0 :       b = subii(b, shifti(at,1));
     826             :     }
     827           0 :     sb = signe(b);
     828           0 :     Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
     829           0 :     if (sa < 0)
     830           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     831           0 :     if (sb < 0)
     832             :     {
     833           0 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     834           0 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     835             :     }
     836           0 :     if (t)
     837             :     {
     838           0 :       gcoeff(U,1,1) = subii( gcoeff(U,1,1), mulii(gcoeff(U,2,1), t));
     839           0 :       gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,2,2), t));
     840             :     }
     841           0 :     W = qfr_redsl2_basecase(Qr, isqrtD);
     842           0 :     Qf = gel(W,1);
     843           0 :     U = ZM2_mul(U,gel(W,2));
     844           0 :     return gc_GEN(av, mkvec2(Qf,U));
     845             :   }
     846             : }
     847             : 
     848             : static GEN
     849     5043354 : qfi_redsl2(GEN Q)
     850             : {
     851     5043354 :   pari_sp av = avma;
     852             :   GEN Qt, U;
     853     5043354 :   if (!qfi_red_fast(Q))
     854     5043048 :     Qt = qfi_redsl2_basecase(Q, &U);
     855             :   else
     856             :   {
     857         311 :     long sb = signe(gel(Q,2));
     858             :     GEN W;
     859         311 :     if (sb < 0) Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
     860         311 :     Q = pqfbred_rec(Q, 0, &U);
     861         311 :     Qt = qfi_redsl2_basecase(Q, &W);
     862         311 :     U = ZM2_mul(U,W);
     863         311 :     if (sb < 0)
     864             :     {
     865         173 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     866         173 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     867             :     }
     868             :   }
     869     5043377 :   return gc_GEN(av, mkvec2(Qt,U));
     870             : }
     871             : 
     872             : GEN
     873     4733039 : redimagsl2(GEN Q, GEN *U)
     874             : {
     875     4733039 :   GEN q = qfi_redsl2(Q);
     876     4733068 :   *U = gel(q,2); return gel(q,1);
     877             : }
     878             : 
     879             : GEN
     880      519892 : qfbredsl2(GEN q, GEN isD)
     881             : {
     882             :   pari_sp av;
     883      519892 :   if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
     884      519892 :   if (qfb_is_qfi(q))
     885             :   {
     886      310310 :     if (isD) pari_err_TYPE("qfbredsl2", isD);
     887      310310 :     return qfi_redsl2(q);
     888             :   }
     889      209582 :   av = avma;
     890      209582 :   if (!isD) isD = sqrti(qfb_disc(q));
     891      208068 :   else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
     892      209575 :   return gc_upto(av, qfr_redsl2(q, isD));
     893             : }
     894             : 
     895             : /* not gc-clean */
     896             : static GEN
     897         427 : qfr_red_i(GEN Q, long flag, GEN isqrtD, GEN sqrtD)
     898             : {
     899         427 :   if (typ(Q) != t_QFB || !qfi_red_fast(Q))
     900         427 :     return qfr_red_basecase_i(Q, flag, isqrtD, sqrtD);
     901             :   else
     902             :   {
     903           0 :     GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     904           0 :     GEN Qr, W, U, t = NULL;
     905           0 :     long sa = signe(a);
     906           0 :     if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
     907           0 :     if (signe(c) < 0)
     908             :     {
     909             :       GEN at;
     910           0 :       t  = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
     911           0 :       at = mulii(a,t);
     912           0 :       c = addii(subii(c, mulii(b, t)), mulii(at, t));
     913           0 :       b = subii(b, shifti(at,1));
     914             :     }
     915           0 :     Qr = pqfbred_rec(mkqfb(a, absi_shallow(b), c, d), 0, &U);
     916           0 :     if (sa < 0)
     917           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     918           0 :     W = qfr_red_basecase_i(Qr, flag, isqrtD, sqrtD);
     919           0 :     return gel(W,1);
     920             :   }
     921             : }
     922             : 
     923             : static GEN
     924          63 : qfr_red_av(pari_sp av, GEN x)
     925          63 : { return gc_GEN(av, qfr_red_i(x,0,NULL,NULL)); }
     926             : GEN
     927           0 : qfr_red(GEN x) { return qfr_red_av(avma, x); }
     928             : 
     929             : static GEN
     930    96284193 : qfi_red_basecase_av(pari_sp av, GEN q)
     931             : {
     932    96284193 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
     933    96284193 :   long cmp, lc = lgefint(c);
     934             : 
     935    96284193 :   if (lgefint(a) == 3 && lc == 3) return qfi_red_1(av, a, b, c, D);
     936      913667 :   cmp = abscmpii(a, b);
     937      919092 :   if (cmp < 0)
     938      436231 :     REDB(a,&b,&c);
     939      482861 :   else if (cmp == 0 && signe(b) < 0)
     940          27 :     b = negi(b);
     941             :   for(;;)
     942             :   {
     943     3126345 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     944     2932679 :     lc = lgefint(a); /* lg(future c): we swap a & c next */
     945     2932679 :     if (lc == 3) return qfi_red_1(av, a, b, c, D);
     946     2207253 :     swap(a,c); b = negi(b); /* apply rho */
     947     2207253 :     REDB(a,&b,&c);
     948     2207253 :     if (gc_needed(av, 2))
     949             :     {
     950           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfi_red, lc = %ld", lc);
     951           0 :       (void)gc_all(av, 3, &a,&b,&c);
     952             :     }
     953             :   }
     954      193666 :   if (cmp == 0 && signe(b) < 0) b = negi(b);
     955      193666 :   return gc_GEN(av, mkqfb(a, b, c, D));
     956             : }
     957             : static GEN
     958    96084300 : qfi_red_av(pari_sp av, GEN Q)
     959             : {
     960    96084300 :   if (qfi_red_fast(Q))
     961             :   {
     962             :     GEN U;
     963           7 :     if (signe(gel(Q,2)) < 0)
     964           0 :       Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
     965           7 :     Q = pqfbred_rec(Q, 0, &U);
     966             :   }
     967    96281050 :   return qfi_red_basecase_av(av, Q);
     968             : }
     969             : 
     970             : GEN
     971    22428261 : qfi_red(GEN q) { return qfi_red_av(avma, q); }
     972             : 
     973             : GEN
     974       54891 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
     975             : {
     976             :   pari_sp av;
     977       54891 :   GEN q = check_qfbext("qfbred",x);
     978       54891 :   if (qfb_is_qfi(q)) return (flag & qf_STEP)? qfi_rho(x): qfi_red(x);
     979         364 :   if (typ(x)==t_QFB) flag |= qf_NOD;
     980          49 :   else               flag &= ~qf_NOD;
     981         364 :   av = avma;
     982         364 :   return gc_GEN(av, qfr_red_i(x,flag,isqrtD,sqrtD));
     983             : }
     984             : /* t_QFB */
     985             : GEN
     986    15202468 : qfbred_i(GEN x) { return qfb_is_qfi(x)? qfi_red(x): qfr_red(x); }
     987             : GEN
     988       53078 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
     989             : /***********************************************************************/
     990             : /**                                                                   **/
     991             : /**                         Composition                               **/
     992             : /**                                                                   **/
     993             : /***********************************************************************/
     994             : 
     995             : static void
     996    33434850 : qfb_sqr(GEN z, GEN x)
     997             : {
     998             :   GEN c, d1, x2, v1, v2, c3, m, p1, r;
     999             : 
    1000    33434850 :   d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
    1001    33434727 :   c = gel(x,3);
    1002    33434727 :   m = mulii(c,x2);
    1003    33434482 :   if (equali1(d1))
    1004    24853205 :     v1 = v2 = gel(x,1);
    1005             :   else
    1006             :   {
    1007     8581330 :     v1 = diviiexact(gel(x,1),d1);
    1008     8581330 :     v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
    1009     8581331 :     c = mulii(c, d1);
    1010             :   }
    1011    33434537 :   togglesign(m);
    1012    33434439 :   r = modii(m,v2);
    1013    33434320 :   p1 = mulii(r, v1);
    1014    33434447 :   c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
    1015    33434421 :   gel(z,1) = mulii(v1,v2);
    1016    33434429 :   gel(z,2) = addii(gel(x,2), shifti(p1,1));
    1017    33434415 :   gel(z,3) = diviiexact(c3,v2);
    1018    33434253 : }
    1019             : /* z <- x * y */
    1020             : static void
    1021    75461366 : qfb_comp(GEN z, GEN x, GEN y)
    1022             : {
    1023             :   GEN n, c, d, y1, v1, v2, c3, m, p1, r;
    1024             : 
    1025    75461366 :   if (x == y) { qfb_sqr(z,x); return; }
    1026    42812319 :   n = shifti(subii(gel(y,2),gel(x,2)), -1);
    1027    42615332 :   v1 = gel(x,1);
    1028    42615332 :   v2 = gel(y,1);
    1029    42615332 :   c  = gel(y,3);
    1030    42615332 :   d = bezout(v2,v1,&y1,NULL);
    1031    42722067 :   if (equali1(d))
    1032    22886491 :     m = mulii(y1,n);
    1033             :   else
    1034             :   {
    1035    19874659 :     GEN s = subii(gel(y,2), n);
    1036    19872771 :     GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
    1037    19875811 :     if (!equali1(d1))
    1038             :     {
    1039    11019120 :       v1 = diviiexact(v1,d1);
    1040    11017837 :       v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
    1041    11017360 :       v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
    1042    11018154 :       c = mulii(c, d1);
    1043             :     }
    1044    19874589 :     m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
    1045             :   }
    1046    42706418 :   togglesign(m);
    1047    42635604 :   r = modii(m, v1);
    1048    42570795 :   p1 = mulii(r, v2);
    1049    42636287 :   c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
    1050    42632853 :   gel(z,1) = mulii(v1,v2);
    1051    42632408 :   gel(z,2) = addii(gel(y,2), shifti(p1,1));
    1052    42623685 :   gel(z,3) = diviiexact(c3,v1);
    1053             : }
    1054             : 
    1055             : /* not meant to be efficient */
    1056             : static GEN
    1057          84 : qfb_comp_gen(GEN x, GEN y)
    1058             : {
    1059          84 :   GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
    1060          84 :   GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
    1061          84 :   GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
    1062          84 :   GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
    1063             : 
    1064          84 :   if (!is_pm1(cx))
    1065             :   {
    1066          14 :     a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
    1067          14 :     c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
    1068             :   }
    1069          84 :   if (!is_pm1(cy))
    1070             :   {
    1071          28 :     a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
    1072          28 :     b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
    1073             :   }
    1074          84 :   D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
    1075         133 :   if (!Z_issquareall(diviiexact(d1, D), &n1) ||
    1076          84 :       !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
    1077          49 :   A = mulii(a1, n2);
    1078          49 :   B = mulii(a2, n1);
    1079          49 :   C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
    1080          49 :   U = ZV_extgcd(mkvec3(A, B, C));
    1081          49 :   m = gel(U,1); U = gmael(U,2,3);
    1082          49 :   A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
    1083          49 :   B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
    1084          49 :   C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
    1085          49 :   C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
    1086          49 :   B = addii(A, addii(B, C));
    1087          49 :   m2 = sqri(m);
    1088          49 :   A = diviiexact(mulii(a1, a2), m2);
    1089          49 :   C = diviiexact(shifti(subii(sqri(B),D), -2), A);
    1090          49 :   cx = mulii(cx, cy);
    1091          49 :   if (!is_pm1(cx))
    1092             :   {
    1093          14 :     A = mulii(A, cx); B = mulii(B, cx);
    1094          14 :     C = mulii(C, cx); D = mulii(D, sqri(cx));
    1095             :   }
    1096          49 :   return mkqfb(A, B, C, D);
    1097             : }
    1098             : 
    1099             : static GEN
    1100    72720126 : qficomp0(GEN x, GEN y, int raw)
    1101             : {
    1102    72720126 :   pari_sp av = avma;
    1103    72720126 :   GEN z = cgetg(5,t_QFB);
    1104    72717218 :   gel(z,4) = gel(x,4);
    1105    72717218 :   qfb_comp(z, x,y);
    1106    72455809 :   if (raw) return gc_GEN(av,z);
    1107    72454031 :   return qfi_red_av(av, z);
    1108             : }
    1109             : static GEN
    1110         441 : qfrcomp0(GEN x, GEN y, int raw)
    1111             : {
    1112         441 :   pari_sp av = avma;
    1113         441 :   GEN dx = NULL, dy = NULL;
    1114         441 :   GEN z = cgetg(5,t_QFB);
    1115         441 :   if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1116         441 :   if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
    1117         441 :   gel(z,4) = gel(x,4);
    1118         441 :   qfb_comp(z, x,y);
    1119         441 :   if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
    1120         441 :   if (raw) return gc_GEN(av, z);
    1121          28 :   return qfr_red_av(av, z);
    1122             : }
    1123             : /* same discriminant, no distance, no checks */
    1124             : GEN
    1125    27418402 : qfbcomp_i(GEN x, GEN y)
    1126    27418402 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
    1127             : GEN
    1128      131057 : qfbcomp(GEN x, GEN y)
    1129             : {
    1130      131057 :   GEN qx = check_qfbext("qfbcomp", x);
    1131      131057 :   GEN qy = check_qfbext("qfbcomp", y);
    1132      131057 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1133             :   {
    1134          63 :     pari_sp av = avma;
    1135          63 :     GEN z = qfb_comp_gen(qx, qy);
    1136          63 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1137           7 :       pari_err_IMPL("Shanks's distance in general composition");
    1138          56 :     if (!z) pari_err_OP("*",x,y);
    1139          21 :     return gc_upto(av, qfbred(z));
    1140             :   }
    1141      130994 :   return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
    1142             : }
    1143             : /* same discriminant, no distance, no checks */
    1144             : GEN
    1145           0 : qfbcompraw_i(GEN x, GEN y)
    1146           0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
    1147             : GEN
    1148        2198 : qfbcompraw(GEN x, GEN y)
    1149             : {
    1150        2198 :   GEN qx = check_qfbext("qfbcompraw", x);
    1151        2198 :   GEN qy = check_qfbext("qfbcompraw", y);
    1152        2198 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1153             :   {
    1154          21 :     pari_sp av = avma;
    1155          21 :     GEN z = qfb_comp_gen(qx, qy);
    1156          21 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1157           0 :       pari_err_IMPL("Shanks's distance in general composition");
    1158          21 :     if (!z) pari_err_OP("qfbcompraw",x,y);
    1159          21 :     return gc_GEN(av, z);
    1160             :   }
    1161        2177 :   if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
    1162        2177 :   return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
    1163             : }
    1164             : 
    1165             : static GEN
    1166      785763 : qfisqr0(GEN x, long raw)
    1167             : {
    1168      785763 :   pari_sp av = avma;
    1169      785763 :   GEN z = cgetg(5,t_QFB);
    1170      785763 :   gel(z,4) = gel(x,4);
    1171      785763 :   qfb_sqr(z,x);
    1172      785763 :   if (raw) return gc_GEN(av,z);
    1173      785763 :   return qfi_red_av(av, z);
    1174             : }
    1175             : static GEN
    1176          35 : qfrsqr0(GEN x, long raw)
    1177             : {
    1178          35 :   pari_sp av = avma;
    1179          35 :   GEN dx = NULL, z = cgetg(5,t_QFB);
    1180          35 :   if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1181          35 :   gel(z,4) = gel(x,4); qfb_sqr(z,x);
    1182          35 :   if (dx) z = mkvec2(z, shiftr(dx,1));
    1183          35 :   if (raw) return gc_GEN(av, z);
    1184          35 :   return qfr_red_av(av, z);
    1185             : }
    1186             : /* same discriminant, no distance, no checks */
    1187             : GEN
    1188      690174 : qfbsqr_i(GEN x)
    1189      690174 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
    1190             : GEN
    1191       95624 : qfbsqr(GEN x)
    1192             : {
    1193       95624 :   GEN qx = check_qfbext("qfbsqr", x);
    1194       95624 :   return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
    1195             : }
    1196             : 
    1197             : static GEN
    1198        6867 : qfr_1_by_disc(GEN D)
    1199             : {
    1200             :   GEN y, r, s;
    1201        6867 :   check_quaddisc_real(D, NULL, "qfr_1_by_disc");
    1202        6867 :   y = cgetg(5,t_QFB);
    1203        6867 :   s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
    1204        6867 :   if (mpodd(r))
    1205             :   {
    1206        3535 :     s = subiu(s,1);
    1207        3535 :     r = subii(r, addiu(shifti(s, 1), 1));
    1208        3535 :     r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
    1209             :   }
    1210             :   else
    1211        3332 :   { r = shifti(r, -2); set_avma((pari_sp)s); }
    1212        6867 :   gel(y,1) = gen_1;
    1213        6867 :   gel(y,2) = s;
    1214        6867 :   gel(y,3) = icopy(r);
    1215        6867 :   gel(y,4) = icopy(D); return y;
    1216             : }
    1217             : 
    1218             : static GEN
    1219          35 : qfr_disc(GEN x)
    1220          35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
    1221             : 
    1222             : static GEN
    1223          35 : qfr_1(GEN x)
    1224          35 : { return qfr_1_by_disc(qfr_disc(x)); }
    1225             : 
    1226             : static void
    1227           0 : qfr_1_fill(GEN y, struct qfr_data *S)
    1228             : {
    1229           0 :   pari_sp av = avma;
    1230           0 :   GEN y2 = S->isqrtD;
    1231           0 :   gel(y,1) = gen_1;
    1232           0 :   if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
    1233           0 :   gel(y,2) = y2; av = avma;
    1234           0 :   gel(y,3) = gc_INT(av, shifti(subii(sqri(y2), S->D),-2));
    1235           0 : }
    1236             : static GEN
    1237           0 : qfr5_1(struct qfr_data *S, long prec)
    1238             : {
    1239           0 :   GEN y = cgetg(6, t_VEC);
    1240           0 :   qfr_1_fill(y, S);
    1241           0 :   gel(y,4) = gen_0;
    1242           0 :   gel(y,5) = real_1(prec); return y;
    1243             : }
    1244             : static GEN
    1245           0 : qfr3_1(struct qfr_data *S)
    1246             : {
    1247           0 :   GEN y = cgetg(4, t_VEC);
    1248           0 :   qfr_1_fill(y, S); return y;
    1249             : }
    1250             : 
    1251             : /* Assume D < 0 is the discriminant of a t_QFB */
    1252             : static GEN
    1253      734976 : qfi_1_by_disc(GEN D)
    1254             : {
    1255      734976 :   GEN b,c, y = cgetg(5,t_QFB);
    1256      734976 :   quadpoly_bc(D, mod2(D), &b,&c);
    1257      734976 :   if (b == gen_m1) b = gen_1;
    1258      734976 :   gel(y,1) = gen_1;
    1259      734976 :   gel(y,2) = b;
    1260      734976 :   gel(y,3) = c;
    1261      734976 :   gel(y,4) = icopy(D); return y;
    1262             : }
    1263             : static GEN
    1264      722997 : qfi_1(GEN x)
    1265             : {
    1266      722997 :   if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
    1267      722997 :   return qfi_1_by_disc(qfb_disc(x));
    1268             : }
    1269             : 
    1270             : GEN
    1271           0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
    1272             : 
    1273             : static GEN
    1274    13610776 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
    1275             : static GEN
    1276    31561364 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
    1277             : static GEN
    1278           7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
    1279             : static GEN
    1280           7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
    1281             : 
    1282             : static GEN
    1283           7 : qfipowraw(GEN x, long n)
    1284             : {
    1285           7 :   pari_sp av = avma;
    1286             :   GEN y;
    1287           7 :   if (!n) return qfi_1(x);
    1288           7 :   if (n== 1) return gcopy(x);
    1289           7 :   if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
    1290           7 :   if (n < 0) x = qfb_inv(x);
    1291           7 :   y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
    1292           7 :   return gc_GEN(av,y);
    1293             : }
    1294             : 
    1295             : static GEN
    1296    15925465 : qfipow(GEN x, GEN n)
    1297             : {
    1298    15925465 :   pari_sp av = avma;
    1299             :   GEN y;
    1300    15925465 :   long s = signe(n);
    1301    15925465 :   if (!s) return qfi_1(x);
    1302    15202468 :   if (s < 0) x = qfb_inv(x);
    1303    15202467 :   y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
    1304    15202471 :   return gc_GEN(av,y);
    1305             : }
    1306             : 
    1307             : static long
    1308      412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
    1309             : {
    1310             :   long z;
    1311      412328 :   *v = gen_0; *v2 = gen_1;
    1312     4351417 :   for (z=0; abscmpii(*v3,L) > 0; z++)
    1313             :   {
    1314     3939089 :     GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
    1315     3939089 :     *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
    1316             :   }
    1317      412328 :   return z;
    1318             : }
    1319             : 
    1320             : /* composition: Shanks' NUCOMP & NUDUPL */
    1321             : /* L = floor((|d|/4)^(1/4)) */
    1322             : GEN
    1323      400722 : nucomp(GEN x, GEN y, GEN L)
    1324             : {
    1325      400722 :   pari_sp av = avma;
    1326             :   long z;
    1327             :   GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
    1328             : 
    1329      400722 :   if (x==y) return nudupl(x,L);
    1330      400680 :   if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
    1331      400680 :   if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
    1332             : 
    1333      400680 :   if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
    1334      400680 :   s = shifti(addii(gel(x,2),gel(y,2)), -1);
    1335      400680 :   n = subii(gel(y,2), s);
    1336      400680 :   a1 = gel(x,1);
    1337      400680 :   a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
    1338      400680 :   if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
    1339      163576 :   else if (dvdii(s,d)) /* d | s */
    1340             :   {
    1341       83503 :     a = negi(mulii(u,n)); d1 = d;
    1342       83503 :     a1 = diviiexact(a1, d1);
    1343       83503 :     a2 = diviiexact(a2, d1);
    1344       83503 :     s = diviiexact(s, d1);
    1345             :   }
    1346             :   else
    1347             :   {
    1348             :     GEN p2, l;
    1349       80073 :     d1 = bezout(s,d,&u1,NULL);
    1350       80073 :     if (!equali1(d1))
    1351             :     {
    1352        2044 :       a1 = diviiexact(a1,d1);
    1353        2044 :       a2 = diviiexact(a2,d1);
    1354        2044 :       s = diviiexact(s,d1);
    1355        2044 :       d = diviiexact(d,d1);
    1356             :     }
    1357       80073 :     p1 = remii(gel(x,3),d);
    1358       80073 :     p2 = remii(gel(y,3),d);
    1359       80073 :     l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
    1360       80073 :     a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
    1361             :   }
    1362      400680 :   a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
    1363      400680 :   d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
    1364      400680 :   Q = cgetg(5,t_QFB);
    1365      400680 :   if (!z) {
    1366       37632 :     g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
    1367       37632 :     b = a2;
    1368       37632 :     b2 = gel(y,2);
    1369       37632 :     v2 = d1;
    1370       37632 :     gel(Q,1) = mulii(d,b);
    1371             :   } else {
    1372             :     GEN e, q3, q4;
    1373      363048 :     if (z&1) { v3 = negi(v3); v2 = negi(v2); }
    1374      363048 :     b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
    1375      363048 :     e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
    1376      363048 :     q3 = mulii(e,v2);
    1377      363048 :     q4 = subii(q3,s);
    1378      363048 :     b2 = addii(q3,q4);
    1379      363048 :     g = diviiexact(q4,v);
    1380      363048 :     if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
    1381      363048 :     gel(Q,1) = addii(mulii(d,b), mulii(e,v));
    1382             :   }
    1383      400680 :   q1 = mulii(b, v3);
    1384      400680 :   q2 = addii(q1,n);
    1385      400680 :   gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
    1386      400680 :   gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
    1387      400680 :   gel(Q,4) = gel(x,4);
    1388      400680 :   return qfi_red_av(av, Q);
    1389             : }
    1390             : 
    1391             : GEN
    1392       11648 : nudupl(GEN x, GEN L)
    1393             : {
    1394       11648 :   pari_sp av = avma;
    1395             :   long z;
    1396             :   GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
    1397             : 
    1398       11648 :   if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
    1399       11648 :   a = gel(x,1);
    1400       11648 :   b = gel(x,2);
    1401       11648 :   d1 = bezout(b,a, &u,NULL);
    1402       11648 :   if (!equali1(d1))
    1403             :   {
    1404        4620 :     a = diviiexact(a, d1);
    1405        4620 :     b = diviiexact(b, d1);
    1406             :   }
    1407       11648 :   c = modii(negi(mulii(u,gel(x,3))), a);
    1408       11648 :   p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
    1409       11648 :   d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
    1410       11648 :   a2 = sqri(d);
    1411       11648 :   c2 = sqri(v3);
    1412       11648 :   Q = cgetg(5,t_QFB);
    1413       11648 :   if (!z) {
    1414        1281 :     g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
    1415        1281 :     b2 = gel(x,2);
    1416        1281 :     v2 = d1;
    1417        1281 :     gel(Q,1) = a2;
    1418             :   } else {
    1419             :     GEN e;
    1420       10367 :     if (z&1) { v = negi(v); d = negi(d); }
    1421       10367 :     e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
    1422       10367 :     g = diviiexact(subii(mulii(e,v2), b), v);
    1423       10367 :     b2 = addii(mulii(e,v2), mulii(v,g));
    1424       10367 :     if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
    1425       10367 :     gel(Q,1) = addii(a2, mulii(e,v));
    1426             :   }
    1427       11648 :   gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
    1428       11648 :   gel(Q,3) = addii(c2, mulii(g,v2));
    1429       11648 :   gel(Q,4) = gel(x,4);
    1430       11648 :   return qfi_red_av(av, Q);
    1431             : }
    1432             : 
    1433             : static GEN
    1434        4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
    1435             : static GEN
    1436       11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
    1437             : GEN
    1438        1008 : nupow(GEN x, GEN n, GEN L)
    1439             : {
    1440             :   pari_sp av;
    1441             :   GEN y, D;
    1442             : 
    1443        1008 :   if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
    1444        1008 :   if (!is_qfi(x)) pari_err_TYPE("nupow",x);
    1445        1008 :   if (gequal1(n)) return gcopy(x);
    1446        1008 :   av = avma;
    1447        1008 :   D = qfb_disc(x);
    1448        1008 :   y = qfi_1_by_disc(D);
    1449        1008 :   if (!signe(n)) return y;
    1450         959 :   if (!L) L = sqrtnint(absi_shallow(D), 4);
    1451         959 :   y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
    1452         959 :   if (signe(n) < 0
    1453          35 :   && !absequalii(gel(y,1),gel(y,2))
    1454          35 :   && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
    1455         959 :   return gc_GEN(av, y);
    1456             : }
    1457             : 
    1458             : /* Not stack-clean */
    1459             : GEN
    1460     1735230 : qfr5_compraw(GEN x, GEN y)
    1461             : {
    1462     1735230 :   GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
    1463     1735230 :   if (x == y)
    1464             :   {
    1465       34552 :     gel(z,4) = shifti(gel(x,4),1);
    1466       34552 :     gel(z,5) = sqrr(gel(x,5));
    1467             :   }
    1468             :   else
    1469             :   {
    1470     1700678 :     gel(z,4) = addii(gel(x,4),gel(y,4));
    1471     1700678 :     gel(z,5) = mulrr(gel(x,5),gel(y,5));
    1472             :   }
    1473     1735230 :   fix_expo(z); return z;
    1474             : }
    1475             : GEN
    1476     1735216 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
    1477     1735216 : { return qfr5_red(qfr5_compraw(x, y), S); }
    1478             : /* Not stack-clean */
    1479             : GEN
    1480     1009397 : qfr3_compraw(GEN x, GEN y)
    1481             : {
    1482     1009397 :   GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
    1483     1009397 :   return z;
    1484             : }
    1485             : GEN
    1486     1009397 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
    1487     1009397 : { return qfr3_red(qfr3_compraw(x,y), S); }
    1488             : 
    1489             : /* m > 0. Not stack-clean */
    1490             : static GEN
    1491           7 : qfr5_powraw(GEN x, long m)
    1492             : {
    1493           7 :   GEN y = NULL;
    1494          14 :   for (; m; m >>= 1)
    1495             :   {
    1496          14 :     if (m&1) y = y? qfr5_compraw(y,x): x;
    1497          14 :     if (m == 1) break;
    1498           7 :     x = qfr5_compraw(x,x);
    1499             :   }
    1500           7 :   return y;
    1501             : }
    1502             : 
    1503             : /* return x^n. Not stack-clean */
    1504             : GEN
    1505          21 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
    1506             : {
    1507          21 :   GEN y = NULL;
    1508          21 :   long i, m, s = signe(n);
    1509          21 :   if (!s) return qfr5_1(S, lg(gel(x,5)));
    1510          21 :   if (s < 0) x = qfb_inv(x);
    1511          42 :   for (i=lgefint(n)-1; i>1; i--)
    1512             :   {
    1513          21 :     m = n[i];
    1514          56 :     for (; m; m>>=1)
    1515             :     {
    1516          56 :       if (m&1) y = y? qfr5_comp(y,x,S): x;
    1517          56 :       if (m == 1 && i == 2) break;
    1518          35 :       x = qfr5_comp(x,x,S);
    1519             :     }
    1520             :   }
    1521          21 :   return y;
    1522             : }
    1523             : /* m > 0; return x^m. Not stack-clean */
    1524             : static GEN
    1525           0 : qfr3_powraw(GEN x, long m)
    1526             : {
    1527           0 :   GEN y = NULL;
    1528           0 :   for (; m; m>>=1)
    1529             :   {
    1530           0 :     if (m&1) y = y? qfr3_compraw(y,x): x;
    1531           0 :     if (m == 1) break;
    1532           0 :     x = qfr3_compraw(x,x);
    1533             :   }
    1534           0 :   return y;
    1535             : }
    1536             : /* return x^n. Not stack-clean */
    1537             : GEN
    1538        4557 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
    1539             : {
    1540        4557 :   GEN y = NULL;
    1541        4557 :   long i, m, s = signe(n);
    1542        4557 :   if (!s) return qfr3_1(S);
    1543        4557 :   if (s < 0) x = qfb_inv(x);
    1544        9130 :   for (i=lgefint(n)-1; i>1; i--)
    1545             :   {
    1546        4573 :     m = n[i];
    1547        5312 :     for (; m; m>>=1)
    1548             :     {
    1549        5292 :       if (m&1) y = y? qfr3_comp(y,x,S): x;
    1550        5292 :       if (m == 1 && i == 2) break;
    1551         739 :       x = qfr3_comp(x,x,S);
    1552             :     }
    1553             :   }
    1554        4557 :   return y;
    1555             : }
    1556             : 
    1557             : static GEN
    1558           7 : qfrinvraw(GEN x)
    1559             : {
    1560           7 :   if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
    1561           7 :  return qfbinv(x);
    1562             : }
    1563             : static GEN
    1564          14 : qfrpowraw(GEN x, long n)
    1565             : {
    1566          14 :   struct qfr_data S = { NULL, NULL, NULL };
    1567          14 :   pari_sp av = avma;
    1568          14 :   if (n==1) return gcopy(x);
    1569          14 :   if (n==-1) return qfrinvraw(x);
    1570           7 :   if (typ(x)==t_QFB)
    1571             :   {
    1572           0 :     GEN D = qfb_disc(x);
    1573           0 :     if (!n) return qfr_1(x);
    1574           0 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1575           0 :     x = qfr3_powraw(x, n);
    1576           0 :     x = qfr3_to_qfr(x, D);
    1577             :   }
    1578             :   else
    1579             :   {
    1580           7 :     GEN d0 = gel(x,2);
    1581           7 :     x = gel(x,1);
    1582           7 :     if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1583           7 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1584           7 :     x = qfr5_init(x, d0, &S);
    1585           7 :     if (labs(n) != 1) x = qfr5_powraw(x, n);
    1586           7 :     x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
    1587             :   }
    1588           7 :   return gc_GEN(av, x);
    1589             : }
    1590             : static GEN
    1591         112 : qfrpow(GEN x, GEN n)
    1592             : {
    1593         112 :   struct qfr_data S = { NULL, NULL, NULL };
    1594         112 :   long s = signe(n);
    1595         112 :   pari_sp av = avma;
    1596         112 :   if (typ(x)==t_QFB)
    1597             :   {
    1598          42 :     if (!s) return qfr_1(x);
    1599          28 :     if (s < 0) x = qfb_inv(x);
    1600          28 :     x = qfr3_init(x, &S);
    1601          28 :     x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
    1602          28 :     x = qfr3_to_qfr(x, S.D);
    1603             :   }
    1604             :   else
    1605             :   {
    1606          70 :     GEN d0 = gel(x,2);
    1607          70 :     x = gel(x,1);
    1608          70 :     if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1609          49 :     if (s < 0) x = qfb_inv(x);
    1610          49 :     x = qfr5_init(x, d0, &S);
    1611          49 :     x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
    1612          49 :     x = qfr5_to_qfr(x, S.D, mulri(d0,n));
    1613             :   }
    1614          77 :   return gc_GEN(av, x);
    1615             : }
    1616             : GEN
    1617          21 : qfbpowraw(GEN x, long n)
    1618             : {
    1619          21 :   GEN q = check_qfbext("qfbpowraw",x);
    1620          21 :   return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
    1621             : }
    1622             : /* same discriminant, no distance, no checks */
    1623             : GEN
    1624    14496436 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
    1625             : GEN
    1626     1429143 : qfbpow(GEN x, GEN n)
    1627             : {
    1628     1429143 :   GEN q = check_qfbext("qfbpow",x);
    1629     1429142 :   return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
    1630             : }
    1631             : GEN
    1632     1330605 : qfbpows(GEN x, long n)
    1633             : {
    1634     1330605 :   long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
    1635     1330605 :   affsi(n, N); return qfbpow(x, N);
    1636             : }
    1637             : 
    1638             : /* Prime forms attached to prime ideals of degree 1 */
    1639             : 
    1640             : /* assume x != 0 a t_INT, p > 0
    1641             :  * Return a t_QFB, but discriminant sign is not checked: can be used for
    1642             :  * real forms as well */
    1643             : GEN
    1644    12306119 : primeform_u(GEN x, ulong p)
    1645             : {
    1646    12306119 :   GEN c, y = cgetg(5, t_QFB);
    1647    12306054 :   pari_sp av = avma;
    1648             :   ulong b;
    1649             :   long s;
    1650             : 
    1651    12306054 :   s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
    1652             :   /* 2 or 3 mod 4 */
    1653    12306104 :   if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
    1654    12306312 :   if (p == 2) {
    1655     3971331 :     switch(s) {
    1656      589600 :       case 0: b = 0; break;
    1657     3030106 :       case 1: b = 1; break;
    1658      351626 :       case 4: b = 2; break;
    1659           0 :       default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1660           0 :                b = 0; /* -Wall */
    1661             :     }
    1662     3971332 :     c = shifti(subsi(s,x), -3);
    1663             :   } else {
    1664     8334981 :     b = Fl_sqrt(umodiu(x,p), p);
    1665     8336254 :     if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1666             :     /* mod(b) != mod2(x) ? */
    1667     8336986 :     if ((b ^ s) & 1) b = p - b;
    1668     8336986 :     c = diviuexact(shifti(subii(sqru(b), x), -2), p);
    1669             :   }
    1670    12291774 :   gel(y,3) = gc_INT(av, c);
    1671    12301367 :   gel(y,4) = icopy(x);
    1672    12304888 :   gel(y,2) = utoi(b);
    1673    12304963 :   gel(y,1) = utoipos(p); return y;
    1674             : }
    1675             : 
    1676             : /* special case: p = 1 return unit form */
    1677             : GEN
    1678      135480 : primeform(GEN x, GEN p)
    1679             : {
    1680      135480 :   const char *f = "primeform";
    1681             :   pari_sp av;
    1682      135480 :   long s, sx = signe(x), sp = signe(p);
    1683             :   GEN y, b, absp;
    1684             : 
    1685      135480 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
    1686      135480 :   if (typ(p) != t_INT) pari_err_TYPE(f,p);
    1687      135480 :   if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
    1688      135480 :   if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
    1689      135480 :   if (lgefint(p) == 3)
    1690             :   {
    1691      135466 :     ulong pp = p[2];
    1692      135466 :     if (pp == 1) {
    1693       17803 :       if (sx < 0) {
    1694             :         long r;
    1695       10971 :         if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1696       10971 :         r = mod4(x);
    1697       10971 :         if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
    1698       10971 :         return qfi_1_by_disc(x);
    1699             :       }
    1700        6832 :       y = qfr_1_by_disc(x);
    1701        6832 :       if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
    1702        6832 :       return y;
    1703             :     }
    1704      117663 :     y = primeform_u(x, pp);
    1705      117656 :     if (sx < 0) {
    1706       89957 :       if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1707       89957 :       return y;
    1708             :     }
    1709       27699 :     if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
    1710       27699 :     return gcopy( qfr3_to_qfr(y, x) );
    1711             :   }
    1712          14 :   s = mod8(x);
    1713          14 :   if (sx < 0)
    1714             :   {
    1715           7 :     if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1716           7 :     if (s) s = 8-s;
    1717             :   }
    1718          14 :   y = cgetg(5, t_QFB);
    1719             :   /* 2 or 3 mod 4 */
    1720          14 :   if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
    1721          14 :   absp = absi_shallow(p); av = avma;
    1722          14 :   b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
    1723          14 :   s &= 1; /* s = x mod 2 */
    1724             :   /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
    1725          14 :   if ((!signe(b) && s) || mod2(b) != s) b = gc_INT(av, subii(absp,b));
    1726             : 
    1727          14 :   av = avma;
    1728          14 :   gel(y,3) = gc_INT(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
    1729          14 :   gel(y,4) = icopy(x);
    1730          14 :   gel(y,2) = b;
    1731          14 :   gel(y,1) = icopy(p);
    1732          14 :   return y;
    1733             : }
    1734             : 
    1735             : static GEN
    1736     2620772 : normforms(GEN D, GEN fa)
    1737             : {
    1738             :   long i, j, k, lB, aN, sa;
    1739             :   GEN a, L, V, B, N, N2;
    1740     2620772 :   int D_odd = mpodd(D);
    1741     2620772 :   a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
    1742     2620772 :   sa = signe(a);
    1743     2620772 :   if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
    1744     1203972 :   V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
    1745     2551766 :            : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
    1746     2551766 :   if (!V) return NULL;
    1747      511966 :   N = gel(V,1); B = gel(V,2); lB = lg(B);
    1748      511966 :   N2 = shifti(N,1);
    1749      511966 :   aN = itou(diviiexact(a, N)); /* |a|/N */
    1750      511965 :   L = cgetg((lB-1)*aN+1, t_VEC);
    1751     2360562 :   for (k = 1, i = 1; i < lB; i++)
    1752             :   {
    1753     1848597 :     GEN b = shifti(gel(B,i), 1), c, C;
    1754     1848594 :     if (D_odd) b = addiu(b , 1);
    1755     1848594 :     c = diviiexact(shifti(subii(sqri(b), D), -2), a);
    1756     1848586 :     for (j = 0;; b = addii(b, N2))
    1757             :     {
    1758     2216660 :       gel(L, k++) = mkqfb(a, b, c, D);
    1759     2216670 :       if (++j == aN) break;
    1760      368072 :       C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
    1761      368074 :       c = sa > 0? addii(c, C): subii(c, C);
    1762             :     }
    1763             :   }
    1764      511965 :   return L;
    1765             : }
    1766             : 
    1767             : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
    1768             : static GEN
    1769      344320 : SL2_div_mul_e1(GEN N, GEN M)
    1770             : {
    1771      344320 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1772      344320 :   GEN A = mulii(gcoeff(N,1,1), d), B = mulii(gcoeff(N,1,2), b);
    1773      344306 :   GEN C = mulii(gcoeff(N,2,1), d), D = mulii(gcoeff(N,2,2), b);
    1774      344311 :   retmkvec2(subii(A,B), subii(C,D));
    1775             : }
    1776             : static GEN
    1777     1445672 : qfisolve_normform(GEN Q, GEN P)
    1778             : {
    1779     1445672 :   GEN a = gel(Q,1), N = gel(Q,2);
    1780     1445672 :   GEN M, b = qfi_redsl2_basecase(P, &M);
    1781     1445681 :   if (!qfb_equal(a,b)) return NULL;
    1782      102127 :   return SL2_div_mul_e1(N,M);
    1783             : }
    1784             : 
    1785             : /* Test equality modulo GL2 of two reduced forms */
    1786             : static int
    1787       61068 : GL2_qfb_equal(GEN a, GEN b)
    1788             : {
    1789       61068 :   return equalii(gel(a,1),gel(b,1))
    1790       11361 :    && absequalii(gel(a,2),gel(b,2))
    1791       72429 :    &&    equalii(gel(a,3),gel(b,3));
    1792             : }
    1793             : 
    1794             : /* Q(u,v) = p; if s < 0 return that solution; else the set of all solutions */
    1795             : static GEN
    1796       48048 : allsols(GEN Q, long s, GEN u, GEN v)
    1797             : {
    1798       48048 :   GEN w = mkvec2(u, v), b;
    1799       48051 :   if (signe(v) < 0) { u = negi(u); v = negi(v); } /* normalize for v >= 0 */
    1800       48051 :   w = mkvec2(u, v); if (s < 0) return w;
    1801       41403 :   if (!s) return mkvec(w);
    1802       38974 :   b = gel(Q,2); /* sum of the 2 solutions (if they exist) is -bv / a */
    1803       38974 :   if (signe(b))
    1804             :   { /* something to check */
    1805             :     GEN r, t;
    1806       13433 :     t = dvmdii(mulii(b, v), gel(Q,1), &r);
    1807       13433 :     if (signe(r)) return mkvec(w);
    1808        1820 :     u = addii(u, t);
    1809             :   }
    1810       27361 :   return mkvec2(w, mkvec2(negi(u), v));
    1811             : }
    1812             : static GEN
    1813      223082 : qfisolvep_all(GEN Q, GEN p, long all)
    1814             : {
    1815      223082 :   GEN R, U, V, M, N, x, q, D = qfb_disc(Q);
    1816      223078 :   long s = kronecker(D, p);
    1817             : 
    1818      223070 :   if (s < 0) return NULL;
    1819      126995 :   if (!all) s = -1; /* to indicate we want a single solution */
    1820             :   /* Solutions iff a class of maximal ideal above p is the class of Q;
    1821             :    * Two solutions iff (s > 0 and the class has order > 2), else one */
    1822      126995 :   if (!signe(gel(Q,2)))
    1823             :   { /* if principal form, use faster cornacchia */
    1824       43674 :     GEN a = gel(Q,1), c = gel(Q,3);
    1825       43674 :     if (equali1(a))
    1826             :     {
    1827       38199 :       if (!cornacchia(c, p, &M,&N)) return NULL;
    1828       33732 :       return allsols(Q, s, M, N);
    1829             :     }
    1830        5474 :     if (equali1(c))
    1831             :     {
    1832        5194 :       if (!cornacchia(a, p, &M,&N)) return NULL;
    1833         721 :       return allsols(Q, s, N, M);
    1834             :     }
    1835             :   }
    1836       83601 :   R = qfi_redsl2_basecase(Q, &U);
    1837       83601 :   if (equali1(gel(R,1)))
    1838             :   { /* principal form */
    1839       22533 :     if (!signe(gel(R,2)))
    1840             :     {
    1841        4396 :       if (!cornacchia(gel(R,3), p, &M,&N)) return NULL;
    1842         812 :       x = mkvec2(M,N);
    1843             :     }
    1844             :     else
    1845             :     { /* x^2 + xy + ((1-D)/4)y^2 = p <==> (2x + y)^2 - D y^2 = 4p */
    1846       18137 :       if (!cornacchia2(negi(D), p, &M, &N)) return NULL;
    1847        2331 :       x = subii(M,N); if (mpodd(x)) return NULL;
    1848        2331 :       x = mkvec2(shifti(x,-1), N);
    1849             :     }
    1850        3143 :     x = ZM_ZC_mul(U, x); x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1851        3143 :     return allsols(Q, s, gel(x,1), gel(x,2));
    1852             :   }
    1853       61068 :   q = qfi_redsl2_basecase(primeform(D, p), &V);
    1854       61068 :   if (!GL2_qfb_equal(R,q)) return NULL;
    1855       10451 :   if (signe(gel(R,2)) != signe(gel(q,2))) gcoeff(V,2,1) = negi(gcoeff(V,2,1));
    1856       10451 :   x = SL2_div_mul_e1(U,V); return allsols(Q, s, gel(x,1), gel(x,2));
    1857             : }
    1858             : GEN
    1859           0 : qfisolvep(GEN Q, GEN p)
    1860             : {
    1861           0 :   pari_sp av = avma;
    1862           0 :   GEN x = qfisolvep_all(Q, p, 0);
    1863           0 :   return x? gc_GEN(av, x): gc_const(av, gen_0);
    1864             : }
    1865             : 
    1866             : static GEN
    1867      770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
    1868             : {
    1869      770126 :   pari_sp av = avma, btop;
    1870      770126 :   GEN M = N, P = qfr_redsl2_basecase(Ps, rd), Q = P;
    1871             : 
    1872      770126 :   btop = avma;
    1873             :   for(;;)
    1874             :   {
    1875     5840681 :     if (qfb_equal(gel(M,1), gel(P,1)))
    1876      154084 :       return gc_upto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
    1877     5686597 :     if (qfb_equal(gel(N,1), gel(Q,1)))
    1878       77658 :       return gc_upto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
    1879     5608939 :     M = qfr_rhosl2(M, rd);
    1880     5608939 :     if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
    1881     5201735 :     Q = qfr_rhosl2(Q, rd);
    1882     5201735 :     if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
    1883     5070555 :     if (gc_needed(btop, 1)) (void)gc_all(btop, 2, &M, &Q);
    1884             :   }
    1885             : }
    1886             : 
    1887             : GEN
    1888           0 : qfrsolvep(GEN Q, GEN p)
    1889             : {
    1890           0 :   pari_sp av = avma;
    1891           0 :   GEN N, x, rd, d = qfb_disc(Q);
    1892             : 
    1893           0 :   if (kronecker(d, p) < 0) return gc_const(av, gen_0);
    1894           0 :   rd = sqrti(d);
    1895           0 :   N = qfr_redsl2(Q, rd);
    1896           0 :   x = qfrsolve_normform(N, primeform(d, p), rd);
    1897           0 :   return x? gc_upto(av, x): gc_const(av, gen_0);
    1898             : }
    1899             : 
    1900             : static GEN
    1901     1862979 : known_prime(GEN v)
    1902             : {
    1903     1862979 :   GEN p, e, fa = check_arith_all(v, "qfbsolve");
    1904     1862974 :   if (!fa) return BPSW_psp(v)? v: NULL;
    1905       42112 :   if (lg(gel(fa,1)) != 2) return NULL;
    1906       29386 :   p = gcoeff(fa,1,1);
    1907       29386 :   e = gcoeff(fa,1,2);
    1908       29386 :   return (equali1(e) && !is_pm1(p) && signe(p) > 0)? p: NULL;
    1909             : }
    1910             : static GEN
    1911     2215798 : qfsolve_normform(GEN Q, GEN f, GEN rd)
    1912     2215798 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
    1913             : static GEN
    1914     2843857 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
    1915             : {
    1916             :   GEN x, W, F, p;
    1917             :   long i, j, l;
    1918     2843857 :   if (!rd && (p = known_prime(fa))) return qfisolvep_all(Q, p, all);
    1919     2620765 :   F = normforms(qfb_disc(Q), fa);
    1920     2620771 :   if (!F) return NULL;
    1921      511965 :   if (!*Qr) *Qr = qfbredsl2(Q, rd);
    1922      511966 :   l = lg(F); W = all? cgetg(l, t_VEC): NULL;
    1923     2727250 :   for (j = i = 1; i < l; i++)
    1924     2215798 :     if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
    1925             :     {
    1926      333859 :       if (!all) return x;
    1927      333348 :       gel(W,j++) = x;
    1928             :     }
    1929      511452 :   if (j == 1) return NULL;
    1930      127453 :   setlg(W,j); return lexsort(W);
    1931             : }
    1932             : 
    1933             : static GEN
    1934     2838554 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
    1935             : static GEN
    1936     2828326 : qfbsolve_primitive(GEN Q, GEN fa, long all)
    1937             : {
    1938     2828326 :   GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
    1939     2828329 :   x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
    1940     2828329 :   if (!x) return cgetg(1, t_VEC);
    1941      174790 :   return x;
    1942             : }
    1943             : 
    1944             : /* f / g^2 */
    1945             : static GEN
    1946        5299 : famat_divsqr(GEN f, GEN g)
    1947        5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
    1948             : static GEN
    1949       10227 : qfbsolve_all(GEN Q, GEN n, long all)
    1950             : {
    1951       10227 :   GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
    1952       10227 :   GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
    1953       10227 :   long i, j, l = lg(D);
    1954       10227 :   W = all? cgetg(l, t_VEC): NULL;
    1955       25151 :   for (i = j = 1; i < l; i++)
    1956             :   {
    1957       15526 :     GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
    1958       15526 :     if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
    1959             :     {
    1960        1218 :       if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
    1961        1218 :       if (!all) return w;
    1962         616 :       gel(W,j++) = w;
    1963             :     }
    1964             :   }
    1965        9625 :   if (j == 1) return cgetg(1, t_VEC);
    1966         525 :   setlg(W,j); return lexsort(shallowconcat1(W));
    1967             : }
    1968             : 
    1969             : GEN
    1970     2838562 : qfbsolve(GEN Q, GEN n, long fl)
    1971             : {
    1972     2838562 :   pari_sp av = avma;
    1973     2838562 :   if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
    1974     2838562 :   if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
    1975     5666885 :   return gc_GEN(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
    1976     2828327 :                                   : qfbsolve_primitive(Q, n, fl & 1));
    1977             : }
    1978             : 
    1979             : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
    1980             :  * Assume d > 0 and p is prime */
    1981             : long
    1982       55273 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
    1983             : {
    1984       55273 :   pari_sp av = avma;
    1985             :   GEN b, c, r;
    1986             : 
    1987       55273 :   *px = *py = gen_0;
    1988       55273 :   b = subii(p, d);
    1989       55223 :   if (signe(b) < 0) return gc_long(av,0);
    1990       55013 :   if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
    1991       55006 :   b = Fp_sqrt(b, p); /* sqrt(-d) */
    1992       55063 :   if (!b) return gc_long(av,0);
    1993       51332 :   b = gmael(halfgcdii(p, b), 2, 2);
    1994       51356 :   c = dvmdii(subii(p, sqri(b)), d, &r);
    1995       51260 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    1996       35509 :   set_avma(av);
    1997       35505 :   *px = icopy(b);
    1998       35496 :   *py = icopy(c); return 1;
    1999             : }
    2000             : 
    2001             : static GEN
    2002     1420527 : lastqi(GEN Q)
    2003             : {
    2004     1420527 :   GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
    2005     1420521 :   if (!signe(q)) return gen_0;
    2006     1420332 :   if (!signe(s)) return p;
    2007     1413948 :   if (is_pm1(q)) return subiu(p,1);
    2008     1413949 :   return divii(p, absi_shallow(q));
    2009             : }
    2010             : 
    2011             : static long
    2012     1420533 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
    2013             : {
    2014             :   GEN M, Q, V, c, r, b2;
    2015     1420533 :   if (!signe(b)) { /* d = p,2p,3p,4p */
    2016          14 :     set_avma(av);
    2017          14 :     if (absequalii(d, px4)){ *py = gen_1; return 1; }
    2018          14 :     if (absequalii(d, p))  { *py = gen_2; return 1; }
    2019           0 :     return 0;
    2020             :   }
    2021     1420519 :   if (mod2(b) != mod2(d)) b = subii(p,b);
    2022     1420497 :   M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
    2023     1420528 :   b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
    2024     1420442 :   b2 = sqri(b);
    2025     1420434 :   if (cmpii(b2,px4) > 0)
    2026             :   {
    2027     1412215 :     b = gel(V,1); b2 = sqri(b);
    2028     1412225 :     if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
    2029             :   }
    2030     1420431 :   c = dvmdii(subii(px4, b2), d, &r);
    2031     1420426 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    2032     1376057 :   set_avma(av);
    2033     1376054 :   *px = icopy(b);
    2034     1376051 :   *py = icopy(c); return 1;
    2035             : }
    2036             : 
    2037             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
    2038             :  * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
    2039             : long
    2040     1381710 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
    2041             : {
    2042     1381710 :   pari_sp av = avma;
    2043     1381710 :   GEN b, p4 = shifti(p,2);
    2044             : 
    2045     1381680 :   *px = *py = gen_0;
    2046     1381680 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2047     1380849 :   if (absequaliu(p, 2))
    2048             :   {
    2049           7 :     set_avma(av);
    2050           7 :     switch (itou_or_0(d)) {
    2051           0 :       case 4: *px = gen_2; break;
    2052           0 :       case 7: *px = gen_1; break;
    2053           7 :       default: return 0;
    2054             :     }
    2055           0 :     *py = gen_1; return 1;
    2056             :   }
    2057     1380847 :   b = Fp_sqrt(negi(d), p);
    2058     1380934 :   if (!b) return gc_long(av,0);
    2059     1380850 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2060             : }
    2061             : 
    2062             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
    2063             : long
    2064       39683 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
    2065             : {
    2066       39683 :   pari_sp av = avma;
    2067       39683 :   GEN p4 = shifti(p,2);
    2068       39683 :   *px = *py = gen_0;
    2069       39683 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2070       39683 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2071             : }
    2072             : 
    2073             : GEN
    2074        7630 : qfbcornacchia(GEN d, GEN p)
    2075             : {
    2076        7630 :   pari_sp av = avma;
    2077             :   GEN x, y;
    2078        7630 :   if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
    2079        7630 :   if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
    2080        7630 :   if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
    2081         287 :     return gc_GEN(av, mkvec2(x, y));
    2082        7343 :   retgc_const(av, cgetg(1, t_VEC));
    2083             : }

Generated by: LCOV version 1.16