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 28574-61f195acfe) Lines: 938 1053 89.1 %
Date: 2023-06-08 07:47:28 Functions: 116 128 90.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      151484 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
      24             : {
      25             :   long r;
      26      151484 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
      27      151470 :   *s = signe(x);
      28      151470 :   if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
      29      151471 :   r = mod4(x); if (*s < 0 && r) r = 4 - r;
      30      151469 :   if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
      31      151455 :   if (pr) *pr = r;
      32      151455 : }
      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      957086 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
      50             : {
      51      957086 :   if (Dodd)
      52             :   {
      53      858659 :     pari_sp av = avma;
      54      858659 :     *b = gen_m1;
      55      858659 :     *c = gerepileuptoint(av, shifti(subui(1,D), -2));
      56             :   }
      57             :   else
      58             :   {
      59       98427 :     *b = gen_0;
      60       98427 :     *c = shifti(D,-2); togglesign(*c);
      61             :   }
      62      957086 : }
      63             : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
      64             : static GEN
      65      243355 : quadpoly_ii(GEN D, long Dmod4)
      66             : {
      67      243355 :   GEN b, c, y = cgetg(5,t_POL);
      68      243355 :   y[1] = evalsigne(1) | evalvarn(0);
      69      243355 :   quadpoly_bc(D, Dmod4, &b,&c);
      70      243355 :   gel(y,2) = c;
      71      243355 :   gel(y,3) = b;
      72      243355 :   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      241318 : 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     1792872 : check_qfbext(const char *fun, GEN x)
     113             : {
     114     1792872 :   long t = typ(x);
     115     1792872 :   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       77161 : qfb3(GEN x, GEN y, GEN z)
     127       77161 : { retmkqfb(icopy(x), icopy(y), icopy(z), qfb_disc3(x,y,z)); }
     128             : 
     129             : static int
     130    24028230 : qfb_equal(GEN x, GEN y)
     131             : {
     132    24028230 :   return equalii(gel(x,1),gel(y,1))
     133     1681743 :       && equalii(gel(x,2),gel(y,2))
     134    25709945 :       && equalii(gel(x,3),gel(y,3));
     135             : }
     136             : 
     137             : /* valid for t_QFB, qfr3, qfr5; shallow */
     138             : static GEN
     139      871764 : qfb_inv(GEN x) {
     140      871764 :   GEN z = shallowcopy(x);
     141      871766 :   gel(z,2) = negi(gel(z,2));
     142      871765 :   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       77189 : Qfb0(GEN a, GEN b, GEN c)
     151             : {
     152             :   GEN q, D;
     153       77189 :   if (!b)
     154             :   {
     155          28 :     if (c) pari_err_TYPE("Qfb",c);
     156          21 :     if (typ(a) == t_VEC && lg(a) == 4)
     157          14 :     { b = gel(a,2); c = gel(a,3); a = gel(a,1); }
     158           7 :     else if (typ(a) == t_POL && degpol(a) == 2)
     159           0 :     { b = gel(a,3); c = gel(a,2); a = gel(a,4); }
     160           7 :     else pari_err_TYPE("Qfb",a);
     161             :   }
     162       77161 :   else if (!c)
     163           7 :     pari_err_TYPE("Qfb",b);
     164       77168 :   if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
     165       77161 :   if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
     166       77161 :   if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
     167       77161 :   q = qfb3(a, b, c); D = qfb_disc(q);
     168       77161 :   if (signe(D) < 0)
     169       42357 :   { if (signe(a) < 0) pari_err_IMPL("negative definite t_QFB"); }
     170       34804 :   else if (Z_issquare(D)) pari_err_DOMAIN("Qfb","issquare(disc)","=", gen_1,q);
     171       77154 :   return q;
     172             : }
     173             : 
     174             : /* composition */
     175             : static void
     176    27071596 : qfb_sqr(GEN z, GEN x)
     177             : {
     178             :   GEN c, d1, x2, v1, v2, c3, m, p1, r;
     179             : 
     180    27071596 :   d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
     181    27071298 :   c = gel(x,3);
     182    27071298 :   m = mulii(c,x2);
     183    27070911 :   if (equali1(d1))
     184    20506219 :     v1 = v2 = gel(x,1);
     185             :   else
     186             :   {
     187     6564896 :     v1 = diviiexact(gel(x,1),d1);
     188     6564898 :     v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
     189     6564901 :     c = mulii(c, d1);
     190             :   }
     191    27071121 :   togglesign(m);
     192    27071176 :   r = modii(m,v2);
     193    27070976 :   p1 = mulii(r, v1);
     194    27070964 :   c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
     195    27070960 :   gel(z,1) = mulii(v1,v2);
     196    27070946 :   gel(z,2) = addii(gel(x,2), shifti(p1,1));
     197    27070966 :   gel(z,3) = diviiexact(c3,v2);
     198    27070916 : }
     199             : /* z <- x * y */
     200             : static void
     201    64761696 : qfb_comp(GEN z, GEN x, GEN y)
     202             : {
     203             :   GEN n, c, d, y1, v1, v2, c3, m, p1, r;
     204             : 
     205    64761696 :   if (x == y) { qfb_sqr(z,x); return; }
     206    38483772 :   n = shifti(subii(gel(y,2),gel(x,2)), -1);
     207    38271535 :   v1 = gel(x,1);
     208    38271535 :   v2 = gel(y,1);
     209    38271535 :   c  = gel(y,3);
     210    38271535 :   d = bezout(v2,v1,&y1,NULL);
     211    38296800 :   if (equali1(d))
     212    22288887 :     m = mulii(y1,n);
     213             :   else
     214             :   {
     215    16075436 :     GEN s = subii(gel(y,2), n);
     216    16074206 :     GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
     217    16077870 :     if (!equali1(d1))
     218             :     {
     219     7856856 :       v1 = diviiexact(v1,d1);
     220     7855909 :       v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
     221     7855618 :       v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
     222     7855520 :       c = mulii(c, d1);
     223             :     }
     224    16076060 :     m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
     225             :   }
     226    38320096 :   togglesign(m);
     227    38269783 :   r = modii(m, v1);
     228    38228353 :   p1 = mulii(r, v2);
     229    38230790 :   c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
     230    38222934 :   gel(z,1) = mulii(v1,v2);
     231    38225029 :   gel(z,2) = addii(gel(y,2), shifti(p1,1));
     232    38226574 :   gel(z,3) = diviiexact(c3,v1);
     233             : }
     234             : 
     235             : /* not meant to be efficient */
     236             : static GEN
     237          84 : qfb_comp_gen(GEN x, GEN y)
     238             : {
     239          84 :   GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
     240          84 :   GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
     241          84 :   GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
     242          84 :   GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
     243             : 
     244          84 :   if (!is_pm1(cx))
     245             :   {
     246          14 :     a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
     247          14 :     c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
     248             :   }
     249          84 :   if (!is_pm1(cy))
     250             :   {
     251          28 :     a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
     252          28 :     b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
     253             :   }
     254          84 :   D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
     255         133 :   if (!Z_issquareall(diviiexact(d1, D), &n1) ||
     256          84 :       !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
     257          49 :   A = mulii(a1, n2);
     258          49 :   B = mulii(a2, n1);
     259          49 :   C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
     260          49 :   U = ZV_extgcd(mkvec3(A, B, C));
     261          49 :   m = gel(U,1); U = gmael(U,2,3);
     262          49 :   A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
     263          49 :   B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
     264          49 :   C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
     265          49 :   C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
     266          49 :   B = addii(A, addii(B, C));
     267          49 :   m2 = sqri(m);
     268          49 :   A = diviiexact(mulii(a1, a2), m2);
     269          49 :   C = diviiexact(shifti(subii(sqri(B),D), -2), A);
     270          49 :   cx = mulii(cx, cy);
     271          49 :   if (!is_pm1(cx))
     272             :   {
     273          14 :     A = mulii(A, cx); B = mulii(B, cx);
     274          14 :     C = mulii(C, cx); D = mulii(D, sqri(cx));
     275             :   }
     276          49 :   return mkqfb(A, B, C, D);
     277             : }
     278             : 
     279             : static GEN redimag_av(pari_sp av, GEN q);
     280             : static GEN
     281    62036963 : qficomp0(GEN x, GEN y, int raw)
     282             : {
     283    62036963 :   pari_sp av = avma;
     284    62036963 :   GEN z = cgetg(5,t_QFB);
     285    62030801 :   gel(z,4) = gel(x,4);
     286    62030801 :   qfb_comp(z, x,y);
     287    61735172 :   if (raw) return gerepilecopy(av,z);
     288    61733401 :   return redimag_av(av, z);
     289             : }
     290             : static GEN redreal(GEN x);
     291             : static GEN
     292         441 : qfrcomp0(GEN x, GEN y, int raw)
     293             : {
     294         441 :   pari_sp av = avma;
     295         441 :   GEN dx = NULL, dy = NULL;
     296         441 :   GEN z = cgetg(5,t_QFB);
     297         441 :   if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
     298         441 :   if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
     299         441 :   gel(z,4) = gel(x,4);
     300         441 :   qfb_comp(z, x,y);
     301         441 :   if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
     302         441 :   if (!raw) z = redreal(z);
     303         441 :   return gerepilecopy(av, z);
     304             : }
     305             : /* same discriminant, no distance, no checks */
     306             : GEN
     307    27193779 : qfbcomp_i(GEN x, GEN y)
     308    27193779 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
     309             : GEN
     310      127606 : qfbcomp(GEN x, GEN y)
     311             : {
     312      127606 :   GEN qx = check_qfbext("qfbcomp", x);
     313      127606 :   GEN qy = check_qfbext("qfbcomp", y);
     314      127606 :   if (!equalii(gel(qx,4),gel(qy,4)))
     315             :   {
     316          63 :     pari_sp av = avma;
     317          63 :     GEN z = qfb_comp_gen(qx, qy);
     318          63 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
     319           7 :       pari_err_IMPL("Shanks's distance in general composition");
     320          56 :     if (!z) pari_err_OP("*",x,y);
     321          21 :     return gerepileupto(av, qfbred(z));
     322             :   }
     323      127543 :   return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
     324             : }
     325             : /* same discriminant, no distance, no checks */
     326             : GEN
     327           0 : qfbcompraw_i(GEN x, GEN y)
     328           0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
     329             : GEN
     330        2191 : qfbcompraw(GEN x, GEN y)
     331             : {
     332        2191 :   GEN qx = check_qfbext("qfbcompraw", x);
     333        2191 :   GEN qy = check_qfbext("qfbcompraw", y);
     334        2191 :   if (!equalii(gel(qx,4),gel(qy,4)))
     335             :   {
     336          21 :     pari_sp av = avma;
     337          21 :     GEN z = qfb_comp_gen(qx, qy);
     338          21 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
     339           0 :       pari_err_IMPL("Shanks's distance in general composition");
     340          21 :     if (!z) pari_err_OP("qfbcompraw",x,y);
     341          21 :     return gerepilecopy(av, z);
     342             :   }
     343        2170 :   if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
     344        2170 :   return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
     345             : }
     346             : 
     347             : static GEN
     348      793666 : qfisqr0(GEN x, long raw)
     349             : {
     350      793666 :   pari_sp av = avma;
     351      793666 :   GEN z = cgetg(5,t_QFB);
     352      793666 :   gel(z,4) = gel(x,4);
     353      793666 :   qfb_sqr(z,x);
     354      793666 :   if (raw) return gerepilecopy(av,z);
     355      793666 :   return redimag_av(av, z);
     356             : }
     357             : static GEN
     358          14 : qfrsqr0(GEN x, long raw)
     359             : {
     360          14 :   pari_sp av = avma;
     361          14 :   GEN dx = NULL, z = cgetg(5,t_QFB);
     362          14 :   if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
     363          14 :   gel(z,4) = gel(x,4); qfb_sqr(z,x);
     364          14 :   if (dx) z = mkvec2(z, shiftr(dx,1));
     365          14 :   if (!raw) z = redreal(z);
     366          14 :   return gerepilecopy(av, z);
     367             : }
     368             : /* same discriminant, no distance, no checks */
     369             : GEN
     370      698063 : qfbsqr_i(GEN x)
     371      698063 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
     372             : GEN
     373       95617 : qfbsqr(GEN x)
     374             : {
     375       95617 :   GEN qx = check_qfbext("qfbsqr", x);
     376       95617 :   return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
     377             : }
     378             : 
     379             : static GEN
     380        6867 : qfr_1_by_disc(GEN D)
     381             : {
     382             :   GEN y, r, s;
     383        6867 :   check_quaddisc_real(D, NULL, "qfr_1_by_disc");
     384        6867 :   y = cgetg(5,t_QFB);
     385        6867 :   s = sqrtremi(D, &r); togglesign(r); /* s^2 - r = D */
     386        6867 :   if (mpodd(r))
     387             :   {
     388        3535 :     s = subiu(s,1);
     389        3535 :     r = subii(r, addiu(shifti(s, 1), 1));
     390        3535 :     r = shifti(r, -2); set_avma((pari_sp)y); s = icopy(s);
     391             :   }
     392             :   else
     393        3332 :   { r = shifti(r, -2); set_avma((pari_sp)s); }
     394        6867 :   gel(y,1) = gen_1;
     395        6867 :   gel(y,2) = s;
     396        6867 :   gel(y,3) = icopy(r);
     397        6867 :   gel(y,4) = icopy(D); return y;
     398             : }
     399             : 
     400             : static GEN
     401          35 : qfr_disc(GEN x)
     402          35 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
     403             : 
     404             : static GEN
     405          35 : qfr_1(GEN x)
     406          35 : { return qfr_1_by_disc(qfr_disc(x)); }
     407             : 
     408             : static void
     409           0 : qfr_1_fill(GEN y, struct qfr_data *S)
     410             : {
     411           0 :   pari_sp av = avma;
     412           0 :   GEN y2 = S->isqrtD;
     413           0 :   gel(y,1) = gen_1;
     414           0 :   if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
     415           0 :   gel(y,2) = y2; av = avma;
     416           0 :   gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(y2), S->D),-2));
     417           0 : }
     418             : static GEN
     419           0 : qfr5_1(struct qfr_data *S, long prec)
     420             : {
     421           0 :   GEN y = cgetg(6, t_VEC);
     422           0 :   qfr_1_fill(y, S);
     423           0 :   gel(y,4) = gen_0;
     424           0 :   gel(y,5) = real_1(prec); return y;
     425             : }
     426             : static GEN
     427           0 : qfr3_1(struct qfr_data *S)
     428             : {
     429           0 :   GEN y = cgetg(4, t_VEC);
     430           0 :   qfr_1_fill(y, S); return y;
     431             : }
     432             : 
     433             : /* Assume D < 0 is the discriminant of a t_QFB */
     434             : static GEN
     435      713731 : qfi_1_by_disc(GEN D)
     436             : {
     437      713731 :   GEN b,c, y = cgetg(5,t_QFB);
     438      713731 :   quadpoly_bc(D, mod2(D), &b,&c);
     439      713731 :   if (b == gen_m1) b = gen_1;
     440      713731 :   gel(y,1) = gen_1;
     441      713731 :   gel(y,2) = b;
     442      713731 :   gel(y,3) = c;
     443      713731 :   gel(y,4) = icopy(D); return y;
     444             : }
     445             : static GEN
     446      701794 : qfi_1(GEN x)
     447             : {
     448      701794 :   if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
     449      701794 :   return qfi_1_by_disc(qfb_disc(x));
     450             : }
     451             : 
     452             : GEN
     453           0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
     454             : 
     455             : static GEN
     456     9506559 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
     457             : static GEN
     458    25212878 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
     459             : static GEN
     460           7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
     461             : static GEN
     462           7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
     463             : 
     464             : static GEN
     465           7 : qfipowraw(GEN x, long n)
     466             : {
     467           7 :   pari_sp av = avma;
     468             :   GEN y;
     469           7 :   if (!n) return qfi_1(x);
     470           7 :   if (n== 1) return gcopy(x);
     471           7 :   if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
     472           7 :   if (n < 0) x = qfb_inv(x);
     473           7 :   y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
     474           7 :   return gerepilecopy(av,y);
     475             : }
     476             : 
     477             : static GEN
     478    11386429 : qfipow(GEN x, GEN n)
     479             : {
     480    11386429 :   pari_sp av = avma;
     481             :   GEN y;
     482    11386429 :   long s = signe(n);
     483    11386429 :   if (!s) return qfi_1(x);
     484    10684635 :   if (s < 0) x = qfb_inv(x);
     485    10684635 :   y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
     486    10684636 :   return gerepilecopy(av,y);
     487             : }
     488             : 
     489             : static long
     490      412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
     491             : {
     492             :   long z;
     493      412328 :   *v = gen_0; *v2 = gen_1;
     494     4351417 :   for (z=0; abscmpii(*v3,L) > 0; z++)
     495             :   {
     496     3939089 :     GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
     497     3939089 :     *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
     498             :   }
     499      412328 :   return z;
     500             : }
     501             : 
     502             : /* composition: Shanks' NUCOMP & NUDUPL */
     503             : /* L = floor((|d|/4)^(1/4)) */
     504             : GEN
     505      400722 : nucomp(GEN x, GEN y, GEN L)
     506             : {
     507      400722 :   pari_sp av = avma;
     508             :   long z;
     509             :   GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
     510             : 
     511      400722 :   if (x==y) return nudupl(x,L);
     512      400680 :   if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
     513      400680 :   if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
     514             : 
     515      400680 :   if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
     516      400680 :   s = shifti(addii(gel(x,2),gel(y,2)), -1);
     517      400680 :   n = subii(gel(y,2), s);
     518      400680 :   a1 = gel(x,1);
     519      400680 :   a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
     520      400680 :   if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
     521      163576 :   else if (dvdii(s,d)) /* d | s */
     522             :   {
     523       83503 :     a = negi(mulii(u,n)); d1 = d;
     524       83503 :     a1 = diviiexact(a1, d1);
     525       83503 :     a2 = diviiexact(a2, d1);
     526       83503 :     s = diviiexact(s, d1);
     527             :   }
     528             :   else
     529             :   {
     530             :     GEN p2, l;
     531       80073 :     d1 = bezout(s,d,&u1,NULL);
     532       80073 :     if (!equali1(d1))
     533             :     {
     534        2044 :       a1 = diviiexact(a1,d1);
     535        2044 :       a2 = diviiexact(a2,d1);
     536        2044 :       s = diviiexact(s,d1);
     537        2044 :       d = diviiexact(d,d1);
     538             :     }
     539       80073 :     p1 = remii(gel(x,3),d);
     540       80073 :     p2 = remii(gel(y,3),d);
     541       80073 :     l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
     542       80073 :     a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
     543             :   }
     544      400680 :   a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
     545      400680 :   d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
     546      400680 :   Q = cgetg(5,t_QFB);
     547      400680 :   if (!z) {
     548       37632 :     g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
     549       37632 :     b = a2;
     550       37632 :     b2 = gel(y,2);
     551       37632 :     v2 = d1;
     552       37632 :     gel(Q,1) = mulii(d,b);
     553             :   } else {
     554             :     GEN e, q3, q4;
     555      363048 :     if (z&1) { v3 = negi(v3); v2 = negi(v2); }
     556      363048 :     b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
     557      363048 :     e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
     558      363048 :     q3 = mulii(e,v2);
     559      363048 :     q4 = subii(q3,s);
     560      363048 :     b2 = addii(q3,q4);
     561      363048 :     g = diviiexact(q4,v);
     562      363048 :     if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
     563      363048 :     gel(Q,1) = addii(mulii(d,b), mulii(e,v));
     564             :   }
     565      400680 :   q1 = mulii(b, v3);
     566      400680 :   q2 = addii(q1,n);
     567      400680 :   gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
     568      400680 :   gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
     569      400680 :   gel(Q,4) = gel(x,4);
     570      400680 :   return redimag_av(av, Q);
     571             : }
     572             : 
     573             : GEN
     574       11648 : nudupl(GEN x, GEN L)
     575             : {
     576       11648 :   pari_sp av = avma;
     577             :   long z;
     578             :   GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
     579             : 
     580       11648 :   if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
     581       11648 :   a = gel(x,1);
     582       11648 :   b = gel(x,2);
     583       11648 :   d1 = bezout(b,a, &u,NULL);
     584       11648 :   if (!equali1(d1))
     585             :   {
     586        4620 :     a = diviiexact(a, d1);
     587        4620 :     b = diviiexact(b, d1);
     588             :   }
     589       11648 :   c = modii(negi(mulii(u,gel(x,3))), a);
     590       11648 :   p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
     591       11648 :   d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
     592       11648 :   a2 = sqri(d);
     593       11648 :   c2 = sqri(v3);
     594       11648 :   Q = cgetg(5,t_QFB);
     595       11648 :   if (!z) {
     596        1281 :     g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
     597        1281 :     b2 = gel(x,2);
     598        1281 :     v2 = d1;
     599        1281 :     gel(Q,1) = a2;
     600             :   } else {
     601             :     GEN e;
     602       10367 :     if (z&1) { v = negi(v); d = negi(d); }
     603       10367 :     e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
     604       10367 :     g = diviiexact(subii(mulii(e,v2), b), v);
     605       10367 :     b2 = addii(mulii(e,v2), mulii(v,g));
     606       10367 :     if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
     607       10367 :     gel(Q,1) = addii(a2, mulii(e,v));
     608             :   }
     609       11648 :   gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
     610       11648 :   gel(Q,3) = addii(c2, mulii(g,v2));
     611       11648 :   gel(Q,4) = gel(x,4);
     612       11648 :   return redimag_av(av, Q);
     613             : }
     614             : 
     615             : static GEN
     616        4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
     617             : static GEN
     618       11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
     619             : GEN
     620        1008 : nupow(GEN x, GEN n, GEN L)
     621             : {
     622             :   pari_sp av;
     623             :   GEN y, D;
     624             : 
     625        1008 :   if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
     626        1008 :   if (!is_qfi(x)) pari_err_TYPE("nupow",x);
     627        1008 :   if (gequal1(n)) return gcopy(x);
     628        1008 :   av = avma;
     629        1008 :   D = qfb_disc(x);
     630        1008 :   y = qfi_1_by_disc(D);
     631        1008 :   if (!signe(n)) return y;
     632         959 :   if (!L) L = sqrtnint(absi_shallow(D), 4);
     633         959 :   y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
     634         959 :   if (signe(n) < 0
     635          35 :   && !absequalii(gel(y,1),gel(y,2))
     636          35 :   && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
     637         959 :   return gerepilecopy(av, y);
     638             : }
     639             : 
     640             : /* Reduction */
     641             : 
     642             : /* assume a > 0. Write b = q*2a + r, with -a < r <= a */
     643             : static GEN
     644     8604871 : dvmdii_round(GEN b, GEN a, GEN *r)
     645             : {
     646     8604871 :   GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
     647     8604573 :   if (signe(b) >= 0) {
     648     4946653 :     if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
     649             :   } else { /* r <= 0 */
     650     3657920 :     if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
     651             :   }
     652     8604228 :   return q;
     653             : }
     654             : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
     655             : static long
     656    96783660 : dvmdsu_round(long b, ulong a, long *r)
     657             : {
     658    96783660 :   ulong a2 = a << 1, q, ub, ur;
     659    96783660 :   if (b >= 0) {
     660    61699387 :     ub = b;
     661    61699387 :     q = ub / a2;
     662    61699387 :     ur = ub % a2;
     663    61699387 :     if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
     664    21743000 :     else *r = (long)ur;
     665    61699387 :     return (long)q;
     666             :   } else { /* r <= 0 */
     667    35084273 :     ub = (ulong)-b; /* |b| */
     668    35084273 :     q = ub / a2;
     669    35084273 :     ur = ub % a2;
     670    35084273 :     if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
     671    19182664 :     else *r = -(long)ur;
     672    35084273 :     return -(long)q;
     673             :   }
     674             : }
     675             : /* reduce b mod 2*a. Update b,c */
     676             : static void
     677     2768824 : REDB(GEN a, GEN *b, GEN *c)
     678             : {
     679     2768824 :   GEN r, q = dvmdii_round(*b, a, &r);
     680     2768824 :   if (!signe(q)) return;
     681     2699896 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     682     2699896 :   *b = r;
     683             : }
     684             : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
     685             : static void
     686    96785095 : sREDB(ulong a, long *b, ulong *c)
     687             : {
     688             :   long r, q;
     689             :   ulong uz;
     690   104969808 :   if (a > LONG_MAX) return; /* b already reduced */
     691    96785095 :   q = dvmdsu_round(*b, a, &r);
     692    96818207 :   if (q == 0) return;
     693             :   /* Final (a,r,c2) satisfies |r| <= |b| hence c2 <= c, c2 = c - q*z,
     694             :    * where z = (b+r) / 2, representable as long, has the same sign as q. */
     695    88633494 :   if (*b < 0)
     696             :   { /* uz = -z >= 0, q < 0 */
     697    31008614 :     if (r >= 0) /* different signs=>no overflow, exact division */
     698    15956188 :       uz = (ulong)-((*b + r)>>1);
     699             :     else
     700             :     {
     701    15052426 :       ulong ub = (ulong)-*b, ur = (ulong)-r;
     702    15052426 :       uz = (ub + ur) >> 1;
     703             :     }
     704    31008614 :     *c -= (-q) * uz; /* c -= qz */
     705             :   }
     706             :   else
     707             :   { /* uz = z >= 0, q > 0 */
     708    57624880 :     if (r <= 0)
     709    40020191 :       uz = (*b + r)>>1;
     710             :     else
     711             :     {
     712    17604689 :       ulong ub = (ulong)*b, ur = (ulong)r;
     713    17604689 :       uz = ((ub + ur) >> 1);
     714             :     }
     715    57624880 :     *c -= q * uz; /* c -= qz */
     716             :   }
     717    88633494 :   *b = r;
     718             : }
     719             : static void
     720     5836059 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
     721             : { /* REDB(a,b,c) */
     722     5836059 :   GEN r, q = dvmdii_round(*b, a, &r);
     723     5835400 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     724     5835336 :   *b = r;
     725     5835336 :   *u2 = subii(*u2, mulii(q, u1));
     726     5835390 : }
     727             : 
     728             : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
     729             : GEN
     730     2160206 : redimagsl2(GEN q, GEN *U)
     731             : {
     732     2160206 :   pari_sp av = avma;
     733             :   GEN z, u1,u2,v1,v2,Q;
     734     2160206 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     735             :   long cmp;
     736     2160206 :   u1 = gen_1; u2 = gen_0;
     737     2160206 :   cmp = abscmpii(a, b);
     738     2160239 :   if (cmp < 0)
     739     1000304 :     REDBU(a,&b,&c, u1,&u2);
     740     1159935 :   else if (cmp == 0 && signe(b) < 0)
     741             :   { /* b = -a */
     742          13 :     b = negi(b);
     743          13 :     u2 = gen_1;
     744             :   }
     745             :   for(;;)
     746             :   {
     747     6995585 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     748     4834898 :     swap(a,c); b = negi(b);
     749     4835767 :     z = u1; u1 = u2; u2 = negi(z);
     750     4835792 :     REDBU(a,&b,&c, u1,&u2);
     751     4835383 :     if (gc_needed(av, 1)) {
     752           0 :       if (DEBUGMEM>1) pari_warn(warnmem, "redimagsl2");
     753           0 :       gerepileall(av, 5, &a,&b,&c, &u1,&u2);
     754             :     }
     755             :   }
     756     2160072 :   if (cmp == 0 && signe(b) < 0)
     757             :   {
     758       16919 :     b = negi(b);
     759       16919 :     z = u1; u1 = u2; u2 = negi(z);
     760             :   }
     761             :   /* Let q = (A,B,C). q o [u1,u2; v1,v2] = Q implies
     762             :    * [v1,v2] = (1/C) [(b-B)/2 u1 - a u2, c u1 - (b+B)/2 u2] */
     763     2160072 :   z = shifti(subii(b, gel(q,2)), -1);
     764     2160004 :   v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
     765     2159948 :   z = subii(z, b);
     766     2159991 :   v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
     767     2159780 :   *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
     768     2160168 :   Q = lg(q)==5 ? mkqfb(a,b,c,gel(q,4)): mkvec3(a,b,c);
     769     2160158 :   return gc_all(av, 2, &Q, U);
     770             : }
     771             : 
     772             : static GEN
     773     1068574 : setq_b0(ulong a, ulong c, GEN D)
     774     1068574 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
     775             : /* assume |sb| = 1 */
     776             : static GEN
     777    79812744 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
     778    79812744 : { retmkqfb(utoipos(a), sb==1? utoipos(b): utoineg(b), utoipos(c), icopy(D)); }
     779             : /* 0 < a, c < 2^BIL, b = 0 */
     780             : static GEN
     781      958893 : redimag_1_b0(ulong a, ulong c, GEN D)
     782      958893 : { return (a <= c)? setq_b0(a, c, D): setq_b0(c, a, D); }
     783             : 
     784             : /* 0 < a, c < 2^BIL: single word affair */
     785             : static GEN
     786    80968642 : redimag_1(pari_sp av, GEN a, GEN b, GEN c, GEN D)
     787             : {
     788             :   ulong ua, ub, uc;
     789             :   long sb;
     790             :   for(;;)
     791      141005 :   { /* at most twice */
     792    80968642 :     long lb = lgefint(b); /* <= 3 after first loop */
     793    80968642 :     if (lb == 2) return redimag_1_b0(a[2],c[2], D);
     794    80009749 :     if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
     795      142046 :     REDB(a,&b,&c);
     796      140142 :     if (uel(a,2) <= uel(c,2))
     797             :     { /* lg(b) <= 3 but may be too large for itos */
     798           0 :       long s = signe(b);
     799           0 :       set_avma(av);
     800           0 :       if (!s) return redimag_1_b0(a[2], c[2], D);
     801           0 :       if (a[2] == c[2]) s = 1;
     802           0 :       return setq(a[2], b[2], c[2], s, D);
     803             :     }
     804      140142 :     swap(a,c); b = negi(b);
     805             :   }
     806             :   /* b != 0 */
     807    79867703 :   set_avma(av);
     808    79877231 :   ua = a[2];
     809    79877231 :   ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
     810    79877231 :   uc = c[2];
     811    79877231 :   if (ua < ub)
     812    28976329 :     sREDB(ua, &sb, &uc);
     813    50900902 :   else if (ua == ub && sb < 0) sb = (long)ub;
     814   147774144 :   while(ua > uc)
     815             :   {
     816    67851774 :     lswap(ua,uc); sb = -sb;
     817    67851774 :     sREDB(ua, &sb, &uc);
     818             :   }
     819    79922370 :   if (!sb) return setq_b0(ua, uc, D);
     820             :   else
     821             :   {
     822    79812688 :     long s = 1;
     823    79812688 :     if (sb < 0)
     824             :     {
     825    29949395 :       sb = -sb;
     826    29949395 :       if (ua != uc) s = -1;
     827             :     }
     828    79812688 :     return setq(ua, sb, uc, s, D);
     829             :   }
     830             : }
     831             : 
     832             : static GEN
     833    80764445 : redimag_av(pari_sp av, GEN q)
     834             : {
     835    80764445 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
     836    80764445 :   long cmp, lc = lgefint(c);
     837             : 
     838    80764445 :   if (lgefint(a) == 3 && lc == 3) return redimag_1(av, a, b, c, D);
     839      906428 :   cmp = abscmpii(a, b);
     840      914075 :   if (cmp < 0)
     841      434190 :     REDB(a,&b,&c);
     842      479885 :   else if (cmp == 0 && signe(b) < 0)
     843          30 :     b = negi(b);
     844             :   for(;;)
     845             :   {
     846     3108567 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     847     2918664 :     lc = lgefint(a); /* lg(future c): we swap a & c next */
     848     2918664 :     if (lc == 3) return redimag_1(av, a, b, c, D);
     849     2194492 :     swap(a,c); b = negi(b); /* apply rho */
     850     2194492 :     REDB(a,&b,&c);
     851             :   }
     852      189903 :   if (cmp == 0 && signe(b) < 0) b = negi(b);
     853      189903 :   return gerepilecopy(av, mkqfb(a, b, c, D));
     854             : }
     855             : static GEN
     856    17820491 : redimag(GEN q) { return redimag_av(avma, q); }
     857             : 
     858             : static GEN
     859           7 : rhoimag(GEN x)
     860             : {
     861           7 :   pari_sp av = avma;
     862           7 :   GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
     863           7 :   int fl = abscmpii(a, c);
     864           7 :   if (fl <= 0)
     865             :   {
     866           7 :     int fg = abscmpii(a, b);
     867           7 :     if (fg >= 0)
     868             :     {
     869           7 :       x = gcopy(x);
     870           7 :       if ((!fl || !fg) && signe(gel(x,2)) < 0) setsigne(gel(x,2), 1);
     871           7 :       return x;
     872             :     }
     873             :   }
     874           0 :   swap(a,c); b = negi(b);
     875           0 :   REDB(a, &b, &c);
     876           0 :   return gerepilecopy(av, mkqfb(a,b,c, qfb_disc(x)));
     877             : }
     878             : 
     879             : /* qfr3 / qfr5 */
     880             : 
     881             : /* t_QFB are unusable: D, sqrtD, isqrtD are recomputed all the time and the
     882             :  * logarithmic Shanks's distance is costly and hard to control.
     883             :  * qfr3 / qfr5 routines take a container of t_INTs (e.g a t_VEC) as argument,
     884             :  * at least 3 (resp. 5) components [it is a feature that they do not check the
     885             :  * precise type or length of the input]. They return a vector of length 3
     886             :  * (resp. 5). A qfr3 [a,b,c] contains the form coeffs, in a qfr5 [a,b,c, e,d]
     887             :  * the t_INT e is a binary exponent, d a t_REAL, coding the distance in
     888             :  * multiplicative form: the true distance is obtained from qfr5_dist.
     889             :  * All other qfr routines are obsolete (inefficient) wrappers */
     890             : 
     891             : /* static functions are not stack-clean. Unless mentionned otherwise, public
     892             :  * functions are. */
     893             : 
     894             : #define EMAX 22
     895             : static void
     896    12022738 : fix_expo(GEN x)
     897             : {
     898    12022738 :   if (expo(gel(x,5)) >= (1L << EMAX)) {
     899           0 :     gel(x,4) = addiu(gel(x,4), 1);
     900           0 :     shiftr_inplace(gel(x,5), - (1L << EMAX));
     901             :   }
     902    12022738 : }
     903             : 
     904             : /* (1/2) log (d * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
     905             : GEN
     906      184254 : qfr5_dist(GEN e, GEN d, long prec)
     907             : {
     908      184254 :   GEN t = logr_abs(d);
     909      184254 :   if (signe(e)) {
     910           0 :     GEN u = mulir(e, mplog2(prec));
     911           0 :     shiftr_inplace(u, EMAX); t = addrr(t, u);
     912             :   }
     913      184254 :   shiftr_inplace(t, -1); return t;
     914             : }
     915             : 
     916             : static void
     917    17290519 : rho_get_BC(GEN *B, GEN *C, GEN b, GEN c, struct qfr_data *S)
     918             : {
     919             :   GEN t, u;
     920    17290519 :   u = shifti(c,1);
     921    17290519 :   t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: c;
     922    17290519 :   u = remii(addii_sign(t,1, b,signe(b)), u);
     923    17290519 :   *B = addii_sign(t, 1, u, -signe(u)); /* |t| - (|t|+b) % |2c| */
     924    17290519 :   if (*B == gen_0)
     925        3413 :   { u = shifti(S->D, -2); setsigne(u, -1); }
     926             :   else
     927    17287106 :     u = shifti(addii_sign(sqri(*B),1, S->D,-1), -2);
     928    17290519 :   *C = diviiexact(u, c); /* = (B^2-D)/4c */
     929    17290519 : }
     930             : /* Not stack-clean */
     931             : GEN
     932     6994478 : qfr3_rho(GEN x, struct qfr_data *S)
     933             : {
     934     6994478 :   GEN B, C, b = gel(x,2), c = gel(x,3);
     935     6994478 :   rho_get_BC(&B, &C, b, c, S);
     936     6994478 :   return mkvec3(c,B,C);
     937             : }
     938             : /* Not stack-clean */
     939             : GEN
     940    10296041 : qfr5_rho(GEN x, struct qfr_data *S)
     941             : {
     942    10296041 :   GEN B, C, y, b = gel(x,2), c = gel(x,3);
     943    10296041 :   long sb = signe(b);
     944             : 
     945    10296041 :   rho_get_BC(&B, &C, b, c, S);
     946    10296041 :   y = mkvec5(c,B,C, gel(x,4), gel(x,5));
     947    10296041 :   if (sb) {
     948    10292324 :     GEN t = subii(sqri(b), S->D);
     949    10292324 :     if (sb < 0)
     950     2299122 :       t = divir(t, sqrr(subir(b,S->sqrtD)));
     951             :     else
     952     7993202 :       t = divri(sqrr(addir(b,S->sqrtD)), t);
     953             :     /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
     954    10292324 :     gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
     955             :   }
     956    10296041 :   return y;
     957             : }
     958             : 
     959             : /* Not stack-clean */
     960             : GEN
     961      217084 : qfr_to_qfr5(GEN x, long prec)
     962      217084 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
     963             : 
     964             : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
     965             : GEN
     966         462 : qfr5_to_qfr(GEN x, GEN D, GEN d0)
     967             : {
     968         462 :   if (d0)
     969             :   {
     970         140 :     GEN n = gel(x,4), d = absr(gel(x,5));
     971         140 :     if (signe(n))
     972             :     {
     973           0 :       n = addis(shifti(n, EMAX), expo(d));
     974           0 :       setexpo(d, 0); d = logr_abs(d);
     975           0 :       if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
     976           0 :       shiftr_inplace(d, -1);
     977           0 :       d0 = addrr(d0, d);
     978             :     }
     979         140 :     else if (!gequal1(d)) /* avoid loss of precision */
     980             :     {
     981          91 :       d = logr_abs(d);
     982          91 :       shiftr_inplace(d, -1);
     983          91 :       d0 = addrr(d0, d);
     984             :     }
     985             :   }
     986         462 :   x = qfr3_to_qfr(x, D);
     987         462 :   return d0 ? mkvec2(x,d0): x;
     988             : }
     989             : 
     990             : /* Not stack-clean */
     991             : GEN
     992       31913 : qfr3_to_qfr(GEN x, GEN d) { retmkqfb(gel(x,1), gel(x,2), gel(x,3), d); }
     993             : 
     994             : static int
     995    20824064 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
     996             : {
     997    20824064 :   if (signe(b) > 0 && abscmpii(b, isqrtD) <= 0)
     998             :   {
     999     5505704 :     GEN t = addii_sign(isqrtD,1, shifti(a,1),-1);
    1000     5505704 :     long l = abscmpii(b, t); /* compare |b| and |floor(sqrt(D)) - |2a|| */
    1001     5505704 :     if (l > 0 || (l == 0 && signe(t) < 0)) return 1;
    1002             :   }
    1003    16727067 :   return 0;
    1004             : }
    1005             : 
    1006             : INLINE int
    1007    17275795 : qfr_isreduced(GEN x, GEN isqrtD)
    1008             : {
    1009    17275795 :   return ab_isreduced(gel(x,1),gel(x,2),isqrtD);
    1010             : }
    1011             : 
    1012             : /* Not stack-clean */
    1013             : GEN
    1014     1947428 : qfr5_red(GEN x, struct qfr_data *S)
    1015             : {
    1016     1947428 :   pari_sp av = avma;
    1017    10247874 :   while (!qfr_isreduced(x, S->isqrtD))
    1018             :   {
    1019     8300446 :     x = qfr5_rho(x, S);
    1020     8300446 :     if (gc_needed(av,2))
    1021             :     {
    1022           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
    1023           0 :       x = gerepilecopy(av, x);
    1024             :     }
    1025             :   }
    1026     1947428 :   return x;
    1027             : }
    1028             : /* Not stack-clean */
    1029             : GEN
    1030     1169914 : qfr3_red(GEN x, struct qfr_data *S)
    1031             : {
    1032     1169914 :   pari_sp av = avma;
    1033     7027921 :   while (!qfr_isreduced(x, S->isqrtD))
    1034             :   {
    1035     5858007 :     x = qfr3_rho(x, S);
    1036     5858007 :     if (gc_needed(av,2))
    1037             :     {
    1038           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
    1039           0 :       x = gerepilecopy(av, x);
    1040             :     }
    1041             :   }
    1042     1169914 :   return x;
    1043             : }
    1044             : 
    1045             : void
    1046        2184 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
    1047             : {
    1048        2184 :   S->D = D;
    1049        2184 :   S->sqrtD = sqrtr(itor(S->D,prec));
    1050        2184 :   S->isqrtD = truncr(S->sqrtD);
    1051        2184 : }
    1052             : 
    1053             : static GEN
    1054         140 : qfr5_init(GEN x, GEN d, struct qfr_data *S)
    1055             : {
    1056         140 :   long prec = realprec(d), l = -expo(d);
    1057         140 :   if (l < BITS_IN_LONG) l = BITS_IN_LONG;
    1058         140 :   prec = maxss(prec, nbits2prec(l));
    1059         140 :   S->D = qfb_disc(x);
    1060         140 :   x = qfr_to_qfr5(x,prec);
    1061         140 :   if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
    1062           0 :   else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
    1063             : 
    1064         140 :   if (!S->isqrtD)
    1065             :   {
    1066         126 :     pari_sp av=avma;
    1067             :     long e;
    1068         126 :     S->isqrtD = gcvtoi(S->sqrtD,&e);
    1069         126 :     if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
    1070             :   }
    1071          14 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
    1072         140 :   return x;
    1073             : }
    1074             : static GEN
    1075         364 : qfr3_init(GEN x, struct qfr_data *S)
    1076             : {
    1077         364 :   S->D = qfb_disc(x);
    1078         364 :   if (!S->isqrtD) S->isqrtD = sqrti(S->D);
    1079         280 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
    1080         364 :   return x;
    1081             : }
    1082             : 
    1083             : #define qf_NOD  2
    1084             : #define qf_STEP 1
    1085             : 
    1086             : static GEN
    1087         399 : redreal_i(GEN x, long flag, GEN isqrtD, GEN sqrtD)
    1088             : {
    1089             :   struct qfr_data S;
    1090         399 :   GEN d = NULL, y;
    1091         399 :   if (typ(x)==t_VEC) { d = gel(x,2); x = gel(x,1); } else flag |= qf_NOD;
    1092         399 :   S.sqrtD = sqrtD;
    1093         399 :   S.isqrtD = isqrtD;
    1094         399 :   y = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, d, &S);
    1095         399 :   switch(flag) {
    1096          56 :     case 0:              y = qfr5_red(y,&S); break;
    1097         315 :     case qf_NOD:         y = qfr3_red(y,&S); break;
    1098          21 :     case qf_STEP:        y = qfr5_rho(y,&S); break;
    1099           7 :     case qf_STEP|qf_NOD: y = qfr3_rho(y,&S); break;
    1100           0 :     default: pari_err_FLAG("qfbred");
    1101             :   }
    1102         399 :   return qfr5_to_qfr(y, qfb_disc(x), d);
    1103             : }
    1104             : static GEN
    1105          42 : redreal(GEN x) { return redreal_i(x,0,NULL,NULL); }
    1106             : 
    1107             : GEN
    1108       54736 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
    1109             : {
    1110             :   pari_sp av;
    1111       54736 :   GEN q = check_qfbext("qfbred",x);
    1112       54736 :   if (qfb_is_qfi(q)) return (flag & qf_STEP)? rhoimag(x): redimag(x);
    1113         357 :   if (typ(x)==t_QFB) flag |= qf_NOD;
    1114          49 :   else               flag &= ~qf_NOD;
    1115         357 :   av = avma;
    1116         357 :   return gerepilecopy(av, redreal_i(x,flag,isqrtD,sqrtD));
    1117             : }
    1118             : /* t_QFB */
    1119             : GEN
    1120    17766121 : qfbred_i(GEN x) { return qfb_is_qfi(x)? redimag(x): redreal(x); }
    1121             : GEN
    1122       52944 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
    1123             : 
    1124             : /* Not stack-clean */
    1125             : GEN
    1126     1730414 : qfr5_compraw(GEN x, GEN y)
    1127             : {
    1128     1730414 :   GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
    1129     1730414 :   if (x == y)
    1130             :   {
    1131       34104 :     gel(z,4) = shifti(gel(x,4),1);
    1132       34104 :     gel(z,5) = sqrr(gel(x,5));
    1133             :   }
    1134             :   else
    1135             :   {
    1136     1696310 :     gel(z,4) = addii(gel(x,4),gel(y,4));
    1137     1696310 :     gel(z,5) = mulrr(gel(x,5),gel(y,5));
    1138             :   }
    1139     1730414 :   fix_expo(z); return z;
    1140             : }
    1141             : GEN
    1142     1730400 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
    1143     1730400 : { return qfr5_red(qfr5_compraw(x, y), S); }
    1144             : /* Not stack-clean */
    1145             : GEN
    1146     1006919 : qfr3_compraw(GEN x, GEN y)
    1147             : {
    1148     1006919 :   GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
    1149     1006919 :   return z;
    1150             : }
    1151             : GEN
    1152     1006919 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
    1153     1006919 : { return qfr3_red(qfr3_compraw(x,y), S); }
    1154             : 
    1155             : /* m > 0. Not stack-clean */
    1156             : static GEN
    1157           7 : qfr5_powraw(GEN x, long m)
    1158             : {
    1159           7 :   GEN y = NULL;
    1160          14 :   for (; m; m >>= 1)
    1161             :   {
    1162          14 :     if (m&1) y = y? qfr5_compraw(y,x): x;
    1163          14 :     if (m == 1) break;
    1164           7 :     x = qfr5_compraw(x,x);
    1165             :   }
    1166           7 :   return y;
    1167             : }
    1168             : 
    1169             : /* return x^n. Not stack-clean */
    1170             : GEN
    1171          28 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
    1172             : {
    1173          28 :   GEN y = NULL;
    1174          28 :   long i, m, s = signe(n);
    1175          28 :   if (!s) return qfr5_1(S, lg(gel(x,5)));
    1176          28 :   if (s < 0) x = qfb_inv(x);
    1177          56 :   for (i=lgefint(n)-1; i>1; i--)
    1178             :   {
    1179          28 :     m = n[i];
    1180          70 :     for (; m; m>>=1)
    1181             :     {
    1182          70 :       if (m&1) y = y? qfr5_comp(y,x,S): x;
    1183          70 :       if (m == 1 && i == 2) break;
    1184          42 :       x = qfr5_comp(x,x,S);
    1185             :     }
    1186             :   }
    1187          28 :   return y;
    1188             : }
    1189             : /* m > 0; return x^m. Not stack-clean */
    1190             : static GEN
    1191           0 : qfr3_powraw(GEN x, long m)
    1192             : {
    1193           0 :   GEN y = NULL;
    1194           0 :   for (; m; m>>=1)
    1195             :   {
    1196           0 :     if (m&1) y = y? qfr3_compraw(y,x): x;
    1197           0 :     if (m == 1) break;
    1198           0 :     x = qfr3_compraw(x,x);
    1199             :   }
    1200           0 :   return y;
    1201             : }
    1202             : /* return x^n. Not stack-clean */
    1203             : GEN
    1204        4529 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
    1205             : {
    1206        4529 :   GEN y = NULL;
    1207        4529 :   long i, m, s = signe(n);
    1208        4529 :   if (!s) return qfr3_1(S);
    1209        4529 :   if (s < 0) x = qfb_inv(x);
    1210        9074 :   for (i=lgefint(n)-1; i>1; i--)
    1211             :   {
    1212        4545 :     m = n[i];
    1213        5137 :     for (; m; m>>=1)
    1214             :     {
    1215        5117 :       if (m&1) y = y? qfr3_comp(y,x,S): x;
    1216        5117 :       if (m == 1 && i == 2) break;
    1217         592 :       x = qfr3_comp(x,x,S);
    1218             :     }
    1219             :   }
    1220        4529 :   return y;
    1221             : }
    1222             : 
    1223             : static GEN
    1224           7 : qfrinvraw(GEN x)
    1225             : {
    1226           7 :   if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
    1227           7 :  return qfbinv(x);
    1228             : }
    1229             : static GEN
    1230          14 : qfrpowraw(GEN x, long n)
    1231             : {
    1232          14 :   struct qfr_data S = { NULL, NULL, NULL };
    1233          14 :   pari_sp av = avma;
    1234          14 :   if (n==1) return gcopy(x);
    1235          14 :   if (n==-1) return qfrinvraw(x);
    1236           7 :   if (typ(x)==t_QFB)
    1237             :   {
    1238           0 :     GEN D = qfb_disc(x);
    1239           0 :     if (!n) return qfr_1(x);
    1240           0 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1241           0 :     x = qfr3_powraw(x, n);
    1242           0 :     x = qfr3_to_qfr(x, D);
    1243             :   }
    1244             :   else
    1245             :   {
    1246           7 :     GEN d0 = gel(x,2);
    1247           7 :     x = gel(x,1);
    1248           7 :     if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1249           7 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1250           7 :     x = qfr5_init(x, d0, &S);
    1251           7 :     if (labs(n) != 1) x = qfr5_powraw(x, n);
    1252           7 :     x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
    1253             :   }
    1254           7 :   return gerepilecopy(av, x);
    1255             : }
    1256             : static GEN
    1257         133 : qfrpow(GEN x, GEN n)
    1258             : {
    1259         133 :   struct qfr_data S = { NULL, NULL, NULL };
    1260         133 :   long s = signe(n);
    1261         133 :   pari_sp av = avma;
    1262         133 :   if (typ(x)==t_QFB)
    1263             :   {
    1264          56 :     if (!s) return qfr_1(x);
    1265          42 :     if (s < 0) x = qfb_inv(x);
    1266          42 :     x = qfr3_init(x, &S);
    1267          42 :     x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
    1268          42 :     x = qfr3_to_qfr(x, S.D);
    1269             :   }
    1270             :   else
    1271             :   {
    1272          77 :     GEN d0 = gel(x,2);
    1273          77 :     x = gel(x,1);
    1274          77 :     if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1275          56 :     if (s < 0) x = qfb_inv(x);
    1276          56 :     x = qfr5_init(x, d0, &S);
    1277          56 :     x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
    1278          56 :     x = qfr5_to_qfr(x, S.D, mulri(d0,n));
    1279             :   }
    1280          98 :   return gerepilecopy(av, x);
    1281             : }
    1282             : GEN
    1283          21 : qfbpowraw(GEN x, long n)
    1284             : {
    1285          21 :   GEN q = check_qfbext("qfbpowraw",x);
    1286          21 :   return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
    1287             : }
    1288             : /* same discriminant, no distance, no checks */
    1289             : GEN
    1290    10003654 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
    1291             : GEN
    1292     1382910 : qfbpow(GEN x, GEN n)
    1293             : {
    1294     1382910 :   GEN q = check_qfbext("qfbpow",x);
    1295     1382909 :   return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
    1296             : }
    1297             : GEN
    1298     1283429 : qfbpows(GEN x, long n)
    1299             : {
    1300     1283429 :   long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
    1301     1283429 :   affsi(n, N); return qfbpow(x, N);
    1302             : }
    1303             : 
    1304             : /* Prime forms attached to prime ideals of degree 1 */
    1305             : 
    1306             : /* assume x != 0 a t_INT, p > 0
    1307             :  * Return a t_QFB, but discriminant sign is not checked: can be used for
    1308             :  * real forms as well */
    1309             : GEN
    1310    12115027 : primeform_u(GEN x, ulong p)
    1311             : {
    1312    12115027 :   GEN c, y = cgetg(5, t_QFB);
    1313    12114936 :   pari_sp av = avma;
    1314             :   ulong b;
    1315             :   long s;
    1316             : 
    1317    12114936 :   s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
    1318             :   /* 2 or 3 mod 4 */
    1319    12114830 :   if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
    1320    12115020 :   if (p == 2) {
    1321     3854291 :     switch(s) {
    1322      589280 :       case 0: b = 0; break;
    1323     2914184 :       case 1: b = 1; break;
    1324      350829 :       case 4: b = 2; break;
    1325           0 :       default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1326           0 :                b = 0; /* -Wall */
    1327             :     }
    1328     3854293 :     c = shifti(subsi(s,x), -3);
    1329             :   } else {
    1330     8260729 :     b = Fl_sqrt(umodiu(x,p), p);
    1331     8262405 :     if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1332             :     /* mod(b) != mod2(x) ? */
    1333     8262840 :     if ((b ^ s) & 1) b = p - b;
    1334     8262840 :     c = diviuexact(shifti(subii(sqru(b), x), -2), p);
    1335             :   }
    1336    12099534 :   gel(y,3) = gerepileuptoint(av, c);
    1337    12114306 :   gel(y,4) = icopy(x);
    1338    12114666 :   gel(y,2) = utoi(b);
    1339    12114611 :   gel(y,1) = utoipos(p); return y;
    1340             : }
    1341             : 
    1342             : /* special case: p = 1 return unit form */
    1343             : GEN
    1344       74370 : primeform(GEN x, GEN p)
    1345             : {
    1346       74370 :   const char *f = "primeform";
    1347             :   pari_sp av;
    1348       74370 :   long s, sx = signe(x), sp = signe(p);
    1349             :   GEN y, b, absp;
    1350             : 
    1351       74370 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
    1352       74370 :   if (typ(p) != t_INT) pari_err_TYPE(f,p);
    1353       74370 :   if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
    1354       74370 :   if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
    1355       74370 :   if (lgefint(p) == 3)
    1356             :   {
    1357       74356 :     ulong pp = p[2];
    1358       74356 :     if (pp == 1) {
    1359       17761 :       if (sx < 0) {
    1360             :         long r;
    1361       10929 :         if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1362       10929 :         r = mod4(x);
    1363       10929 :         if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
    1364       10929 :         return qfi_1_by_disc(x);
    1365             :       }
    1366        6832 :       y = qfr_1_by_disc(x);
    1367        6832 :       if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
    1368        6832 :       return y;
    1369             :     }
    1370       56595 :     y = primeform_u(x, pp);
    1371       56588 :     if (sx < 0) {
    1372       28889 :       if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1373       28889 :       return y;
    1374             :     }
    1375       27699 :     if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
    1376       27699 :     return gcopy( qfr3_to_qfr(y, x) );
    1377             :   }
    1378          14 :   s = mod8(x);
    1379          14 :   if (sx < 0)
    1380             :   {
    1381           7 :     if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1382           7 :     if (s) s = 8-s;
    1383             :   }
    1384          14 :   y = cgetg(5, t_QFB);
    1385             :   /* 2 or 3 mod 4 */
    1386          14 :   if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
    1387          14 :   absp = absi_shallow(p); av = avma;
    1388          14 :   b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
    1389          14 :   s &= 1; /* s = x mod 2 */
    1390             :   /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
    1391          14 :   if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(absp,b));
    1392             : 
    1393          14 :   av = avma;
    1394          14 :   gel(y,3) = gerepileuptoint(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
    1395          14 :   gel(y,4) = icopy(x);
    1396          14 :   gel(y,2) = b;
    1397          14 :   gel(y,1) = icopy(p);
    1398          14 :   return y;
    1399             : }
    1400             : 
    1401             : static GEN
    1402     2843871 : normforms(GEN D, GEN fa)
    1403             : {
    1404             :   long i, j, k, lB, aN, sa;
    1405             :   GEN a, L, V, B, N, N2;
    1406     2843871 :   int D_odd = mpodd(D);
    1407     2843871 :   a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
    1408     2843871 :   sa = signe(a);
    1409     2843871 :   if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
    1410     1299046 :   V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
    1411     2774865 :            : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
    1412     2774859 :   if (!V) return NULL;
    1413      638970 :   N = gel(V,1); B = gel(V,2); lB = lg(B);
    1414      638970 :   N2 = shifti(N,1);
    1415      638930 :   aN = itou(diviiexact(a, N)); /* |a|/N */
    1416      638912 :   L = cgetg((lB-1)*aN+1, t_VEC);
    1417     2738687 :   for (k = 1, i = 1; i < lB; i++)
    1418             :   {
    1419     2099726 :     GEN b = shifti(gel(B,i), 1), c, C;
    1420     2099622 :     if (D_odd) b = addiu(b , 1);
    1421     2099622 :     c = diviiexact(shifti(subii(sqri(b), D), -2), a);
    1422     2099543 :     for (j = 0;; b = addii(b, N2))
    1423             :     {
    1424     2467652 :       gel(L, k++) = mkqfb(a, b, c, D);
    1425     2467803 :       if (++j == aN) break;
    1426      368064 :       C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
    1427      368109 :       c = sa > 0? addii(c, C): subii(c, C);
    1428             :     }
    1429             :   }
    1430      638961 :   return L;
    1431             : }
    1432             : 
    1433             : static void
    1434      412843 : SL2_div_init(GEN N, GEN M, GEN *A, GEN *B, GEN *C, GEN *D)
    1435             : {
    1436      412843 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1437      412843 :   *A = mulii(gcoeff(N,1,1), d); *B = mulii(gcoeff(N,1,2), b);
    1438      412683 :   *C = mulii(gcoeff(N,2,1), d); *D = mulii(gcoeff(N,2,2), b);
    1439      412707 : }
    1440             : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
    1441             : static GEN
    1442      412840 : SL2_div_mul_e1(GEN N, GEN M)
    1443             : {
    1444      412840 :   GEN A, B, C, D; SL2_div_init(N, M, &A, &B, &C, &D);
    1445      412700 :   retmkvec2(subii(A,B), subii(C,D));
    1446             : }
    1447             : /* Let M and N in SL2(Z), return (N*[1,0;0,-1]*M^-1)[,1] */
    1448             : static GEN
    1449           0 : SL2_swap_div_mul_e1(GEN N, GEN M)
    1450             : {
    1451           0 :   GEN A, B, C, D; SL2_div_init(N, M, &A, &B, &C, &D);
    1452           0 :   retmkvec2(addii(A,B), addii(C,D));
    1453             : }
    1454             : 
    1455             : /* Test equality modulo GL2 of two reduced forms */
    1456             : static int
    1457           0 : GL2_qfb_equal(GEN a, GEN b)
    1458             : {
    1459           0 :   return equalii(gel(a,1),gel(b,1))
    1460           0 :    && absequalii(gel(a,2),gel(b,2))
    1461           0 :    &&    equalii(gel(a,3),gel(b,3));
    1462             : }
    1463             : 
    1464             : static GEN
    1465     1690208 : qfisolve_normform(GEN Q, GEN P)
    1466             : {
    1467     1690208 :   GEN a = gel(Q,1), N = gel(Q,2);
    1468     1690208 :   GEN M, b = redimagsl2(P, &M);
    1469     1690289 :   if (!qfb_equal(a,b)) return NULL;
    1470      181099 :   return SL2_div_mul_e1(N,M);
    1471             : }
    1472             : 
    1473             : static GEN
    1474           0 : qfbsolve_cornacchia(GEN c, GEN p)
    1475             : {
    1476             :   GEN M, N;
    1477           0 :   if (kronecker(negi(c), p) < 0 || !cornacchia(c, p, &M,&N)) return NULL;
    1478           0 :   return mkvec2(M, N);
    1479             : }
    1480             : 
    1481             : GEN
    1482           0 : qfisolvep(GEN Q, GEN p)
    1483             : {
    1484             :   GEN M, N, x, y, a, c, d, q;
    1485           0 :   pari_sp av = avma;
    1486           0 :   if (!signe(gel(Q,2)))
    1487             :   {
    1488           0 :     a = gel(Q,1);
    1489           0 :     c = gel(Q,3); /* if principal form, use faster cornacchia */
    1490           0 :     if (equali1(a))
    1491             :     {
    1492           0 :       if (!(x = qfbsolve_cornacchia(c, p))) return gc_const(av, gen_0);
    1493           0 :       return gerepilecopy(av, x);
    1494             :     }
    1495           0 :     if (equali1(c))
    1496             :     {
    1497           0 :       if (!(x = qfbsolve_cornacchia(a, p))) return gc_const(av, gen_0);
    1498           0 :       swap(gel(x,1), gel(x,2)); return gerepilecopy(av, x);
    1499             :     }
    1500             :   }
    1501           0 :   d = qfb_disc(Q); if (kronecker(d,p) < 0) return gen_0;
    1502           0 :   Q = redimagsl2(Q, &N);
    1503           0 :   if (equali1(gel(Q,1))) /* principal form */
    1504             :   {
    1505           0 :     if (!signe(gel(Q,2)))
    1506           0 :     { if (!(x = qfbsolve_cornacchia(gel(Q,3), p))) return gc_const(av, gen_0); }
    1507             :     else
    1508             :     { /* x^2 + xy + ((1-d)/4)y^2 = p <==> (2x + y)^2 - d y^2 = 4p */
    1509           0 :       if (!cornacchia2(negi(d), p, &x, &y)) return gc_const(av, gen_0);
    1510           0 :       x = subii(x,y); if (mpodd(x)) return gc_const(av, gen_0);
    1511           0 :       x = mkvec2(shifti(x,-1), y);
    1512             :     }
    1513           0 :     x = ZM_ZC_mul(N, x);
    1514           0 :     x[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1515           0 :     return gerepileupto(av, x);
    1516             :   }
    1517           0 :   q = redimagsl2(primeform(d, p), &M);
    1518           0 :   if (!GL2_qfb_equal(Q,q)) return gc_const(av, gen_0);
    1519           0 :   if (signe(gel(Q,2))==signe(gel(q,2)))
    1520           0 :     x = SL2_div_mul_e1(N,M);
    1521             :   else
    1522           0 :     x = SL2_swap_div_mul_e1(N,M);
    1523           0 :   return gerepileupto(av, x);
    1524             : }
    1525             : 
    1526             : static void
    1527    13379288 : _rhorealsl2(GEN *pa, GEN *pb, GEN *pc, GEN *pu1, GEN *pu2, GEN *pv1,
    1528             :             GEN *pv2, GEN d, GEN rd)
    1529             : {
    1530    13379288 :   GEN C = mpabs_shallow(*pc), t = addii(*pb, gmax_shallow(rd,C));
    1531    13379288 :   GEN r, q = truedvmdii(t, shifti(C,1), &r);
    1532    13379288 :   *pb = subii(t, addii(r, *pb));
    1533    13379288 :   *pa = *pc;
    1534    13379288 :   *pc = diviiexact(subii(sqri(*pb), d), shifti(*pa, 2));
    1535    13379288 :   if (signe(*pa) < 0) togglesign(q);
    1536    13379288 :   r = *pu1; *pu1 = *pv1; *pv1 = subii(mulii(q, *pv1), r);
    1537    13379288 :   r = *pu2; *pu2 = *pv2; *pv2 = subii(mulii(q, *pv2), r);
    1538    13379288 : }
    1539             : 
    1540             : static GEN
    1541    10810674 : rhorealsl2(GEN A, GEN rd)
    1542             : {
    1543    10810674 :   GEN V = gel(A,1), M = gel(A,2);
    1544    10810674 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
    1545    10810674 :   GEN u1 = gcoeff(M,1,1), v1 = gcoeff(M,1,2);
    1546    10810674 :   GEN u2 = gcoeff(M,2,1), v2 = gcoeff(M,2,2);
    1547    10810674 :   _rhorealsl2(&a,&b,&c, &u1,&u2,&v1,&v2, d, rd);
    1548    10810674 :   return mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2));
    1549             : }
    1550             : 
    1551             : static GEN
    1552      979655 : redrealsl2(GEN V, GEN rd)
    1553             : {
    1554      979655 :   pari_sp av = avma;
    1555      979655 :   GEN u1 = gen_1, u2 = gen_0, v1 = gen_0, v2 = gen_1;
    1556      979655 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
    1557     3548269 :   while (!ab_isreduced(a,b,rd))
    1558             :   {
    1559     2568614 :     _rhorealsl2(&a,&b,&c, &u1,&u2,&v1,&v2, d, rd);
    1560     2568614 :     if (gc_needed(av, 1))
    1561             :     {
    1562           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"redrealsl2");
    1563           0 :       gerepileall(av, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
    1564             :     }
    1565             :   }
    1566      979655 :   return gerepilecopy(av, mkvec2(mkqfb(a,b,c,d), mkmat22(u1,v1,u2,v2)));
    1567             : }
    1568             : 
    1569             : GEN
    1570      646543 : qfbredsl2(GEN q, GEN isD)
    1571             : {
    1572             :   pari_sp av;
    1573      646543 :   if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
    1574      646543 :   if (qfb_is_qfi(q))
    1575             :   {
    1576      437007 :     GEN v = cgetg(3,t_VEC);
    1577      437003 :     if (isD) pari_err_TYPE("qfbredsl2", isD);
    1578      437003 :     gel(v,1) = redimagsl2(q, &gel(v,2)); return v;
    1579             :   }
    1580      209536 :   av = avma;
    1581      209536 :   if (!isD) isD = sqrti(qfb_disc(q));
    1582      208068 :   else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
    1583      209529 :   return gerepileupto(av, redrealsl2(q, isD));
    1584             : }
    1585             : 
    1586             : static GEN
    1587      770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
    1588             : {
    1589      770126 :   pari_sp av = avma, btop;
    1590      770126 :   GEN M = N, P = redrealsl2(Ps, rd), Q = P;
    1591             : 
    1592      770126 :   btop = avma;
    1593             :   for(;;)
    1594             :   {
    1595     5840681 :     if (qfb_equal(gel(M,1), gel(P,1)))
    1596      154084 :       return gerepileupto(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
    1597     5686597 :     if (qfb_equal(gel(N,1), gel(Q,1)))
    1598       77658 :       return gerepileupto(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
    1599     5608939 :     M = rhorealsl2(M, rd);
    1600     5608939 :     if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
    1601     5201735 :     Q = rhorealsl2(Q, rd);
    1602     5201735 :     if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
    1603     5070555 :     if (gc_needed(btop, 1)) gerepileall(btop, 2, &M, &Q);
    1604             :   }
    1605             : }
    1606             : 
    1607             : GEN
    1608           0 : qfrsolvep(GEN Q, GEN p)
    1609             : {
    1610           0 :   pari_sp av = avma;
    1611           0 :   GEN N, x, rd, d = qfb_disc(Q);
    1612             : 
    1613           0 :   if (kronecker(d, p) < 0) return gc_const(av, gen_0);
    1614           0 :   rd = sqrti(d);
    1615           0 :   N = redrealsl2(Q, rd);
    1616           0 :   x = qfrsolve_normform(N, primeform(d, p), rd);
    1617           0 :   return x? gerepileupto(av, x): gc_const(av, gen_0);
    1618             : }
    1619             : 
    1620             : static GEN
    1621     2460334 : qfsolve_normform(GEN Q, GEN f, GEN rd)
    1622     2460334 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
    1623             : static GEN
    1624     2843871 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
    1625             : {
    1626     2843871 :   GEN x, W, F = normforms(qfb_disc(Q), fa);
    1627             :   long i, j, l;
    1628     2843855 :   if (!F) return NULL;
    1629      638960 :   if (!*Qr) *Qr = qfbredsl2(Q, rd);
    1630      638975 :   l = lg(F); W = all? cgetg(l, t_VEC): NULL;
    1631     3092068 :   for (j = i = 1; i < l; i++)
    1632     2460336 :     if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
    1633             :     {
    1634      412717 :       if (!all) return x;
    1635      405570 :       gel(W,j++) = x;
    1636             :     }
    1637      631732 :   if (j == 1) return NULL;
    1638      168780 :   setlg(W,j); return W;
    1639             : }
    1640             : 
    1641             : static GEN
    1642     2838571 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
    1643             : static GEN
    1644     2828344 : qfbsolve_primitive(GEN Q, GEN fa, long all)
    1645             : {
    1646     2828344 :   GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
    1647     2828345 :   x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
    1648     2828306 :   if (!x) return cgetg(1, t_VEC);
    1649      174767 :   return x;
    1650             : }
    1651             : 
    1652             : /* f / g^2 */
    1653             : static GEN
    1654        5299 : famat_divsqr(GEN f, GEN g)
    1655        5299 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
    1656             : static GEN
    1657       10227 : qfbsolve_all(GEN Q, GEN n, long all)
    1658             : {
    1659       10227 :   GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
    1660       10227 :   GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
    1661       10227 :   long i, j, l = lg(D);
    1662       10227 :   W = all? cgetg(l, t_VEC): NULL;
    1663       25151 :   for (i = j = 1; i < l; i++)
    1664             :   {
    1665       15526 :     GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
    1666       15526 :     if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
    1667             :     {
    1668        1218 :       if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
    1669        1218 :       if (!all) return w;
    1670         616 :       gel(W,j++) = w;
    1671             :     }
    1672             :   }
    1673        9625 :   if (j == 1) return cgetg(1, t_VEC);
    1674         525 :   setlg(W,j); return shallowconcat1(W);
    1675             : }
    1676             : 
    1677             : GEN
    1678     2838577 : qfbsolve(GEN Q, GEN n, long fl)
    1679             : {
    1680     2838577 :   pari_sp av = avma;
    1681     2838577 :   if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
    1682     2838577 :   if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
    1683     5666876 :   return gerepilecopy(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
    1684     2828342 :                                   : qfbsolve_primitive(Q, n, fl & 1));
    1685             : }
    1686             : 
    1687             : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
    1688             :  * Assume d > 0 and p is prime */
    1689             : long
    1690        7483 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
    1691             : {
    1692        7483 :   pari_sp av = avma;
    1693             :   GEN b, c, r;
    1694             : 
    1695        7483 :   *px = *py = gen_0;
    1696        7483 :   b = subii(p, d);
    1697        7483 :   if (signe(b) < 0) return gc_long(av,0);
    1698        7483 :   if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
    1699        7476 :   b = Fp_sqrt(b, p); /* sqrt(-d) */
    1700        7476 :   if (!b) return gc_long(av,0);
    1701        3745 :   b = gmael(halfgcdii(p, b), 2, 2);
    1702        3745 :   c = dvmdii(subii(p, sqri(b)), d, &r);
    1703        3745 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    1704         231 :   set_avma(av);
    1705         231 :   *px = icopy(b);
    1706         231 :   *py = icopy(c); return 1;
    1707             : }
    1708             : 
    1709             : static GEN
    1710     1111228 : lastqi(GEN Q)
    1711             : {
    1712     1111228 :   GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
    1713     1111226 :   if (!signe(q)) return gen_0;
    1714     1111219 :   if (!signe(s)) return p;
    1715     1106340 :   if (is_pm1(q)) return subiu(p,1);
    1716     1106339 :   return divii(p, absi_shallow(q));
    1717             : }
    1718             : 
    1719             : static long
    1720     1111233 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
    1721             : {
    1722             :   GEN M, Q, V, c, r, b2;
    1723     1111233 :   if (!signe(b)) { /* d = p,2p,3p,4p */
    1724           7 :     set_avma(av);
    1725           7 :     if (absequalii(d, px4)){ *py = gen_1; return 1; }
    1726           7 :     if (absequalii(d, p))  { *py = gen_2; return 1; }
    1727           0 :     return 0;
    1728             :   }
    1729     1111226 :   if (mod2(b) != mod2(d)) b = subii(p,b);
    1730     1111200 :   M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
    1731     1111232 :   b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
    1732     1111111 :   b2 = sqri(b);
    1733     1111097 :   if (cmpii(b2,px4) > 0)
    1734             :   {
    1735     1104371 :     b = gel(V,1); b2 = sqri(b);
    1736     1104342 :     if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
    1737             :   }
    1738     1111137 :   c = dvmdii(subii(px4, b2), d, &r);
    1739     1111130 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    1740     1081403 :   set_avma(av);
    1741     1081400 :   *px = icopy(b);
    1742     1081399 :   *py = icopy(c); return 1;
    1743             : }
    1744             : 
    1745             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
    1746             :  * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
    1747             : long
    1748     1070632 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
    1749             : {
    1750     1070632 :   pari_sp av = avma;
    1751     1070632 :   GEN b, p4 = shifti(p,2);
    1752             : 
    1753     1070573 :   *px = *py = gen_0;
    1754     1070573 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    1755     1070616 :   if (absequaliu(p, 2))
    1756             :   {
    1757           7 :     set_avma(av);
    1758           7 :     switch (itou_or_0(d)) {
    1759           0 :       case 4: *px = gen_2; break;
    1760           0 :       case 7: *px = gen_1; break;
    1761           7 :       default: return 0;
    1762             :     }
    1763           0 :     *py = gen_1; return 1;
    1764             :   }
    1765     1070610 :   b = Fp_sqrt(negi(d), p);
    1766     1070638 :   if (!b) return gc_long(av,0);
    1767     1070554 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    1768             : }
    1769             : 
    1770             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
    1771             : long
    1772       40677 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
    1773             : {
    1774       40677 :   pari_sp av = avma;
    1775       40677 :   GEN p4 = shifti(p,2);
    1776       40677 :   *px = *py = gen_0;
    1777       40677 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    1778       40677 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    1779             : }
    1780             : 
    1781             : GEN
    1782        7630 : qfbcornacchia(GEN d, GEN p)
    1783             : {
    1784        7630 :   pari_sp av = avma;
    1785             :   GEN x, y;
    1786        7630 :   if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
    1787        7630 :   if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
    1788        7630 :   if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
    1789         287 :     return gerepilecopy(av, mkvec2(x, y));
    1790        7343 :   set_avma(av); return cgetg(1, t_VEC);
    1791             : }

Generated by: LCOV version 1.14