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.16.2 lcov report (development 29115-f22e516b23) Lines: 1102 1224 90.0 %
Date: 2024-03-18 08:03:28 Functions: 131 140 93.6 %
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      151514 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
      24             : {
      25             :   long r;
      26      151514 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
      27      151500 :   *s = signe(x);
      28      151500 :   if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
      29      151500 :   r = mod4(x); if (*s < 0 && r) r = 4 - r;
      30      151500 :   if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
      31      151486 :   if (pr) *pr = r;
      32      151486 : }
      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        2170 : check_quaddisc_imag(GEN x, long *r, const char *f)
      41             : {
      42        2170 :   long sx; check_quaddisc(x, &sx, r, f);
      43        2163 :   if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
      44        2163 : }
      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      963477 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
      50             : {
      51      963477 :   if (Dodd)
      52             :   {
      53      865043 :     pari_sp av = avma;
      54      865043 :     *b = gen_m1;
      55      865043 :     *c = gerepileuptoint(av, shifti(subui(1,D), -2));
      56             :   }
      57             :   else
      58             :   {
      59       98434 :     *b = gen_0;
      60       98434 :     *c = shifti(D,-2); togglesign(*c);
      61             :   }
      62      963477 : }
      63             : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
      64             : static GEN
      65      243362 : quadpoly_ii(GEN D, long Dmod4)
      66             : {
      67      243362 :   GEN b, c, y = cgetg(5,t_POL);
      68      243362 :   y[1] = evalsigne(1) | evalvarn(0);
      69      243362 :   quadpoly_bc(D, Dmod4, &b,&c);
      70      243362 :   gel(y,2) = c;
      71      243362 :   gel(y,3) = b;
      72      243362 :   gel(y,4) = gen_1; return y;
      73             : }
      74             : GEN
      75        2044 : quadpoly(GEN D)
      76             : {
      77             :   long s, Dmod4;
      78        2044 :   check_quaddisc(D, &s, &Dmod4, "quadpoly");
      79        2037 :   return quadpoly_ii(D, Dmod4);
      80             : }
      81             : GEN /* no checks */
      82      241325 : 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         560 : quadgen0(GEN x, long v)
      98             : {
      99         560 :   if (v==-1) v = fetch_user_var("w");
     100         560 :   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     1806998 : check_qfbext(const char *fun, GEN x)
     113             : {
     114     1806998 :   long t = typ(x);
     115     1806998 :   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      144368 : qfb3(GEN x, GEN y, GEN z)
     127      144368 : { retmkqfb(icopy(x), icopy(y), icopy(z), qfb_disc3(x,y,z)); }
     128             : 
     129             : static int
     130    23783633 : qfb_equal(GEN x, GEN y)
     131             : {
     132    23783633 :   return equalii(gel(x,1),gel(y,1))
     133     1592911 :       && equalii(gel(x,2),gel(y,2))
     134    25376538 :       && equalii(gel(x,3),gel(y,3));
     135             : }
     136             : 
     137             : /* valid for t_QFB, qfr3, qfr5; shallow */
     138             : static GEN
     139      878189 : qfb_inv(GEN x) {
     140      878189 :   GEN z = shallowcopy(x);
     141      878187 :   gel(z,2) = negi(gel(z,2));
     142      878187 :   return z;
     143             : }
     144             : /* valid for t_QFB, gerepile-safe */
     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       77224 : Qfb0(GEN a, GEN b, GEN c)
     151             : {
     152             :   GEN q, D;
     153       77224 :   if (!b)
     154             :   {
     155          42 :     if (c) pari_err_TYPE("Qfb",c);
     156          35 :     if (typ(a) == t_VEC && lg(a) == 4)
     157          14 :     { 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       77203 :   if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
     171       77196 :   if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
     172       77196 :   if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
     173       77196 :   q = qfb3(a, b, c); D = qfb_disc(q);
     174       77196 :   if (signe(D) < 0)
     175       42385 :   { 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       77189 :   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    16163740 : dvmdii_round(GEN b, GEN a, GEN *r)
     189             : {
     190    16163740 :   GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
     191    16163681 :   if (signe(b) >= 0) {
     192     8899574 :     if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
     193             :   } else { /* r <= 0 */
     194     7264107 :     if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
     195             :   }
     196    16163599 :   return q;
     197             : }
     198             : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
     199             : static long
     200    97204535 : dvmdsu_round(long b, ulong a, long *r)
     201             : {
     202    97204535 :   ulong a2 = a << 1, q, ub, ur;
     203    97204535 :   if (b >= 0) {
     204    61994121 :     ub = b;
     205    61994121 :     q = ub / a2;
     206    61994121 :     ur = ub % a2;
     207    61994121 :     if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
     208    21828779 :     else *r = (long)ur;
     209    61994121 :     return (long)q;
     210             :   } else { /* r <= 0 */
     211    35210414 :     ub = (ulong)-b; /* |b| */
     212    35210414 :     q = ub / a2;
     213    35210414 :     ur = ub % a2;
     214    35210414 :     if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
     215    19261101 :     else *r = -(long)ur;
     216    35210414 :     return -(long)q;
     217             :   }
     218             : }
     219             : /* reduce b mod 2*a. Update b,c */
     220             : static void
     221     2778583 : REDB(GEN a, GEN *b, GEN *c)
     222             : {
     223     2778583 :   GEN r, q = dvmdii_round(*b, a, &r);
     224     2778583 :   if (!signe(q)) return;
     225     2709447 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     226     2709447 :   *b = r;
     227             : }
     228             : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
     229             : static void
     230    97198670 : sREDB(ulong a, long *b, ulong *c)
     231             : {
     232             :   long r, q;
     233             :   ulong uz;
     234   105425002 :   if (a > LONG_MAX) return; /* b already reduced */
     235    97198670 :   q = dvmdsu_round(*b, a, &r);
     236    97213023 :   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    88986691 :   if (*b < 0)
     240             :   { /* uz = -z >= 0, q < 0 */
     241    31100502 :     if (r >= 0) /* different signs=>no overflow, exact division */
     242    16003534 :       uz = (ulong)-((*b + r)>>1);
     243             :     else
     244             :     {
     245    15096968 :       ulong ub = (ulong)-*b, ur = (ulong)-r;
     246    15096968 :       uz = (ub + ur) >> 1;
     247             :     }
     248    31100502 :     *c -= (-q) * uz; /* c -= qz */
     249             :   }
     250             :   else
     251             :   { /* uz = z >= 0, q > 0 */
     252    57886189 :     if (r <= 0)
     253    40220980 :       uz = (*b + r)>>1;
     254             :     else
     255             :     {
     256    17665209 :       ulong ub = (ulong)*b, ur = (ulong)r;
     257    17665209 :       uz = ((ub + ur) >> 1);
     258             :     }
     259    57886189 :     *c -= q * uz; /* c -= qz */
     260             :   }
     261    88986691 :   *b = r;
     262             : }
     263             : static void
     264    13385157 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
     265             : { /* REDB(a,b,c) */
     266    13385157 :   GEN r, q = dvmdii_round(*b, a, &r);
     267    13385013 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     268    13385020 :   *b = r;
     269    13385020 :   *u2 = subii(*u2, mulii(q, u1));
     270    13385034 : }
     271             : 
     272             : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
     273             : static GEN
     274     6619350 : qfbredsl2_imag_basecase(GEN q, GEN *U)
     275             : {
     276     6619350 :   pari_sp av = avma;
     277             :   GEN z, u1,u2,v1,v2,Q;
     278     6619350 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     279             :   long cmp;
     280     6619350 :   u1 = gen_1; u2 = gen_0;
     281     6619350 :   cmp = abscmpii(a, b);
     282     6619349 :   if (cmp < 0)
     283     2084245 :     REDBU(a,&b,&c, u1,&u2);
     284     4535104 :   else if (cmp == 0 && signe(b) < 0)
     285             :   { /* b = -a */
     286       11919 :     b = negi(b);
     287       11919 :     u2 = gen_1;
     288             :   }
     289             :   for(;;)
     290             :   {
     291    17920184 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     292    11300718 :     swap(a,c); b = negi(b);
     293    11300912 :     z = u1; u1 = u2; u2 = negi(z);
     294    11300920 :     REDBU(a,&b,&c, u1,&u2);
     295    11300840 :     if (gc_needed(av, 1)) {
     296           7 :       if (DEBUGMEM>1) pari_warn(warnmem, "qfbredsl2");
     297           7 :       gerepileall(av, 5, &a,&b,&c, &u1,&u2);
     298             :     }
     299             :   }
     300     6619334 :   if (cmp == 0 && signe(b) < 0)
     301             :   {
     302       19138 :     b = negi(b);
     303       19138 :     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     6619334 :   z = shifti(subii(b, gel(q,2)), -1);
     308     6619337 :   v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
     309     6619316 :   z = subii(z, b);
     310     6619321 :   v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
     311     6619297 :   *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
     312     6619358 :   Q = mkqfb(a,b,c,gel(q,4));
     313     6619358 :   return gc_all(av, 2, &Q, U);
     314             : }
     315             : 
     316             : static GEN
     317     1068236 : setq_b0(ulong a, ulong c, GEN D)
     318     1068236 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
     319             : /* assume |sb| = 1 */
     320             : static GEN
     321    80283607 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
     322    80283607 : { 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      958720 : qfbred_imag_1_b0(ulong a, ulong c, GEN D)
     326      958720 : { 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    81448930 : qfbred_imag_1(pari_sp av, GEN a, GEN b, GEN c, GEN D)
     331             : {
     332             :   ulong ua, ub, uc;
     333             :   long sb;
     334             :   for(;;)
     335      142538 :   { /* at most twice */
     336    81448930 :     long lb = lgefint(b); /* <= 3 after first loop */
     337    81448930 :     if (lb == 2) return qfbred_imag_1_b0(a[2],c[2], D);
     338    80490210 :     if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
     339      141227 :     REDB(a,&b,&c);
     340      141047 :     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 qfbred_imag_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      141047 :     swap(a,c); b = negi(b);
     349             :   }
     350             :   /* b != 0 */
     351    80348983 :   set_avma(av);
     352    80351088 :   ua = a[2];
     353    80351088 :   ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
     354    80351088 :   uc = c[2];
     355    80351088 :   if (ua < ub)
     356    29148523 :     sREDB(ua, &sb, &uc);
     357    51202565 :   else if (ua == ub && sb < 0) sb = (long)ub;
     358   148460389 :   while(ua > uc)
     359             :   {
     360    68072023 :     lswap(ua,uc); sb = -sb;
     361    68072023 :     sREDB(ua, &sb, &uc);
     362             :   }
     363    80388366 :   if (!sb) return setq_b0(ua, uc, D);
     364             :   else
     365             :   {
     366    80278848 :     long s = 1;
     367    80278848 :     if (sb < 0)
     368             :     {
     369    30125751 :       sb = -sb;
     370    30125751 :       if (ua != uc) s = -1;
     371             :     }
     372    80278848 :     return setq(ua, sb, uc, s, D);
     373             :   }
     374             : }
     375             : 
     376             : static GEN
     377           7 : rhoimag(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 gerepilecopy(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    10217732 : fix_expo(GEN x)
     415             : {
     416    10217732 :   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    10217732 : }
     421             : 
     422             : /* (1/2) log (d * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
     423             : GEN
     424      184660 : qfr5_dist(GEN e, GEN d, long prec)
     425             : {
     426      184660 :   GEN t = logr_abs(d);
     427      184660 :   if (signe(e)) {
     428           0 :     GEN u = mulir(e, mplog2(prec));
     429           0 :     shiftr_inplace(u, EMAX); t = addrr(t, u);
     430             :   }
     431      184660 :   shiftr_inplace(t, -1); return t;
     432             : }
     433             : 
     434             : static void
     435    14152279 : rho_get_BC(GEN *B, GEN *C, GEN a, GEN b, GEN c, struct qfr_data *S)
     436             : {
     437             :   GEN t, u, q;
     438    14152279 :   t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: absi_shallow(c);
     439    14152279 :   q = truedvmdii(addii(t, b), shifti(c,1), &u);
     440    14152279 :   *B = subii(t, u); /* |t| - ((|t|+b) % 2c) */
     441    14152279 :   *C = subii(a, mulii(q, subii(b, mulii(q,c))));
     442    14152279 : }
     443             : /* Not stack-clean */
     444             : GEN
     445     1139439 : qfr3_rho(GEN x, struct qfr_data *S)
     446             : {
     447     1139439 :   GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3);
     448     1139439 :   rho_get_BC(&B, &C, a, b, c, S);
     449     1139439 :   return mkvec3(c, B, C);
     450             : }
     451             : 
     452             : /* Not stack-clean */
     453             : GEN
     454     8486548 : qfr5_rho(GEN x, struct qfr_data *S)
     455             : {
     456     8486548 :   GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3), y;
     457     8486548 :   long sb = signe(b);
     458     8486548 :   rho_get_BC(&B, &C, a, b, c, S);
     459     8486548 :   y = mkvec5(c, B, C, gel(x,4), gel(x,5));
     460     8486548 :   if (sb) {
     461     8482565 :     GEN t = subii(sqri(b), S->D);
     462     8482565 :     if (sb < 0)
     463     2509423 :       t = divir(t, sqrr(subir(b,S->sqrtD)));
     464             :     else
     465     5973142 :       t = divri(sqrr(addir(b,S->sqrtD)), t);
     466             :     /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
     467     8482565 :     gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
     468        3983 :   } else gel(y,5) = negr(gel(y,5));
     469     8486548 :   return y;
     470             : }
     471             : 
     472             : /* Not stack-clean */
     473             : GEN
     474      217700 : qfr_to_qfr5(GEN x, long prec)
     475      217700 : { 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    17712399 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
     509             : {
     510             :   GEN t;
     511    17712399 :   if (signe(b) <= 0 || abscmpii(b, isqrtD) > 0) return 0;
     512     5271321 :   t = addii_sign(isqrtD,1, shifti(a,1),-1); /* floor(sqrt(D)) - |2a| */
     513     1099083 :   return signe(t) < 0 ? abscmpii(b, t) >= 0
     514     6370404 :                       : abscmpii(b, t) > 0;
     515             : }
     516             : 
     517             : /* Not stack-clean */
     518             : GEN
     519     1952804 : qfr5_red(GEN x, struct qfr_data *S)
     520             : {
     521     1952804 :   pari_sp av = avma;
     522     8465121 :   while (!ab_isreduced(gel(x,1), gel(x,2), S->isqrtD))
     523             :   {
     524     6512317 :     x = qfr5_rho(x, S);
     525     6512317 :     if (gc_needed(av,2))
     526             :     {
     527           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
     528           0 :       x = gerepilecopy(av, x);
     529             :     }
     530             :   }
     531     1952804 :   return x;
     532             : }
     533             : /* Not stack-clean */
     534             : GEN
     535     1172756 : qfr3_red(GEN x, struct qfr_data *S)
     536             : {
     537     1172756 :   pari_sp av = avma;
     538     1172756 :   GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
     539     5699048 :   while (!ab_isreduced(a, b, S->isqrtD))
     540             :   {
     541             :     GEN B, C;
     542     4526292 :     rho_get_BC(&B, &C, a, b, c, S);
     543     4526292 :     a = c; b = B; c = C;
     544     4526292 :     if (gc_needed(av,2))
     545             :     {
     546           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
     547           0 :       gerepileall(av, 3, &a, &b, &c);
     548             :     }
     549             :   }
     550     1172756 :   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 : qfbred_real_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    13379252 : _rhorealsl2(GEN *pa, GEN *pb, GEN *pc, GEN *pu1, GEN *pu2, GEN *pv1,
     615             :             GEN *pv2, GEN rd)
     616             : {
     617    13379252 :   GEN C = mpabs_shallow(*pc), t = addii(*pb, gmax_shallow(rd,C));
     618    13379252 :   GEN r, q = truedvmdii(t, shifti(C,1), &r);
     619    13379252 :   GEN a = *pa, b= *pb, c = *pc;
     620    13379252 :   if (signe(c) < 0) togglesign(q);
     621    13379252 :   *pa = *pc;
     622    13379252 :   *pb = subii(t, addii(r, *pb));
     623    13379252 :   *pc = subii(a,mulii(q,subii(b, mulii(q,c))));
     624    13379252 :   r = *pu1; *pu1 = *pv1; *pv1 = subii(mulii(q, *pv1), r);
     625    13379252 :   r = *pu2; *pu2 = *pv2; *pv2 = subii(mulii(q, *pv2), r);
     626    13379252 : }
     627             : 
     628             : static GEN
     629    10810674 : rhorealsl2(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 :   _rhorealsl2(&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      979652 : qfbredsl2_real_basecase(GEN V, GEN rd)
     641             : {
     642      979652 :   pari_sp av = avma;
     643      979652 :   GEN u1 = gen_1, u2 = gen_0, v1 = gen_0, v2 = gen_1;
     644      979652 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
     645     3548230 :   while (!ab_isreduced(a,b,rd))
     646             :   {
     647     2568578 :     _rhorealsl2(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
     648     2568578 :     if (gc_needed(av, 1))
     649             :     {
     650           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfbredsl2");
     651           0 :       gerepileall(av, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
     652             :     }
     653             :   }
     654      979652 :   return gerepilecopy(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             : #define QFBRED_LIMIT 9000
     667             : 
     668             : /* pqfb: qf with positive coefficients */
     669             : 
     670             : static int
     671     5314748 : lti2n(GEN a, long m)
     672             : {
     673     5314748 :   return signe(a) < 0 || expi(a) < m ? 1 : 0;
     674             : }
     675             : 
     676             : static GEN
     677     2065206 : pqfbred_1(GEN Q, long m, GEN U)
     678             : {
     679     2065206 :   GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     680     2065206 :   if (abscmpii(a, c) < 0)
     681             :   {
     682             :     GEN t, at, r;
     683     1032327 :     GEN r2 = addii(shifti(a, m + 2), d);
     684     1032327 :     long e2 = expi(r2);
     685     1032327 :     r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
     686     1032327 :     t = truedivii(subii(b, r), shifti(a,1));
     687     1032327 :     if (signe(t)==0) pari_err_BUG("pqfbred_1");
     688     1032327 :     at = mulii(a,t);
     689     1032327 :     c = addii(subii(c, mulii(b, t)), mulii(at, t));
     690     1032327 :     b = subii(b, shifti(at,1));
     691     1032327 :     gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,1,1), t));
     692     1032327 :     gcoeff(U,2,2) = subii( gcoeff(U,2,2), mulii(gcoeff(U,2,1), t));
     693             :   } else
     694             :   {
     695             :     GEN t, ct, r;
     696     1032879 :     GEN r2 = addii(shifti(c, m + 2), d);
     697     1032879 :     long e2 = expi(r2);
     698     1032879 :     r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
     699     1032879 :     t = truedivii(subii(b, r), shifti(c,1));
     700     1032879 :     if (signe(t)==0) pari_err_BUG("pqfbred_1");
     701     1032879 :     ct = mulii(c, t);
     702     1032879 :     a = addii(subii(a, mulii(b, t)), mulii(ct, t));
     703     1032879 :     b = subii(b, shifti(ct, 1));
     704     1032879 :     gcoeff(U,1,1) = subii(gcoeff(U,1,1), mulii(gcoeff(U,1,2), t));
     705     1032879 :     gcoeff(U,2,1) = subii(gcoeff(U,2,1), mulii(gcoeff(U,2,2), t));
     706             :   }
     707     2065206 :   return mkqfb(a,b,c,d);
     708             : }
     709             : 
     710             : static int
     711     2201585 : is_minimal(GEN Q, long m)
     712             : {
     713     2201585 :   pari_sp av = avma;
     714     2201585 :   GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3);
     715     5314748 :   return gc_bool(av, lti2n(addii(subii(a,b), c), m)
     716     2072745 :                  || (lti2n(subii(b, shifti(a,1)), m+1)
     717     1040418 :                      && lti2n(subii(b, shifti(c,1)), m+1)));
     718             : }
     719             : 
     720             : static GEN
     721      134615 : pqfbred_iter_1(GEN Q, ulong m, GEN U)
     722             : {
     723      134615 :   pari_sp av = avma;
     724     2068245 :   while (!is_minimal(Q,m))
     725             :   {
     726     1933630 :     Q = pqfbred_1(Q, m, U);
     727     1933630 :     if (gc_needed(av, 1))
     728             :     {
     729           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"pqfbred_iter_1, lc = %ld", expi(gel(Q,3)));
     730           0 :       gerepileall(av, 3, &Q, &gel(U,1), &gel(U,2));
     731             :     }
     732             :   }
     733      134615 :   return Q;
     734             : }
     735             : 
     736             : static GEN
     737       66501 : pqfbred_basecase(GEN Q, ulong m, GEN *pt_U)
     738             : {
     739       66501 :   pari_sp av = avma;
     740       66501 :   GEN  U = matid(2);
     741       66501 :   Q = pqfbred_iter_1(Q, m, U);
     742       66501 :   *pt_U = U;
     743       66501 :   return gc_all(av, 2, &Q, pt_U);
     744             : }
     745             : 
     746             : static long
     747    86884283 : qfb_maxexpi(GEN Q)
     748    86884283 : { return 1+maxss(expi(gel(Q,1)), maxss(expi(gel(Q,2)), expi(gel(Q,3)))); }
     749             : 
     750             : static long
     751      136228 : qfb_minexpi(GEN Q)
     752             : {
     753      136228 :   long m =  minss(expi(gel(Q,1)), minss(expi(gel(Q,2)), expi(gel(Q,3))));
     754      136228 :   return m < 0 ? 0: m;
     755             : }
     756             : 
     757             : static GEN
     758      134615 : pqfbred_rec(GEN Q, long m, GEN *pt_U)
     759             : {
     760      134615 :   pari_sp av = avma;
     761             :   GEN Q0, Q1, QR;
     762      134615 :   GEN d = qfb_disc(Q);
     763      134615 :   long n = qfb_maxexpi(Q) - m;
     764             :   long h;
     765      134615 :   int going_to_r8 = 0;
     766             :   GEN U;
     767             : 
     768      134615 :   if (n < 170)
     769       66501 :     return pqfbred_basecase(Q, m, pt_U);
     770             : 
     771       68114 :   if (qfb_minexpi(Q) <= m + 2)
     772             :   {
     773           0 :     U = matid(2);
     774           0 :     QR = Q;
     775             :   }
     776             :   else
     777             :   {
     778             :     long p, mm;
     779       68114 :     if (m <= n)
     780             :     {
     781         942 :       mm = m;
     782         942 :       p = 0;
     783         942 :       Q1 = Q;
     784             :     } else
     785             :     {
     786       67172 :       mm = n;
     787       67172 :       p = m + 1 - n;
     788       67172 :       Q0 = mkvec3(remi2n(gel(Q,1), p), remi2n(gel(Q,2), p), remi2n(gel(Q,3), p));
     789       67172 :       Q1 = qfb3(shifti(gel(Q,1),-p), shifti(gel(Q,2),-p), shifti(gel(Q,3),-p));
     790             :     }
     791       68114 :     h = mm + (n>>1);
     792       68114 :     if (qfb_minexpi(Q1) <= h)
     793             :     {
     794         288 :       U = matid(2);
     795         288 :       QR = Q1;
     796             :     }
     797             :     else
     798       67826 :       QR = pqfbred_rec(Q1, h, &U);
     799      199690 :     while (qfb_maxexpi(QR) > h)
     800             :     {
     801      133340 :       if (is_minimal(QR, mm))
     802             :       {
     803        1764 :         going_to_r8 = 1;
     804        1764 :         break;
     805             :       }
     806      131576 :       QR = pqfbred_1(QR, mm, U);
     807             :     }
     808       68114 :     if (!going_to_r8)
     809             :     {
     810             :       GEN V;
     811       66350 :       QR = pqfbred_rec(QR, mm, &V);
     812       66350 :       U = ZM2_mul(U,V);
     813             :     }
     814       68114 :     if (p > 0)
     815             :     {
     816       67172 :       GEN Q0U = qfb_ZM_apply(Q0,U);
     817      134344 :       QR = mkqfb(addii(shifti(gel(QR,1), p), gel(Q0U,1)),
     818       67172 :                  addii(shifti(gel(QR,2), p), gel(Q0U,2)),
     819       67172 :                  addii(shifti(gel(QR,3), p), gel(Q0U,3)), d);
     820             :     }
     821             :   }
     822       68114 :   QR = pqfbred_iter_1(QR, m, U);
     823       68114 :   *pt_U = U;
     824       68114 :   return gc_all(av, 2, &QR, pt_U);
     825             : }
     826             : 
     827             : static GEN
     828      209526 : qfbredsl2_real(GEN Q, GEN isqrtD)
     829             : {
     830      209526 :   pari_sp av = avma;
     831      209526 :   if (2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     832      209526 :     return qfbredsl2_real_basecase(Q, isqrtD);
     833             :   else
     834             :   {
     835           0 :     GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     836           0 :     GEN Qf, Qr, W, U, t = NULL;
     837           0 :     long sa = signe(a), sb;
     838           0 :     if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
     839           0 :     if (signe(c) < 0)
     840             :     {
     841             :       GEN at;
     842           0 :       t  = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
     843           0 :       at = mulii(a,t);
     844           0 :       c = addii(subii(c, mulii(b, t)), mulii(at, t));
     845           0 :       b = subii(b, shifti(at,1));
     846             :     }
     847           0 :     sb = signe(b);
     848           0 :     Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
     849           0 :     if (sa < 0)
     850           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     851           0 :     if (sb < 0)
     852             :     {
     853           0 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     854           0 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     855             :     }
     856           0 :     if (t)
     857             :     {
     858           0 :       gcoeff(U,1,1) = subii( gcoeff(U,1,1), mulii(gcoeff(U,2,1), t));
     859           0 :       gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,2,2), t));
     860             :     }
     861           0 :     W = qfbredsl2_real_basecase(Qr, isqrtD);
     862           0 :     Qf = gel(W,1);
     863           0 :     U = ZM2_mul(U,gel(W,2));
     864           0 :     return gerepilecopy(av, mkvec2(Qf,U));
     865             :   }
     866             : }
     867             : 
     868             : static GEN
     869     5029010 : qfbredsl2_imag(GEN Q)
     870             : {
     871     5029010 :   pari_sp av = avma;
     872             :   GEN Qt, U;
     873     5029010 :   if (2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     874     5028582 :     Qt = qfbredsl2_imag_basecase(Q, &U);
     875             :   else
     876             :   {
     877         432 :     long sb = signe(gel(Q,2));
     878             :     GEN Qr, W;
     879         432 :     Qr = pqfbred_rec(sb < 0 ? mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4)): Q, 0, &U);
     880         432 :     Qt = qfbredsl2_imag_basecase(Qr,&W);
     881         432 :     U = ZM2_mul(U,W);
     882         432 :     if (sb < 0)
     883             :     {
     884         230 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     885         230 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     886             :     }
     887             :   }
     888     5029028 :   return gerepilecopy(av, mkvec2(Qt,U));
     889             : }
     890             : 
     891             : GEN
     892     4718767 : redimagsl2(GEN Q, GEN *U)
     893             : {
     894     4718767 :   GEN q = qfbredsl2_imag(Q);
     895     4718794 :   *U = gel(q,2); return gel(q,1);
     896             : }
     897             : 
     898             : GEN
     899      519767 : qfbredsl2(GEN q, GEN isD)
     900             : {
     901             :   pari_sp av;
     902      519767 :   if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
     903      519767 :   if (qfb_is_qfi(q))
     904             :   {
     905      310234 :     if (isD) pari_err_TYPE("qfbredsl2", isD);
     906      310234 :     return qfbredsl2_imag(q);
     907             :   }
     908      209533 :   av = avma;
     909      209533 :   if (!isD) isD = sqrti(qfb_disc(q));
     910      208068 :   else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
     911      209526 :   return gerepileupto(av, qfbredsl2_real(q, isD));
     912             : }
     913             : 
     914             : static GEN
     915         427 : qfbred_real_i(GEN Q, long flag, GEN isqrtD, GEN sqrtD)
     916             : {
     917         427 :   if (typ(Q)!=t_QFB || 2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     918         427 :     return qfbred_real_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), sb;
     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 :     sb = signe(b);
     934           0 :     Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
     935           0 :     if (sa < 0)
     936           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     937           0 :     W = qfbred_real_basecase_i(Qr, flag, isqrtD, sqrtD);
     938           0 :     return gel(W,1);
     939             :   }
     940             : }
     941             : 
     942             : static GEN
     943          63 : qfbred_real(GEN x) { return qfbred_real_i(x,0,NULL,NULL); }
     944             : 
     945             : static GEN
     946    81481824 : qfbred_imag_basecase_av(pari_sp av, GEN q)
     947             : {
     948    81481824 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
     949    81481824 :   long cmp, lc = lgefint(c);
     950             : 
     951    81481824 :   if (lgefint(a) == 3 && lc == 3) return qfbred_imag_1(av, a, b, c, D);
     952      912199 :   cmp = abscmpii(a, b);
     953      917155 :   if (cmp < 0)
     954      436395 :     REDB(a,&b,&c);
     955      480760 :   else if (cmp == 0 && signe(b) < 0)
     956          33 :     b = negi(b);
     957             :   for(;;)
     958             :   {
     959     3118296 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     960     2927388 :     lc = lgefint(a); /* lg(future c): we swap a & c next */
     961     2927388 :     if (lc == 3) return qfbred_imag_1(av, a, b, c, D);
     962     2201141 :     swap(a,c); b = negi(b); /* apply rho */
     963     2201141 :     REDB(a,&b,&c);
     964     2201141 :     if (gc_needed(av, 2))
     965             :     {
     966           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfbred_imag, lc = %ld", lc);
     967           0 :       gerepileall(av, 3, &a,&b,&c);
     968             :     }
     969             :   }
     970      190908 :   if (cmp == 0 && signe(b) < 0) b = negi(b);
     971      190908 :   return gerepilecopy(av, mkqfb(a, b, c, D));
     972             : }
     973             : static GEN
     974    81313533 : qfbred_imag_av(pari_sp av, GEN Q)
     975             : {
     976    81313533 :   if (2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     977    81478856 :     return qfbred_imag_basecase_av(av, Q);
     978             :   else
     979             :   {
     980           7 :     long sb = signe(gel(Q,2));
     981           7 :     GEN U, Qr = pqfbred_rec(sb < 0 ? mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4)): Q, 0, &U);
     982           7 :     return qfbred_imag_basecase_av(av, Qr);
     983             :   }
     984             : }
     985             : 
     986             : static GEN
     987    17936045 : qfbred_imag(GEN q) { return qfbred_imag_av(avma, q); }
     988             : 
     989             : GEN
     990       54877 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
     991             : {
     992             :   pari_sp av;
     993       54877 :   GEN q = check_qfbext("qfbred",x);
     994       54877 :   if (qfb_is_qfi(q)) return (flag & qf_STEP)? rhoimag(x): qfbred_imag(x);
     995         364 :   if (typ(x)==t_QFB) flag |= qf_NOD;
     996          49 :   else               flag &= ~qf_NOD;
     997         364 :   av = avma;
     998         364 :   return gerepilecopy(av, qfbred_real_i(x,flag,isqrtD,sqrtD));
     999             : }
    1000             : /* t_QFB */
    1001             : GEN
    1002    17881542 : qfbred_i(GEN x) { return qfb_is_qfi(x)? qfbred_imag(x): qfbred_real(x); }
    1003             : GEN
    1004       53064 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
    1005             : /***********************************************************************/
    1006             : /**                                                                   **/
    1007             : /**                         Composition                               **/
    1008             : /**                                                                   **/
    1009             : /***********************************************************************/
    1010             : 
    1011             : static void
    1012    27275279 : qfb_sqr(GEN z, GEN x)
    1013             : {
    1014             :   GEN c, d1, x2, v1, v2, c3, m, p1, r;
    1015             : 
    1016    27275279 :   d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
    1017    27275159 :   c = gel(x,3);
    1018    27275159 :   m = mulii(c,x2);
    1019    27274916 :   if (equali1(d1))
    1020    20652776 :     v1 = v2 = gel(x,1);
    1021             :   else
    1022             :   {
    1023     6622180 :     v1 = diviiexact(gel(x,1),d1);
    1024     6622180 :     v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
    1025     6622178 :     c = mulii(c, d1);
    1026             :   }
    1027    27274952 :   togglesign(m);
    1028    27274848 :   r = modii(m,v2);
    1029    27274713 :   p1 = mulii(r, v1);
    1030    27274835 :   c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
    1031    27274875 :   gel(z,1) = mulii(v1,v2);
    1032    27274878 :   gel(z,2) = addii(gel(x,2), shifti(p1,1));
    1033    27274901 :   gel(z,3) = diviiexact(c3,v2);
    1034    27274780 : }
    1035             : /* z <- x * y */
    1036             : static void
    1037    65148797 : qfb_comp(GEN z, GEN x, GEN y)
    1038             : {
    1039             :   GEN n, c, d, y1, v1, v2, c3, m, p1, r;
    1040             : 
    1041    65148797 :   if (x == y) { qfb_sqr(z,x); return; }
    1042    38667216 :   n = shifti(subii(gel(y,2),gel(x,2)), -1);
    1043    38496675 :   v1 = gel(x,1);
    1044    38496675 :   v2 = gel(y,1);
    1045    38496675 :   c  = gel(y,3);
    1046    38496675 :   d = bezout(v2,v1,&y1,NULL);
    1047    38583561 :   if (equali1(d))
    1048    22429507 :     m = mulii(y1,n);
    1049             :   else
    1050             :   {
    1051    16191903 :     GEN s = subii(gel(y,2), n);
    1052    16190151 :     GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
    1053    16193157 :     if (!equali1(d1))
    1054             :     {
    1055     7912511 :       v1 = diviiexact(v1,d1);
    1056     7911378 :       v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
    1057     7910902 :       v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
    1058     7911482 :       c = mulii(c, d1);
    1059             :     }
    1060    16191810 :     m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
    1061             :   }
    1062    38565241 :   togglesign(m);
    1063    38497778 :   r = modii(m, v1);
    1064    38447163 :   p1 = mulii(r, v2);
    1065    38496589 :   c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
    1066    38494906 :   gel(z,1) = mulii(v1,v2);
    1067    38490977 :   gel(z,2) = addii(gel(y,2), shifti(p1,1));
    1068    38493759 :   gel(z,3) = diviiexact(c3,v1);
    1069             : }
    1070             : 
    1071             : /* not meant to be efficient */
    1072             : static GEN
    1073          84 : qfb_comp_gen(GEN x, GEN y)
    1074             : {
    1075          84 :   GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
    1076          84 :   GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
    1077          84 :   GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
    1078          84 :   GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
    1079             : 
    1080          84 :   if (!is_pm1(cx))
    1081             :   {
    1082          14 :     a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
    1083          14 :     c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
    1084             :   }
    1085          84 :   if (!is_pm1(cy))
    1086             :   {
    1087          28 :     a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
    1088          28 :     b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
    1089             :   }
    1090          84 :   D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
    1091         133 :   if (!Z_issquareall(diviiexact(d1, D), &n1) ||
    1092          84 :       !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
    1093          49 :   A = mulii(a1, n2);
    1094          49 :   B = mulii(a2, n1);
    1095          49 :   C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
    1096          49 :   U = ZV_extgcd(mkvec3(A, B, C));
    1097          49 :   m = gel(U,1); U = gmael(U,2,3);
    1098          49 :   A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
    1099          49 :   B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
    1100          49 :   C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
    1101          49 :   C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
    1102          49 :   B = addii(A, addii(B, C));
    1103          49 :   m2 = sqri(m);
    1104          49 :   A = diviiexact(mulii(a1, a2), m2);
    1105          49 :   C = diviiexact(shifti(subii(sqri(B),D), -2), A);
    1106          49 :   cx = mulii(cx, cy);
    1107          49 :   if (!is_pm1(cx))
    1108             :   {
    1109          14 :     A = mulii(A, cx); B = mulii(B, cx);
    1110          14 :     C = mulii(C, cx); D = mulii(D, sqri(cx));
    1111             :   }
    1112          49 :   return mkqfb(A, B, C, D);
    1113             : }
    1114             : 
    1115             : static GEN
    1116    62408380 : qficomp0(GEN x, GEN y, int raw)
    1117             : {
    1118    62408380 :   pari_sp av = avma;
    1119    62408380 :   GEN z = cgetg(5,t_QFB);
    1120    62404543 :   gel(z,4) = gel(x,4);
    1121    62404543 :   qfb_comp(z, x,y);
    1122    62166817 :   if (raw) return gerepilecopy(av,z);
    1123    62165039 :   return qfbred_imag_av(av, z);
    1124             : }
    1125             : static GEN
    1126         441 : qfrcomp0(GEN x, GEN y, int raw)
    1127             : {
    1128         441 :   pari_sp av = avma;
    1129         441 :   GEN dx = NULL, dy = NULL;
    1130         441 :   GEN z = cgetg(5,t_QFB);
    1131         441 :   if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1132         441 :   if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
    1133         441 :   gel(z,4) = gel(x,4);
    1134         441 :   qfb_comp(z, x,y);
    1135         441 :   if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
    1136         441 :   if (!raw) z = qfbred_real(z);
    1137         441 :   return gerepilecopy(av, z);
    1138             : }
    1139             : /* same discriminant, no distance, no checks */
    1140             : GEN
    1141    27285482 : qfbcomp_i(GEN x, GEN y)
    1142    27285482 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
    1143             : GEN
    1144      129398 : qfbcomp(GEN x, GEN y)
    1145             : {
    1146      129398 :   GEN qx = check_qfbext("qfbcomp", x);
    1147      129398 :   GEN qy = check_qfbext("qfbcomp", y);
    1148      129398 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1149             :   {
    1150          63 :     pari_sp av = avma;
    1151          63 :     GEN z = qfb_comp_gen(qx, qy);
    1152          63 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1153           7 :       pari_err_IMPL("Shanks's distance in general composition");
    1154          56 :     if (!z) pari_err_OP("*",x,y);
    1155          21 :     return gerepileupto(av, qfbred(z));
    1156             :   }
    1157      129335 :   return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
    1158             : }
    1159             : /* same discriminant, no distance, no checks */
    1160             : GEN
    1161           0 : qfbcompraw_i(GEN x, GEN y)
    1162           0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
    1163             : GEN
    1164        2198 : qfbcompraw(GEN x, GEN y)
    1165             : {
    1166        2198 :   GEN qx = check_qfbext("qfbcompraw", x);
    1167        2198 :   GEN qy = check_qfbext("qfbcompraw", y);
    1168        2198 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1169             :   {
    1170          21 :     pari_sp av = avma;
    1171          21 :     GEN z = qfb_comp_gen(qx, qy);
    1172          21 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1173           0 :       pari_err_IMPL("Shanks's distance in general composition");
    1174          21 :     if (!z) pari_err_OP("qfbcompraw",x,y);
    1175          21 :     return gerepilecopy(av, z);
    1176             :   }
    1177        2177 :   if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
    1178        2177 :   return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
    1179             : }
    1180             : 
    1181             : static GEN
    1182      793666 : qfisqr0(GEN x, long raw)
    1183             : {
    1184      793666 :   pari_sp av = avma;
    1185      793666 :   GEN z = cgetg(5,t_QFB);
    1186      793666 :   gel(z,4) = gel(x,4);
    1187      793666 :   qfb_sqr(z,x);
    1188      793666 :   if (raw) return gerepilecopy(av,z);
    1189      793666 :   return qfbred_imag_av(av, z);
    1190             : }
    1191             : static GEN
    1192          35 : qfrsqr0(GEN x, long raw)
    1193             : {
    1194          35 :   pari_sp av = avma;
    1195          35 :   GEN dx = NULL, z = cgetg(5,t_QFB);
    1196          35 :   if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1197          35 :   gel(z,4) = gel(x,4); qfb_sqr(z,x);
    1198          35 :   if (dx) z = mkvec2(z, shiftr(dx,1));
    1199          35 :   if (!raw) z = qfbred_real(z);
    1200          35 :   return gerepilecopy(av, z);
    1201             : }
    1202             : /* same discriminant, no distance, no checks */
    1203             : GEN
    1204      698063 : qfbsqr_i(GEN x)
    1205      698063 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
    1206             : GEN
    1207       95638 : qfbsqr(GEN x)
    1208             : {
    1209       95638 :   GEN qx = check_qfbext("qfbsqr", x);
    1210       95638 :   return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
    1211             : }
    1212             : 
    1213             : static GEN
    1214        6867 : qfr_1_by_disc(GEN D)
    1215             : {
    1216             :   GEN y, r, s;
    1217        6867 :   check_quaddisc_real(D, NULL, "qfr_1_by_disc");
    1218        6867 :   y = cgetg(5,t_QFB);
    1219        6867 :   s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
    1220        6867 :   if (mpodd(r))
    1221             :   {
    1222        3535 :     s = subiu(s,1);
    1223        3535 :     r = subii(r, addiu(shifti(s, 1), 1));
    1224        3535 :     r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
    1225             :   }
    1226             :   else
    1227        3332 :   { r = shifti(r, -2); set_avma((pari_sp)s); }
    1228        6867 :   gel(y,1) = gen_1;
    1229        6867 :   gel(y,2) = s;
    1230        6867 :   gel(y,3) = icopy(r);
    1231        6867 :   gel(y,4) = icopy(D); return y;
    1232             : }
    1233             : 
    1234             : static GEN
    1235          35 : qfr_disc(GEN x)
    1236          35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
    1237             : 
    1238             : static GEN
    1239          35 : qfr_1(GEN x)
    1240          35 : { return qfr_1_by_disc(qfr_disc(x)); }
    1241             : 
    1242             : static void
    1243           0 : qfr_1_fill(GEN y, struct qfr_data *S)
    1244             : {
    1245           0 :   pari_sp av = avma;
    1246           0 :   GEN y2 = S->isqrtD;
    1247           0 :   gel(y,1) = gen_1;
    1248           0 :   if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
    1249           0 :   gel(y,2) = y2; av = avma;
    1250           0 :   gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(y2), S->D),-2));
    1251           0 : }
    1252             : static GEN
    1253           0 : qfr5_1(struct qfr_data *S, long prec)
    1254             : {
    1255           0 :   GEN y = cgetg(6, t_VEC);
    1256           0 :   qfr_1_fill(y, S);
    1257           0 :   gel(y,4) = gen_0;
    1258           0 :   gel(y,5) = real_1(prec); return y;
    1259             : }
    1260             : static GEN
    1261           0 : qfr3_1(struct qfr_data *S)
    1262             : {
    1263           0 :   GEN y = cgetg(4, t_VEC);
    1264           0 :   qfr_1_fill(y, S); return y;
    1265             : }
    1266             : 
    1267             : /* Assume D < 0 is the discriminant of a t_QFB */
    1268             : static GEN
    1269      720115 : qfi_1_by_disc(GEN D)
    1270             : {
    1271      720115 :   GEN b,c, y = cgetg(5,t_QFB);
    1272      720115 :   quadpoly_bc(D, mod2(D), &b,&c);
    1273      720115 :   if (b == gen_m1) b = gen_1;
    1274      720115 :   gel(y,1) = gen_1;
    1275      720115 :   gel(y,2) = b;
    1276      720115 :   gel(y,3) = c;
    1277      720115 :   gel(y,4) = icopy(D); return y;
    1278             : }
    1279             : static GEN
    1280      708157 : qfi_1(GEN x)
    1281             : {
    1282      708157 :   if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
    1283      708157 :   return qfi_1_by_disc(qfb_disc(x));
    1284             : }
    1285             : 
    1286             : GEN
    1287           0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
    1288             : 
    1289             : static GEN
    1290     9583809 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
    1291             : static GEN
    1292    25410012 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
    1293             : static GEN
    1294           7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
    1295             : static GEN
    1296           7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
    1297             : 
    1298             : static GEN
    1299           7 : qfipowraw(GEN x, long n)
    1300             : {
    1301           7 :   pari_sp av = avma;
    1302             :   GEN y;
    1303           7 :   if (!n) return qfi_1(x);
    1304           7 :   if (n== 1) return gcopy(x);
    1305           7 :   if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
    1306           7 :   if (n < 0) x = qfb_inv(x);
    1307           7 :   y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
    1308           7 :   return gerepilecopy(av,y);
    1309             : }
    1310             : 
    1311             : static GEN
    1312    11475733 : qfipow(GEN x, GEN n)
    1313             : {
    1314    11475733 :   pari_sp av = avma;
    1315             :   GEN y;
    1316    11475733 :   long s = signe(n);
    1317    11475733 :   if (!s) return qfi_1(x);
    1318    10767576 :   if (s < 0) x = qfb_inv(x);
    1319    10767574 :   y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
    1320    10767575 :   return gerepilecopy(av,y);
    1321             : }
    1322             : 
    1323             : static long
    1324      412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
    1325             : {
    1326             :   long z;
    1327      412328 :   *v = gen_0; *v2 = gen_1;
    1328     4351417 :   for (z=0; abscmpii(*v3,L) > 0; z++)
    1329             :   {
    1330     3939089 :     GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
    1331     3939089 :     *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
    1332             :   }
    1333      412328 :   return z;
    1334             : }
    1335             : 
    1336             : /* composition: Shanks' NUCOMP & NUDUPL */
    1337             : /* L = floor((|d|/4)^(1/4)) */
    1338             : GEN
    1339      400722 : nucomp(GEN x, GEN y, GEN L)
    1340             : {
    1341      400722 :   pari_sp av = avma;
    1342             :   long z;
    1343             :   GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
    1344             : 
    1345      400722 :   if (x==y) return nudupl(x,L);
    1346      400680 :   if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
    1347      400680 :   if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
    1348             : 
    1349      400680 :   if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
    1350      400680 :   s = shifti(addii(gel(x,2),gel(y,2)), -1);
    1351      400680 :   n = subii(gel(y,2), s);
    1352      400680 :   a1 = gel(x,1);
    1353      400680 :   a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
    1354      400680 :   if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
    1355      163576 :   else if (dvdii(s,d)) /* d | s */
    1356             :   {
    1357       83503 :     a = negi(mulii(u,n)); d1 = d;
    1358       83503 :     a1 = diviiexact(a1, d1);
    1359       83503 :     a2 = diviiexact(a2, d1);
    1360       83503 :     s = diviiexact(s, d1);
    1361             :   }
    1362             :   else
    1363             :   {
    1364             :     GEN p2, l;
    1365       80073 :     d1 = bezout(s,d,&u1,NULL);
    1366       80073 :     if (!equali1(d1))
    1367             :     {
    1368        2044 :       a1 = diviiexact(a1,d1);
    1369        2044 :       a2 = diviiexact(a2,d1);
    1370        2044 :       s = diviiexact(s,d1);
    1371        2044 :       d = diviiexact(d,d1);
    1372             :     }
    1373       80073 :     p1 = remii(gel(x,3),d);
    1374       80073 :     p2 = remii(gel(y,3),d);
    1375       80073 :     l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
    1376       80073 :     a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
    1377             :   }
    1378      400680 :   a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
    1379      400680 :   d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
    1380      400680 :   Q = cgetg(5,t_QFB);
    1381      400680 :   if (!z) {
    1382       37632 :     g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
    1383       37632 :     b = a2;
    1384       37632 :     b2 = gel(y,2);
    1385       37632 :     v2 = d1;
    1386       37632 :     gel(Q,1) = mulii(d,b);
    1387             :   } else {
    1388             :     GEN e, q3, q4;
    1389      363048 :     if (z&1) { v3 = negi(v3); v2 = negi(v2); }
    1390      363048 :     b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
    1391      363048 :     e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
    1392      363048 :     q3 = mulii(e,v2);
    1393      363048 :     q4 = subii(q3,s);
    1394      363048 :     b2 = addii(q3,q4);
    1395      363048 :     g = diviiexact(q4,v);
    1396      363048 :     if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
    1397      363048 :     gel(Q,1) = addii(mulii(d,b), mulii(e,v));
    1398             :   }
    1399      400680 :   q1 = mulii(b, v3);
    1400      400680 :   q2 = addii(q1,n);
    1401      400680 :   gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
    1402      400680 :   gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
    1403      400680 :   gel(Q,4) = gel(x,4);
    1404      400680 :   return qfbred_imag_av(av, Q);
    1405             : }
    1406             : 
    1407             : GEN
    1408       11648 : nudupl(GEN x, GEN L)
    1409             : {
    1410       11648 :   pari_sp av = avma;
    1411             :   long z;
    1412             :   GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
    1413             : 
    1414       11648 :   if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
    1415       11648 :   a = gel(x,1);
    1416       11648 :   b = gel(x,2);
    1417       11648 :   d1 = bezout(b,a, &u,NULL);
    1418       11648 :   if (!equali1(d1))
    1419             :   {
    1420        4620 :     a = diviiexact(a, d1);
    1421        4620 :     b = diviiexact(b, d1);
    1422             :   }
    1423       11648 :   c = modii(negi(mulii(u,gel(x,3))), a);
    1424       11648 :   p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
    1425       11648 :   d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
    1426       11648 :   a2 = sqri(d);
    1427       11648 :   c2 = sqri(v3);
    1428       11648 :   Q = cgetg(5,t_QFB);
    1429       11648 :   if (!z) {
    1430        1281 :     g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
    1431        1281 :     b2 = gel(x,2);
    1432        1281 :     v2 = d1;
    1433        1281 :     gel(Q,1) = a2;
    1434             :   } else {
    1435             :     GEN e;
    1436       10367 :     if (z&1) { v = negi(v); d = negi(d); }
    1437       10367 :     e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
    1438       10367 :     g = diviiexact(subii(mulii(e,v2), b), v);
    1439       10367 :     b2 = addii(mulii(e,v2), mulii(v,g));
    1440       10367 :     if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
    1441       10367 :     gel(Q,1) = addii(a2, mulii(e,v));
    1442             :   }
    1443       11648 :   gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
    1444       11648 :   gel(Q,3) = addii(c2, mulii(g,v2));
    1445       11648 :   gel(Q,4) = gel(x,4);
    1446       11648 :   return qfbred_imag_av(av, Q);
    1447             : }
    1448             : 
    1449             : static GEN
    1450        4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
    1451             : static GEN
    1452       11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
    1453             : GEN
    1454        1008 : nupow(GEN x, GEN n, GEN L)
    1455             : {
    1456             :   pari_sp av;
    1457             :   GEN y, D;
    1458             : 
    1459        1008 :   if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
    1460        1008 :   if (!is_qfi(x)) pari_err_TYPE("nupow",x);
    1461        1008 :   if (gequal1(n)) return gcopy(x);
    1462        1008 :   av = avma;
    1463        1008 :   D = qfb_disc(x);
    1464        1008 :   y = qfi_1_by_disc(D);
    1465        1008 :   if (!signe(n)) return y;
    1466         959 :   if (!L) L = sqrtnint(absi_shallow(D), 4);
    1467         959 :   y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
    1468         959 :   if (signe(n) < 0
    1469          35 :   && !absequalii(gel(y,1),gel(y,2))
    1470          35 :   && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
    1471         959 :   return gerepilecopy(av, y);
    1472             : }
    1473             : 
    1474             : /* Not stack-clean */
    1475             : GEN
    1476     1735167 : qfr5_compraw(GEN x, GEN y)
    1477             : {
    1478     1735167 :   GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
    1479     1735167 :   if (x == y)
    1480             :   {
    1481       34552 :     gel(z,4) = shifti(gel(x,4),1);
    1482       34552 :     gel(z,5) = sqrr(gel(x,5));
    1483             :   }
    1484             :   else
    1485             :   {
    1486     1700615 :     gel(z,4) = addii(gel(x,4),gel(y,4));
    1487     1700615 :     gel(z,5) = mulrr(gel(x,5),gel(y,5));
    1488             :   }
    1489     1735167 :   fix_expo(z); return z;
    1490             : }
    1491             : GEN
    1492     1735153 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
    1493     1735153 : { return qfr5_red(qfr5_compraw(x, y), S); }
    1494             : /* Not stack-clean */
    1495             : GEN
    1496     1009334 : qfr3_compraw(GEN x, GEN y)
    1497             : {
    1498     1009334 :   GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
    1499     1009334 :   return z;
    1500             : }
    1501             : GEN
    1502     1009334 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
    1503     1009334 : { return qfr3_red(qfr3_compraw(x,y), S); }
    1504             : 
    1505             : /* m > 0. Not stack-clean */
    1506             : static GEN
    1507           7 : qfr5_powraw(GEN x, long m)
    1508             : {
    1509           7 :   GEN y = NULL;
    1510          14 :   for (; m; m >>= 1)
    1511             :   {
    1512          14 :     if (m&1) y = y? qfr5_compraw(y,x): x;
    1513          14 :     if (m == 1) break;
    1514           7 :     x = qfr5_compraw(x,x);
    1515             :   }
    1516           7 :   return y;
    1517             : }
    1518             : 
    1519             : /* return x^n. Not stack-clean */
    1520             : GEN
    1521          21 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
    1522             : {
    1523          21 :   GEN y = NULL;
    1524          21 :   long i, m, s = signe(n);
    1525          21 :   if (!s) return qfr5_1(S, lg(gel(x,5)));
    1526          21 :   if (s < 0) x = qfb_inv(x);
    1527          42 :   for (i=lgefint(n)-1; i>1; i--)
    1528             :   {
    1529          21 :     m = n[i];
    1530          56 :     for (; m; m>>=1)
    1531             :     {
    1532          56 :       if (m&1) y = y? qfr5_comp(y,x,S): x;
    1533          56 :       if (m == 1 && i == 2) break;
    1534          35 :       x = qfr5_comp(x,x,S);
    1535             :     }
    1536             :   }
    1537          21 :   return y;
    1538             : }
    1539             : /* m > 0; return x^m. Not stack-clean */
    1540             : static GEN
    1541           0 : qfr3_powraw(GEN x, long m)
    1542             : {
    1543           0 :   GEN y = NULL;
    1544           0 :   for (; m; m>>=1)
    1545             :   {
    1546           0 :     if (m&1) y = y? qfr3_compraw(y,x): x;
    1547           0 :     if (m == 1) break;
    1548           0 :     x = qfr3_compraw(x,x);
    1549             :   }
    1550           0 :   return y;
    1551             : }
    1552             : /* return x^n. Not stack-clean */
    1553             : GEN
    1554        4557 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
    1555             : {
    1556        4557 :   GEN y = NULL;
    1557        4557 :   long i, m, s = signe(n);
    1558        4557 :   if (!s) return qfr3_1(S);
    1559        4557 :   if (s < 0) x = qfb_inv(x);
    1560        9130 :   for (i=lgefint(n)-1; i>1; i--)
    1561             :   {
    1562        4573 :     m = n[i];
    1563        5312 :     for (; m; m>>=1)
    1564             :     {
    1565        5292 :       if (m&1) y = y? qfr3_comp(y,x,S): x;
    1566        5292 :       if (m == 1 && i == 2) break;
    1567         739 :       x = qfr3_comp(x,x,S);
    1568             :     }
    1569             :   }
    1570        4557 :   return y;
    1571             : }
    1572             : 
    1573             : static GEN
    1574           7 : qfrinvraw(GEN x)
    1575             : {
    1576           7 :   if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
    1577           7 :  return qfbinv(x);
    1578             : }
    1579             : static GEN
    1580          14 : qfrpowraw(GEN x, long n)
    1581             : {
    1582          14 :   struct qfr_data S = { NULL, NULL, NULL };
    1583          14 :   pari_sp av = avma;
    1584          14 :   if (n==1) return gcopy(x);
    1585          14 :   if (n==-1) return qfrinvraw(x);
    1586           7 :   if (typ(x)==t_QFB)
    1587             :   {
    1588           0 :     GEN D = qfb_disc(x);
    1589           0 :     if (!n) return qfr_1(x);
    1590           0 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1591           0 :     x = qfr3_powraw(x, n);
    1592           0 :     x = qfr3_to_qfr(x, D);
    1593             :   }
    1594             :   else
    1595             :   {
    1596           7 :     GEN d0 = gel(x,2);
    1597           7 :     x = gel(x,1);
    1598           7 :     if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1599           7 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1600           7 :     x = qfr5_init(x, d0, &S);
    1601           7 :     if (labs(n) != 1) x = qfr5_powraw(x, n);
    1602           7 :     x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
    1603             :   }
    1604           7 :   return gerepilecopy(av, x);
    1605             : }
    1606             : static GEN
    1607         112 : qfrpow(GEN x, GEN n)
    1608             : {
    1609         112 :   struct qfr_data S = { NULL, NULL, NULL };
    1610         112 :   long s = signe(n);
    1611         112 :   pari_sp av = avma;
    1612         112 :   if (typ(x)==t_QFB)
    1613             :   {
    1614          42 :     if (!s) return qfr_1(x);
    1615          28 :     if (s < 0) x = qfb_inv(x);
    1616          28 :     x = qfr3_init(x, &S);
    1617          28 :     x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
    1618          28 :     x = qfr3_to_qfr(x, S.D);
    1619             :   }
    1620             :   else
    1621             :   {
    1622          70 :     GEN d0 = gel(x,2);
    1623          70 :     x = gel(x,1);
    1624          70 :     if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1625          49 :     if (s < 0) x = qfb_inv(x);
    1626          49 :     x = qfr5_init(x, d0, &S);
    1627          49 :     x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
    1628          49 :     x = qfr5_to_qfr(x, S.D, mulri(d0,n));
    1629             :   }
    1630          77 :   return gerepilecopy(av, x);
    1631             : }
    1632             : GEN
    1633          21 : qfbpowraw(GEN x, long n)
    1634             : {
    1635          21 :   GEN q = check_qfbext("qfbpowraw",x);
    1636          21 :   return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
    1637             : }
    1638             : /* same discriminant, no distance, no checks */
    1639             : GEN
    1640    10082572 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
    1641             : GEN
    1642     1393273 : qfbpow(GEN x, GEN n)
    1643             : {
    1644     1393273 :   GEN q = check_qfbext("qfbpow",x);
    1645     1393272 :   return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
    1646             : }
    1647             : GEN
    1648     1293725 : qfbpows(GEN x, long n)
    1649             : {
    1650     1293725 :   long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
    1651     1293725 :   affsi(n, N); return qfbpow(x, N);
    1652             : }
    1653             : 
    1654             : /* Prime forms attached to prime ideals of degree 1 */
    1655             : 
    1656             : /* assume x != 0 a t_INT, p > 0
    1657             :  * Return a t_QFB, but discriminant sign is not checked: can be used for
    1658             :  * real forms as well */
    1659             : GEN
    1660    12217259 : primeform_u(GEN x, ulong p)
    1661             : {
    1662    12217259 :   GEN c, y = cgetg(5, t_QFB);
    1663    12217018 :   pari_sp av = avma;
    1664             :   ulong b;
    1665             :   long s;
    1666             : 
    1667    12217018 :   s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
    1668             :   /* 2 or 3 mod 4 */
    1669    12216982 :   if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
    1670    12217191 :   if (p == 2) {
    1671     3881439 :     switch(s) {
    1672      589604 :       case 0: b = 0; break;
    1673     2939923 :       case 1: b = 1; break;
    1674      351914 :       case 4: b = 2; break;
    1675           0 :       default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1676           0 :                b = 0; /* -Wall */
    1677             :     }
    1678     3881441 :     c = shifti(subsi(s,x), -3);
    1679             :   } else {
    1680     8335752 :     b = Fl_sqrt(umodiu(x,p), p);
    1681     8337051 :     if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1682             :     /* mod(b) != mod2(x) ? */
    1683     8337913 :     if ((b ^ s) & 1) b = p - b;
    1684     8337913 :     c = diviuexact(shifti(subii(sqru(b), x), -2), p);
    1685             :   }
    1686    12204400 :   gel(y,3) = gerepileuptoint(av, c);
    1687    12212868 :   gel(y,4) = icopy(x);
    1688    12215536 :   gel(y,2) = utoi(b);
    1689    12215644 :   gel(y,1) = utoipos(p); return y;
    1690             : }
    1691             : 
    1692             : /* special case: p = 1 return unit form */
    1693             : GEN
    1694      135459 : primeform(GEN x, GEN p)
    1695             : {
    1696      135459 :   const char *f = "primeform";
    1697             :   pari_sp av;
    1698      135459 :   long s, sx = signe(x), sp = signe(p);
    1699             :   GEN y, b, absp;
    1700             : 
    1701      135459 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
    1702      135459 :   if (typ(p) != t_INT) pari_err_TYPE(f,p);
    1703      135459 :   if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
    1704      135459 :   if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
    1705      135459 :   if (lgefint(p) == 3)
    1706             :   {
    1707      135445 :     ulong pp = p[2];
    1708      135445 :     if (pp == 1) {
    1709       17782 :       if (sx < 0) {
    1710             :         long r;
    1711       10950 :         if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1712       10950 :         r = mod4(x);
    1713       10950 :         if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
    1714       10950 :         return qfi_1_by_disc(x);
    1715             :       }
    1716        6832 :       y = qfr_1_by_disc(x);
    1717        6832 :       if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
    1718        6832 :       return y;
    1719             :     }
    1720      117663 :     y = primeform_u(x, pp);
    1721      117656 :     if (sx < 0) {
    1722       89957 :       if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1723       89957 :       return y;
    1724             :     }
    1725       27699 :     if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
    1726       27699 :     return gcopy( qfr3_to_qfr(y, x) );
    1727             :   }
    1728          14 :   s = mod8(x);
    1729          14 :   if (sx < 0)
    1730             :   {
    1731           7 :     if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1732           7 :     if (s) s = 8-s;
    1733             :   }
    1734          14 :   y = cgetg(5, t_QFB);
    1735             :   /* 2 or 3 mod 4 */
    1736          14 :   if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
    1737          14 :   absp = absi_shallow(p); av = avma;
    1738          14 :   b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
    1739          14 :   s &= 1; /* s = x mod 2 */
    1740             :   /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
    1741          14 :   if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(absp,b));
    1742             : 
    1743          14 :   av = avma;
    1744          14 :   gel(y,3) = gerepileuptoint(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
    1745          14 :   gel(y,4) = icopy(x);
    1746          14 :   gel(y,2) = b;
    1747          14 :   gel(y,1) = icopy(p);
    1748          14 :   return y;
    1749             : }
    1750             : 
    1751             : static GEN
    1752     2620772 : normforms(GEN D, GEN fa)
    1753             : {
    1754             :   long i, j, k, lB, aN, sa;
    1755             :   GEN a, L, V, B, N, N2;
    1756     2620772 :   int D_odd = mpodd(D);
    1757     2620772 :   a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
    1758     2620772 :   sa = signe(a);
    1759     2620772 :   if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
    1760     1203972 :   V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
    1761     2551766 :            : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
    1762     2551766 :   if (!V) return NULL;
    1763      511966 :   N = gel(V,1); B = gel(V,2); lB = lg(B);
    1764      511966 :   N2 = shifti(N,1);
    1765      511966 :   aN = itou(diviiexact(a, N)); /* |a|/N */
    1766      511966 :   L = cgetg((lB-1)*aN+1, t_VEC);
    1767     2360564 :   for (k = 1, i = 1; i < lB; i++)
    1768             :   {
    1769     1848598 :     GEN b = shifti(gel(B,i), 1), c, C;
    1770     1848593 :     if (D_odd) b = addiu(b , 1);
    1771     1848593 :     c = diviiexact(shifti(subii(sqri(b), D), -2), a);
    1772     1848582 :     for (j = 0;; b = addii(b, N2))
    1773             :     {
    1774     2216656 :       gel(L, k++) = mkqfb(a, b, c, D);
    1775     2216664 :       if (++j == aN) break;
    1776      368069 :       C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
    1777      368074 :       c = sa > 0? addii(c, C): subii(c, C);
    1778             :     }
    1779             :   }
    1780      511966 :   return L;
    1781             : }
    1782             : 
    1783             : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
    1784             : static GEN
    1785      344315 : SL2_div_mul_e1(GEN N, GEN M)
    1786             : {
    1787      344315 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1788      344315 :   GEN A = mulii(gcoeff(N,1,1), d), B = mulii(gcoeff(N,1,2), b);
    1789      344307 :   GEN C = mulii(gcoeff(N,2,1), d), D = mulii(gcoeff(N,2,2), b);
    1790      344305 :   retmkvec2(subii(A,B), subii(C,D));
    1791             : }
    1792             : static GEN
    1793     1445665 : qfisolve_normform(GEN Q, GEN P)
    1794             : {
    1795     1445665 :   GEN a = gel(Q,1), N = gel(Q,2);
    1796     1445665 :   GEN M, b = qfbredsl2_imag_basecase(P, &M);
    1797     1445681 :   if (!qfb_equal(a,b)) return NULL;
    1798      102122 :   return SL2_div_mul_e1(N,M);
    1799             : }
    1800             : 
    1801             : /* Test equality modulo GL2 of two reduced forms */
    1802             : static int
    1803       61068 : GL2_qfb_equal(GEN a, GEN b)
    1804             : {
    1805       61068 :   return equalii(gel(a,1),gel(b,1))
    1806       11361 :    && absequalii(gel(a,2),gel(b,2))
    1807       72429 :    &&    equalii(gel(a,3),gel(b,3));
    1808             : }
    1809             : 
    1810             : /* Q(u,v) = p; if s < 0 return that solution; else the set of all solutions */
    1811             : static GEN
    1812       48019 : allsols(GEN Q, long s, GEN u, GEN v)
    1813             : {
    1814       48019 :   GEN w = mkvec2(u, v), b;
    1815       48010 :   if (signe(v) < 0) { u = negi(u); v = negi(v); } /* normalize for v >= 0 */
    1816       48010 :   w = mkvec2(u, v); if (s < 0) return w;
    1817       41364 :   if (!s) return mkvec(w);
    1818       38949 :   b = gel(Q,2); /* sum of the 2 solutions (if they exist) is -bv / a */
    1819       38949 :   if (signe(b))
    1820             :   { /* something to check */
    1821             :     GEN r, t;
    1822       13433 :     t = dvmdii(mulii(b, v), gel(Q,1), &r);
    1823       13433 :     if (signe(r)) return mkvec(w);
    1824        1820 :     u = addii(u, t);
    1825             :   }
    1826       27336 :   return mkvec2(w, mkvec2(negi(u), v));
    1827             : }
    1828             : static GEN
    1829      223046 : qfisolvep_all(GEN Q, GEN p, long all)
    1830             : {
    1831      223046 :   GEN R, U, V, M, N, x, q, D = qfb_disc(Q);
    1832      223045 :   long s = kronecker(D, p);
    1833             : 
    1834      223039 :   if (s < 0) return NULL;
    1835      126964 :   if (!all) s = -1; /* to indicate we want a single solution */
    1836             :   /* Solutions iff a class of maximal ideal above p is the class of Q;
    1837             :    * Two solutions iff (s > 0 and the class has order > 2), else one */
    1838      126964 :   if (!signe(gel(Q,2)))
    1839             :   { /* if principal form, use faster cornacchia */
    1840       43643 :     GEN a = gel(Q,1), c = gel(Q,3);
    1841       43643 :     if (equali1(a))
    1842             :     {
    1843       38166 :       if (!cornacchia(c, p, &M,&N)) return NULL;
    1844       33702 :       return allsols(Q, s, M, N);
    1845             :     }
    1846        5474 :     if (equali1(c))
    1847             :     {
    1848        5194 :       if (!cornacchia(a, p, &M,&N)) return NULL;
    1849         721 :       return allsols(Q, s, N, M);
    1850             :     }
    1851             :   }
    1852       83601 :   R = qfbredsl2_imag_basecase(Q, &U);
    1853       83601 :   if (equali1(gel(R,1)))
    1854             :   { /* principal form */
    1855       22533 :     if (!signe(gel(R,2)))
    1856             :     {
    1857        4396 :       if (!cornacchia(gel(R,3), p, &M,&N)) return NULL;
    1858         812 :       x = mkvec2(M,N);
    1859             :     }
    1860             :     else
    1861             :     { /* x^2 + xy + ((1-D)/4)y^2 = p <==> (2x + y)^2 - D y^2 = 4p */
    1862       18137 :       if (!cornacchia2(negi(D), p, &M, &N)) return NULL;
    1863        2331 :       x = subii(M,N); if (mpodd(x)) return NULL;
    1864        2331 :       x = mkvec2(shifti(x,-1), N);
    1865             :     }
    1866        3143 :     x = ZM_ZC_mul(U, x); x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1867        3143 :     return allsols(Q, s, gel(x,1), gel(x,2));
    1868             :   }
    1869       61068 :   q = qfbredsl2_imag_basecase(primeform(D, p), &V);
    1870       61068 :   if (!GL2_qfb_equal(R,q)) return NULL;
    1871       10451 :   if (signe(gel(R,2)) != signe(gel(q,2))) gcoeff(V,2,1) = negi(gcoeff(V,2,1));
    1872       10451 :   x = SL2_div_mul_e1(U,V); return allsols(Q, s, gel(x,1), gel(x,2));
    1873             : }
    1874             : GEN
    1875           0 : qfisolvep(GEN Q, GEN p)
    1876             : {
    1877           0 :   pari_sp av = avma;
    1878           0 :   GEN x = qfisolvep_all(Q, p, 0);
    1879           0 :   return x? gerepilecopy(av, x): gc_const(av, gen_0);
    1880             : }
    1881             : 
    1882             : static GEN
    1883      770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
    1884             : {
    1885      770126 :   pari_sp av = avma, btop;
    1886      770126 :   GEN M = N, P = qfbredsl2_real_basecase(Ps, rd), Q = P;
    1887             : 
    1888      770126 :   btop = avma;
    1889             :   for(;;)
    1890             :   {
    1891     5840681 :     if (qfb_equal(gel(M,1), gel(P,1)))
    1892      154084 :       return gerepileupto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
    1893     5686597 :     if (qfb_equal(gel(N,1), gel(Q,1)))
    1894       77658 :       return gerepileupto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
    1895     5608939 :     M = rhorealsl2(M, rd);
    1896     5608939 :     if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
    1897     5201735 :     Q = rhorealsl2(Q, rd);
    1898     5201735 :     if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
    1899     5070555 :     if (gc_needed(btop, 1)) gerepileall(btop, 2, &M, &Q);
    1900             :   }
    1901             : }
    1902             : 
    1903             : GEN
    1904           0 : qfrsolvep(GEN Q, GEN p)
    1905             : {
    1906           0 :   pari_sp av = avma;
    1907           0 :   GEN N, x, rd, d = qfb_disc(Q);
    1908             : 
    1909           0 :   if (kronecker(d, p) < 0) return gc_const(av, gen_0);
    1910           0 :   rd = sqrti(d);
    1911           0 :   N = qfbredsl2_real(Q, rd);
    1912           0 :   x = qfrsolve_normform(N, primeform(d, p), rd);
    1913           0 :   return x? gerepileupto(av, x): gc_const(av, gen_0);
    1914             : }
    1915             : 
    1916             : static GEN
    1917     1862936 : known_prime(GEN v)
    1918             : {
    1919     1862936 :   GEN p, e, fa = check_arith_all(v, "qfbsolve");
    1920     1862941 :   if (!fa) return BPSW_psp(v)? v: NULL;
    1921       42083 :   if (lg(gel(fa,1)) != 2) return NULL;
    1922       29357 :   p = gcoeff(fa,1,1);
    1923       29357 :   e = gcoeff(fa,1,2);
    1924       29357 :   return (equali1(e) && !is_pm1(p) && signe(p) > 0)? p: NULL;
    1925             : }
    1926             : static GEN
    1927     2215791 : qfsolve_normform(GEN Q, GEN f, GEN rd)
    1928     2215791 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
    1929             : static GEN
    1930     2843813 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
    1931             : {
    1932             :   GEN x, W, F, p;
    1933             :   long i, j, l;
    1934     2843813 :   if (!rd && (p = known_prime(fa))) return qfisolvep_all(Q, p, all);
    1935     2620760 :   F = normforms(qfb_disc(Q), fa);
    1936     2620772 :   if (!F) return NULL;
    1937      511966 :   if (!*Qr) *Qr = qfbredsl2(Q, rd);
    1938      511966 :   l = lg(F); W = all? cgetg(l, t_VEC): NULL;
    1939     2727245 :   for (j = i = 1; i < l; i++)
    1940     2215791 :     if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
    1941             :     {
    1942      333855 :       if (!all) return x;
    1943      333344 :       gel(W,j++) = x;
    1944             :     }
    1945      511454 :   if (j == 1) return NULL;
    1946      127455 :   setlg(W,j); return lexsort(W);
    1947             : }
    1948             : 
    1949             : static GEN
    1950     2838513 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
    1951             : static GEN
    1952     2828288 : qfbsolve_primitive(GEN Q, GEN fa, long all)
    1953             : {
    1954     2828288 :   GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
    1955     2828285 :   x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
    1956     2828286 :   if (!x) return cgetg(1, t_VEC);
    1957      174747 :   return x;
    1958             : }
    1959             : 
    1960             : /* f / g^2 */
    1961             : static GEN
    1962        5299 : famat_divsqr(GEN f, GEN g)
    1963        5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
    1964             : static GEN
    1965       10227 : qfbsolve_all(GEN Q, GEN n, long all)
    1966             : {
    1967       10227 :   GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
    1968       10227 :   GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
    1969       10227 :   long i, j, l = lg(D);
    1970       10227 :   W = all? cgetg(l, t_VEC): NULL;
    1971       25151 :   for (i = j = 1; i < l; i++)
    1972             :   {
    1973       15526 :     GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
    1974       15526 :     if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
    1975             :     {
    1976        1218 :       if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
    1977        1218 :       if (!all) return w;
    1978         616 :       gel(W,j++) = w;
    1979             :     }
    1980             :   }
    1981        9625 :   if (j == 1) return cgetg(1, t_VEC);
    1982         525 :   setlg(W,j); return lexsort(shallowconcat1(W));
    1983             : }
    1984             : 
    1985             : GEN
    1986     2838519 : qfbsolve(GEN Q, GEN n, long fl)
    1987             : {
    1988     2838519 :   pari_sp av = avma;
    1989     2838519 :   if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
    1990     2838519 :   if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
    1991     5666798 :   return gerepilecopy(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
    1992     2828285 :                                   : qfbsolve_primitive(Q, n, fl & 1));
    1993             : }
    1994             : 
    1995             : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
    1996             :  * Assume d > 0 and p is prime */
    1997             : long
    1998       55242 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
    1999             : {
    2000       55242 :   pari_sp av = avma;
    2001             :   GEN b, c, r;
    2002             : 
    2003       55242 :   *px = *py = gen_0;
    2004       55242 :   b = subii(p, d);
    2005       55190 :   if (signe(b) < 0) return gc_long(av,0);
    2006       54980 :   if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
    2007       54973 :   b = Fp_sqrt(b, p); /* sqrt(-d) */
    2008       55037 :   if (!b) return gc_long(av,0);
    2009       51306 :   b = gmael(halfgcdii(p, b), 2, 2);
    2010       51339 :   c = dvmdii(subii(p, sqri(b)), d, &r);
    2011       51209 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    2012       35488 :   set_avma(av);
    2013       35489 :   *px = icopy(b);
    2014       35472 :   *py = icopy(c); return 1;
    2015             : }
    2016             : 
    2017             : static GEN
    2018     1128546 : lastqi(GEN Q)
    2019             : {
    2020     1128546 :   GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
    2021     1128535 :   if (!signe(q)) return gen_0;
    2022     1128346 :   if (!signe(s)) return p;
    2023     1122522 :   if (is_pm1(q)) return subiu(p,1);
    2024     1122525 :   return divii(p, absi_shallow(q));
    2025             : }
    2026             : 
    2027             : static long
    2028     1128536 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
    2029             : {
    2030             :   GEN M, Q, V, c, r, b2;
    2031     1128536 :   if (!signe(b)) { /* d = p,2p,3p,4p */
    2032          14 :     set_avma(av);
    2033          14 :     if (absequalii(d, px4)){ *py = gen_1; return 1; }
    2034          14 :     if (absequalii(d, p))  { *py = gen_2; return 1; }
    2035           0 :     return 0;
    2036             :   }
    2037     1128522 :   if (mod2(b) != mod2(d)) b = subii(p,b);
    2038     1128512 :   M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
    2039     1128546 :   b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
    2040     1128466 :   b2 = sqri(b);
    2041     1128464 :   if (cmpii(b2,px4) > 0)
    2042             :   {
    2043     1120750 :     b = gel(V,1); b2 = sqri(b);
    2044     1120763 :     if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
    2045             :   }
    2046     1128465 :   c = dvmdii(subii(px4, b2), d, &r);
    2047     1128474 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    2048     1083727 :   set_avma(av);
    2049     1083726 :   *px = icopy(b);
    2050     1083722 :   *py = icopy(c); return 1;
    2051             : }
    2052             : 
    2053             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
    2054             :  * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
    2055             : long
    2056     1088737 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
    2057             : {
    2058     1088737 :   pari_sp av = avma;
    2059     1088737 :   GEN b, p4 = shifti(p,2);
    2060             : 
    2061     1088718 :   *px = *py = gen_0;
    2062     1088718 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2063     1087894 :   if (absequaliu(p, 2))
    2064             :   {
    2065           7 :     set_avma(av);
    2066           7 :     switch (itou_or_0(d)) {
    2067           0 :       case 4: *px = gen_2; break;
    2068           0 :       case 7: *px = gen_1; break;
    2069           7 :       default: return 0;
    2070             :     }
    2071           0 :     *py = gen_1; return 1;
    2072             :   }
    2073     1087888 :   b = Fp_sqrt(negi(d), p);
    2074     1087942 :   if (!b) return gc_long(av,0);
    2075     1087858 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2076             : }
    2077             : 
    2078             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
    2079             : long
    2080       40677 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
    2081             : {
    2082       40677 :   pari_sp av = avma;
    2083       40677 :   GEN p4 = shifti(p,2);
    2084       40677 :   *px = *py = gen_0;
    2085       40677 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2086       40677 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2087             : }
    2088             : 
    2089             : GEN
    2090        7630 : qfbcornacchia(GEN d, GEN p)
    2091             : {
    2092        7630 :   pari_sp av = avma;
    2093             :   GEN x, y;
    2094        7630 :   if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
    2095        7630 :   if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
    2096        7630 :   if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
    2097         287 :     return gerepilecopy(av, mkvec2(x, y));
    2098        7343 :   set_avma(av); return cgetg(1, t_VEC);
    2099             : }

Generated by: LCOV version 1.14