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 30371-bd27ad7662) Lines: 1108 1229 90.2 %
Date: 2025-07-06 09:22:54 Functions: 133 143 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      151918 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
      24             : {
      25             :   long r;
      26      151918 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
      27      151904 :   *s = signe(x);
      28      151904 :   if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
      29      151905 :   r = mod4(x); if (*s < 0 && r) r = 4 - r;
      30      151905 :   if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
      31      151891 :   if (pr) *pr = r;
      32      151891 : }
      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     1019309 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
      50             : {
      51     1019309 :   if (Dodd)
      52             :   {
      53      920945 :     pari_sp av = avma;
      54      920945 :     *b = gen_m1;
      55      920945 :     *c = gc_INT(av, shifti(subui(1,D), -2));
      56             :   }
      57             :   else
      58             :   {
      59       98364 :     *b = gen_0;
      60       98364 :     *c = shifti(D,-2); togglesign(*c);
      61             :   }
      62     1019309 : }
      63             : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
      64             : static GEN
      65      243348 : quadpoly_ii(GEN D, long Dmod4)
      66             : {
      67      243348 :   GEN b, c, y = cgetg(5,t_POL);
      68      243348 :   y[1] = evalsigne(1) | evalvarn(0);
      69      243348 :   quadpoly_bc(D, Dmod4, &b,&c);
      70      243348 :   gel(y,2) = c;
      71      243348 :   gel(y,3) = b;
      72      243348 :   gel(y,4) = gen_1; return y;
      73             : }
      74             : GEN
      75        2051 : quadpoly(GEN D)
      76             : {
      77             :   long s, Dmod4;
      78        2051 :   check_quaddisc(D, &s, &Dmod4, "quadpoly");
      79        2044 :   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        1036 : quadpoly0(GEN x, long v)
      86             : {
      87        1036 :   GEN T = quadpoly(x);
      88        1029 :   if (v > 0) setvarn(T, v);
      89        1029 :   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     1959730 : check_qfbext(const char *fun, GEN x)
     113             : {
     114     1959730 :   long t = typ(x);
     115     1959730 :   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    23783634 : qfb_equal(GEN x, GEN y)
     131             : {
     132    23783634 :   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      933620 : qfb_inv(GEN x) {
     140      933620 :   GEN z = shallowcopy(x);
     141      933620 :   gel(z,2) = negi(gel(z,2));
     142      933619 :   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    16909031 : dvmdii_round(GEN b, GEN a, GEN *r)
     189             : {
     190    16909031 :   GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
     191    16908635 :   if (signe(b) >= 0) {
     192     9290420 :     if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
     193             :   } else { /* r <= 0 */
     194     7618215 :     if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
     195             :   }
     196    16908279 :   return q;
     197             : }
     198             : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
     199             : static long
     200   114712099 : dvmdsu_round(long b, ulong a, long *r)
     201             : {
     202   114712099 :   ulong a2 = a << 1, q, ub, ur;
     203   114712099 :   if (b >= 0) {
     204    73308620 :     ub = b;
     205    73308620 :     q = ub / a2;
     206    73308620 :     ur = ub % a2;
     207    73308620 :     if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
     208    26744957 :     else *r = (long)ur;
     209    73308620 :     return (long)q;
     210             :   } else { /* r <= 0 */
     211    41403479 :     ub = (ulong)-b; /* |b| */
     212    41403479 :     q = ub / a2;
     213    41403479 :     ur = ub % a2;
     214    41403479 :     if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
     215    22030149 :     else *r = -(long)ur;
     216    41403479 :     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   114708315 : sREDB(ulong a, long *b, ulong *c)
     231             : {
     232             :   long r, q;
     233             :   ulong uz;
     234   124561586 :   if (a > LONG_MAX) return; /* b already reduced */
     235   114708315 :   q = dvmdsu_round(*b, a, &r);
     236   114727082 :   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   104873811 :   if (*b < 0)
     240             :   { /* uz = -z >= 0, q < 0 */
     241    36431786 :     if (r >= 0) /* different signs=>no overflow, exact division */
     242    19429412 :       uz = (ulong)-((*b + r)>>1);
     243             :     else
     244             :     {
     245    17002374 :       ulong ub = (ulong)-*b, ur = (ulong)-r;
     246    17002374 :       uz = (ub + ur) >> 1;
     247             :     }
     248    36431786 :     *c -= (-q) * uz; /* c -= qz */
     249             :   }
     250             :   else
     251             :   { /* uz = z >= 0, q > 0 */
     252    68442025 :     if (r <= 0)
     253    46624929 :       uz = (*b + r)>>1;
     254             :     else
     255             :     {
     256    21817096 :       ulong ub = (ulong)*b, ur = (ulong)r;
     257    21817096 :       uz = ((ub + ur) >> 1);
     258             :     }
     259    68442025 :     *c -= q * uz; /* c -= qz */
     260             :   }
     261   104873811 :   *b = r;
     262             : }
     263             : static void
     264    14125047 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
     265             : { /* REDB(a,b,c) */
     266    14125047 :   GEN r, q = dvmdii_round(*b, a, &r);
     267    14124273 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     268    14124186 :   *b = r;
     269    14124186 :   *u2 = subii(*u2, mulii(q, u1));
     270    14124282 : }
     271             : 
     272             : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
     273             : static GEN
     274     6773428 : qfi_redsl2_basecase(GEN q, GEN *U)
     275             : {
     276     6773428 :   pari_sp av = avma;
     277             :   GEN z, u1,u2,v1,v2,Q;
     278     6773428 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     279             :   long cmp;
     280     6773428 :   u1 = gen_1; u2 = gen_0;
     281     6773428 :   cmp = abscmpii(a, b);
     282     6773434 :   if (cmp < 0)
     283     2195253 :     REDBU(a,&b,&c, u1,&u2);
     284     4578181 :   else if (cmp == 0 && signe(b) < 0)
     285             :   { /* b = -a */
     286       11784 :     b = negi(b);
     287       11784 :     u2 = gen_1;
     288             :   }
     289             :   for(;;)
     290             :   {
     291    18702625 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     292    11928832 :     swap(a,c); b = negi(b);
     293    11929755 :     z = u1; u1 = u2; u2 = negi(z);
     294    11929783 :     REDBU(a,&b,&c, u1,&u2);
     295    11929206 :     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     6773384 :   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     6773384 :   z = shifti(subii(b, gel(q,2)), -1);
     308     6773379 :   v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
     309     6773397 :   z = subii(z, b);
     310     6773383 :   v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
     311     6773370 :   *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
     312     6773448 :   Q = mkqfb(a,b,c,gel(q,4));
     313     6773444 :   return gc_all(av, 2, &Q, U);
     314             : }
     315             : 
     316             : static GEN
     317     1061596 : setq_b0(ulong a, ulong c, GEN D)
     318     1061596 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
     319             : /* assume |sb| = 1 */
     320             : static GEN
     321    99178319 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
     322    99178319 : { 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      947373 : qfi_red_1_b0(ulong a, ulong c, GEN D)
     326      947373 : { 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   100344566 : 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      141086 :   { /* at most twice */
     336   100344566 :     long lb = lgefint(b); /* <= 3 after first loop */
     337   100344566 :     if (lb == 2) return qfi_red_1_b0(a[2],c[2], D);
     338    99397193 :     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    99258168 :   set_avma(av);
     352    99257381 :   ua = a[2];
     353    99257381 :   ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
     354    99257381 :   uc = c[2];
     355    99257381 :   if (ua < ub)
     356    37333064 :     sREDB(ua, &sb, &uc);
     357    61924317 :   else if (ua == ub && sb < 0) sb = (long)ub;
     358   176695093 :   while(ua > uc)
     359             :   {
     360    77405109 :     lswap(ua,uc); sb = -sb;
     361    77405109 :     sREDB(ua, &sb, &uc);
     362             :   }
     363    99289984 :   if (!sb) return setq_b0(ua, uc, D);
     364             :   else
     365             :   {
     366    99175761 :     long s = 1;
     367    99175761 :     if (sb < 0)
     368             :     {
     369    37249751 :       sb = -sb;
     370    37249751 :       if (ua != uc) s = -1;
     371             :     }
     372    99175761 :     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   105969619 : qfb_maxexpi(GEN Q)
     743   105969619 : { 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   105661344 : qfi_red_fast(GEN Q)
     748             : {
     749   105661344 :   const long QFBRED_LIMIT = 9000;
     750   105661344 :   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             : GEN
     761       61913 : qfb3_SL2_apply(GEN q, GEN M)
     762             : {
     763       61913 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     764       61913 :   GEN x = gcoeff(M,1,1), y = gcoeff(M,2,1);
     765       61913 :   GEN z = gcoeff(M,1,2), t = gcoeff(M,2,2);
     766       61913 :   GEN by = mulii(b,y), bt = mulii(b,t), bz  = mulii(b,z);
     767       61913 :   GEN a2 = shifti(a,1), c2 = shifti(c,1);
     768             : 
     769       61913 :   GEN A1 = mulii(x, addii(mulii(a,x), by));
     770       61913 :   GEN A2 = mulii(c, sqri(y));
     771       61913 :   GEN B1 = mulii(x, addii(mulii(a2,z), bt));
     772       61913 :   GEN B2 = mulii(y, addii(mulii(c2,t), bz));
     773       61913 :   GEN C1 = mulii(z, addii(mulii(a,z), bt));
     774       61913 :   GEN C2 = mulii(c, sqri(t));
     775       61913 :   retmkvec3(addii(A1,A2), addii(B1,B2), addii(C1, C2));
     776             : }
     777             : 
     778             : static GEN
     779      124020 : pqfbred_rec(GEN Q, long m, GEN *pt_U)
     780             : {
     781      124020 :   pari_sp av = avma;
     782      124020 :   GEN U, Q0, Q1, QR, d = qfb_disc(Q);
     783      124020 :   long h, n = qfb_maxexpi(Q) - m;
     784      124020 :   int going_to_r8 = 0;
     785             : 
     786      124020 :   if (n < 170) return pqfbred_basecase(Q, m, pt_U);
     787       62528 :   if (qfb_minexpi(Q) <= m + 2) { U = matid(2); QR = Q; }
     788             :   else
     789             :   {
     790             :     long p, mm;
     791       62528 :     if (m <= n) { mm = m; p = 0; Q1 = Q; }
     792             :     else
     793             :     {
     794       61878 :       mm = n; p = m + 1 - n;
     795       61878 :       Q0 = mkvec3(remi2n(gel(Q,1),p), remi2n(gel(Q,2),p), remi2n(gel(Q,3),p));
     796       61878 :       Q1 = qfb3(shifti(gel(Q,1),-p), shifti(gel(Q,2),-p), shifti(gel(Q,3),-p));
     797             :     }
     798       62528 :     h = mm + (n>>1);
     799       62528 :     if (qfb_minexpi(Q1) <= h) { U = matid(2); QR = Q1; }
     800             :     else
     801       62321 :       QR = pqfbred_rec(Q1, h, &U);
     802      184249 :     while (qfb_maxexpi(QR) > h)
     803             :     {
     804      122868 :       if (is_minimal(QR, mm)) { going_to_r8 = 1; break; }
     805      121721 :       QR = pqfbred_1(QR, mm, U);
     806             :     }
     807       62528 :     if (!going_to_r8)
     808             :     {
     809             :       GEN V;
     810       61381 :       QR = pqfbred_rec(QR, mm, &V);
     811       61381 :       U = ZM2_mul(U,V);
     812             :     }
     813       62528 :     if (p > 0)
     814             :     {
     815       61878 :       GEN Q0U = qfb3_SL2_apply(Q0,U);
     816      123756 :       QR = mkqfb(addii(shifti(gel(QR,1), p), gel(Q0U,1)),
     817       61878 :                  addii(shifti(gel(QR,2), p), gel(Q0U,2)),
     818       61878 :                  addii(shifti(gel(QR,3), p), gel(Q0U,3)), d);
     819             :     }
     820             :   }
     821       62528 :   QR = pqfbred_iter_1(QR, m, U);
     822       62528 :   *pt_U = U; return gc_all(av, 2, &QR, pt_U);
     823             : }
     824             : 
     825             : static GEN
     826      209575 : qfr_redsl2(GEN Q, GEN isqrtD)
     827             : {
     828      209575 :   pari_sp av = avma;
     829      209575 :   if (!qfi_red_fast(Q))
     830      209575 :     return qfr_redsl2_basecase(Q, isqrtD);
     831             :   else
     832             :   {
     833           0 :     GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     834           0 :     GEN Qf, Qr, W, U, t = NULL;
     835           0 :     long sa = signe(a), sb;
     836           0 :     if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
     837           0 :     if (signe(c) < 0)
     838             :     {
     839             :       GEN at;
     840           0 :       t  = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
     841           0 :       at = mulii(a,t);
     842           0 :       c = addii(subii(c, mulii(b, t)), mulii(at, t));
     843           0 :       b = subii(b, shifti(at,1));
     844             :     }
     845           0 :     sb = signe(b);
     846           0 :     Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
     847           0 :     if (sa < 0)
     848           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     849           0 :     if (sb < 0)
     850             :     {
     851           0 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     852           0 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     853             :     }
     854           0 :     if (t)
     855             :     {
     856           0 :       gcoeff(U,1,1) = subii( gcoeff(U,1,1), mulii(gcoeff(U,2,1), t));
     857           0 :       gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,2,2), t));
     858             :     }
     859           0 :     W = qfr_redsl2_basecase(Qr, isqrtD);
     860           0 :     Qf = gel(W,1);
     861           0 :     U = ZM2_mul(U,gel(W,2));
     862           0 :     return gc_GEN(av, mkvec2(Qf,U));
     863             :   }
     864             : }
     865             : 
     866             : static GEN
     867     5183080 : qfi_redsl2(GEN Q)
     868             : {
     869     5183080 :   pari_sp av = avma;
     870             :   GEN Qt, U;
     871     5183080 :   if (!qfi_red_fast(Q))
     872     5182773 :     Qt = qfi_redsl2_basecase(Q, &U);
     873             :   else
     874             :   {
     875         311 :     long sb = signe(gel(Q,2));
     876             :     GEN W;
     877         311 :     if (sb < 0) Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
     878         311 :     Q = pqfbred_rec(Q, 0, &U);
     879         311 :     Qt = qfi_redsl2_basecase(Q, &W);
     880         311 :     U = ZM2_mul(U,W);
     881         311 :     if (sb < 0)
     882             :     {
     883         173 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     884         173 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     885             :     }
     886             :   }
     887     5183108 :   return gc_GEN(av, mkvec2(Qt,U));
     888             : }
     889             : 
     890             : GEN
     891     4872765 : redimagsl2(GEN Q, GEN *U)
     892             : {
     893     4872765 :   GEN q = qfi_redsl2(Q);
     894     4872804 :   *U = gel(q,2); return gel(q,1);
     895             : }
     896             : 
     897             : GEN
     898      519893 : qfbredsl2(GEN q, GEN isD)
     899             : {
     900             :   pari_sp av;
     901      519893 :   if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
     902      519893 :   if (qfb_is_qfi(q))
     903             :   {
     904      310311 :     if (isD) pari_err_TYPE("qfbredsl2", isD);
     905      310311 :     return qfi_redsl2(q);
     906             :   }
     907      209582 :   av = avma;
     908      209582 :   if (!isD) isD = sqrti(qfb_disc(q));
     909      208068 :   else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
     910      209575 :   return gc_upto(av, qfr_redsl2(q, isD));
     911             : }
     912             : 
     913             : /* not gc-clean */
     914             : static GEN
     915         427 : qfr_red_i(GEN Q, long flag, GEN isqrtD, GEN sqrtD)
     916             : {
     917         427 :   if (typ(Q) != t_QFB || !qfi_red_fast(Q))
     918         427 :     return qfr_red_basecase_i(Q, flag, isqrtD, sqrtD);
     919             :   else
     920             :   {
     921           0 :     GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     922           0 :     GEN Qr, W, U, t = NULL;
     923           0 :     long sa = signe(a);
     924           0 :     if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
     925           0 :     if (signe(c) < 0)
     926             :     {
     927             :       GEN at;
     928           0 :       t  = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
     929           0 :       at = mulii(a,t);
     930           0 :       c = addii(subii(c, mulii(b, t)), mulii(at, t));
     931           0 :       b = subii(b, shifti(at,1));
     932             :     }
     933           0 :     Qr = pqfbred_rec(mkqfb(a, absi_shallow(b), c, d), 0, &U);
     934           0 :     if (sa < 0)
     935           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     936           0 :     W = qfr_red_basecase_i(Qr, flag, isqrtD, sqrtD);
     937           0 :     return gel(W,1);
     938             :   }
     939             : }
     940             : 
     941             : static GEN
     942          63 : qfr_red_av(pari_sp av, GEN x)
     943          63 : { return gc_GEN(av, qfr_red_i(x,0,NULL,NULL)); }
     944             : GEN
     945           0 : qfr_red(GEN x) { return qfr_red_av(avma, x); }
     946             : 
     947             : static GEN
     948   100385388 : qfi_red_basecase_av(pari_sp av, GEN q)
     949             : {
     950   100385388 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
     951   100385388 :   long cmp, lc = lgefint(c);
     952             : 
     953   100385388 :   if (lgefint(a) == 3 && lc == 3) return qfi_red_1(av, a, b, c, D);
     954      914299 :   cmp = abscmpii(a, b);
     955      919092 :   if (cmp < 0)
     956      436231 :     REDB(a,&b,&c);
     957      482861 :   else if (cmp == 0 && signe(b) < 0)
     958          27 :     b = negi(b);
     959             :   for(;;)
     960             :   {
     961     3126345 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     962     2932679 :     lc = lgefint(a); /* lg(future c): we swap a & c next */
     963     2932679 :     if (lc == 3) return qfi_red_1(av, a, b, c, D);
     964     2207253 :     swap(a,c); b = negi(b); /* apply rho */
     965     2207253 :     REDB(a,&b,&c);
     966     2207253 :     if (gc_needed(av, 2))
     967             :     {
     968           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfi_red, lc = %ld", lc);
     969           0 :       (void)gc_all(av, 3, &a,&b,&c);
     970             :     }
     971             :   }
     972      193666 :   if (cmp == 0 && signe(b) < 0) b = negi(b);
     973      193666 :   return gc_GEN(av, mkqfb(a, b, c, D));
     974             : }
     975             : static GEN
     976   100270796 : qfi_red_av(pari_sp av, GEN Q)
     977             : {
     978   100270796 :   if (qfi_red_fast(Q))
     979             :   {
     980             :     GEN U;
     981           7 :     if (signe(gel(Q,2)) < 0)
     982           0 :       Q = mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4));
     983           7 :     Q = pqfbred_rec(Q, 0, &U);
     984             :   }
     985   100383019 :   return qfi_red_basecase_av(av, Q);
     986             : }
     987             : 
     988             : GEN
     989    23628873 : qfi_red(GEN q) { return qfi_red_av(avma, q); }
     990             : 
     991             : GEN
     992       55042 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
     993             : {
     994             :   pari_sp av;
     995       55042 :   GEN q = check_qfbext("qfbred",x);
     996       55042 :   if (qfb_is_qfi(q)) return (flag & qf_STEP)? qfi_rho(x): qfi_red(x);
     997         364 :   if (typ(x)==t_QFB) flag |= qf_NOD;
     998          49 :   else               flag &= ~qf_NOD;
     999         364 :   av = avma;
    1000         364 :   return gc_GEN(av, qfr_red_i(x,flag,isqrtD,sqrtD));
    1001             : }
    1002             : /* t_QFB */
    1003             : GEN
    1004    16026978 : qfbred_i(GEN x) { return qfb_is_qfi(x)? qfi_red(x): qfr_red(x); }
    1005             : GEN
    1006       53229 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
    1007             : /***********************************************************************/
    1008             : /**                                                                   **/
    1009             : /**                         Composition                               **/
    1010             : /**                                                                   **/
    1011             : /***********************************************************************/
    1012             : 
    1013             : static void
    1014    35228436 : qfb_sqr(GEN z, GEN x)
    1015             : {
    1016             :   GEN c, d1, x2, v1, v2, c3, m, p1, r;
    1017             : 
    1018    35228436 :   d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
    1019    35228353 :   c = gel(x,3);
    1020    35228353 :   m = mulii(c,x2);
    1021    35228133 :   if (equali1(d1))
    1022    26193180 :     v1 = v2 = gel(x,1);
    1023             :   else
    1024             :   {
    1025     9035001 :     v1 = diviiexact(gel(x,1),d1);
    1026     9035002 :     v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
    1027     9035002 :     c = mulii(c, d1);
    1028             :   }
    1029    35228183 :   togglesign(m);
    1030    35228096 :   r = modii(m,v2);
    1031    35227963 :   p1 = mulii(r, v1);
    1032    35228099 :   c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
    1033    35228085 :   gel(z,1) = mulii(v1,v2);
    1034    35228078 :   gel(z,2) = addii(gel(x,2), shifti(p1,1));
    1035    35228117 :   gel(z,3) = diviiexact(c3,v2);
    1036    35228074 : }
    1037             : /* z <- x * y */
    1038             : static void
    1039    78316575 : qfb_comp(GEN z, GEN x, GEN y)
    1040             : {
    1041             :   GEN n, c, d, y1, v1, v2, c3, m, p1, r;
    1042             : 
    1043    78316575 :   if (x == y) { qfb_sqr(z,x); return; }
    1044    43916078 :   n = shifti(subii(gel(y,2),gel(x,2)), -1);
    1045    43774989 :   v1 = gel(x,1);
    1046    43774989 :   v2 = gel(y,1);
    1047    43774989 :   c  = gel(y,3);
    1048    43774989 :   d = bezout(v2,v1,&y1,NULL);
    1049    43842264 :   if (equali1(d))
    1050    23154661 :     m = mulii(y1,n);
    1051             :   else
    1052             :   {
    1053    20723877 :     GEN s = subii(gel(y,2), n);
    1054    20723143 :     GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
    1055    20725311 :     if (!equali1(d1))
    1056             :     {
    1057    11456845 :       v1 = diviiexact(v1,d1);
    1058    11455763 :       v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
    1059    11455419 :       v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
    1060    11455984 :       c = mulii(c, d1);
    1061             :     }
    1062    20724110 :     m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
    1063             :   }
    1064    43828209 :   togglesign(m);
    1065    43765007 :   r = modii(m, v1);
    1066    43710771 :   p1 = mulii(r, v2);
    1067    43767453 :   c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
    1068    43759096 :   gel(z,1) = mulii(v1,v2);
    1069    43767073 :   gel(z,2) = addii(gel(y,2), shifti(p1,1));
    1070    43769605 :   gel(z,3) = diviiexact(c3,v1);
    1071             : }
    1072             : 
    1073             : /* not meant to be efficient */
    1074             : static GEN
    1075          84 : qfb_comp_gen(GEN x, GEN y)
    1076             : {
    1077          84 :   GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
    1078          84 :   GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
    1079          84 :   GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
    1080          84 :   GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
    1081             : 
    1082          84 :   if (!is_pm1(cx))
    1083             :   {
    1084          14 :     a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
    1085          14 :     c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
    1086             :   }
    1087          84 :   if (!is_pm1(cy))
    1088             :   {
    1089          28 :     a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
    1090          28 :     b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
    1091             :   }
    1092          84 :   D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
    1093         133 :   if (!Z_issquareall(diviiexact(d1, D), &n1) ||
    1094          84 :       !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
    1095          49 :   A = mulii(a1, n2);
    1096          49 :   B = mulii(a2, n1);
    1097          49 :   C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
    1098          49 :   U = ZV_extgcd(mkvec3(A, B, C));
    1099          49 :   m = gel(U,1); U = gmael(U,2,3);
    1100          49 :   A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
    1101          49 :   B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
    1102          49 :   C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
    1103          49 :   C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
    1104          49 :   B = addii(A, addii(B, C));
    1105          49 :   m2 = sqri(m);
    1106          49 :   A = diviiexact(mulii(a1, a2), m2);
    1107          49 :   C = diviiexact(shifti(subii(sqri(B),D), -2), A);
    1108          49 :   cx = mulii(cx, cy);
    1109          49 :   if (!is_pm1(cx))
    1110             :   {
    1111          14 :     A = mulii(A, cx); B = mulii(B, cx);
    1112          14 :     C = mulii(C, cx); D = mulii(D, sqri(cx));
    1113             :   }
    1114          49 :   return mkqfb(A, B, C, D);
    1115             : }
    1116             : 
    1117             : static GEN
    1118    75574581 : qficomp0(GEN x, GEN y, int raw)
    1119             : {
    1120    75574581 :   pari_sp av = avma;
    1121    75574581 :   GEN z = cgetg(5,t_QFB);
    1122    75572850 :   gel(z,4) = gel(x,4);
    1123    75572850 :   qfb_comp(z, x,y);
    1124    75402009 :   if (raw) return gc_GEN(av,z);
    1125    75400231 :   return qfi_red_av(av, z);
    1126             : }
    1127             : static GEN
    1128         441 : qfrcomp0(GEN x, GEN y, int raw)
    1129             : {
    1130         441 :   pari_sp av = avma;
    1131         441 :   GEN dx = NULL, dy = NULL;
    1132         441 :   GEN z = cgetg(5,t_QFB);
    1133         441 :   if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1134         441 :   if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
    1135         441 :   gel(z,4) = gel(x,4);
    1136         441 :   qfb_comp(z, x,y);
    1137         441 :   if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
    1138         441 :   if (raw) return gc_GEN(av, z);
    1139          28 :   return qfr_red_av(av, z);
    1140             : }
    1141             : /* same discriminant, no distance, no checks */
    1142             : GEN
    1143    27847633 : qfbcomp_i(GEN x, GEN y)
    1144    27847633 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
    1145             : GEN
    1146      136887 : qfbcomp(GEN x, GEN y)
    1147             : {
    1148      136887 :   GEN qx = check_qfbext("qfbcomp", x);
    1149      136887 :   GEN qy = check_qfbext("qfbcomp", y);
    1150      136887 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1151             :   {
    1152          63 :     pari_sp av = avma;
    1153          63 :     GEN z = qfb_comp_gen(qx, qy);
    1154          63 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1155           7 :       pari_err_IMPL("Shanks's distance in general composition");
    1156          56 :     if (!z) pari_err_OP("*",x,y);
    1157          21 :     return gc_upto(av, qfbred(z));
    1158             :   }
    1159      136824 :   return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
    1160             : }
    1161             : /* same discriminant, no distance, no checks */
    1162             : GEN
    1163           0 : qfbcompraw_i(GEN x, GEN y)
    1164           0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
    1165             : GEN
    1166        2198 : qfbcompraw(GEN x, GEN y)
    1167             : {
    1168        2198 :   GEN qx = check_qfbext("qfbcompraw", x);
    1169        2198 :   GEN qy = check_qfbext("qfbcompraw", y);
    1170        2198 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1171             :   {
    1172          21 :     pari_sp av = avma;
    1173          21 :     GEN z = qfb_comp_gen(qx, qy);
    1174          21 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1175           0 :       pari_err_IMPL("Shanks's distance in general composition");
    1176          21 :     if (!z) pari_err_OP("qfbcompraw",x,y);
    1177          21 :     return gc_GEN(av, z);
    1178             :   }
    1179        2177 :   if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
    1180        2177 :   return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
    1181             : }
    1182             : 
    1183             : static GEN
    1184      827908 : qfisqr0(GEN x, long raw)
    1185             : {
    1186      827908 :   pari_sp av = avma;
    1187      827908 :   GEN z = cgetg(5,t_QFB);
    1188      827908 :   gel(z,4) = gel(x,4);
    1189      827908 :   qfb_sqr(z,x);
    1190      827908 :   if (raw) return gc_GEN(av,z);
    1191      827908 :   return qfi_red_av(av, z);
    1192             : }
    1193             : static GEN
    1194          35 : qfrsqr0(GEN x, long raw)
    1195             : {
    1196          35 :   pari_sp av = avma;
    1197          35 :   GEN dx = NULL, z = cgetg(5,t_QFB);
    1198          35 :   if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1199          35 :   gel(z,4) = gel(x,4); qfb_sqr(z,x);
    1200          35 :   if (dx) z = mkvec2(z, shiftr(dx,1));
    1201          35 :   if (raw) return gc_GEN(av, z);
    1202          35 :   return qfr_red_av(av, z);
    1203             : }
    1204             : /* same discriminant, no distance, no checks */
    1205             : GEN
    1206      697020 : qfbsqr_i(GEN x)
    1207      697020 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
    1208             : GEN
    1209      130923 : qfbsqr(GEN x)
    1210             : {
    1211      130923 :   GEN qx = check_qfbext("qfbsqr", x);
    1212      130923 :   return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
    1213             : }
    1214             : 
    1215             : static GEN
    1216        6867 : qfr_1_by_disc(GEN D)
    1217             : {
    1218             :   GEN y, r, s;
    1219        6867 :   check_quaddisc_real(D, NULL, "qfr_1_by_disc");
    1220        6867 :   y = cgetg(5,t_QFB);
    1221        6867 :   s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
    1222        6867 :   if (mpodd(r))
    1223             :   {
    1224        3535 :     s = subiu(s,1);
    1225        3535 :     r = subii(r, addiu(shifti(s, 1), 1));
    1226        3535 :     r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
    1227             :   }
    1228             :   else
    1229        3332 :   { r = shifti(r, -2); set_avma((pari_sp)s); }
    1230        6867 :   gel(y,1) = gen_1;
    1231        6867 :   gel(y,2) = s;
    1232        6867 :   gel(y,3) = icopy(r);
    1233        6867 :   gel(y,4) = icopy(D); return y;
    1234             : }
    1235             : 
    1236             : static GEN
    1237          35 : qfr_disc(GEN x)
    1238          35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
    1239             : 
    1240             : static GEN
    1241          35 : qfr_1(GEN x)
    1242          35 : { return qfr_1_by_disc(qfr_disc(x)); }
    1243             : 
    1244             : static void
    1245           0 : qfr_1_fill(GEN y, struct qfr_data *S)
    1246             : {
    1247           0 :   pari_sp av = avma;
    1248           0 :   GEN y2 = S->isqrtD;
    1249           0 :   gel(y,1) = gen_1;
    1250           0 :   if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
    1251           0 :   gel(y,2) = y2; av = avma;
    1252           0 :   gel(y,3) = gc_INT(av, shifti(subii(sqri(y2), S->D),-2));
    1253           0 : }
    1254             : static GEN
    1255           0 : qfr5_1(struct qfr_data *S, long prec)
    1256             : {
    1257           0 :   GEN y = cgetg(6, t_VEC);
    1258           0 :   qfr_1_fill(y, S);
    1259           0 :   gel(y,4) = gen_0;
    1260           0 :   gel(y,5) = real_1(prec); return y;
    1261             : }
    1262             : static GEN
    1263           0 : qfr3_1(struct qfr_data *S)
    1264             : {
    1265           0 :   GEN y = cgetg(4, t_VEC);
    1266           0 :   qfr_1_fill(y, S); return y;
    1267             : }
    1268             : 
    1269             : /* Assume D < 0 is the discriminant of a t_QFB */
    1270             : static GEN
    1271      775961 : qfi_1_by_disc(GEN D)
    1272             : {
    1273      775961 :   GEN b,c, y = cgetg(5,t_QFB);
    1274      775961 :   quadpoly_bc(D, mod2(D), &b,&c);
    1275      775961 :   if (b == gen_m1) b = gen_1;
    1276      775961 :   gel(y,1) = gen_1;
    1277      775961 :   gel(y,2) = b;
    1278      775961 :   gel(y,3) = c;
    1279      775961 :   gel(y,4) = icopy(D); return y;
    1280             : }
    1281             : static GEN
    1282      763695 : qfi_1(GEN x)
    1283             : {
    1284      763695 :   if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
    1285      763695 :   return qfi_1_by_disc(qfb_disc(x));
    1286             : }
    1287             : 
    1288             : GEN
    1289           0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
    1290             : 
    1291             : static GEN
    1292    14319672 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
    1293             : static GEN
    1294    33271548 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
    1295             : static GEN
    1296           7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
    1297             : static GEN
    1298           7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
    1299             : 
    1300             : static GEN
    1301           7 : qfipowraw(GEN x, long n)
    1302             : {
    1303           7 :   pari_sp av = avma;
    1304             :   GEN y;
    1305           7 :   if (!n) return qfi_1(x);
    1306           7 :   if (n== 1) return gcopy(x);
    1307           7 :   if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
    1308           7 :   if (n < 0) x = qfb_inv(x);
    1309           7 :   y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
    1310           7 :   return gc_GEN(av,y);
    1311             : }
    1312             : 
    1313             : static GEN
    1314    16790674 : qfipow(GEN x, GEN n)
    1315             : {
    1316    16790674 :   pari_sp av = avma;
    1317             :   GEN y;
    1318    16790674 :   long s = signe(n);
    1319    16790674 :   if (!s) return qfi_1(x);
    1320    16026979 :   if (s < 0) x = qfb_inv(x);
    1321    16026978 :   y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
    1322    16026985 :   return gc_GEN(av,y);
    1323             : }
    1324             : 
    1325             : static long
    1326      412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
    1327             : {
    1328             :   long z;
    1329      412328 :   *v = gen_0; *v2 = gen_1;
    1330     4351417 :   for (z=0; abscmpii(*v3,L) > 0; z++)
    1331             :   {
    1332     3939089 :     GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
    1333     3939089 :     *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
    1334             :   }
    1335      412328 :   return z;
    1336             : }
    1337             : 
    1338             : /* composition: Shanks' NUCOMP & NUDUPL */
    1339             : /* L = floor((|d|/4)^(1/4)) */
    1340             : GEN
    1341      400722 : nucomp(GEN x, GEN y, GEN L)
    1342             : {
    1343      400722 :   pari_sp av = avma;
    1344             :   long z;
    1345             :   GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
    1346             : 
    1347      400722 :   if (x==y) return nudupl(x,L);
    1348      400680 :   if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
    1349      400680 :   if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
    1350             : 
    1351      400680 :   if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
    1352      400680 :   s = shifti(addii(gel(x,2),gel(y,2)), -1);
    1353      400680 :   n = subii(gel(y,2), s);
    1354      400680 :   a1 = gel(x,1);
    1355      400680 :   a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
    1356      400680 :   if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
    1357      163576 :   else if (dvdii(s,d)) /* d | s */
    1358             :   {
    1359       83503 :     a = negi(mulii(u,n)); d1 = d;
    1360       83503 :     a1 = diviiexact(a1, d1);
    1361       83503 :     a2 = diviiexact(a2, d1);
    1362       83503 :     s = diviiexact(s, d1);
    1363             :   }
    1364             :   else
    1365             :   {
    1366             :     GEN p2, l;
    1367       80073 :     d1 = bezout(s,d,&u1,NULL);
    1368       80073 :     if (!equali1(d1))
    1369             :     {
    1370        2044 :       a1 = diviiexact(a1,d1);
    1371        2044 :       a2 = diviiexact(a2,d1);
    1372        2044 :       s = diviiexact(s,d1);
    1373        2044 :       d = diviiexact(d,d1);
    1374             :     }
    1375       80073 :     p1 = remii(gel(x,3),d);
    1376       80073 :     p2 = remii(gel(y,3),d);
    1377       80073 :     l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
    1378       80073 :     a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
    1379             :   }
    1380      400680 :   a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
    1381      400680 :   d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
    1382      400680 :   Q = cgetg(5,t_QFB);
    1383      400680 :   if (!z) {
    1384       37632 :     g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
    1385       37632 :     b = a2;
    1386       37632 :     b2 = gel(y,2);
    1387       37632 :     v2 = d1;
    1388       37632 :     gel(Q,1) = mulii(d,b);
    1389             :   } else {
    1390             :     GEN e, q3, q4;
    1391      363048 :     if (z&1) { v3 = negi(v3); v2 = negi(v2); }
    1392      363048 :     b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
    1393      363048 :     e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
    1394      363048 :     q3 = mulii(e,v2);
    1395      363048 :     q4 = subii(q3,s);
    1396      363048 :     b2 = addii(q3,q4);
    1397      363048 :     g = diviiexact(q4,v);
    1398      363048 :     if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
    1399      363048 :     gel(Q,1) = addii(mulii(d,b), mulii(e,v));
    1400             :   }
    1401      400680 :   q1 = mulii(b, v3);
    1402      400680 :   q2 = addii(q1,n);
    1403      400680 :   gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
    1404      400680 :   gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
    1405      400680 :   gel(Q,4) = gel(x,4);
    1406      400680 :   return qfi_red_av(av, Q);
    1407             : }
    1408             : 
    1409             : GEN
    1410       11648 : nudupl(GEN x, GEN L)
    1411             : {
    1412       11648 :   pari_sp av = avma;
    1413             :   long z;
    1414             :   GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
    1415             : 
    1416       11648 :   if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
    1417       11648 :   a = gel(x,1);
    1418       11648 :   b = gel(x,2);
    1419       11648 :   d1 = bezout(b,a, &u,NULL);
    1420       11648 :   if (!equali1(d1))
    1421             :   {
    1422        4620 :     a = diviiexact(a, d1);
    1423        4620 :     b = diviiexact(b, d1);
    1424             :   }
    1425       11648 :   c = modii(negi(mulii(u,gel(x,3))), a);
    1426       11648 :   p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
    1427       11648 :   d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
    1428       11648 :   a2 = sqri(d);
    1429       11648 :   c2 = sqri(v3);
    1430       11648 :   Q = cgetg(5,t_QFB);
    1431       11648 :   if (!z) {
    1432        1281 :     g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
    1433        1281 :     b2 = gel(x,2);
    1434        1281 :     v2 = d1;
    1435        1281 :     gel(Q,1) = a2;
    1436             :   } else {
    1437             :     GEN e;
    1438       10367 :     if (z&1) { v = negi(v); d = negi(d); }
    1439       10367 :     e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
    1440       10367 :     g = diviiexact(subii(mulii(e,v2), b), v);
    1441       10367 :     b2 = addii(mulii(e,v2), mulii(v,g));
    1442       10367 :     if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
    1443       10367 :     gel(Q,1) = addii(a2, mulii(e,v));
    1444             :   }
    1445       11648 :   gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
    1446       11648 :   gel(Q,3) = addii(c2, mulii(g,v2));
    1447       11648 :   gel(Q,4) = gel(x,4);
    1448       11648 :   return qfi_red_av(av, Q);
    1449             : }
    1450             : 
    1451             : static GEN
    1452        4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
    1453             : static GEN
    1454       11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
    1455             : GEN
    1456        1008 : nupow(GEN x, GEN n, GEN L)
    1457             : {
    1458             :   pari_sp av;
    1459             :   GEN y, D;
    1460             : 
    1461        1008 :   if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
    1462        1008 :   if (!is_qfi(x)) pari_err_TYPE("nupow",x);
    1463        1008 :   if (gequal1(n)) return gcopy(x);
    1464        1008 :   av = avma;
    1465        1008 :   D = qfb_disc(x);
    1466        1008 :   y = qfi_1_by_disc(D);
    1467        1008 :   if (!signe(n)) return y;
    1468         959 :   if (!L) L = sqrtnint(absi_shallow(D), 4);
    1469         959 :   y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
    1470         959 :   if (signe(n) < 0
    1471          35 :   && !absequalii(gel(y,1),gel(y,2))
    1472          35 :   && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
    1473         959 :   return gc_GEN(av, y);
    1474             : }
    1475             : 
    1476             : /* Not stack-clean */
    1477             : GEN
    1478     1735230 : qfr5_compraw(GEN x, GEN y)
    1479             : {
    1480     1735230 :   GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
    1481     1735230 :   if (x == y)
    1482             :   {
    1483       34552 :     gel(z,4) = shifti(gel(x,4),1);
    1484       34552 :     gel(z,5) = sqrr(gel(x,5));
    1485             :   }
    1486             :   else
    1487             :   {
    1488     1700678 :     gel(z,4) = addii(gel(x,4),gel(y,4));
    1489     1700678 :     gel(z,5) = mulrr(gel(x,5),gel(y,5));
    1490             :   }
    1491     1735230 :   fix_expo(z); return z;
    1492             : }
    1493             : GEN
    1494     1735216 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
    1495     1735216 : { return qfr5_red(qfr5_compraw(x, y), S); }
    1496             : /* Not stack-clean */
    1497             : GEN
    1498     1009397 : qfr3_compraw(GEN x, GEN y)
    1499             : {
    1500     1009397 :   GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
    1501     1009397 :   return z;
    1502             : }
    1503             : GEN
    1504     1009397 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
    1505     1009397 : { return qfr3_red(qfr3_compraw(x,y), S); }
    1506             : 
    1507             : /* m > 0. Not stack-clean */
    1508             : static GEN
    1509           7 : qfr5_powraw(GEN x, long m)
    1510             : {
    1511           7 :   GEN y = NULL;
    1512          14 :   for (; m; m >>= 1)
    1513             :   {
    1514          14 :     if (m&1) y = y? qfr5_compraw(y,x): x;
    1515          14 :     if (m == 1) break;
    1516           7 :     x = qfr5_compraw(x,x);
    1517             :   }
    1518           7 :   return y;
    1519             : }
    1520             : 
    1521             : /* return x^n. Not stack-clean */
    1522             : GEN
    1523          21 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
    1524             : {
    1525          21 :   GEN y = NULL;
    1526          21 :   long i, m, s = signe(n);
    1527          21 :   if (!s) return qfr5_1(S, lg(gel(x,5)));
    1528          21 :   if (s < 0) x = qfb_inv(x);
    1529          42 :   for (i=lgefint(n)-1; i>1; i--)
    1530             :   {
    1531          21 :     m = n[i];
    1532          56 :     for (; m; m>>=1)
    1533             :     {
    1534          56 :       if (m&1) y = y? qfr5_comp(y,x,S): x;
    1535          56 :       if (m == 1 && i == 2) break;
    1536          35 :       x = qfr5_comp(x,x,S);
    1537             :     }
    1538             :   }
    1539          21 :   return y;
    1540             : }
    1541             : /* m > 0; return x^m. Not stack-clean */
    1542             : static GEN
    1543           0 : qfr3_powraw(GEN x, long m)
    1544             : {
    1545           0 :   GEN y = NULL;
    1546           0 :   for (; m; m>>=1)
    1547             :   {
    1548           0 :     if (m&1) y = y? qfr3_compraw(y,x): x;
    1549           0 :     if (m == 1) break;
    1550           0 :     x = qfr3_compraw(x,x);
    1551             :   }
    1552           0 :   return y;
    1553             : }
    1554             : /* return x^n. Not stack-clean */
    1555             : GEN
    1556        4557 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
    1557             : {
    1558        4557 :   GEN y = NULL;
    1559        4557 :   long i, m, s = signe(n);
    1560        4557 :   if (!s) return qfr3_1(S);
    1561        4557 :   if (s < 0) x = qfb_inv(x);
    1562        9130 :   for (i=lgefint(n)-1; i>1; i--)
    1563             :   {
    1564        4573 :     m = n[i];
    1565        5312 :     for (; m; m>>=1)
    1566             :     {
    1567        5292 :       if (m&1) y = y? qfr3_comp(y,x,S): x;
    1568        5292 :       if (m == 1 && i == 2) break;
    1569         739 :       x = qfr3_comp(x,x,S);
    1570             :     }
    1571             :   }
    1572        4557 :   return y;
    1573             : }
    1574             : 
    1575             : static GEN
    1576           7 : qfrinvraw(GEN x)
    1577             : {
    1578           7 :   if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
    1579           7 :  return qfbinv(x);
    1580             : }
    1581             : static GEN
    1582          14 : qfrpowraw(GEN x, long n)
    1583             : {
    1584          14 :   struct qfr_data S = { NULL, NULL, NULL };
    1585          14 :   pari_sp av = avma;
    1586          14 :   if (n==1) return gcopy(x);
    1587          14 :   if (n==-1) return qfrinvraw(x);
    1588           7 :   if (typ(x)==t_QFB)
    1589             :   {
    1590           0 :     GEN D = qfb_disc(x);
    1591           0 :     if (!n) return qfr_1(x);
    1592           0 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1593           0 :     x = qfr3_powraw(x, n);
    1594           0 :     x = qfr3_to_qfr(x, D);
    1595             :   }
    1596             :   else
    1597             :   {
    1598           7 :     GEN d0 = gel(x,2);
    1599           7 :     x = gel(x,1);
    1600           7 :     if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1601           7 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1602           7 :     x = qfr5_init(x, d0, &S);
    1603           7 :     if (labs(n) != 1) x = qfr5_powraw(x, n);
    1604           7 :     x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
    1605             :   }
    1606           7 :   return gc_GEN(av, x);
    1607             : }
    1608             : static GEN
    1609         112 : qfrpow(GEN x, GEN n)
    1610             : {
    1611         112 :   struct qfr_data S = { NULL, NULL, NULL };
    1612         112 :   long s = signe(n);
    1613         112 :   pari_sp av = avma;
    1614         112 :   if (typ(x)==t_QFB)
    1615             :   {
    1616          42 :     if (!s) return qfr_1(x);
    1617          28 :     if (s < 0) x = qfb_inv(x);
    1618          28 :     x = qfr3_init(x, &S);
    1619          28 :     x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
    1620          28 :     x = qfr3_to_qfr(x, S.D);
    1621             :   }
    1622             :   else
    1623             :   {
    1624          70 :     GEN d0 = gel(x,2);
    1625          70 :     x = gel(x,1);
    1626          70 :     if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1627          49 :     if (s < 0) x = qfb_inv(x);
    1628          49 :     x = qfr5_init(x, d0, &S);
    1629          49 :     x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
    1630          49 :     x = qfr5_to_qfr(x, S.D, mulri(d0,n));
    1631             :   }
    1632          77 :   return gc_GEN(av, x);
    1633             : }
    1634             : GEN
    1635          21 : qfbpowraw(GEN x, long n)
    1636             : {
    1637          21 :   GEN q = check_qfbext("qfbpowraw",x);
    1638          21 :   return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
    1639             : }
    1640             : /* same discriminant, no distance, no checks */
    1641             : GEN
    1642    15295208 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
    1643             : GEN
    1644     1495578 : qfbpow(GEN x, GEN n)
    1645             : {
    1646     1495578 :   GEN q = check_qfbext("qfbpow",x);
    1647     1495578 :   return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
    1648             : }
    1649             : GEN
    1650     1396403 : qfbpows(GEN x, long n)
    1651             : {
    1652     1396403 :   long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
    1653     1396403 :   affsi(n, N); return qfbpow(x, N);
    1654             : }
    1655             : 
    1656             : /* Prime forms attached to prime ideals of degree 1 */
    1657             : 
    1658             : /* assume x != 0 a t_INT, p > 0
    1659             :  * Return a t_QFB, but discriminant sign is not checked: can be used for
    1660             :  * real forms as well */
    1661             : GEN
    1662    12609651 : primeform_u(GEN x, ulong p)
    1663             : {
    1664    12609651 :   GEN c, y = cgetg(5, t_QFB);
    1665    12609571 :   pari_sp av = avma;
    1666             :   ulong b;
    1667             :   long s;
    1668             : 
    1669    12609571 :   s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
    1670             :   /* 2 or 3 mod 4 */
    1671    12609637 :   if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
    1672    12609851 :   if (p == 2) {
    1673     4163174 :     switch(s) {
    1674      589544 :       case 0: b = 0; break;
    1675     3222158 :       case 1: b = 1; break;
    1676      351472 :       case 4: b = 2; break;
    1677           0 :       default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1678           0 :                b = 0; /* -Wall */
    1679             :     }
    1680     4163174 :     c = shifti(subsi(s,x), -3);
    1681             :   } else {
    1682     8446677 :     b = Fl_sqrt(umodiu(x,p), p);
    1683     8448538 :     if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1684             :     /* mod(b) != mod2(x) ? */
    1685     8449380 :     if ((b ^ s) & 1) b = p - b;
    1686     8449380 :     c = diviuexact(shifti(subii(sqru(b), x), -2), p);
    1687             :   }
    1688    12599614 :   gel(y,3) = gc_INT(av, c);
    1689    12606912 :   gel(y,4) = icopy(x);
    1690    12609740 :   gel(y,2) = utoi(b);
    1691    12609884 :   gel(y,1) = utoipos(p); return y;
    1692             : }
    1693             : 
    1694             : /* special case: p = 1 return unit form */
    1695             : GEN
    1696      135767 : primeform(GEN x, GEN p)
    1697             : {
    1698      135767 :   const char *f = "primeform";
    1699             :   pari_sp av;
    1700      135767 :   long s, sx = signe(x), sp = signe(p);
    1701             :   GEN y, b, absp;
    1702             : 
    1703      135767 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
    1704      135767 :   if (typ(p) != t_INT) pari_err_TYPE(f,p);
    1705      135767 :   if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
    1706      135767 :   if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
    1707      135767 :   if (lgefint(p) == 3)
    1708             :   {
    1709      135753 :     ulong pp = p[2];
    1710      135753 :     if (pp == 1) {
    1711       18090 :       if (sx < 0) {
    1712             :         long r;
    1713       11258 :         if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1714       11258 :         r = mod4(x);
    1715       11258 :         if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
    1716       11258 :         return qfi_1_by_disc(x);
    1717             :       }
    1718        6832 :       y = qfr_1_by_disc(x);
    1719        6832 :       if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
    1720        6832 :       return y;
    1721             :     }
    1722      117663 :     y = primeform_u(x, pp);
    1723      117656 :     if (sx < 0) {
    1724       89957 :       if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1725       89957 :       return y;
    1726             :     }
    1727       27699 :     if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
    1728       27699 :     return gcopy( qfr3_to_qfr(y, x) );
    1729             :   }
    1730          14 :   s = mod8(x);
    1731          14 :   if (sx < 0)
    1732             :   {
    1733           7 :     if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1734           7 :     if (s) s = 8-s;
    1735             :   }
    1736          14 :   y = cgetg(5, t_QFB);
    1737             :   /* 2 or 3 mod 4 */
    1738          14 :   if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
    1739          14 :   absp = absi_shallow(p); av = avma;
    1740          14 :   b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
    1741          14 :   s &= 1; /* s = x mod 2 */
    1742             :   /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
    1743          14 :   if ((!signe(b) && s) || mod2(b) != s) b = gc_INT(av, subii(absp,b));
    1744             : 
    1745          14 :   av = avma;
    1746          14 :   gel(y,3) = gc_INT(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
    1747          14 :   gel(y,4) = icopy(x);
    1748          14 :   gel(y,2) = b;
    1749          14 :   gel(y,1) = icopy(p);
    1750          14 :   return y;
    1751             : }
    1752             : 
    1753             : static GEN
    1754     2620772 : normforms(GEN D, GEN fa)
    1755             : {
    1756             :   long i, j, k, lB, aN, sa;
    1757             :   GEN a, L, V, B, N, N2;
    1758     2620772 :   int D_odd = mpodd(D);
    1759     2620772 :   a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
    1760     2620772 :   sa = signe(a);
    1761     2620772 :   if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
    1762     1203972 :   V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
    1763     2551766 :            : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
    1764     2551766 :   if (!V) return NULL;
    1765      511966 :   N = gel(V,1); B = gel(V,2); lB = lg(B);
    1766      511966 :   N2 = shifti(N,1);
    1767      511965 :   aN = itou(diviiexact(a, N)); /* |a|/N */
    1768      511964 :   L = cgetg((lB-1)*aN+1, t_VEC);
    1769     2360563 :   for (k = 1, i = 1; i < lB; i++)
    1770             :   {
    1771     1848597 :     GEN b = shifti(gel(B,i), 1), c, C;
    1772     1848593 :     if (D_odd) b = addiu(b , 1);
    1773     1848593 :     c = diviiexact(shifti(subii(sqri(b), D), -2), a);
    1774     1848585 :     for (j = 0;; b = addii(b, N2))
    1775             :     {
    1776     2216659 :       gel(L, k++) = mkqfb(a, b, c, D);
    1777     2216672 :       if (++j == aN) break;
    1778      368074 :       C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
    1779      368074 :       c = sa > 0? addii(c, C): subii(c, C);
    1780             :     }
    1781             :   }
    1782      511966 :   return L;
    1783             : }
    1784             : 
    1785             : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
    1786             : static GEN
    1787      344321 : SL2_div_mul_e1(GEN N, GEN M)
    1788             : {
    1789      344321 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1790      344321 :   GEN A = mulii(gcoeff(N,1,1), d), B = mulii(gcoeff(N,1,2), b);
    1791      344313 :   GEN C = mulii(gcoeff(N,2,1), d), D = mulii(gcoeff(N,2,2), b);
    1792      344310 :   retmkvec2(subii(A,B), subii(C,D));
    1793             : }
    1794             : static GEN
    1795     1445674 : qfisolve_normform(GEN Q, GEN P)
    1796             : {
    1797     1445674 :   GEN a = gel(Q,1), N = gel(Q,2);
    1798     1445674 :   GEN M, b = qfi_redsl2_basecase(P, &M);
    1799     1445682 :   if (!qfb_equal(a,b)) return NULL;
    1800      102128 :   return SL2_div_mul_e1(N,M);
    1801             : }
    1802             : 
    1803             : /* Test equality modulo GL2 of two reduced forms */
    1804             : static int
    1805       61068 : GL2_qfb_equal(GEN a, GEN b)
    1806             : {
    1807       61068 :   return equalii(gel(a,1),gel(b,1))
    1808       11361 :    && absequalii(gel(a,2),gel(b,2))
    1809       72429 :    &&    equalii(gel(a,3),gel(b,3));
    1810             : }
    1811             : 
    1812             : /* Q(u,v) = p; if s < 0 return that solution; else the set of all solutions */
    1813             : static GEN
    1814       48020 : allsols(GEN Q, long s, GEN u, GEN v)
    1815             : {
    1816       48020 :   GEN w = mkvec2(u, v), b;
    1817       48016 :   if (signe(v) < 0) { u = negi(u); v = negi(v); } /* normalize for v >= 0 */
    1818       48016 :   w = mkvec2(u, v); if (s < 0) return w;
    1819       41363 :   if (!s) return mkvec(w);
    1820       38934 :   b = gel(Q,2); /* sum of the 2 solutions (if they exist) is -bv / a */
    1821       38934 :   if (signe(b))
    1822             :   { /* something to check */
    1823             :     GEN r, t;
    1824       13433 :     t = dvmdii(mulii(b, v), gel(Q,1), &r);
    1825       13433 :     if (signe(r)) return mkvec(w);
    1826        1820 :     u = addii(u, t);
    1827             :   }
    1828       27321 :   return mkvec2(w, mkvec2(negi(u), v));
    1829             : }
    1830             : static GEN
    1831      223057 : qfisolvep_all(GEN Q, GEN p, long all)
    1832             : {
    1833      223057 :   GEN R, U, V, M, N, x, q, D = qfb_disc(Q);
    1834      223053 :   long s = kronecker(D, p);
    1835             : 
    1836      223047 :   if (s < 0) return NULL;
    1837      126972 :   if (!all) s = -1; /* to indicate we want a single solution */
    1838             :   /* Solutions iff a class of maximal ideal above p is the class of Q;
    1839             :    * Two solutions iff (s > 0 and the class has order > 2), else one */
    1840      126972 :   if (!signe(gel(Q,2)))
    1841             :   { /* if principal form, use faster cornacchia */
    1842       43651 :     GEN a = gel(Q,1), c = gel(Q,3);
    1843       43651 :     if (equali1(a))
    1844             :     {
    1845       38173 :       if (!cornacchia(c, p, &M,&N)) return NULL;
    1846       33703 :       return allsols(Q, s, M, N);
    1847             :     }
    1848        5474 :     if (equali1(c))
    1849             :     {
    1850        5194 :       if (!cornacchia(a, p, &M,&N)) return NULL;
    1851         721 :       return allsols(Q, s, N, M);
    1852             :     }
    1853             :   }
    1854       83601 :   R = qfi_redsl2_basecase(Q, &U);
    1855       83601 :   if (equali1(gel(R,1)))
    1856             :   { /* principal form */
    1857       22533 :     if (!signe(gel(R,2)))
    1858             :     {
    1859        4396 :       if (!cornacchia(gel(R,3), p, &M,&N)) return NULL;
    1860         812 :       x = mkvec2(M,N);
    1861             :     }
    1862             :     else
    1863             :     { /* x^2 + xy + ((1-D)/4)y^2 = p <==> (2x + y)^2 - D y^2 = 4p */
    1864       18137 :       if (!cornacchia2(negi(D), p, &M, &N)) return NULL;
    1865        2331 :       x = subii(M,N); if (mpodd(x)) return NULL;
    1866        2331 :       x = mkvec2(shifti(x,-1), N);
    1867             :     }
    1868        3143 :     x = ZM_ZC_mul(U, x); x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1869        3143 :     return allsols(Q, s, gel(x,1), gel(x,2));
    1870             :   }
    1871       61068 :   q = qfi_redsl2_basecase(primeform(D, p), &V);
    1872       61068 :   if (!GL2_qfb_equal(R,q)) return NULL;
    1873       10451 :   if (signe(gel(R,2)) != signe(gel(q,2))) gcoeff(V,2,1) = negi(gcoeff(V,2,1));
    1874       10451 :   x = SL2_div_mul_e1(U,V); return allsols(Q, s, gel(x,1), gel(x,2));
    1875             : }
    1876             : GEN
    1877           0 : qfisolvep(GEN Q, GEN p)
    1878             : {
    1879           0 :   pari_sp av = avma;
    1880           0 :   GEN x = qfisolvep_all(Q, p, 0);
    1881           0 :   return x? gc_GEN(av, x): gc_const(av, gen_0);
    1882             : }
    1883             : 
    1884             : static GEN
    1885      770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
    1886             : {
    1887      770126 :   pari_sp av = avma, btop;
    1888      770126 :   GEN M = N, P = qfr_redsl2_basecase(Ps, rd), Q = P;
    1889             : 
    1890      770126 :   btop = avma;
    1891             :   for(;;)
    1892             :   {
    1893     5840681 :     if (qfb_equal(gel(M,1), gel(P,1)))
    1894      154084 :       return gc_upto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
    1895     5686597 :     if (qfb_equal(gel(N,1), gel(Q,1)))
    1896       77658 :       return gc_upto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
    1897     5608939 :     M = qfr_rhosl2(M, rd);
    1898     5608939 :     if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
    1899     5201735 :     Q = qfr_rhosl2(Q, rd);
    1900     5201735 :     if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
    1901     5070555 :     if (gc_needed(btop, 1)) (void)gc_all(btop, 2, &M, &Q);
    1902             :   }
    1903             : }
    1904             : 
    1905             : GEN
    1906           0 : qfrsolvep(GEN Q, GEN p)
    1907             : {
    1908           0 :   pari_sp av = avma;
    1909           0 :   GEN N, x, rd, d = qfb_disc(Q);
    1910             : 
    1911           0 :   if (kronecker(d, p) < 0) return gc_const(av, gen_0);
    1912           0 :   rd = sqrti(d);
    1913           0 :   N = qfr_redsl2(Q, rd);
    1914           0 :   x = qfrsolve_normform(N, primeform(d, p), rd);
    1915           0 :   return x? gc_upto(av, x): gc_const(av, gen_0);
    1916             : }
    1917             : 
    1918             : static GEN
    1919     1862957 : known_prime(GEN v)
    1920             : {
    1921     1862957 :   GEN p, e, fa = check_arith_all(v, "qfbsolve");
    1922     1862955 :   if (!fa) return BPSW_psp(v)? v: NULL;
    1923       42092 :   if (lg(gel(fa,1)) != 2) return NULL;
    1924       29366 :   p = gcoeff(fa,1,1);
    1925       29366 :   e = gcoeff(fa,1,2);
    1926       29366 :   return (equali1(e) && !is_pm1(p) && signe(p) > 0)? p: NULL;
    1927             : }
    1928             : static GEN
    1929     2215798 : qfsolve_normform(GEN Q, GEN f, GEN rd)
    1930     2215798 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
    1931             : static GEN
    1932     2843832 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
    1933             : {
    1934             :   GEN x, W, F, p;
    1935             :   long i, j, l;
    1936     2843832 :   if (!rd && (p = known_prime(fa))) return qfisolvep_all(Q, p, all);
    1937     2620767 :   F = normforms(qfb_disc(Q), fa);
    1938     2620772 :   if (!F) return NULL;
    1939      511966 :   if (!*Qr) *Qr = qfbredsl2(Q, rd);
    1940      511965 :   l = lg(F); W = all? cgetg(l, t_VEC): NULL;
    1941     2727253 :   for (j = i = 1; i < l; i++)
    1942     2215798 :     if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
    1943             :     {
    1944      333863 :       if (!all) return x;
    1945      333352 :       gel(W,j++) = x;
    1946             :     }
    1947      511455 :   if (j == 1) return NULL;
    1948      127456 :   setlg(W,j); return lexsort(W);
    1949             : }
    1950             : 
    1951             : static GEN
    1952     2838529 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
    1953             : static GEN
    1954     2828304 : qfbsolve_primitive(GEN Q, GEN fa, long all)
    1955             : {
    1956     2828304 :   GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
    1957     2828301 :   x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
    1958     2828286 :   if (!x) return cgetg(1, t_VEC);
    1959      174747 :   return x;
    1960             : }
    1961             : 
    1962             : /* f / g^2 */
    1963             : static GEN
    1964        5299 : famat_divsqr(GEN f, GEN g)
    1965        5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
    1966             : static GEN
    1967       10227 : qfbsolve_all(GEN Q, GEN n, long all)
    1968             : {
    1969       10227 :   GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
    1970       10227 :   GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
    1971       10227 :   long i, j, l = lg(D);
    1972       10227 :   W = all? cgetg(l, t_VEC): NULL;
    1973       25151 :   for (i = j = 1; i < l; i++)
    1974             :   {
    1975       15526 :     GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
    1976       15526 :     if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
    1977             :     {
    1978        1218 :       if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
    1979        1218 :       if (!all) return w;
    1980         616 :       gel(W,j++) = w;
    1981             :     }
    1982             :   }
    1983        9625 :   if (j == 1) return cgetg(1, t_VEC);
    1984         525 :   setlg(W,j); return lexsort(shallowconcat1(W));
    1985             : }
    1986             : 
    1987             : GEN
    1988     2838536 : qfbsolve(GEN Q, GEN n, long fl)
    1989             : {
    1990     2838536 :   pari_sp av = avma;
    1991     2838536 :   if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
    1992     2838536 :   if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
    1993     5666815 :   return gc_GEN(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
    1994     2828303 :                                   : qfbsolve_primitive(Q, n, fl & 1));
    1995             : }
    1996             : 
    1997             : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
    1998             :  * Assume d > 0 and p is prime */
    1999             : long
    2000       55247 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
    2001             : {
    2002       55247 :   pari_sp av = avma;
    2003             :   GEN b, c, r;
    2004             : 
    2005       55247 :   *px = *py = gen_0;
    2006       55247 :   b = subii(p, d);
    2007       55195 :   if (signe(b) < 0) return gc_long(av,0);
    2008       54984 :   if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
    2009       54977 :   b = Fp_sqrt(b, p); /* sqrt(-d) */
    2010       55045 :   if (!b) return gc_long(av,0);
    2011       51314 :   b = gmael(halfgcdii(p, b), 2, 2);
    2012       51354 :   c = dvmdii(subii(p, sqri(b)), d, &r);
    2013       51211 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    2014       35489 :   set_avma(av);
    2015       35483 :   *px = icopy(b);
    2016       35470 :   *py = icopy(c); return 1;
    2017             : }
    2018             : 
    2019             : static GEN
    2020     1423558 : lastqi(GEN Q)
    2021             : {
    2022     1423558 :   GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
    2023     1423553 :   if (!signe(q)) return gen_0;
    2024     1423364 :   if (!signe(s)) return p;
    2025     1416980 :   if (is_pm1(q)) return subiu(p,1);
    2026     1416984 :   return divii(p, absi_shallow(q));
    2027             : }
    2028             : 
    2029             : static long
    2030     1423563 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
    2031             : {
    2032             :   GEN M, Q, V, c, r, b2;
    2033     1423563 :   if (!signe(b)) { /* d = p,2p,3p,4p */
    2034          14 :     set_avma(av);
    2035          14 :     if (absequalii(d, px4)){ *py = gen_1; return 1; }
    2036          14 :     if (absequalii(d, p))  { *py = gen_2; return 1; }
    2037           0 :     return 0;
    2038             :   }
    2039     1423549 :   if (mod2(b) != mod2(d)) b = subii(p,b);
    2040     1423521 :   M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
    2041     1423561 :   b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
    2042     1423518 :   b2 = sqri(b);
    2043     1423513 :   if (cmpii(b2,px4) > 0)
    2044             :   {
    2045     1415137 :     b = gel(V,1); b2 = sqri(b);
    2046     1415152 :     if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
    2047             :   }
    2048     1423522 :   c = dvmdii(subii(px4, b2), d, &r);
    2049     1423508 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    2050     1378133 :   set_avma(av);
    2051     1378131 :   *px = icopy(b);
    2052     1378126 :   *py = icopy(c); return 1;
    2053             : }
    2054             : 
    2055             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
    2056             :  * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
    2057             : long
    2058     1381740 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
    2059             : {
    2060     1381740 :   pari_sp av = avma;
    2061     1381740 :   GEN b, p4 = shifti(p,2);
    2062             : 
    2063     1381719 :   *px = *py = gen_0;
    2064     1381719 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2065     1380897 :   if (absequaliu(p, 2))
    2066             :   {
    2067           7 :     set_avma(av);
    2068           7 :     switch (itou_or_0(d)) {
    2069           0 :       case 4: *px = gen_2; break;
    2070           0 :       case 7: *px = gen_1; break;
    2071           7 :       default: return 0;
    2072             :     }
    2073           0 :     *py = gen_1; return 1;
    2074             :   }
    2075     1380890 :   b = Fp_sqrt(negi(d), p);
    2076     1380934 :   if (!b) return gc_long(av,0);
    2077     1380850 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2078             : }
    2079             : 
    2080             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
    2081             : long
    2082       42714 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
    2083             : {
    2084       42714 :   pari_sp av = avma;
    2085       42714 :   GEN p4 = shifti(p,2);
    2086       42714 :   *px = *py = gen_0;
    2087       42714 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2088       42714 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2089             : }
    2090             : 
    2091             : GEN
    2092        7630 : qfbcornacchia(GEN d, GEN p)
    2093             : {
    2094        7630 :   pari_sp av = avma;
    2095             :   GEN x, y;
    2096        7630 :   if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
    2097        7630 :   if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
    2098        7630 :   if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
    2099         287 :     return gc_GEN(av, mkvec2(x, y));
    2100        7343 :   retgc_const(av, cgetg(1, t_VEC));
    2101             : }

Generated by: LCOV version 1.16