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.1 lcov report (development 28904-c3aa21e911) Lines: 1103 1225 90.0 %
Date: 2023-12-04 07:51:13 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      151510 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
      24             : {
      25             :   long r;
      26      151510 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
      27      151496 :   *s = signe(x);
      28      151496 :   if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
      29      151498 :   r = mod4(x); if (*s < 0 && r) r = 4 - r;
      30      151499 :   if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
      31      151485 :   if (pr) *pr = r;
      32      151485 : }
      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     1806960 : check_qfbext(const char *fun, GEN x)
     113             : {
     114     1806960 :   long t = typ(x);
     115     1806960 :   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      144361 : qfb3(GEN x, GEN y, GEN z)
     127      144361 : { retmkqfb(icopy(x), icopy(y), icopy(z), qfb_disc3(x,y,z)); }
     128             : 
     129             : static int
     130    23783631 : qfb_equal(GEN x, GEN y)
     131             : {
     132    23783631 :   return equalii(gel(x,1),gel(y,1))
     133     1592909 :       && equalii(gel(x,2),gel(y,2))
     134    25376539 :       && equalii(gel(x,3),gel(y,3));
     135             : }
     136             : 
     137             : /* valid for t_QFB, qfr3, qfr5; shallow */
     138             : static GEN
     139      878188 : qfb_inv(GEN x) {
     140      878188 :   GEN z = shallowcopy(x);
     141      878191 :   gel(z,2) = negi(gel(z,2));
     142      878188 :   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       77217 : Qfb0(GEN a, GEN b, GEN c)
     151             : {
     152             :   GEN q, D;
     153       77217 :   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       77175 :   else if (!c)
     169           7 :     pari_err_TYPE("Qfb",b);
     170       77196 :   if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
     171       77189 :   if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
     172       77189 :   if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
     173       77189 :   q = qfb3(a, b, c); D = qfb_disc(q);
     174       77189 :   if (signe(D) < 0)
     175       42385 :   { if (signe(a) < 0) pari_err_IMPL("negative definite t_QFB"); }
     176       34804 :   else if (Z_issquare(D)) pari_err_DOMAIN("Qfb","issquare(disc)","=", gen_1,q);
     177       77182 :   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    16146573 : dvmdii_round(GEN b, GEN a, GEN *r)
     189             : {
     190    16146573 :   GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
     191    16146540 :   if (signe(b) >= 0) {
     192     8887883 :     if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
     193             :   } else { /* r <= 0 */
     194     7258657 :     if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
     195             :   }
     196    16146506 :   return q;
     197             : }
     198             : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
     199             : static long
     200    97107912 : dvmdsu_round(long b, ulong a, long *r)
     201             : {
     202    97107912 :   ulong a2 = a << 1, q, ub, ur;
     203    97107912 :   if (b >= 0) {
     204    61928683 :     ub = b;
     205    61928683 :     q = ub / a2;
     206    61928683 :     ur = ub % a2;
     207    61928683 :     if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
     208    21805584 :     else *r = (long)ur;
     209    61928683 :     return (long)q;
     210             :   } else { /* r <= 0 */
     211    35179229 :     ub = (ulong)-b; /* |b| */
     212    35179229 :     q = ub / a2;
     213    35179229 :     ur = ub % a2;
     214    35179229 :     if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
     215    19235776 :     else *r = -(long)ur;
     216    35179229 :     return -(long)q;
     217             :   }
     218             : }
     219             : /* reduce b mod 2*a. Update b,c */
     220             : static void
     221     2768839 : REDB(GEN a, GEN *b, GEN *c)
     222             : {
     223     2768839 :   GEN r, q = dvmdii_round(*b, a, &r);
     224     2768839 :   if (!signe(q)) return;
     225     2699909 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     226     2699909 :   *b = r;
     227             : }
     228             : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
     229             : static void
     230    97110107 : sREDB(ulong a, long *b, ulong *c)
     231             : {
     232             :   long r, q;
     233             :   ulong uz;
     234   105325724 :   if (a > LONG_MAX) return; /* b already reduced */
     235    97110107 :   q = dvmdsu_round(*b, a, &r);
     236    97126932 :   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    88911315 :   if (*b < 0)
     240             :   { /* uz = -z >= 0, q < 0 */
     241    31079125 :     if (r >= 0) /* different signs=>no overflow, exact division */
     242    15997768 :       uz = (ulong)-((*b + r)>>1);
     243             :     else
     244             :     {
     245    15081357 :       ulong ub = (ulong)-*b, ur = (ulong)-r;
     246    15081357 :       uz = (ub + ur) >> 1;
     247             :     }
     248    31079125 :     *c -= (-q) * uz; /* c -= qz */
     249             :   }
     250             :   else
     251             :   { /* uz = z >= 0, q > 0 */
     252    57832190 :     if (r <= 0)
     253    40182088 :       uz = (*b + r)>>1;
     254             :     else
     255             :     {
     256    17650102 :       ulong ub = (ulong)*b, ur = (ulong)r;
     257    17650102 :       uz = ((ub + ur) >> 1);
     258             :     }
     259    57832190 :     *c -= q * uz; /* c -= qz */
     260             :   }
     261    88911315 :   *b = r;
     262             : }
     263             : static void
     264    13377734 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
     265             : { /* REDB(a,b,c) */
     266    13377734 :   GEN r, q = dvmdii_round(*b, a, &r);
     267    13377666 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     268    13377646 :   *b = r;
     269    13377646 :   *u2 = subii(*u2, mulii(q, u1));
     270    13377656 : }
     271             : 
     272             : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
     273             : static GEN
     274     6616345 : qfbredsl2_imag_basecase(GEN q, GEN *U)
     275             : {
     276     6616345 :   pari_sp av = avma;
     277             :   GEN z, u1,u2,v1,v2,Q;
     278     6616345 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     279             :   long cmp;
     280     6616345 :   u1 = gen_1; u2 = gen_0;
     281     6616345 :   cmp = abscmpii(a, b);
     282     6616347 :   if (cmp < 0)
     283     2083510 :     REDBU(a,&b,&c, u1,&u2);
     284     4532837 :   else if (cmp == 0 && signe(b) < 0)
     285             :   { /* b = -a */
     286       11904 :     b = negi(b);
     287       11904 :     u2 = gen_1;
     288             :   }
     289             :   for(;;)
     290             :   {
     291    17910523 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     292    11294102 :     swap(a,c); b = negi(b);
     293    11294214 :     z = u1; u1 = u2; u2 = negi(z);
     294    11294227 :     REDBU(a,&b,&c, u1,&u2);
     295    11294177 :     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     6616331 :   if (cmp == 0 && signe(b) < 0)
     301             :   {
     302       19116 :     b = negi(b);
     303       19116 :     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     6616331 :   z = shifti(subii(b, gel(q,2)), -1);
     308     6616311 :   v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
     309     6616320 :   z = subii(z, b);
     310     6616306 :   v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
     311     6616316 :   *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
     312     6616350 :   Q = mkqfb(a,b,c,gel(q,4));
     313     6616345 :   return gc_all(av, 2, &Q, U);
     314             : }
     315             : 
     316             : static GEN
     317     1068563 : setq_b0(ulong a, ulong c, GEN D)
     318     1068563 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
     319             : /* assume |sb| = 1 */
     320             : static GEN
     321    80263514 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
     322    80263514 : { 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      958983 : qfbred_imag_1_b0(ulong a, ulong c, GEN D)
     326      958983 : { 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    81442043 : 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      141105 :   { /* at most twice */
     336    81442043 :     long lb = lgefint(b); /* <= 3 after first loop */
     337    81442043 :     if (lb == 2) return qfbred_imag_1_b0(a[2],c[2], D);
     338    80483060 :     if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
     339      139160 :     REDB(a,&b,&c);
     340      140144 :     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      140144 :     swap(a,c); b = negi(b);
     349             :   }
     350             :   /* b != 0 */
     351    80346000 :   set_avma(av);
     352    80354251 :   ua = a[2];
     353    80354251 :   ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
     354    80354251 :   uc = c[2];
     355    80354251 :   if (ua < ub)
     356    29135675 :     sREDB(ua, &sb, &uc);
     357    51218576 :   else if (ua == ub && sb < 0) sb = (long)ub;
     358   148390223 :   while(ua > uc)
     359             :   {
     360    68014275 :     lswap(ua,uc); sb = -sb;
     361    68014275 :     sREDB(ua, &sb, &uc);
     362             :   }
     363    80375948 :   if (!sb) return setq_b0(ua, uc, D);
     364             :   else
     365             :   {
     366    80266368 :     long s = 1;
     367    80266368 :     if (sb < 0)
     368             :     {
     369    30107317 :       sb = -sb;
     370    30107317 :       if (ua != uc) s = -1;
     371             :     }
     372    80266368 :     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    12022731 : fix_expo(GEN x)
     415             : {
     416    12022731 :   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    12022731 : }
     421             : 
     422             : /* (1/2) log (d * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
     423             : GEN
     424      184254 : qfr5_dist(GEN e, GEN d, long prec)
     425             : {
     426      184254 :   GEN t = logr_abs(d);
     427      184254 :   if (signe(e)) {
     428           0 :     GEN u = mulir(e, mplog2(prec));
     429           0 :     shiftr_inplace(u, EMAX); t = addrr(t, u);
     430             :   }
     431      184254 :   shiftr_inplace(t, -1); return t;
     432             : }
     433             : 
     434             : static void
     435    17290519 : rho_get_BC(GEN *B, GEN *C, GEN a, GEN b, GEN c, struct qfr_data *S)
     436             : {
     437             :   GEN t, u, q;
     438    17290519 :   u = shifti(c,1);
     439    17290519 :   t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: c;
     440    17290519 :   q = dvmdii(addii_sign(t,1, b,signe(b)), u, &u);
     441    17290519 :   *B = addii_sign(t, 1, u, -signe(u)); /* |t| - (|t|+b) % |2c| */
     442    17290519 :   *C = subii(a, mulii(q, subii(b, mulii(q,c))));
     443    17290519 : }
     444             : /* Not stack-clean */
     445             : GEN
     446     1136471 : qfr3_rho(GEN x, struct qfr_data *S)
     447             : {
     448     1136471 :   GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3);
     449     1136471 :   rho_get_BC(&B, &C, a, b, c, S);
     450     1136471 :   return mkvec3(c, B, C);
     451             : }
     452             : 
     453             : /* Not stack-clean */
     454             : GEN
     455    10296041 : qfr5_rho(GEN x, struct qfr_data *S)
     456             : {
     457    10296041 :   GEN B, C, a = gel(x,1), b = gel(x,2), c = gel(x,3), y;
     458    10296041 :   long sb = signe(b);
     459    10296041 :   rho_get_BC(&B, &C, a, b, c, S);
     460    10296041 :   y = mkvec5(c, B, C, gel(x,4), gel(x,5));
     461    10296041 :   if (sb) {
     462    10292324 :     GEN t = subii(sqri(b), S->D);
     463    10292324 :     if (sb < 0)
     464     2299122 :       t = divir(t, sqrr(subir(b,S->sqrtD)));
     465             :     else
     466     7993202 :       t = divri(sqrr(addir(b,S->sqrtD)), t);
     467             :     /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
     468    10292324 :     gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
     469        3717 :   } else gel(y,5) = negr(gel(y,5));
     470    10296041 :   return y;
     471             : }
     472             : 
     473             : /* Not stack-clean */
     474             : GEN
     475      217084 : qfr_to_qfr5(GEN x, long prec)
     476      217084 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
     477             : 
     478             : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
     479             : GEN
     480         476 : qfr5_to_qfr(GEN x, GEN D, GEN d0)
     481             : {
     482         476 :   if (d0)
     483             :   {
     484         140 :     GEN n = gel(x,4), d = absr(gel(x,5));
     485         140 :     if (signe(n))
     486             :     {
     487           0 :       n = addis(shifti(n, EMAX), expo(d));
     488           0 :       setexpo(d, 0); d = logr_abs(d);
     489           0 :       if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
     490           0 :       shiftr_inplace(d, -1);
     491           0 :       d0 = addrr(d0, d);
     492             :     }
     493         140 :     else if (!gequal1(d)) /* avoid loss of precision */
     494             :     {
     495          91 :       d = logr_abs(d);
     496          91 :       shiftr_inplace(d, -1);
     497          91 :       d0 = addrr(d0, d);
     498             :     }
     499             :   }
     500         476 :   x = qfr3_to_qfr(x, D);
     501         476 :   return d0 ? mkvec2(x,d0): x;
     502             : }
     503             : 
     504             : /* Not stack-clean */
     505             : GEN
     506       31913 : qfr3_to_qfr(GEN x, GEN d) { retmkqfb(gel(x,1), gel(x,2), gel(x,3), d); }
     507             : 
     508             : static int
     509    20824028 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
     510             : {
     511             :   GEN t;
     512    20824028 :   if (signe(b) <= 0 || abscmpii(b, isqrtD) > 0) return 0;
     513     5505679 :   t = addii_sign(isqrtD,1, shifti(a,1),-1); /* floor(sqrt(D)) - |2a| */
     514     1312351 :   return signe(t) < 0 ? abscmpii(b, t) >= 0
     515     6818030 :                       : abscmpii(b, t) > 0;
     516             : }
     517             : 
     518             : /* Not stack-clean */
     519             : GEN
     520     1947428 : qfr5_red(GEN x, struct qfr_data *S)
     521             : {
     522     1947428 :   pari_sp av = avma;
     523    10247874 :   while (!ab_isreduced(gel(x,1), gel(x,2), S->isqrtD))
     524             :   {
     525     8300446 :     x = qfr5_rho(x, S);
     526     8300446 :     if (gc_needed(av,2))
     527             :     {
     528           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
     529           0 :       x = gerepilecopy(av, x);
     530             :     }
     531             :   }
     532     1947428 :   return x;
     533             : }
     534             : /* Not stack-clean */
     535             : GEN
     536     1169914 : qfr3_red(GEN x, struct qfr_data *S)
     537             : {
     538     1169914 :   pari_sp av = avma;
     539     1169914 :   GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
     540     7027921 :   while (!ab_isreduced(a, b, S->isqrtD))
     541             :   {
     542             :     GEN B, C;
     543     5858007 :     rho_get_BC(&B, &C, a, b, c, S);
     544     5858007 :     a = c; b = B; c = C;
     545     5858007 :     if (gc_needed(av,2))
     546             :     {
     547           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
     548           0 :       gerepileall(av, 3, &a, &b, &c);
     549             :     }
     550             :   }
     551     1169914 :   return mkvec3(a, b, c);
     552             : }
     553             : 
     554             : void
     555        2170 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
     556             : {
     557        2170 :   S->D = D;
     558        2170 :   S->sqrtD = sqrtr(itor(S->D,prec));
     559        2170 :   S->isqrtD = truncr(S->sqrtD);
     560        2170 : }
     561             : 
     562             : static GEN
     563         140 : qfr5_init(GEN x, GEN d, struct qfr_data *S)
     564             : {
     565         140 :   long prec = realprec(d), l = -expo(d);
     566         140 :   if (l < BITS_IN_LONG) l = BITS_IN_LONG;
     567         140 :   prec = maxss(prec, nbits2prec(l));
     568         140 :   S->D = qfb_disc(x);
     569         140 :   x = qfr_to_qfr5(x,prec);
     570         140 :   if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
     571           0 :   else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
     572             : 
     573         140 :   if (!S->isqrtD)
     574             :   {
     575         126 :     pari_sp av=avma;
     576             :     long e;
     577         126 :     S->isqrtD = gcvtoi(S->sqrtD,&e);
     578         126 :     if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
     579             :   }
     580          14 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
     581         140 :   return x;
     582             : }
     583             : static GEN
     584         364 : qfr3_init(GEN x, struct qfr_data *S)
     585             : {
     586         364 :   S->D = qfb_disc(x);
     587         364 :   if (!S->isqrtD) S->isqrtD = sqrti(S->D);
     588         280 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
     589         364 :   return x;
     590             : }
     591             : 
     592             : #define qf_NOD  2
     593             : #define qf_STEP 1
     594             : 
     595             : static GEN
     596         420 : qfbred_real_basecase_i(GEN x, long flag, GEN isqrtD, GEN sqrtD)
     597             : {
     598             :   struct qfr_data S;
     599         420 :   GEN d = NULL, y;
     600         420 :   if (typ(x)==t_VEC) { d = gel(x,2); x = gel(x,1); } else flag |= qf_NOD;
     601         420 :   S.sqrtD = sqrtD;
     602         420 :   S.isqrtD = isqrtD;
     603         420 :   y = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, d, &S);
     604         420 :   switch(flag) {
     605          63 :     case 0:              y = qfr5_red(y,&S); break;
     606         329 :     case qf_NOD:         y = qfr3_red(y,&S); break;
     607          21 :     case qf_STEP:        y = qfr5_rho(y,&S); break;
     608           7 :     case qf_STEP|qf_NOD: y = qfr3_rho(y,&S); break;
     609           0 :     default: pari_err_FLAG("qfbred");
     610             :   }
     611         420 :   return qfr5_to_qfr(y, qfb_disc(x), d);
     612             : }
     613             : 
     614             : static void
     615    13379256 : _rhorealsl2(GEN *pa, GEN *pb, GEN *pc, GEN *pu1, GEN *pu2, GEN *pv1,
     616             :             GEN *pv2, GEN rd)
     617             : {
     618    13379256 :   GEN C = mpabs_shallow(*pc), t = addii(*pb, gmax_shallow(rd,C));
     619    13379256 :   GEN r, q = truedvmdii(t, shifti(C,1), &r);
     620    13379256 :   GEN a = *pa, b= *pb, c = *pc;
     621    13379256 :   if (signe(c) < 0) togglesign(q);
     622    13379256 :   *pa = *pc;
     623    13379256 :   *pb = subii(t, addii(r, *pb));
     624    13379256 :   *pc = subii(a,mulii(q,subii(b, mulii(q,c))));
     625    13379256 :   r = *pu1; *pu1 = *pv1; *pv1 = subii(mulii(q, *pv1), r);
     626    13379256 :   r = *pu2; *pu2 = *pv2; *pv2 = subii(mulii(q, *pv2), r);
     627    13379256 : }
     628             : 
     629             : static GEN
     630    10810674 : rhorealsl2(GEN A, GEN rd)
     631             : {
     632    10810674 :   GEN V = gel(A,1), M = gel(A,2);
     633    10810674 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
     634    10810674 :   GEN u1 = gcoeff(M,1,1), v1 = gcoeff(M,1,2);
     635    10810674 :   GEN u2 = gcoeff(M,2,1), v2 = gcoeff(M,2,2);
     636    10810674 :   _rhorealsl2(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
     637    10810674 :   return mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2));
     638             : }
     639             : 
     640             : static GEN
     641      979651 : qfbredsl2_real_basecase(GEN V, GEN rd)
     642             : {
     643      979651 :   pari_sp av = avma;
     644      979651 :   GEN u1 = gen_1, u2 = gen_0, v1 = gen_0, v2 = gen_1;
     645      979651 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
     646     3548233 :   while (!ab_isreduced(a,b,rd))
     647             :   {
     648     2568582 :     _rhorealsl2(&a,&b,&c, &u1,&u2,&v1,&v2, rd);
     649     2568582 :     if (gc_needed(av, 1))
     650             :     {
     651           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfbredsl2");
     652           0 :       gerepileall(av, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
     653             :     }
     654             :   }
     655      979651 :   return gerepilecopy(av, mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2)));
     656             : }
     657             : 
     658             : /* fast reduction of qfb with positive coefficients, based on
     659             : Arnold Schoenhage, Fast reduction and composition of binary quadratic forms,
     660             : Proc. of Intern. Symp. on Symbolic and Algebraic Computation (Bonn) (S. M.
     661             : Watt, ed.), ACM Press, 1991, pp. 128-133.
     662             : <https://dl.acm.org/doi/pdf/10.1145/120694.120711>
     663             : With thanks to Keegan Ryan
     664             : BA20230927
     665             : */
     666             : 
     667             : #define QFBRED_LIMIT 9000
     668             : 
     669             : /* pqfb: qf with positive coefficients */
     670             : 
     671             : static int
     672     5314748 : lti2n(GEN a, long m)
     673             : {
     674     5314748 :   return signe(a) < 0 || expi(a) < m ? 1 : 0;
     675             : }
     676             : 
     677             : static GEN
     678     2065206 : pqfbred_1(GEN Q, long m, GEN U)
     679             : {
     680     2065206 :   GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     681     2065206 :   if (abscmpii(a, c) < 0)
     682             :   {
     683             :     GEN t, at, r;
     684     1032327 :     GEN r2 = addii(shifti(a, m + 2), d);
     685     1032327 :     long e2 = expi(r2);
     686     1032327 :     r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
     687     1032327 :     t = truedivii(subii(b, r), shifti(a,1));
     688     1032327 :     if (signe(t)==0) pari_err_BUG("pqfbred_1");
     689     1032327 :     at = mulii(a,t);
     690     1032327 :     c = addii(subii(c, mulii(b, t)), mulii(at, t));
     691     1032327 :     b = subii(b, shifti(at,1));
     692     1032327 :     gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,1,1), t));
     693     1032327 :     gcoeff(U,2,2) = subii( gcoeff(U,2,2), mulii(gcoeff(U,2,1), t));
     694             :   } else
     695             :   {
     696             :     GEN t, ct, r;
     697     1032879 :     GEN r2 = addii(shifti(c, m + 2), d);
     698     1032879 :     long e2 = expi(r2);
     699     1032879 :     r = int2n(signe(r2) < 0 || e2 < 2*m+2 ? m+1 : e2>>1);
     700     1032879 :     t = truedivii(subii(b, r), shifti(c,1));
     701     1032879 :     if (signe(t)==0) pari_err_BUG("pqfbred_1");
     702     1032879 :     ct = mulii(c, t);
     703     1032879 :     a = addii(subii(a, mulii(b, t)), mulii(ct, t));
     704     1032879 :     b = subii(b, shifti(ct, 1));
     705     1032879 :     gcoeff(U,1,1) = subii(gcoeff(U,1,1), mulii(gcoeff(U,1,2), t));
     706     1032879 :     gcoeff(U,2,1) = subii(gcoeff(U,2,1), mulii(gcoeff(U,2,2), t));
     707             :   }
     708     2065206 :   return mkqfb(a,b,c,d);
     709             : }
     710             : 
     711             : static int
     712     2201585 : is_minimal(GEN Q, long m)
     713             : {
     714     2201585 :   pari_sp av = avma;
     715     2201585 :   GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3);
     716     5314748 :   return gc_bool(av, lti2n(addii(subii(a,b), c), m)
     717     2072745 :                  || (lti2n(subii(b, shifti(a,1)), m+1)
     718     1040418 :                      && lti2n(subii(b, shifti(c,1)), m+1)));
     719             : }
     720             : 
     721             : static GEN
     722      134615 : pqfbred_iter_1(GEN Q, ulong m, GEN U)
     723             : {
     724      134615 :   pari_sp av = avma;
     725     2068245 :   while (!is_minimal(Q,m))
     726             :   {
     727     1933630 :     Q = pqfbred_1(Q, m, U);
     728     1933630 :     if (gc_needed(av, 1))
     729             :     {
     730           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"pqfbred_iter_1, lc = %ld", expi(gel(Q,3)));
     731           0 :       gerepileall(av, 3, &Q, &gel(U,1), &gel(U,2));
     732             :     }
     733             :   }
     734      134615 :   return Q;
     735             : }
     736             : 
     737             : static GEN
     738       66501 : pqfbred_basecase(GEN Q, ulong m, GEN *pt_U)
     739             : {
     740       66501 :   pari_sp av = avma;
     741       66501 :   GEN  U = matid(2);
     742       66501 :   Q = pqfbred_iter_1(Q, m, U);
     743       66501 :   *pt_U = U;
     744       66501 :   return gc_all(av, 2, &Q, pt_U);
     745             : }
     746             : 
     747             : static long
     748    86874095 : qfb_maxexpi(GEN Q)
     749    86874095 : { return 1+maxss(expi(gel(Q,1)), maxss(expi(gel(Q,2)), expi(gel(Q,3)))); }
     750             : 
     751             : static long
     752      136228 : qfb_minexpi(GEN Q)
     753             : {
     754      136228 :   long m =  minss(expi(gel(Q,1)), minss(expi(gel(Q,2)), expi(gel(Q,3))));
     755      136228 :   return m < 0 ? 0: m;
     756             : }
     757             : 
     758             : static GEN
     759      134615 : pqfbred_rec(GEN Q, long m, GEN *pt_U)
     760             : {
     761      134615 :   pari_sp av = avma;
     762             :   GEN Q0, Q1, QR;
     763      134615 :   GEN d = qfb_disc(Q);
     764      134615 :   long n = qfb_maxexpi(Q) - m;
     765             :   long h;
     766      134615 :   int going_to_r8 = 0;
     767             :   GEN U;
     768             : 
     769      134615 :   if (n < 170)
     770       66501 :     return pqfbred_basecase(Q, m, pt_U);
     771             : 
     772       68114 :   if (qfb_minexpi(Q) <= m + 2)
     773             :   {
     774           0 :     U = matid(2);
     775           0 :     QR = Q;
     776             :   }
     777             :   else
     778             :   {
     779             :     long p, mm;
     780       68114 :     if (m <= n)
     781             :     {
     782         942 :       mm = m;
     783         942 :       p = 0;
     784         942 :       Q1 = Q;
     785             :     } else
     786             :     {
     787       67172 :       mm = n;
     788       67172 :       p = m + 1 - n;
     789       67172 :       Q0 = mkvec3(remi2n(gel(Q,1), p), remi2n(gel(Q,2), p), remi2n(gel(Q,3), p));
     790       67172 :       Q1 = qfb3(shifti(gel(Q,1),-p), shifti(gel(Q,2),-p), shifti(gel(Q,3),-p));
     791             :     }
     792       68114 :     h = mm + (n>>1);
     793       68114 :     if (qfb_minexpi(Q1) <= h)
     794             :     {
     795         288 :       U = matid(2);
     796         288 :       QR = Q1;
     797             :     }
     798             :     else
     799       67826 :       QR = pqfbred_rec(Q1, h, &U);
     800      199690 :     while (qfb_maxexpi(QR) > h)
     801             :     {
     802      133340 :       if (is_minimal(QR, mm))
     803             :       {
     804        1764 :         going_to_r8 = 1;
     805        1764 :         break;
     806             :       }
     807      131576 :       QR = pqfbred_1(QR, mm, U);
     808             :     }
     809       68114 :     if (!going_to_r8)
     810             :     {
     811             :       GEN V;
     812       66350 :       QR = pqfbred_rec(QR, mm, &V);
     813       66350 :       U = ZM2_mul(U,V);
     814             :     }
     815       68114 :     if (p > 0)
     816             :     {
     817       67172 :       GEN Q0U = qfb_ZM_apply(Q0,U);
     818      134344 :       QR = mkqfb(addii(shifti(gel(QR,1), p), gel(Q0U,1)),
     819       67172 :                  addii(shifti(gel(QR,2), p), gel(Q0U,2)),
     820       67172 :                  addii(shifti(gel(QR,3), p), gel(Q0U,3)), d);
     821             :     }
     822             :   }
     823       68114 :   QR = pqfbred_iter_1(QR, m, U);
     824       68114 :   *pt_U = U;
     825       68114 :   return gc_all(av, 2, &QR, pt_U);
     826             : }
     827             : 
     828             : static GEN
     829      209525 : qfbredsl2_real(GEN Q, GEN isqrtD)
     830             : {
     831      209525 :   pari_sp av = avma;
     832      209525 :   if (2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     833      209525 :     return qfbredsl2_real_basecase(Q, isqrtD);
     834             :   else
     835             :   {
     836           0 :     GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     837           0 :     GEN Qf, Qr, W, U, t = NULL;
     838           0 :     long sa = signe(a), sb;
     839           0 :     if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
     840           0 :     if (signe(c) < 0)
     841             :     {
     842             :       GEN at;
     843           0 :       t  = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
     844           0 :       at = mulii(a,t);
     845           0 :       c = addii(subii(c, mulii(b, t)), mulii(at, t));
     846           0 :       b = subii(b, shifti(at,1));
     847             :     }
     848           0 :     sb = signe(b);
     849           0 :     Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
     850           0 :     if (sa < 0)
     851           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     852           0 :     if (sb < 0)
     853             :     {
     854           0 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     855           0 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     856             :     }
     857           0 :     if (t)
     858             :     {
     859           0 :       gcoeff(U,1,1) = subii( gcoeff(U,1,1), mulii(gcoeff(U,2,1), t));
     860           0 :       gcoeff(U,1,2) = subii( gcoeff(U,1,2), mulii(gcoeff(U,2,2), t));
     861             :     }
     862           0 :     W = qfbredsl2_real_basecase(Qr, isqrtD);
     863           0 :     Qf = gel(W,1);
     864           0 :     U = ZM2_mul(U,gel(W,2));
     865           0 :     return gerepilecopy(av, mkvec2(Qf,U));
     866             :   }
     867             : }
     868             : 
     869             : static GEN
     870     5025994 : qfbredsl2_imag(GEN Q)
     871             : {
     872     5025994 :   pari_sp av = avma;
     873             :   GEN Qt, U;
     874     5025994 :   if (2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     875     5025562 :     Qt = qfbredsl2_imag_basecase(Q, &U);
     876             :   else
     877             :   {
     878         432 :     long sb = signe(gel(Q,2));
     879             :     GEN Qr, W;
     880         432 :     Qr = pqfbred_rec(sb < 0 ? mkqfb(gel(Q,1), negi(gel(Q,2)), gel(Q,3), gel(Q,4)): Q, 0, &U);
     881         432 :     Qt = qfbredsl2_imag_basecase(Qr,&W);
     882         432 :     U = ZM2_mul(U,W);
     883         432 :     if (sb < 0)
     884             :     {
     885         230 :       gcoeff(U,2,1) = negi(gcoeff(U,2,1));
     886         230 :       gcoeff(U,2,2) = negi(gcoeff(U,2,2));
     887             :     }
     888             :   }
     889     5026014 :   return gerepilecopy(av, mkvec2(Qt,U));
     890             : }
     891             : 
     892             : GEN
     893     4715771 : redimagsl2(GEN Q, GEN *U)
     894             : {
     895     4715771 :   GEN q = qfbredsl2_imag(Q);
     896     4715791 :   *U = gel(q,2); return gel(q,1);
     897             : }
     898             : 
     899             : GEN
     900      519752 : qfbredsl2(GEN q, GEN isD)
     901             : {
     902             :   pari_sp av;
     903      519752 :   if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
     904      519752 :   if (qfb_is_qfi(q))
     905             :   {
     906      310220 :     if (isD) pari_err_TYPE("qfbredsl2", isD);
     907      310220 :     return qfbredsl2_imag(q);
     908             :   }
     909      209532 :   av = avma;
     910      209532 :   if (!isD) isD = sqrti(qfb_disc(q));
     911      208068 :   else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
     912      209525 :   return gerepileupto(av, qfbredsl2_real(q, isD));
     913             : }
     914             : 
     915             : static GEN
     916         420 : qfbred_real_i(GEN Q, long flag, GEN isqrtD, GEN sqrtD)
     917             : {
     918         420 :   if (typ(Q)!=t_QFB || 2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     919         420 :     return qfbred_real_basecase_i(Q, flag, isqrtD, sqrtD);
     920             :   else
     921             :   {
     922           0 :     GEN a = gel(Q,1), b = gel(Q,2), c = gel(Q,3), d = gel(Q,4);
     923           0 :     GEN Qr, W, U, t = NULL;
     924           0 :     long sa = signe(a), sb;
     925           0 :     if (sa < 0) { a = negi(a); b = negi(b); c = negi(c); }
     926           0 :     if (signe(c) < 0)
     927             :     {
     928             :       GEN at;
     929           0 :       t  = addiu(truedivii(subii(isqrtD,b),shifti(a,1)),1);
     930           0 :       at = mulii(a,t);
     931           0 :       c = addii(subii(c, mulii(b, t)), mulii(at, t));
     932           0 :       b = subii(b, shifti(at,1));
     933             :     }
     934           0 :     sb = signe(b);
     935           0 :     Qr = pqfbred_rec(mkqfb(a, sb < 0 ? negi(b): b, c, d), 0, &U);
     936           0 :     if (sa < 0)
     937           0 :       Qr = mkqfb(negi(gel(Qr,1)), negi(gel(Qr,2)), negi(gel(Qr,3)), gel(Qr,4));
     938           0 :     W = qfbred_real_basecase_i(Qr, flag, isqrtD, sqrtD);
     939           0 :     return gel(W,1);
     940             :   }
     941             : }
     942             : 
     943             : static GEN
     944          63 : qfbred_real(GEN x) { return qfbred_real_i(x,0,NULL,NULL); }
     945             : 
     946             : static GEN
     947    81453933 : qfbred_imag_basecase_av(pari_sp av, GEN q)
     948             : {
     949    81453933 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
     950    81453933 :   long cmp, lc = lgefint(c);
     951             : 
     952    81453933 :   if (lgefint(a) == 3 && lc == 3) return qfbred_imag_1(av, a, b, c, D);
     953      906434 :   cmp = abscmpii(a, b);
     954      914082 :   if (cmp < 0)
     955      434198 :     REDB(a,&b,&c);
     956      479884 :   else if (cmp == 0 && signe(b) < 0)
     957          30 :     b = negi(b);
     958             :   for(;;)
     959             :   {
     960     3108579 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     961     2918668 :     lc = lgefint(a); /* lg(future c): we swap a & c next */
     962     2918668 :     if (lc == 3) return qfbred_imag_1(av, a, b, c, D);
     963     2194497 :     swap(a,c); b = negi(b); /* apply rho */
     964     2194497 :     REDB(a,&b,&c);
     965     2194497 :     if (gc_needed(av, 2))
     966             :     {
     967           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfbred_imag, lc = %ld", lc);
     968           0 :       gerepileall(av, 3, &a,&b,&c);
     969             :     }
     970             :   }
     971      189911 :   if (cmp == 0 && signe(b) < 0) b = negi(b);
     972      189911 :   return gerepilecopy(av, mkqfb(a, b, c, D));
     973             : }
     974             : static GEN
     975    81250364 : qfbred_imag_av(pari_sp av, GEN Q)
     976             : {
     977    81250364 :   if (2*qfb_maxexpi(Q)-expi(gel(Q,4)) <= QFBRED_LIMIT)
     978    81458578 :     return qfbred_imag_basecase_av(av, Q);
     979             :   else
     980             :   {
     981           7 :     long sb = signe(gel(Q,2));
     982           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);
     983           7 :     return qfbred_imag_basecase_av(av, Qr);
     984             :   }
     985             : }
     986             : 
     987             : static GEN
     988    17936011 : qfbred_imag(GEN q) { return qfbred_imag_av(avma, q); }
     989             : 
     990             : GEN
     991       54782 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
     992             : {
     993             :   pari_sp av;
     994       54782 :   GEN q = check_qfbext("qfbred",x);
     995       54782 :   if (qfb_is_qfi(q)) return (flag & qf_STEP)? rhoimag(x): qfbred_imag(x);
     996         357 :   if (typ(x)==t_QFB) flag |= qf_NOD;
     997          49 :   else               flag &= ~qf_NOD;
     998         357 :   av = avma;
     999         357 :   return gerepilecopy(av, qfbred_real_i(x,flag,isqrtD,sqrtD));
    1000             : }
    1001             : /* t_QFB */
    1002             : GEN
    1003    17881596 : qfbred_i(GEN x) { return qfb_is_qfi(x)? qfbred_imag(x): qfbred_real(x); }
    1004             : GEN
    1005       52976 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
    1006             : /***********************************************************************/
    1007             : /**                                                                   **/
    1008             : /**                         Composition                               **/
    1009             : /**                                                                   **/
    1010             : /***********************************************************************/
    1011             : 
    1012             : static void
    1013    27276531 : qfb_sqr(GEN z, GEN x)
    1014             : {
    1015             :   GEN c, d1, x2, v1, v2, c3, m, p1, r;
    1016             : 
    1017    27276531 :   d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
    1018    27276280 :   c = gel(x,3);
    1019    27276280 :   m = mulii(c,x2);
    1020    27275893 :   if (equali1(d1))
    1021    20653626 :     v1 = v2 = gel(x,1);
    1022             :   else
    1023             :   {
    1024     6622449 :     v1 = diviiexact(gel(x,1),d1);
    1025     6622444 :     v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
    1026     6622448 :     c = mulii(c, d1);
    1027             :   }
    1028    27276074 :   togglesign(m);
    1029    27276125 :   r = modii(m,v2);
    1030    27275952 :   p1 = mulii(r, v1);
    1031    27275930 :   c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
    1032    27275907 :   gel(z,1) = mulii(v1,v2);
    1033    27275870 :   gel(z,2) = addii(gel(x,2), shifti(p1,1));
    1034    27275958 :   gel(z,3) = diviiexact(c3,v2);
    1035    27275923 : }
    1036             : /* z <- x * y */
    1037             : static void
    1038    65092564 : qfb_comp(GEN z, GEN x, GEN y)
    1039             : {
    1040             :   GEN n, c, d, y1, v1, v2, c3, m, p1, r;
    1041             : 
    1042    65092564 :   if (x == y) { qfb_sqr(z,x); return; }
    1043    38609716 :   n = shifti(subii(gel(y,2),gel(x,2)), -1);
    1044    38403411 :   v1 = gel(x,1);
    1045    38403411 :   v2 = gel(y,1);
    1046    38403411 :   c  = gel(y,3);
    1047    38403411 :   d = bezout(v2,v1,&y1,NULL);
    1048    38434985 :   if (equali1(d))
    1049    22325701 :     m = mulii(y1,n);
    1050             :   else
    1051             :   {
    1052    16175380 :     GEN s = subii(gel(y,2), n);
    1053    16173996 :     GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
    1054    16177672 :     if (!equali1(d1))
    1055             :     {
    1056     7903875 :       v1 = diviiexact(v1,d1);
    1057     7902976 :       v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
    1058     7902568 :       v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
    1059     7902472 :       c = mulii(c, d1);
    1060             :     }
    1061    16176042 :     m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
    1062             :   }
    1063    38454399 :   togglesign(m);
    1064    38398397 :   r = modii(m, v1);
    1065    38368318 :   p1 = mulii(r, v2);
    1066    38363259 :   c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
    1067    38353912 :   gel(z,1) = mulii(v1,v2);
    1068    38358629 :   gel(z,2) = addii(gel(y,2), shifti(p1,1));
    1069    38360765 :   gel(z,3) = diviiexact(c3,v1);
    1070             : }
    1071             : 
    1072             : /* not meant to be efficient */
    1073             : static GEN
    1074          84 : qfb_comp_gen(GEN x, GEN y)
    1075             : {
    1076          84 :   GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
    1077          84 :   GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
    1078          84 :   GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
    1079          84 :   GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
    1080             : 
    1081          84 :   if (!is_pm1(cx))
    1082             :   {
    1083          14 :     a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
    1084          14 :     c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
    1085             :   }
    1086          84 :   if (!is_pm1(cy))
    1087             :   {
    1088          28 :     a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
    1089          28 :     b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
    1090             :   }
    1091          84 :   D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
    1092         133 :   if (!Z_issquareall(diviiexact(d1, D), &n1) ||
    1093          84 :       !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
    1094          49 :   A = mulii(a1, n2);
    1095          49 :   B = mulii(a2, n1);
    1096          49 :   C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
    1097          49 :   U = ZV_extgcd(mkvec3(A, B, C));
    1098          49 :   m = gel(U,1); U = gmael(U,2,3);
    1099          49 :   A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
    1100          49 :   B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
    1101          49 :   C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
    1102          49 :   C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
    1103          49 :   B = addii(A, addii(B, C));
    1104          49 :   m2 = sqri(m);
    1105          49 :   A = diviiexact(mulii(a1, a2), m2);
    1106          49 :   C = diviiexact(shifti(subii(sqri(B),D), -2), A);
    1107          49 :   cx = mulii(cx, cy);
    1108          49 :   if (!is_pm1(cx))
    1109             :   {
    1110          14 :     A = mulii(A, cx); B = mulii(B, cx);
    1111          14 :     C = mulii(C, cx); D = mulii(D, sqri(cx));
    1112             :   }
    1113          49 :   return mkqfb(A, B, C, D);
    1114             : }
    1115             : 
    1116             : static GEN
    1117    62369211 : qficomp0(GEN x, GEN y, int raw)
    1118             : {
    1119    62369211 :   pari_sp av = avma;
    1120    62369211 :   GEN z = cgetg(5,t_QFB);
    1121    62361782 :   gel(z,4) = gel(x,4);
    1122    62361782 :   qfb_comp(z, x,y);
    1123    62072476 :   if (raw) return gerepilecopy(av,z);
    1124    62070698 :   return qfbred_imag_av(av, z);
    1125             : }
    1126             : static GEN
    1127         441 : qfrcomp0(GEN x, GEN y, int raw)
    1128             : {
    1129         441 :   pari_sp av = avma;
    1130         441 :   GEN dx = NULL, dy = NULL;
    1131         441 :   GEN z = cgetg(5,t_QFB);
    1132         441 :   if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1133         441 :   if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
    1134         441 :   gel(z,4) = gel(x,4);
    1135         441 :   qfb_comp(z, x,y);
    1136         441 :   if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
    1137         441 :   if (!raw) z = qfbred_real(z);
    1138         441 :   return gerepilecopy(av, z);
    1139             : }
    1140             : /* same discriminant, no distance, no checks */
    1141             : GEN
    1142    27245488 : qfbcomp_i(GEN x, GEN y)
    1143    27245488 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
    1144             : GEN
    1145      129398 : qfbcomp(GEN x, GEN y)
    1146             : {
    1147      129398 :   GEN qx = check_qfbext("qfbcomp", x);
    1148      129398 :   GEN qy = check_qfbext("qfbcomp", y);
    1149      129398 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1150             :   {
    1151          63 :     pari_sp av = avma;
    1152          63 :     GEN z = qfb_comp_gen(qx, qy);
    1153          63 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1154           7 :       pari_err_IMPL("Shanks's distance in general composition");
    1155          56 :     if (!z) pari_err_OP("*",x,y);
    1156          21 :     return gerepileupto(av, qfbred(z));
    1157             :   }
    1158      129335 :   return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
    1159             : }
    1160             : /* same discriminant, no distance, no checks */
    1161             : GEN
    1162           0 : qfbcompraw_i(GEN x, GEN y)
    1163           0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
    1164             : GEN
    1165        2198 : qfbcompraw(GEN x, GEN y)
    1166             : {
    1167        2198 :   GEN qx = check_qfbext("qfbcompraw", x);
    1168        2198 :   GEN qy = check_qfbext("qfbcompraw", y);
    1169        2198 :   if (!equalii(gel(qx,4),gel(qy,4)))
    1170             :   {
    1171          21 :     pari_sp av = avma;
    1172          21 :     GEN z = qfb_comp_gen(qx, qy);
    1173          21 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
    1174           0 :       pari_err_IMPL("Shanks's distance in general composition");
    1175          21 :     if (!z) pari_err_OP("qfbcompraw",x,y);
    1176          21 :     return gerepilecopy(av, z);
    1177             :   }
    1178        2177 :   if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
    1179        2177 :   return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
    1180             : }
    1181             : 
    1182             : static GEN
    1183      793666 : qfisqr0(GEN x, long raw)
    1184             : {
    1185      793666 :   pari_sp av = avma;
    1186      793666 :   GEN z = cgetg(5,t_QFB);
    1187      793666 :   gel(z,4) = gel(x,4);
    1188      793666 :   qfb_sqr(z,x);
    1189      793666 :   if (raw) return gerepilecopy(av,z);
    1190      793666 :   return qfbred_imag_av(av, z);
    1191             : }
    1192             : static GEN
    1193          35 : qfrsqr0(GEN x, long raw)
    1194             : {
    1195          35 :   pari_sp av = avma;
    1196          35 :   GEN dx = NULL, z = cgetg(5,t_QFB);
    1197          35 :   if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
    1198          35 :   gel(z,4) = gel(x,4); qfb_sqr(z,x);
    1199          35 :   if (dx) z = mkvec2(z, shiftr(dx,1));
    1200          35 :   if (!raw) z = qfbred_real(z);
    1201          35 :   return gerepilecopy(av, z);
    1202             : }
    1203             : /* same discriminant, no distance, no checks */
    1204             : GEN
    1205      698063 : qfbsqr_i(GEN x)
    1206      698063 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
    1207             : GEN
    1208       95638 : qfbsqr(GEN x)
    1209             : {
    1210       95638 :   GEN qx = check_qfbext("qfbsqr", x);
    1211       95638 :   return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
    1212             : }
    1213             : 
    1214             : static GEN
    1215        6867 : qfr_1_by_disc(GEN D)
    1216             : {
    1217             :   GEN y, r, s;
    1218        6867 :   check_quaddisc_real(D, NULL, "qfr_1_by_disc");
    1219        6867 :   y = cgetg(5,t_QFB);
    1220        6867 :   s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
    1221        6867 :   if (mpodd(r))
    1222             :   {
    1223        3535 :     s = subiu(s,1);
    1224        3535 :     r = subii(r, addiu(shifti(s, 1), 1));
    1225        3535 :     r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
    1226             :   }
    1227             :   else
    1228        3332 :   { r = shifti(r, -2); set_avma((pari_sp)s); }
    1229        6867 :   gel(y,1) = gen_1;
    1230        6867 :   gel(y,2) = s;
    1231        6867 :   gel(y,3) = icopy(r);
    1232        6867 :   gel(y,4) = icopy(D); return y;
    1233             : }
    1234             : 
    1235             : static GEN
    1236          35 : qfr_disc(GEN x)
    1237          35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
    1238             : 
    1239             : static GEN
    1240          35 : qfr_1(GEN x)
    1241          35 : { return qfr_1_by_disc(qfr_disc(x)); }
    1242             : 
    1243             : static void
    1244           0 : qfr_1_fill(GEN y, struct qfr_data *S)
    1245             : {
    1246           0 :   pari_sp av = avma;
    1247           0 :   GEN y2 = S->isqrtD;
    1248           0 :   gel(y,1) = gen_1;
    1249           0 :   if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
    1250           0 :   gel(y,2) = y2; av = avma;
    1251           0 :   gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(y2), S->D),-2));
    1252           0 : }
    1253             : static GEN
    1254           0 : qfr5_1(struct qfr_data *S, long prec)
    1255             : {
    1256           0 :   GEN y = cgetg(6, t_VEC);
    1257           0 :   qfr_1_fill(y, S);
    1258           0 :   gel(y,4) = gen_0;
    1259           0 :   gel(y,5) = real_1(prec); return y;
    1260             : }
    1261             : static GEN
    1262           0 : qfr3_1(struct qfr_data *S)
    1263             : {
    1264           0 :   GEN y = cgetg(4, t_VEC);
    1265           0 :   qfr_1_fill(y, S); return y;
    1266             : }
    1267             : 
    1268             : /* Assume D < 0 is the discriminant of a t_QFB */
    1269             : static GEN
    1270      720115 : qfi_1_by_disc(GEN D)
    1271             : {
    1272      720115 :   GEN b,c, y = cgetg(5,t_QFB);
    1273      720115 :   quadpoly_bc(D, mod2(D), &b,&c);
    1274      720115 :   if (b == gen_m1) b = gen_1;
    1275      720115 :   gel(y,1) = gen_1;
    1276      720115 :   gel(y,2) = b;
    1277      720115 :   gel(y,3) = c;
    1278      720115 :   gel(y,4) = icopy(D); return y;
    1279             : }
    1280             : static GEN
    1281      708157 : qfi_1(GEN x)
    1282             : {
    1283      708157 :   if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
    1284      708157 :   return qfi_1_by_disc(qfb_disc(x));
    1285             : }
    1286             : 
    1287             : GEN
    1288           0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
    1289             : 
    1290             : static GEN
    1291     9584177 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
    1292             : static GEN
    1293    25411372 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
    1294             : static GEN
    1295           7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
    1296             : static GEN
    1297           7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
    1298             : 
    1299             : static GEN
    1300           7 : qfipowraw(GEN x, long n)
    1301             : {
    1302           7 :   pari_sp av = avma;
    1303             :   GEN y;
    1304           7 :   if (!n) return qfi_1(x);
    1305           7 :   if (n== 1) return gcopy(x);
    1306           7 :   if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
    1307           7 :   if (n < 0) x = qfb_inv(x);
    1308           7 :   y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
    1309           7 :   return gerepilecopy(av,y);
    1310             : }
    1311             : 
    1312             : static GEN
    1313    11475786 : qfipow(GEN x, GEN n)
    1314             : {
    1315    11475786 :   pari_sp av = avma;
    1316             :   GEN y;
    1317    11475786 :   long s = signe(n);
    1318    11475786 :   if (!s) return qfi_1(x);
    1319    10767629 :   if (s < 0) x = qfb_inv(x);
    1320    10767628 :   y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
    1321    10767634 :   return gerepilecopy(av,y);
    1322             : }
    1323             : 
    1324             : static long
    1325      412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
    1326             : {
    1327             :   long z;
    1328      412328 :   *v = gen_0; *v2 = gen_1;
    1329     4351417 :   for (z=0; abscmpii(*v3,L) > 0; z++)
    1330             :   {
    1331     3939089 :     GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
    1332     3939089 :     *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
    1333             :   }
    1334      412328 :   return z;
    1335             : }
    1336             : 
    1337             : /* composition: Shanks' NUCOMP & NUDUPL */
    1338             : /* L = floor((|d|/4)^(1/4)) */
    1339             : GEN
    1340      400722 : nucomp(GEN x, GEN y, GEN L)
    1341             : {
    1342      400722 :   pari_sp av = avma;
    1343             :   long z;
    1344             :   GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
    1345             : 
    1346      400722 :   if (x==y) return nudupl(x,L);
    1347      400680 :   if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
    1348      400680 :   if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
    1349             : 
    1350      400680 :   if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
    1351      400680 :   s = shifti(addii(gel(x,2),gel(y,2)), -1);
    1352      400680 :   n = subii(gel(y,2), s);
    1353      400680 :   a1 = gel(x,1);
    1354      400680 :   a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
    1355      400680 :   if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
    1356      163576 :   else if (dvdii(s,d)) /* d | s */
    1357             :   {
    1358       83503 :     a = negi(mulii(u,n)); d1 = d;
    1359       83503 :     a1 = diviiexact(a1, d1);
    1360       83503 :     a2 = diviiexact(a2, d1);
    1361       83503 :     s = diviiexact(s, d1);
    1362             :   }
    1363             :   else
    1364             :   {
    1365             :     GEN p2, l;
    1366       80073 :     d1 = bezout(s,d,&u1,NULL);
    1367       80073 :     if (!equali1(d1))
    1368             :     {
    1369        2044 :       a1 = diviiexact(a1,d1);
    1370        2044 :       a2 = diviiexact(a2,d1);
    1371        2044 :       s = diviiexact(s,d1);
    1372        2044 :       d = diviiexact(d,d1);
    1373             :     }
    1374       80073 :     p1 = remii(gel(x,3),d);
    1375       80073 :     p2 = remii(gel(y,3),d);
    1376       80073 :     l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
    1377       80073 :     a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
    1378             :   }
    1379      400680 :   a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
    1380      400680 :   d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
    1381      400680 :   Q = cgetg(5,t_QFB);
    1382      400680 :   if (!z) {
    1383       37632 :     g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
    1384       37632 :     b = a2;
    1385       37632 :     b2 = gel(y,2);
    1386       37632 :     v2 = d1;
    1387       37632 :     gel(Q,1) = mulii(d,b);
    1388             :   } else {
    1389             :     GEN e, q3, q4;
    1390      363048 :     if (z&1) { v3 = negi(v3); v2 = negi(v2); }
    1391      363048 :     b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
    1392      363048 :     e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
    1393      363048 :     q3 = mulii(e,v2);
    1394      363048 :     q4 = subii(q3,s);
    1395      363048 :     b2 = addii(q3,q4);
    1396      363048 :     g = diviiexact(q4,v);
    1397      363048 :     if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
    1398      363048 :     gel(Q,1) = addii(mulii(d,b), mulii(e,v));
    1399             :   }
    1400      400680 :   q1 = mulii(b, v3);
    1401      400680 :   q2 = addii(q1,n);
    1402      400680 :   gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
    1403      400680 :   gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
    1404      400680 :   gel(Q,4) = gel(x,4);
    1405      400680 :   return qfbred_imag_av(av, Q);
    1406             : }
    1407             : 
    1408             : GEN
    1409       11648 : nudupl(GEN x, GEN L)
    1410             : {
    1411       11648 :   pari_sp av = avma;
    1412             :   long z;
    1413             :   GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
    1414             : 
    1415       11648 :   if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
    1416       11648 :   a = gel(x,1);
    1417       11648 :   b = gel(x,2);
    1418       11648 :   d1 = bezout(b,a, &u,NULL);
    1419       11648 :   if (!equali1(d1))
    1420             :   {
    1421        4620 :     a = diviiexact(a, d1);
    1422        4620 :     b = diviiexact(b, d1);
    1423             :   }
    1424       11648 :   c = modii(negi(mulii(u,gel(x,3))), a);
    1425       11648 :   p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
    1426       11648 :   d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
    1427       11648 :   a2 = sqri(d);
    1428       11648 :   c2 = sqri(v3);
    1429       11648 :   Q = cgetg(5,t_QFB);
    1430       11648 :   if (!z) {
    1431        1281 :     g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
    1432        1281 :     b2 = gel(x,2);
    1433        1281 :     v2 = d1;
    1434        1281 :     gel(Q,1) = a2;
    1435             :   } else {
    1436             :     GEN e;
    1437       10367 :     if (z&1) { v = negi(v); d = negi(d); }
    1438       10367 :     e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
    1439       10367 :     g = diviiexact(subii(mulii(e,v2), b), v);
    1440       10367 :     b2 = addii(mulii(e,v2), mulii(v,g));
    1441       10367 :     if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
    1442       10367 :     gel(Q,1) = addii(a2, mulii(e,v));
    1443             :   }
    1444       11648 :   gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
    1445       11648 :   gel(Q,3) = addii(c2, mulii(g,v2));
    1446       11648 :   gel(Q,4) = gel(x,4);
    1447       11648 :   return qfbred_imag_av(av, Q);
    1448             : }
    1449             : 
    1450             : static GEN
    1451        4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
    1452             : static GEN
    1453       11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
    1454             : GEN
    1455        1008 : nupow(GEN x, GEN n, GEN L)
    1456             : {
    1457             :   pari_sp av;
    1458             :   GEN y, D;
    1459             : 
    1460        1008 :   if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
    1461        1008 :   if (!is_qfi(x)) pari_err_TYPE("nupow",x);
    1462        1008 :   if (gequal1(n)) return gcopy(x);
    1463        1008 :   av = avma;
    1464        1008 :   D = qfb_disc(x);
    1465        1008 :   y = qfi_1_by_disc(D);
    1466        1008 :   if (!signe(n)) return y;
    1467         959 :   if (!L) L = sqrtnint(absi_shallow(D), 4);
    1468         959 :   y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
    1469         959 :   if (signe(n) < 0
    1470          35 :   && !absequalii(gel(y,1),gel(y,2))
    1471          35 :   && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
    1472         959 :   return gerepilecopy(av, y);
    1473             : }
    1474             : 
    1475             : /* Not stack-clean */
    1476             : GEN
    1477     1730407 : qfr5_compraw(GEN x, GEN y)
    1478             : {
    1479     1730407 :   GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
    1480     1730407 :   if (x == y)
    1481             :   {
    1482       34097 :     gel(z,4) = shifti(gel(x,4),1);
    1483       34097 :     gel(z,5) = sqrr(gel(x,5));
    1484             :   }
    1485             :   else
    1486             :   {
    1487     1696310 :     gel(z,4) = addii(gel(x,4),gel(y,4));
    1488     1696310 :     gel(z,5) = mulrr(gel(x,5),gel(y,5));
    1489             :   }
    1490     1730407 :   fix_expo(z); return z;
    1491             : }
    1492             : GEN
    1493     1730393 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
    1494     1730393 : { return qfr5_red(qfr5_compraw(x, y), S); }
    1495             : /* Not stack-clean */
    1496             : GEN
    1497     1006905 : qfr3_compraw(GEN x, GEN y)
    1498             : {
    1499     1006905 :   GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
    1500     1006905 :   return z;
    1501             : }
    1502             : GEN
    1503     1006905 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
    1504     1006905 : { return qfr3_red(qfr3_compraw(x,y), S); }
    1505             : 
    1506             : /* m > 0. Not stack-clean */
    1507             : static GEN
    1508           7 : qfr5_powraw(GEN x, long m)
    1509             : {
    1510           7 :   GEN y = NULL;
    1511          14 :   for (; m; m >>= 1)
    1512             :   {
    1513          14 :     if (m&1) y = y? qfr5_compraw(y,x): x;
    1514          14 :     if (m == 1) break;
    1515           7 :     x = qfr5_compraw(x,x);
    1516             :   }
    1517           7 :   return y;
    1518             : }
    1519             : 
    1520             : /* return x^n. Not stack-clean */
    1521             : GEN
    1522          21 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
    1523             : {
    1524          21 :   GEN y = NULL;
    1525          21 :   long i, m, s = signe(n);
    1526          21 :   if (!s) return qfr5_1(S, lg(gel(x,5)));
    1527          21 :   if (s < 0) x = qfb_inv(x);
    1528          42 :   for (i=lgefint(n)-1; i>1; i--)
    1529             :   {
    1530          21 :     m = n[i];
    1531          56 :     for (; m; m>>=1)
    1532             :     {
    1533          56 :       if (m&1) y = y? qfr5_comp(y,x,S): x;
    1534          56 :       if (m == 1 && i == 2) break;
    1535          35 :       x = qfr5_comp(x,x,S);
    1536             :     }
    1537             :   }
    1538          21 :   return y;
    1539             : }
    1540             : /* m > 0; return x^m. Not stack-clean */
    1541             : static GEN
    1542           0 : qfr3_powraw(GEN x, long m)
    1543             : {
    1544           0 :   GEN y = NULL;
    1545           0 :   for (; m; m>>=1)
    1546             :   {
    1547           0 :     if (m&1) y = y? qfr3_compraw(y,x): x;
    1548           0 :     if (m == 1) break;
    1549           0 :     x = qfr3_compraw(x,x);
    1550             :   }
    1551           0 :   return y;
    1552             : }
    1553             : /* return x^n. Not stack-clean */
    1554             : GEN
    1555        4515 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
    1556             : {
    1557        4515 :   GEN y = NULL;
    1558        4515 :   long i, m, s = signe(n);
    1559        4515 :   if (!s) return qfr3_1(S);
    1560        4515 :   if (s < 0) x = qfb_inv(x);
    1561        9046 :   for (i=lgefint(n)-1; i>1; i--)
    1562             :   {
    1563        4531 :     m = n[i];
    1564        5109 :     for (; m; m>>=1)
    1565             :     {
    1566        5089 :       if (m&1) y = y? qfr3_comp(y,x,S): x;
    1567        5089 :       if (m == 1 && i == 2) break;
    1568         578 :       x = qfr3_comp(x,x,S);
    1569             :     }
    1570             :   }
    1571        4515 :   return y;
    1572             : }
    1573             : 
    1574             : static GEN
    1575           7 : qfrinvraw(GEN x)
    1576             : {
    1577           7 :   if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
    1578           7 :  return qfbinv(x);
    1579             : }
    1580             : static GEN
    1581          14 : qfrpowraw(GEN x, long n)
    1582             : {
    1583          14 :   struct qfr_data S = { NULL, NULL, NULL };
    1584          14 :   pari_sp av = avma;
    1585          14 :   if (n==1) return gcopy(x);
    1586          14 :   if (n==-1) return qfrinvraw(x);
    1587           7 :   if (typ(x)==t_QFB)
    1588             :   {
    1589           0 :     GEN D = qfb_disc(x);
    1590           0 :     if (!n) return qfr_1(x);
    1591           0 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1592           0 :     x = qfr3_powraw(x, n);
    1593           0 :     x = qfr3_to_qfr(x, D);
    1594             :   }
    1595             :   else
    1596             :   {
    1597           7 :     GEN d0 = gel(x,2);
    1598           7 :     x = gel(x,1);
    1599           7 :     if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1600           7 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1601           7 :     x = qfr5_init(x, d0, &S);
    1602           7 :     if (labs(n) != 1) x = qfr5_powraw(x, n);
    1603           7 :     x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
    1604             :   }
    1605           7 :   return gerepilecopy(av, x);
    1606             : }
    1607             : static GEN
    1608         112 : qfrpow(GEN x, GEN n)
    1609             : {
    1610         112 :   struct qfr_data S = { NULL, NULL, NULL };
    1611         112 :   long s = signe(n);
    1612         112 :   pari_sp av = avma;
    1613         112 :   if (typ(x)==t_QFB)
    1614             :   {
    1615          42 :     if (!s) return qfr_1(x);
    1616          28 :     if (s < 0) x = qfb_inv(x);
    1617          28 :     x = qfr3_init(x, &S);
    1618          28 :     x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
    1619          28 :     x = qfr3_to_qfr(x, S.D);
    1620             :   }
    1621             :   else
    1622             :   {
    1623          70 :     GEN d0 = gel(x,2);
    1624          70 :     x = gel(x,1);
    1625          70 :     if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1626          49 :     if (s < 0) x = qfb_inv(x);
    1627          49 :     x = qfr5_init(x, d0, &S);
    1628          49 :     x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
    1629          49 :     x = qfr5_to_qfr(x, S.D, mulri(d0,n));
    1630             :   }
    1631          77 :   return gerepilecopy(av, x);
    1632             : }
    1633             : GEN
    1634          21 : qfbpowraw(GEN x, long n)
    1635             : {
    1636          21 :   GEN q = check_qfbext("qfbpowraw",x);
    1637          21 :   return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
    1638             : }
    1639             : /* same discriminant, no distance, no checks */
    1640             : GEN
    1641    10082572 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
    1642             : GEN
    1643     1393328 : qfbpow(GEN x, GEN n)
    1644             : {
    1645     1393328 :   GEN q = check_qfbext("qfbpow",x);
    1646     1393328 :   return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
    1647             : }
    1648             : GEN
    1649     1293818 : qfbpows(GEN x, long n)
    1650             : {
    1651     1293818 :   long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
    1652     1293818 :   affsi(n, N); return qfbpow(x, N);
    1653             : }
    1654             : 
    1655             : /* Prime forms attached to prime ideals of degree 1 */
    1656             : 
    1657             : /* assume x != 0 a t_INT, p > 0
    1658             :  * Return a t_QFB, but discriminant sign is not checked: can be used for
    1659             :  * real forms as well */
    1660             : GEN
    1661    12208764 : primeform_u(GEN x, ulong p)
    1662             : {
    1663    12208764 :   GEN c, y = cgetg(5, t_QFB);
    1664    12208644 :   pari_sp av = avma;
    1665             :   ulong b;
    1666             :   long s;
    1667             : 
    1668    12208644 :   s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
    1669             :   /* 2 or 3 mod 4 */
    1670    12208411 :   if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
    1671    12208724 :   if (p == 2) {
    1672     3880189 :     switch(s) {
    1673      589280 :       case 0: b = 0; break;
    1674     2940031 :       case 1: b = 1; break;
    1675      350878 :       case 4: b = 2; break;
    1676           0 :       default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1677           0 :                b = 0; /* -Wall */
    1678             :     }
    1679     3880189 :     c = shifti(subsi(s,x), -3);
    1680             :   } else {
    1681     8328535 :     b = Fl_sqrt(umodiu(x,p), p);
    1682     8330832 :     if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1683             :     /* mod(b) != mod2(x) ? */
    1684     8331362 :     if ((b ^ s) & 1) b = p - b;
    1685     8331362 :     c = diviuexact(shifti(subii(sqru(b), x), -2), p);
    1686             :   }
    1687    12193200 :   gel(y,3) = gerepileuptoint(av, c);
    1688    12208316 :   gel(y,4) = icopy(x);
    1689    12207961 :   gel(y,2) = utoi(b);
    1690    12208871 :   gel(y,1) = utoipos(p); return y;
    1691             : }
    1692             : 
    1693             : /* special case: p = 1 return unit form */
    1694             : GEN
    1695      135459 : primeform(GEN x, GEN p)
    1696             : {
    1697      135459 :   const char *f = "primeform";
    1698             :   pari_sp av;
    1699      135459 :   long s, sx = signe(x), sp = signe(p);
    1700             :   GEN y, b, absp;
    1701             : 
    1702      135459 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
    1703      135459 :   if (typ(p) != t_INT) pari_err_TYPE(f,p);
    1704      135459 :   if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
    1705      135459 :   if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
    1706      135459 :   if (lgefint(p) == 3)
    1707             :   {
    1708      135445 :     ulong pp = p[2];
    1709      135445 :     if (pp == 1) {
    1710       17782 :       if (sx < 0) {
    1711             :         long r;
    1712       10950 :         if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1713       10950 :         r = mod4(x);
    1714       10950 :         if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
    1715       10950 :         return qfi_1_by_disc(x);
    1716             :       }
    1717        6832 :       y = qfr_1_by_disc(x);
    1718        6832 :       if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
    1719        6832 :       return y;
    1720             :     }
    1721      117663 :     y = primeform_u(x, pp);
    1722      117656 :     if (sx < 0) {
    1723       89957 :       if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1724       89957 :       return y;
    1725             :     }
    1726       27699 :     if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
    1727       27699 :     return gcopy( qfr3_to_qfr(y, x) );
    1728             :   }
    1729          14 :   s = mod8(x);
    1730          14 :   if (sx < 0)
    1731             :   {
    1732           7 :     if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1733           7 :     if (s) s = 8-s;
    1734             :   }
    1735          14 :   y = cgetg(5, t_QFB);
    1736             :   /* 2 or 3 mod 4 */
    1737          14 :   if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
    1738          14 :   absp = absi_shallow(p); av = avma;
    1739          14 :   b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
    1740          14 :   s &= 1; /* s = x mod 2 */
    1741             :   /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
    1742          14 :   if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(absp,b));
    1743             : 
    1744          14 :   av = avma;
    1745          14 :   gel(y,3) = gerepileuptoint(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
    1746          14 :   gel(y,4) = icopy(x);
    1747          14 :   gel(y,2) = b;
    1748          14 :   gel(y,1) = icopy(p);
    1749          14 :   return y;
    1750             : }
    1751             : 
    1752             : static GEN
    1753     2620772 : normforms(GEN D, GEN fa)
    1754             : {
    1755             :   long i, j, k, lB, aN, sa;
    1756             :   GEN a, L, V, B, N, N2;
    1757     2620772 :   int D_odd = mpodd(D);
    1758     2620772 :   a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
    1759     2620772 :   sa = signe(a);
    1760     2620772 :   if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
    1761     1203972 :   V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
    1762     2551766 :            : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
    1763     2551766 :   if (!V) return NULL;
    1764      511966 :   N = gel(V,1); B = gel(V,2); lB = lg(B);
    1765      511966 :   N2 = shifti(N,1);
    1766      511966 :   aN = itou(diviiexact(a, N)); /* |a|/N */
    1767      511965 :   L = cgetg((lB-1)*aN+1, t_VEC);
    1768     2360561 :   for (k = 1, i = 1; i < lB; i++)
    1769             :   {
    1770     1848596 :     GEN b = shifti(gel(B,i), 1), c, C;
    1771     1848593 :     if (D_odd) b = addiu(b , 1);
    1772     1848593 :     c = diviiexact(shifti(subii(sqri(b), D), -2), a);
    1773     1848585 :     for (j = 0;; b = addii(b, N2))
    1774             :     {
    1775     2216659 :       gel(L, k++) = mkqfb(a, b, c, D);
    1776     2216665 :       if (++j == aN) break;
    1777      368071 :       C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
    1778      368074 :       c = sa > 0? addii(c, C): subii(c, C);
    1779             :     }
    1780             :   }
    1781      511965 :   return L;
    1782             : }
    1783             : 
    1784             : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
    1785             : static GEN
    1786      344319 : SL2_div_mul_e1(GEN N, GEN M)
    1787             : {
    1788      344319 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1789      344319 :   GEN A = mulii(gcoeff(N,1,1), d), B = mulii(gcoeff(N,1,2), b);
    1790      344316 :   GEN C = mulii(gcoeff(N,2,1), d), D = mulii(gcoeff(N,2,2), b);
    1791      344315 :   retmkvec2(subii(A,B), subii(C,D));
    1792             : }
    1793             : static GEN
    1794     1445677 : qfisolve_normform(GEN Q, GEN P)
    1795             : {
    1796     1445677 :   GEN a = gel(Q,1), N = gel(Q,2);
    1797     1445677 :   GEN M, b = qfbredsl2_imag_basecase(P, &M);
    1798     1445679 :   if (!qfb_equal(a,b)) return NULL;
    1799      102126 :   return SL2_div_mul_e1(N,M);
    1800             : }
    1801             : 
    1802             : /* Test equality modulo GL2 of two reduced forms */
    1803             : static int
    1804       61068 : GL2_qfb_equal(GEN a, GEN b)
    1805             : {
    1806       61068 :   return equalii(gel(a,1),gel(b,1))
    1807       11361 :    && absequalii(gel(a,2),gel(b,2))
    1808       72429 :    &&    equalii(gel(a,3),gel(b,3));
    1809             : }
    1810             : 
    1811             : /* Q(u,v) = p; if s < 0 return that solution; else the set of all solutions */
    1812             : static GEN
    1813       48012 : allsols(GEN Q, long s, GEN u, GEN v)
    1814             : {
    1815       48012 :   GEN w = mkvec2(u, v), b;
    1816       48010 :   if (signe(v) < 0) { u = negi(u); v = negi(v); } /* normalize for v >= 0 */
    1817       48010 :   w = mkvec2(u, v); if (s < 0) return w;
    1818       41369 :   if (!s) return mkvec(w);
    1819       38954 :   b = gel(Q,2); /* sum of the 2 solutions (if they exist) is -bv / a */
    1820       38954 :   if (signe(b))
    1821             :   { /* something to check */
    1822             :     GEN r, t;
    1823       13433 :     t = dvmdii(mulii(b, v), gel(Q,1), &r);
    1824       13433 :     if (signe(r)) return mkvec(w);
    1825        1820 :     u = addii(u, t);
    1826             :   }
    1827       27341 :   return mkvec2(w, mkvec2(negi(u), v));
    1828             : }
    1829             : static GEN
    1830      223031 : qfisolvep_all(GEN Q, GEN p, long all)
    1831             : {
    1832      223031 :   GEN R, U, V, M, N, x, q, D = qfb_disc(Q);
    1833      223027 :   long s = kronecker(D, p);
    1834             : 
    1835      223038 :   if (s < 0) return NULL;
    1836      126963 :   if (!all) s = -1; /* to indicate we want a single solution */
    1837             :   /* Solutions iff a class of maximal ideal above p is the class of Q;
    1838             :    * Two solutions iff (s > 0 and the class has order > 2), else one */
    1839      126963 :   if (!signe(gel(Q,2)))
    1840             :   { /* if principal form, use faster cornacchia */
    1841       43643 :     GEN a = gel(Q,1), c = gel(Q,3);
    1842       43643 :     if (equali1(a))
    1843             :     {
    1844       38167 :       if (!cornacchia(c, p, &M,&N)) return NULL;
    1845       33697 :       return allsols(Q, s, M, N);
    1846             :     }
    1847        5474 :     if (equali1(c))
    1848             :     {
    1849        5194 :       if (!cornacchia(a, p, &M,&N)) return NULL;
    1850         721 :       return allsols(Q, s, N, M);
    1851             :     }
    1852             :   }
    1853       83600 :   R = qfbredsl2_imag_basecase(Q, &U);
    1854       83601 :   if (equali1(gel(R,1)))
    1855             :   { /* principal form */
    1856       22533 :     if (!signe(gel(R,2)))
    1857             :     {
    1858        4396 :       if (!cornacchia(gel(R,3), p, &M,&N)) return NULL;
    1859         812 :       x = mkvec2(M,N);
    1860             :     }
    1861             :     else
    1862             :     { /* x^2 + xy + ((1-D)/4)y^2 = p <==> (2x + y)^2 - D y^2 = 4p */
    1863       18137 :       if (!cornacchia2(negi(D), p, &M, &N)) return NULL;
    1864        2331 :       x = subii(M,N); if (mpodd(x)) return NULL;
    1865        2331 :       x = mkvec2(shifti(x,-1), N);
    1866             :     }
    1867        3143 :     x = ZM_ZC_mul(U, x); x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1868        3143 :     return allsols(Q, s, gel(x,1), gel(x,2));
    1869             :   }
    1870       61068 :   q = qfbredsl2_imag_basecase(primeform(D, p), &V);
    1871       61068 :   if (!GL2_qfb_equal(R,q)) return NULL;
    1872       10451 :   if (signe(gel(R,2)) != signe(gel(q,2))) gcoeff(V,2,1) = negi(gcoeff(V,2,1));
    1873       10451 :   x = SL2_div_mul_e1(U,V); return allsols(Q, s, gel(x,1), gel(x,2));
    1874             : }
    1875             : GEN
    1876           0 : qfisolvep(GEN Q, GEN p)
    1877             : {
    1878           0 :   pari_sp av = avma;
    1879           0 :   GEN x = qfisolvep_all(Q, p, 0);
    1880           0 :   return x? gerepilecopy(av, x): gc_const(av, gen_0);
    1881             : }
    1882             : 
    1883             : static GEN
    1884      770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
    1885             : {
    1886      770126 :   pari_sp av = avma, btop;
    1887      770126 :   GEN M = N, P = qfbredsl2_real_basecase(Ps, rd), Q = P;
    1888             : 
    1889      770126 :   btop = avma;
    1890             :   for(;;)
    1891             :   {
    1892     5840681 :     if (qfb_equal(gel(M,1), gel(P,1)))
    1893      154084 :       return gerepileupto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
    1894     5686597 :     if (qfb_equal(gel(N,1), gel(Q,1)))
    1895       77658 :       return gerepileupto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
    1896     5608939 :     M = rhorealsl2(M, rd);
    1897     5608939 :     if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
    1898     5201735 :     Q = rhorealsl2(Q, rd);
    1899     5201735 :     if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
    1900     5070555 :     if (gc_needed(btop, 1)) gerepileall(btop, 2, &M, &Q);
    1901             :   }
    1902             : }
    1903             : 
    1904             : GEN
    1905           0 : qfrsolvep(GEN Q, GEN p)
    1906             : {
    1907           0 :   pari_sp av = avma;
    1908           0 :   GEN N, x, rd, d = qfb_disc(Q);
    1909             : 
    1910           0 :   if (kronecker(d, p) < 0) return gc_const(av, gen_0);
    1911           0 :   rd = sqrti(d);
    1912           0 :   N = qfbredsl2_real(Q, rd);
    1913           0 :   x = qfrsolve_normform(N, primeform(d, p), rd);
    1914           0 :   return x? gerepileupto(av, x): gc_const(av, gen_0);
    1915             : }
    1916             : 
    1917             : static GEN
    1918     1862930 : known_prime(GEN v)
    1919             : {
    1920     1862930 :   GEN p, e, fa = check_arith_all(v, "qfbsolve");
    1921     1862928 :   if (!fa) return BPSW_psp(v)? v: NULL;
    1922       42066 :   if (lg(gel(fa,1)) != 2) return NULL;
    1923       29340 :   p = gcoeff(fa,1,1);
    1924       29340 :   e = gcoeff(fa,1,2);
    1925       29340 :   return (equali1(e) && !is_pm1(p) && signe(p) > 0)? p: NULL;
    1926             : }
    1927             : static GEN
    1928     2215803 : qfsolve_normform(GEN Q, GEN f, GEN rd)
    1929     2215803 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
    1930             : static GEN
    1931     2843810 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
    1932             : {
    1933             :   GEN x, W, F, p;
    1934             :   long i, j, l;
    1935     2843810 :   if (!rd && (p = known_prime(fa))) return qfisolvep_all(Q, p, all);
    1936     2620765 :   F = normforms(qfb_disc(Q), fa);
    1937     2620771 :   if (!F) return NULL;
    1938      511965 :   if (!*Qr) *Qr = qfbredsl2(Q, rd);
    1939      511966 :   l = lg(F); W = all? cgetg(l, t_VEC): NULL;
    1940     2727258 :   for (j = i = 1; i < l; i++)
    1941     2215803 :     if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
    1942             :     {
    1943      333867 :       if (!all) return x;
    1944      333356 :       gel(W,j++) = x;
    1945             :     }
    1946      511455 :   if (j == 1) return NULL;
    1947      127456 :   setlg(W,j); return lexsort(W);
    1948             : }
    1949             : 
    1950             : static GEN
    1951     2838516 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
    1952             : static GEN
    1953     2828289 : qfbsolve_primitive(GEN Q, GEN fa, long all)
    1954             : {
    1955     2828289 :   GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
    1956     2828283 :   x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
    1957     2828289 :   if (!x) return cgetg(1, t_VEC);
    1958      174750 :   return x;
    1959             : }
    1960             : 
    1961             : /* f / g^2 */
    1962             : static GEN
    1963        5299 : famat_divsqr(GEN f, GEN g)
    1964        5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
    1965             : static GEN
    1966       10227 : qfbsolve_all(GEN Q, GEN n, long all)
    1967             : {
    1968       10227 :   GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
    1969       10227 :   GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
    1970       10227 :   long i, j, l = lg(D);
    1971       10227 :   W = all? cgetg(l, t_VEC): NULL;
    1972       25151 :   for (i = j = 1; i < l; i++)
    1973             :   {
    1974       15526 :     GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
    1975       15526 :     if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
    1976             :     {
    1977        1218 :       if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
    1978        1218 :       if (!all) return w;
    1979         616 :       gel(W,j++) = w;
    1980             :     }
    1981             :   }
    1982        9625 :   if (j == 1) return cgetg(1, t_VEC);
    1983         525 :   setlg(W,j); return lexsort(shallowconcat1(W));
    1984             : }
    1985             : 
    1986             : GEN
    1987     2838520 : qfbsolve(GEN Q, GEN n, long fl)
    1988             : {
    1989     2838520 :   pari_sp av = avma;
    1990     2838520 :   if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
    1991     2838520 :   if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
    1992     5666803 :   return gerepilecopy(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
    1993     2828288 :                                   : qfbsolve_primitive(Q, n, fl & 1));
    1994             : }
    1995             : 
    1996             : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
    1997             :  * Assume d > 0 and p is prime */
    1998             : long
    1999       55233 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
    2000             : {
    2001       55233 :   pari_sp av = avma;
    2002             :   GEN b, c, r;
    2003             : 
    2004       55233 :   *px = *py = gen_0;
    2005       55233 :   b = subii(p, d);
    2006       55177 :   if (signe(b) < 0) return gc_long(av,0);
    2007       54967 :   if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
    2008       54960 :   b = Fp_sqrt(b, p); /* sqrt(-d) */
    2009       55037 :   if (!b) return gc_long(av,0);
    2010       51306 :   b = gmael(halfgcdii(p, b), 2, 2);
    2011       51341 :   c = dvmdii(subii(p, sqri(b)), d, &r);
    2012       51220 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    2013       35484 :   set_avma(av);
    2014       35479 :   *px = icopy(b);
    2015       35467 :   *py = icopy(c); return 1;
    2016             : }
    2017             : 
    2018             : static GEN
    2019     1128535 : lastqi(GEN Q)
    2020             : {
    2021     1128535 :   GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
    2022     1128540 :   if (!signe(q)) return gen_0;
    2023     1128351 :   if (!signe(s)) return p;
    2024     1122527 :   if (is_pm1(q)) return subiu(p,1);
    2025     1122527 :   return divii(p, absi_shallow(q));
    2026             : }
    2027             : 
    2028             : static long
    2029     1128545 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
    2030             : {
    2031             :   GEN M, Q, V, c, r, b2;
    2032     1128545 :   if (!signe(b)) { /* d = p,2p,3p,4p */
    2033          14 :     set_avma(av);
    2034          14 :     if (absequalii(d, px4)){ *py = gen_1; return 1; }
    2035          14 :     if (absequalii(d, p))  { *py = gen_2; return 1; }
    2036           0 :     return 0;
    2037             :   }
    2038     1128531 :   if (mod2(b) != mod2(d)) b = subii(p,b);
    2039     1128512 :   M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
    2040     1128536 :   b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
    2041     1128428 :   b2 = sqri(b);
    2042     1128434 :   if (cmpii(b2,px4) > 0)
    2043             :   {
    2044     1120758 :     b = gel(V,1); b2 = sqri(b);
    2045     1120752 :     if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
    2046             :   }
    2047     1128461 :   c = dvmdii(subii(px4, b2), d, &r);
    2048     1128466 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    2049     1083725 :   set_avma(av);
    2050     1083723 :   *px = icopy(b);
    2051     1083728 :   *py = icopy(c); return 1;
    2052             : }
    2053             : 
    2054             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
    2055             :  * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
    2056             : long
    2057     1088755 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
    2058             : {
    2059     1088755 :   pari_sp av = avma;
    2060     1088755 :   GEN b, p4 = shifti(p,2);
    2061             : 
    2062     1088699 :   *px = *py = gen_0;
    2063     1088699 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2064     1087921 :   if (absequaliu(p, 2))
    2065             :   {
    2066           7 :     set_avma(av);
    2067           7 :     switch (itou_or_0(d)) {
    2068           0 :       case 4: *px = gen_2; break;
    2069           0 :       case 7: *px = gen_1; break;
    2070           7 :       default: return 0;
    2071             :     }
    2072           0 :     *py = gen_1; return 1;
    2073             :   }
    2074     1087915 :   b = Fp_sqrt(negi(d), p);
    2075     1087950 :   if (!b) return gc_long(av,0);
    2076     1087866 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2077             : }
    2078             : 
    2079             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
    2080             : long
    2081       40677 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
    2082             : {
    2083       40677 :   pari_sp av = avma;
    2084       40677 :   GEN p4 = shifti(p,2);
    2085       40677 :   *px = *py = gen_0;
    2086       40677 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    2087       40677 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    2088             : }
    2089             : 
    2090             : GEN
    2091        7630 : qfbcornacchia(GEN d, GEN p)
    2092             : {
    2093        7630 :   pari_sp av = avma;
    2094             :   GEN x, y;
    2095        7630 :   if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
    2096        7630 :   if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
    2097        7630 :   if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
    2098         287 :     return gerepilecopy(av, mkvec2(x, y));
    2099        7343 :   set_avma(av); return cgetg(1, t_VEC);
    2100             : }

Generated by: LCOV version 1.14