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.14.0 lcov report (development 27775-aca467eab2) Lines: 943 1057 89.2 %
Date: 2022-07-03 07:33:15 Functions: 115 126 91.3 %
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      217720 : check_quaddisc(GEN x, long *s, long *pr, const char *f)
      24             : {
      25             :   long r;
      26      217720 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
      27      217699 :   *s = signe(x);
      28      217699 :   if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
      29      217699 :   r = mod4(x); if (*s < 0 && r) r = 4 - r;
      30      217699 :   if (r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
      31      217685 :   if (pr) *pr = r;
      32      217685 : }
      33             : void
      34       75950 : check_quaddisc_real(GEN x, long *r, const char *f)
      35             : {
      36       75950 :   long sx; check_quaddisc(x, &sx, r, f);
      37       75943 :   if (sx < 0) pari_err_DOMAIN(f, "disc","<",gen_0,x);
      38       75943 : }
      39             : void
      40         623 : check_quaddisc_imag(GEN x, long *r, const char *f)
      41             : {
      42         623 :   long sx; check_quaddisc(x, &sx, r, f);
      43         616 :   if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
      44         616 : }
      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      871771 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
      50             : {
      51      871771 :   if (Dodd)
      52             :   {
      53      783221 :     pari_sp av = avma;
      54      783221 :     *b = gen_m1;
      55      783221 :     *c = gerepileuptoint(av, shifti(subui(1,D), -2));
      56             :   }
      57             :   else
      58             :   {
      59       88550 :     *b = gen_0;
      60       88550 :     *c = shifti(D,-2); togglesign(*c);
      61             :   }
      62      871771 : }
      63             : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
      64             : static GEN
      65      216874 : quadpoly_ii(GEN D, long Dmod4)
      66             : {
      67      216874 :   GEN b, c, y = cgetg(5,t_POL);
      68      216874 :   y[1] = evalsigne(1) | evalvarn(0);
      69      216874 :   quadpoly_bc(D, Dmod4, &b,&c);
      70      216874 :   gel(y,2) = c;
      71      216874 :   gel(y,3) = b;
      72      216874 :   gel(y,4) = gen_1; return y;
      73             : }
      74             : GEN
      75        2023 : quadpoly(GEN D)
      76             : {
      77             :   long s, Dmod4;
      78        2023 :   check_quaddisc(D, &s, &Dmod4, "quadpoly");
      79        2016 :   return quadpoly_ii(D, Dmod4);
      80             : }
      81             : GEN /* no checks */
      82      214858 : quadpoly_i(GEN D) { return quadpoly_ii(D, Mod4(D)); }
      83             : 
      84             : GEN
      85        1022 : quadpoly0(GEN x, long v)
      86             : {
      87        1022 :   GEN T = quadpoly(x);
      88        1015 :   if (v > 0) setvarn(T, v);
      89        1015 :   return T;
      90             : }
      91             : 
      92             : GEN
      93           0 : quadgen(GEN x)
      94           0 : { retmkquad(quadpoly(x), gen_0, gen_1); }
      95             : 
      96             : GEN
      97         546 : quadgen0(GEN x, long v)
      98             : {
      99         546 :   if (v==-1) v = fetch_user_var("w");
     100         546 :   retmkquad(quadpoly0(x, v), gen_0, gen_1);
     101             : }
     102             : 
     103             : /***********************************************************************/
     104             : /**                                                                   **/
     105             : /**                      BINARY QUADRATIC FORMS                       **/
     106             : /**                                                                   **/
     107             : /***********************************************************************/
     108             : static int
     109      814191 : is_qfi(GEN q) { return typ(q)==t_QFB && qfb_is_qfi(q); }
     110             : 
     111             : static GEN
     112     1646617 : check_qfbext(const char *fun, GEN x)
     113             : {
     114     1646617 :   long t = typ(x);
     115     1646617 :   if (t == t_QFB) return x;
     116         175 :   if (t == t_VEC && lg(x)==3)
     117             :   {
     118         175 :     GEN q = gel(x,1);
     119         175 :     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       77133 : qfb3(GEN x, GEN y, GEN z)
     127       77133 : { retmkqfb(icopy(x), icopy(y), icopy(z), qfb_disc3(x,y,z)); }
     128             : 
     129             : static int
     130    22337952 : qfb_equal(GEN x, GEN y)
     131             : {
     132    22337952 :   return equalii(gel(x,1),gel(y,1))
     133     1426670 :       && equalii(gel(x,2),gel(y,2))
     134    23764622 :       && equalii(gel(x,3),gel(y,3));
     135             : }
     136             : 
     137             : /* valid for t_QFB, qfr3, qfr5; shallow */
     138             : static GEN
     139      817088 : qfb_inv(GEN x) {
     140      817088 :   GEN z = shallowcopy(x);
     141      817087 :   gel(z,2) = negi(gel(z,2));
     142      817087 :   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       77161 : Qfb0(GEN a, GEN b, GEN c)
     151             : {
     152             :   GEN q, D;
     153       77161 :   if (!b)
     154             :   {
     155          28 :     if (c) pari_err_TYPE("Qfb",c);
     156          21 :     if (typ(a) != t_VEC || lg(a) != 4) pari_err_TYPE("Qfb",a);
     157          14 :     b = gel(a,2);
     158          14 :     c = gel(a,3);
     159          14 :     a = gel(a,1);
     160             :   }
     161       77133 :   else if (!c)
     162           7 :     pari_err_TYPE("Qfb",b);
     163       77140 :   if (typ(a)!=t_INT) pari_err_TYPE("Qfb",a);
     164       77133 :   if (typ(b)!=t_INT) pari_err_TYPE("Qfb",b);
     165       77133 :   if (typ(c)!=t_INT) pari_err_TYPE("Qfb",c);
     166       77133 :   q = qfb3(a, b, c); D = qfb_disc(q);
     167       77133 :   if (signe(D) < 0)
     168       42329 :   { if (signe(a) < 0) pari_err_IMPL("negative definite t_QFB"); }
     169       34804 :   else if (Z_issquare(D)) pari_err_DOMAIN("Qfb","issquare(disc)","=", gen_1,q);
     170       77126 :   return q;
     171             : }
     172             : 
     173             : /* composition */
     174             : static void
     175    24774998 : qfb_sqr(GEN z, GEN x)
     176             : {
     177             :   GEN c, d1, x2, v1, v2, c3, m, p1, r;
     178             : 
     179    24774998 :   d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
     180    24774694 :   c = gel(x,3);
     181    24774694 :   m = mulii(c,x2);
     182    24774339 :   if (equali1(d1))
     183    18904647 :     v1 = v2 = gel(x,1);
     184             :   else
     185             :   {
     186     5869893 :     v1 = diviiexact(gel(x,1),d1);
     187     5869891 :     v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
     188     5869890 :     c = mulii(c, d1);
     189             :   }
     190    24774536 :   togglesign(m);
     191    24774618 :   r = modii(m,v2);
     192    24774372 :   p1 = mulii(r, v1);
     193    24774341 :   c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
     194    24774374 :   gel(z,1) = mulii(v1,v2);
     195    24774334 :   gel(z,2) = addii(gel(x,2), shifti(p1,1));
     196    24774352 :   gel(z,3) = diviiexact(c3,v2);
     197    24774347 : }
     198             : /* z <- x * y */
     199             : static void
     200    61145573 : qfb_comp(GEN z, GEN x, GEN y)
     201             : {
     202             :   GEN n, c, d, y1, v1, v2, c3, m, p1, r;
     203             : 
     204    61145573 :   if (x == y) { qfb_sqr(z,x); return; }
     205    37159763 :   n = shifti(subii(gel(y,2),gel(x,2)), -1);
     206    36943069 :   v1 = gel(x,1);
     207    36943069 :   v2 = gel(y,1);
     208    36943069 :   c  = gel(y,3);
     209    36943069 :   d = bezout(v2,v1,&y1,NULL);
     210    36972406 :   if (equali1(d))
     211    21979019 :     m = mulii(y1,n);
     212             :   else
     213             :   {
     214    15063880 :     GEN s = subii(gel(y,2), n);
     215    15062240 :     GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
     216    15065830 :     if (!equali1(d1))
     217             :     {
     218     7388201 :       v1 = diviiexact(v1,d1);
     219     7387314 :       v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
     220     7387013 :       v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
     221     7386845 :       c = mulii(c, d1);
     222             :     }
     223    15064504 :     m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
     224             :   }
     225    36997378 :   togglesign(m);
     226    36954799 :   r = modii(m, v1);
     227    36912241 :   p1 = mulii(r, v2);
     228    36911458 :   c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
     229    36902165 :   gel(z,1) = mulii(v1,v2);
     230    36902055 :   gel(z,2) = addii(gel(y,2), shifti(p1,1));
     231    36906165 :   gel(z,3) = diviiexact(c3,v1);
     232             : }
     233             : 
     234             : /* not meant to be efficient */
     235             : static GEN
     236          84 : qfb_comp_gen(GEN x, GEN y)
     237             : {
     238          84 :   GEN d1 = qfb_disc(x), d2 = qfb_disc(y);
     239          84 :   GEN a1 = gel(x,1), b1 = gel(x,2), c1 = gel(x,3), n1;
     240          84 :   GEN a2 = gel(y,1), b2 = gel(y,2), c2 = gel(y,3), n2;
     241          84 :   GEN cx = content(x), cy = content(y), A, B, C, D, U, m, m2;
     242             : 
     243          84 :   if (!is_pm1(cx))
     244             :   {
     245          14 :     a1 = diviiexact(a1, cx); b1 = diviiexact(b1, cx);
     246          14 :     c1 = diviiexact(c1, cx); d1 = diviiexact(d1, sqri(cx));
     247             :   }
     248          84 :   if (!is_pm1(cy))
     249             :   {
     250          28 :     a2 = diviiexact(a2, cy); c2 = diviiexact(c2, cy);
     251          28 :     b2 = diviiexact(b2, cy); d2 = diviiexact(d2, sqri(cy));
     252             :   }
     253          84 :   D = gcdii(d1, d2); if (signe(d1) < 0) setsigne(D, -1);
     254         133 :   if (!Z_issquareall(diviiexact(d1, D), &n1) ||
     255          84 :       !Z_issquareall(diviiexact(d2, D), &n2)) return NULL;
     256          49 :   A = mulii(a1, n2);
     257          49 :   B = mulii(a2, n1);
     258          49 :   C = shifti(addii(mulii(b1, n2), mulii(b2, n1)), -1);
     259          49 :   U = ZV_extgcd(mkvec3(A, B, C));
     260          49 :   m = gel(U,1); U = gmael(U,2,3);
     261          49 :   A = mulii(diviiexact(mulii(a1,b2),m), gel(U,1));
     262          49 :   B = mulii(diviiexact(mulii(a2,b1),m), gel(U,2));
     263          49 :   C = addii(mulii(b1,b2), mulii(D, mulii(n1,n2)));
     264          49 :   C = mulii(diviiexact(shifti(C,-1), m), gel(U,3));
     265          49 :   B = addii(A, addii(B, C));
     266          49 :   m2 = sqri(m);
     267          49 :   A = diviiexact(mulii(a1, a2), m2);
     268          49 :   C = diviiexact(shifti(subii(sqri(B),D), -2), A);
     269          49 :   cx = mulii(cx, cy);
     270          49 :   if (!is_pm1(cx))
     271             :   {
     272          14 :     A = mulii(A, cx); B = mulii(B, cx);
     273          14 :     C = mulii(C, cx); D = mulii(D, sqri(cx));
     274             :   }
     275          49 :   return mkqfb(A, B, C, D);
     276             : }
     277             : 
     278             : static GEN redimag_av(pari_sp av, GEN q);
     279             : static GEN
     280    58419873 : qficomp0(GEN x, GEN y, int raw)
     281             : {
     282    58419873 :   pari_sp av = avma;
     283    58419873 :   GEN z = cgetg(5,t_QFB);
     284    58414552 :   gel(z,4) = gel(x,4);
     285    58414552 :   qfb_comp(z, x,y);
     286    58128775 :   if (raw) return gerepilecopy(av,z);
     287    58127004 :   return redimag_av(av, z);
     288             : }
     289             : static GEN redreal(GEN x);
     290             : static GEN
     291         441 : qfrcomp0(GEN x, GEN y, int raw)
     292             : {
     293         441 :   pari_sp av = avma;
     294         441 :   GEN dx = NULL, dy = NULL;
     295         441 :   GEN z = cgetg(5,t_QFB);
     296         441 :   if (typ(x)==t_VEC) { dx = gel(x,2); x = gel(x,1); }
     297         441 :   if (typ(y)==t_VEC) { dy = gel(y,2); y = gel(y,1); }
     298         441 :   gel(z,4) = gel(x,4);
     299         441 :   qfb_comp(z, x,y);
     300         441 :   if (dx) z = mkvec2(z, dy? addrr(dx, dy): dx); else if (dy) z = mkvec2(z, dy);
     301         441 :   if (!raw) z = redreal(z);
     302         441 :   return gerepilecopy(av, z);
     303             : }
     304             : /* same discriminant, no distance, no checks */
     305             : GEN
     306    26611526 : qfbcomp_i(GEN x, GEN y)
     307    26611526 : { return qfb_is_qfi(x)? qficomp0(x,y,0): qfrcomp0(x,y,0); }
     308             : GEN
     309      121123 : qfbcomp(GEN x, GEN y)
     310             : {
     311      121123 :   GEN qx = check_qfbext("qfbcomp", x);
     312      121123 :   GEN qy = check_qfbext("qfbcomp", y);
     313      121123 :   if (!equalii(gel(qx,4),gel(qy,4)))
     314             :   {
     315          63 :     pari_sp av = avma;
     316          63 :     GEN z = qfb_comp_gen(qx, qy);
     317          63 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
     318           7 :       pari_err_IMPL("Shanks's distance in general composition");
     319          56 :     if (!z) pari_err_OP("*",x,y);
     320          21 :     return gerepileupto(av, qfbred(z));
     321             :   }
     322      121060 :   return qfb_is_qfi(qx)? qficomp0(x,y,0): qfrcomp0(x,y,0);
     323             : }
     324             : /* same discriminant, no distance, no checks */
     325             : GEN
     326           0 : qfbcompraw_i(GEN x, GEN y)
     327           0 : { return qfb_is_qfi(x)? qficomp0(x,y,1): qfrcomp0(x,y,1); }
     328             : GEN
     329        2191 : qfbcompraw(GEN x, GEN y)
     330             : {
     331        2191 :   GEN qx = check_qfbext("qfbcompraw", x);
     332        2191 :   GEN qy = check_qfbext("qfbcompraw", y);
     333        2191 :   if (!equalii(gel(qx,4),gel(qy,4)))
     334             :   {
     335          21 :     pari_sp av = avma;
     336          21 :     GEN z = qfb_comp_gen(qx, qy);
     337          21 :     if (typ(x) == t_VEC || typ(y) == t_VEC)
     338           0 :       pari_err_IMPL("Shanks's distance in general composition");
     339          21 :     if (!z) pari_err_OP("qfbcompraw",x,y);
     340          21 :     return gerepilecopy(av, z);
     341             :   }
     342        2170 :   if (!equalii(gel(qx,4),gel(qy,4))) pari_err_OP("qfbcompraw",x,y);
     343        2170 :   return qfb_is_qfi(qx)? qficomp0(x,y,1): qfrcomp0(x,y,1);
     344             : }
     345             : 
     346             : static GEN
     347      789179 : qfisqr0(GEN x, long raw)
     348             : {
     349      789179 :   pari_sp av = avma;
     350      789179 :   GEN z = cgetg(5,t_QFB);
     351      789179 :   gel(z,4) = gel(x,4);
     352      789179 :   qfb_sqr(z,x);
     353      789179 :   if (raw) return gerepilecopy(av,z);
     354      789179 :   return redimag_av(av, z);
     355             : }
     356             : static GEN
     357          14 : qfrsqr0(GEN x, long raw)
     358             : {
     359          14 :   pari_sp av = avma;
     360          14 :   GEN dx = NULL, z = cgetg(5,t_QFB);
     361          14 :   if (typ(x) == t_VEC) { dx = gel(x,2); x = gel(x,1); }
     362          14 :   gel(z,4) = gel(x,4); qfb_sqr(z,x);
     363          14 :   if (dx) z = mkvec2(z, shiftr(dx,1));
     364          14 :   if (!raw) z = redreal(z);
     365          14 :   return gerepilecopy(av, z);
     366             : }
     367             : /* same discriminant, no distance, no checks */
     368             : GEN
     369      693576 : qfbsqr_i(GEN x)
     370      693576 : { return qfb_is_qfi(x)? qfisqr0(x,0): qfrsqr0(x,0); }
     371             : GEN
     372       95617 : qfbsqr(GEN x)
     373             : {
     374       95617 :   GEN qx = check_qfbext("qfbsqr", x);
     375       95617 :   return qfb_is_qfi(qx)? qfisqr0(x,0): qfrsqr0(x,0);
     376             : }
     377             : 
     378             : static GEN
     379        6860 : qfr_1_by_disc(GEN D)
     380             : {
     381        6860 :   GEN y = cgetg(5,t_QFB), isqrtD;
     382        6860 :   pari_sp av = avma;
     383             :   long r;
     384             : 
     385        6860 :   check_quaddisc_real(D, &r, "qfr_1_by_disc");
     386        6860 :   gel(y,1) = gen_1; isqrtD = sqrti(D);
     387        6860 :   if ((r & 1) != mod2(isqrtD)) /* we know isqrtD > 0 */
     388        3535 :     isqrtD = gerepileuptoint(av, subiu(isqrtD,1));
     389        6860 :   gel(y,2) = isqrtD; av = avma;
     390        6860 :   gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(isqrtD), D),-2));
     391        6860 :   gel(y,4) = icopy(D); return y;
     392             : }
     393             : 
     394             : static GEN
     395          28 : qfr_disc(GEN x)
     396          28 : { return qfb_disc(typ(x)==t_VEC ? gel(x,1): x); }
     397             : 
     398             : static GEN
     399          28 : qfr_1(GEN x)
     400          28 : { return qfr_1_by_disc(qfr_disc(x)); }
     401             : 
     402             : static void
     403           0 : qfr_1_fill(GEN y, struct qfr_data *S)
     404             : {
     405           0 :   pari_sp av = avma;
     406           0 :   GEN y2 = S->isqrtD;
     407           0 :   gel(y,1) = gen_1;
     408           0 :   if (mod2(S->D) != mod2(y2)) y2 = subiu(y,1);
     409           0 :   gel(y,2) = y2; av = avma;
     410           0 :   gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(y2), S->D),-2));
     411           0 : }
     412             : static GEN
     413           0 : qfr5_1(struct qfr_data *S, long prec)
     414             : {
     415           0 :   GEN y = cgetg(6, t_VEC);
     416           0 :   qfr_1_fill(y, S);
     417           0 :   gel(y,4) = gen_0;
     418           0 :   gel(y,5) = real_1(prec); return y;
     419             : }
     420             : static GEN
     421           0 : qfr3_1(struct qfr_data *S)
     422             : {
     423           0 :   GEN y = cgetg(4, t_VEC);
     424           0 :   qfr_1_fill(y, S); return y;
     425             : }
     426             : 
     427             : /* Assume D < 0 is the discriminant of a t_QFB */
     428             : static GEN
     429      654897 : qfi_1_by_disc(GEN D)
     430             : {
     431      654897 :   GEN b,c, y = cgetg(5,t_QFB);
     432      654897 :   quadpoly_bc(D, mod2(D), &b,&c);
     433      654897 :   if (b == gen_m1) b = gen_1;
     434      654897 :   gel(y,1) = gen_1;
     435      654897 :   gel(y,2) = b;
     436      654897 :   gel(y,3) = c;
     437      654897 :   gel(y,4) = icopy(D); return y;
     438             : }
     439             : static GEN
     440      647288 : qfi_1(GEN x)
     441             : {
     442      647288 :   if (typ(x) != t_QFB) pari_err_TYPE("qfi_1",x);
     443      647288 :   return qfi_1_by_disc(qfb_disc(x));
     444             : }
     445             : 
     446             : GEN
     447           0 : qfb_1(GEN x) { return qfb_is_qfi(x) ? qfi_1(x): qfr_1(x); }
     448             : 
     449             : static GEN
     450     8711765 : _qfimul(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,0); }
     451             : static GEN
     452    22976882 : _qfisqr(void *E, GEN x) { (void) E; return qficomp0(x,x,0); }
     453             : static GEN
     454           7 : _qfimulraw(void *E, GEN x, GEN y) { (void) E; return qficomp0(x,y,1); }
     455             : static GEN
     456           7 : _qfisqrraw(void *E, GEN x) { (void) E; return qficomp0(x,x,1); }
     457             : 
     458             : static GEN
     459           7 : qfipowraw(GEN x, long n)
     460             : {
     461           7 :   pari_sp av = avma;
     462             :   GEN y;
     463           7 :   if (!n) return qfi_1(x);
     464           7 :   if (n== 1) return gcopy(x);
     465           7 :   if (n==-1) { x = gcopy(x); togglesign(gel(x,2)); return x; }
     466           7 :   if (n < 0) x = qfb_inv(x);
     467           7 :   y = gen_powu(x, labs(n), NULL, &_qfisqrraw, &_qfimulraw);
     468           7 :   return gerepilecopy(av,y);
     469             : }
     470             : 
     471             : static GEN
     472    10460310 : qfipow(GEN x, GEN n)
     473             : {
     474    10460310 :   pari_sp av = avma;
     475             :   GEN y;
     476    10460310 :   long s = signe(n);
     477    10460310 :   if (!s) return qfi_1(x);
     478     9813022 :   if (s < 0) x = qfb_inv(x);
     479     9813020 :   y = gen_pow(qfbred_i(x), n, NULL, &_qfisqr, &_qfimul);
     480     9813026 :   return gerepilecopy(av,y);
     481             : }
     482             : 
     483             : static long
     484      412328 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
     485             : {
     486             :   long z;
     487      412328 :   *v = gen_0; *v2 = gen_1;
     488     4351417 :   for (z=0; abscmpii(*v3,L) > 0; z++)
     489             :   {
     490     3939089 :     GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
     491     3939089 :     *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
     492             :   }
     493      412328 :   return z;
     494             : }
     495             : 
     496             : /* composition: Shanks' NUCOMP & NUDUPL */
     497             : /* L = floor((|d|/4)^(1/4)) */
     498             : GEN
     499      400722 : nucomp(GEN x, GEN y, GEN L)
     500             : {
     501      400722 :   pari_sp av = avma;
     502             :   long z;
     503             :   GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
     504             : 
     505      400722 :   if (x==y) return nudupl(x,L);
     506      400680 :   if (!is_qfi(x)) pari_err_TYPE("nucomp",x);
     507      400680 :   if (!is_qfi(y)) pari_err_TYPE("nucomp",y);
     508             : 
     509      400680 :   if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
     510      400680 :   s = shifti(addii(gel(x,2),gel(y,2)), -1);
     511      400680 :   n = subii(gel(y,2), s);
     512      400680 :   a1 = gel(x,1);
     513      400680 :   a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
     514      400680 :   if (equali1(d)) { a = negi(mulii(u,n)); d1 = d; }
     515      163576 :   else if (dvdii(s,d)) /* d | s */
     516             :   {
     517       83503 :     a = negi(mulii(u,n)); d1 = d;
     518       83503 :     a1 = diviiexact(a1, d1);
     519       83503 :     a2 = diviiexact(a2, d1);
     520       83503 :     s = diviiexact(s, d1);
     521             :   }
     522             :   else
     523             :   {
     524             :     GEN p2, l;
     525       80073 :     d1 = bezout(s,d,&u1,NULL);
     526       80073 :     if (!equali1(d1))
     527             :     {
     528        2044 :       a1 = diviiexact(a1,d1);
     529        2044 :       a2 = diviiexact(a2,d1);
     530        2044 :       s = diviiexact(s,d1);
     531        2044 :       d = diviiexact(d,d1);
     532             :     }
     533       80073 :     p1 = remii(gel(x,3),d);
     534       80073 :     p2 = remii(gel(y,3),d);
     535       80073 :     l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
     536       80073 :     a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
     537             :   }
     538      400680 :   a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
     539      400680 :   d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
     540      400680 :   Q = cgetg(5,t_QFB);
     541      400680 :   if (!z) {
     542       37632 :     g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
     543       37632 :     b = a2;
     544       37632 :     b2 = gel(y,2);
     545       37632 :     v2 = d1;
     546       37632 :     gel(Q,1) = mulii(d,b);
     547             :   } else {
     548             :     GEN e, q3, q4;
     549      363048 :     if (z&1) { v3 = negi(v3); v2 = negi(v2); }
     550      363048 :     b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
     551      363048 :     e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
     552      363048 :     q3 = mulii(e,v2);
     553      363048 :     q4 = subii(q3,s);
     554      363048 :     b2 = addii(q3,q4);
     555      363048 :     g = diviiexact(q4,v);
     556      363048 :     if (!equali1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
     557      363048 :     gel(Q,1) = addii(mulii(d,b), mulii(e,v));
     558             :   }
     559      400680 :   q1 = mulii(b, v3);
     560      400680 :   q2 = addii(q1,n);
     561      400680 :   gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
     562      400680 :   gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
     563      400680 :   gel(Q,4) = gel(x,4);
     564      400680 :   return redimag_av(av, Q);
     565             : }
     566             : 
     567             : GEN
     568       11648 : nudupl(GEN x, GEN L)
     569             : {
     570       11648 :   pari_sp av = avma;
     571             :   long z;
     572             :   GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
     573             : 
     574       11648 :   if (!is_qfi(x)) pari_err_TYPE("nudupl",x);
     575       11648 :   a = gel(x,1);
     576       11648 :   b = gel(x,2);
     577       11648 :   d1 = bezout(b,a, &u,NULL);
     578       11648 :   if (!equali1(d1))
     579             :   {
     580        4620 :     a = diviiexact(a, d1);
     581        4620 :     b = diviiexact(b, d1);
     582             :   }
     583       11648 :   c = modii(negi(mulii(u,gel(x,3))), a);
     584       11648 :   p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
     585       11648 :   d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
     586       11648 :   a2 = sqri(d);
     587       11648 :   c2 = sqri(v3);
     588       11648 :   Q = cgetg(5,t_QFB);
     589       11648 :   if (!z) {
     590        1281 :     g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
     591        1281 :     b2 = gel(x,2);
     592        1281 :     v2 = d1;
     593        1281 :     gel(Q,1) = a2;
     594             :   } else {
     595             :     GEN e;
     596       10367 :     if (z&1) { v = negi(v); d = negi(d); }
     597       10367 :     e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
     598       10367 :     g = diviiexact(subii(mulii(e,v2), b), v);
     599       10367 :     b2 = addii(mulii(e,v2), mulii(v,g));
     600       10367 :     if (!equali1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
     601       10367 :     gel(Q,1) = addii(a2, mulii(e,v));
     602             :   }
     603       11648 :   gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
     604       11648 :   gel(Q,3) = addii(c2, mulii(g,v2));
     605       11648 :   gel(Q,4) = gel(x,4);
     606       11648 :   return redimag_av(av, Q);
     607             : }
     608             : 
     609             : static GEN
     610        4739 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
     611             : static GEN
     612       11606 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
     613             : GEN
     614        1008 : nupow(GEN x, GEN n, GEN L)
     615             : {
     616             :   pari_sp av;
     617             :   GEN y, D;
     618             : 
     619        1008 :   if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
     620        1008 :   if (!is_qfi(x)) pari_err_TYPE("nupow",x);
     621        1008 :   if (gequal1(n)) return gcopy(x);
     622        1008 :   av = avma;
     623        1008 :   D = qfb_disc(x);
     624        1008 :   y = qfi_1_by_disc(D);
     625        1008 :   if (!signe(n)) return y;
     626         959 :   if (!L) L = sqrtnint(absi_shallow(D), 4);
     627         959 :   y = gen_pow_i(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
     628         959 :   if (signe(n) < 0
     629          35 :   && !absequalii(gel(y,1),gel(y,2))
     630          35 :   && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
     631         959 :   return gerepilecopy(av, y);
     632             : }
     633             : 
     634             : /* Reduction */
     635             : 
     636             : /* assume a > 0. Write b = q*2a + r, with -a < r <= a */
     637             : static GEN
     638     8602616 : dvmdii_round(GEN b, GEN a, GEN *r)
     639             : {
     640     8602616 :   GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
     641     8602405 :   if (signe(b) >= 0) {
     642     4945559 :     if (abscmpii(*r, a) > 0) { q = addiu(q, 1); *r = subii(*r, a2); }
     643             :   } else { /* r <= 0 */
     644     3656846 :     if (abscmpii(*r, a) >= 0){ q = subiu(q, 1); *r = addii(*r, a2); }
     645             :   }
     646     8602074 :   return q;
     647             : }
     648             : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
     649             : static long
     650    93622805 : dvmdsu_round(long b, ulong a, long *r)
     651             : {
     652    93622805 :   ulong a2 = a << 1, q, ub, ur;
     653    93622805 :   if (b >= 0) {
     654    59405319 :     ub = b;
     655    59405319 :     q = ub / a2;
     656    59405319 :     ur = ub % a2;
     657    59405319 :     if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
     658    21045070 :     else *r = (long)ur;
     659    59405319 :     return (long)q;
     660             :   } else { /* r <= 0 */
     661    34217486 :     ub = (ulong)-b; /* |b| */
     662    34217486 :     q = ub / a2;
     663    34217486 :     ur = ub % a2;
     664    34217486 :     if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
     665    18752810 :     else *r = -(long)ur;
     666    34217486 :     return -(long)q;
     667             :   }
     668             : }
     669             : /* reduce b mod 2*a. Update b,c */
     670             : static void
     671     2768824 : REDB(GEN a, GEN *b, GEN *c)
     672             : {
     673     2768824 :   GEN r, q = dvmdii_round(*b, a, &r);
     674     2768824 :   if (!signe(q)) return;
     675     2699896 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     676     2699896 :   *b = r;
     677             : }
     678             : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
     679             : static void
     680    93622851 : sREDB(ulong a, long *b, ulong *c)
     681             : {
     682             :   long r, q;
     683             :   ulong uz;
     684   101509618 :   if (a > LONG_MAX) return; /* b already reduced */
     685    93622851 :   q = dvmdsu_round(*b, a, &r);
     686    93641225 :   if (q == 0) return;
     687             :   /* Final (a,r,c2) satisfies |r| <= |b| hence c2 <= c, c2 = c - q*z,
     688             :    * where z = (b+r) / 2, representable as long, has the same sign as q. */
     689    85754458 :   if (*b < 0)
     690             :   { /* uz = -z >= 0, q < 0 */
     691    30276763 :     if (r >= 0) /* different signs=>no overflow, exact division */
     692    15511840 :       uz = (ulong)-((*b + r)>>1);
     693             :     else
     694             :     {
     695    14764923 :       ulong ub = (ulong)-*b, ur = (ulong)-r;
     696    14764923 :       uz = (ub + ur) >> 1;
     697             :     }
     698    30276763 :     *c -= (-q) * uz; /* c -= qz */
     699             :   }
     700             :   else
     701             :   { /* uz = z >= 0, q > 0 */
     702    55477695 :     if (r <= 0)
     703    38418692 :       uz = (*b + r)>>1;
     704             :     else
     705             :     {
     706    17059003 :       ulong ub = (ulong)*b, ur = (ulong)r;
     707    17059003 :       uz = ((ub + ur) >> 1);
     708             :     }
     709    55477695 :     *c -= q * uz; /* c -= qz */
     710             :   }
     711    85754458 :   *b = r;
     712             : }
     713             : static void
     714     5833802 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
     715             : { /* REDB(a,b,c) */
     716     5833802 :   GEN r, q = dvmdii_round(*b, a, &r);
     717     5833251 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     718     5833157 :   *b = r;
     719     5833157 :   *u2 = subii(*u2, mulii(q, u1));
     720     5833269 : }
     721             : 
     722             : /* q t_QFB, return reduced representative and set base change U in Sl2(Z) */
     723             : GEN
     724     2159964 : redimagsl2(GEN q, GEN *U)
     725             : {
     726     2159964 :   pari_sp av = avma;
     727             :   GEN z, u1,u2,v1,v2,Q;
     728     2159964 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     729             :   long cmp;
     730     2159964 :   u1 = gen_1; u2 = gen_0;
     731     2159964 :   cmp = abscmpii(a, b);
     732     2159973 :   if (cmp < 0)
     733     1000255 :     REDBU(a,&b,&c, u1,&u2);
     734     1159718 :   else if (cmp == 0 && signe(b) < 0)
     735             :   { /* b = -a */
     736          13 :     b = negi(b);
     737          13 :     u2 = gen_1;
     738             :   }
     739             :   for(;;)
     740             :   {
     741     6993282 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     742     4832895 :     swap(a,c); b = negi(b);
     743     4833602 :     z = u1; u1 = u2; u2 = negi(z);
     744     4833599 :     REDBU(a,&b,&c, u1,&u2);
     745     4833321 :     if (gc_needed(av, 1)) {
     746           0 :       if (DEBUGMEM>1) pari_warn(warnmem, "redimagsl2");
     747           0 :       gerepileall(av, 5, &a,&b,&c, &u1,&u2);
     748             :     }
     749             :   }
     750     2159871 :   if (cmp == 0 && signe(b) < 0)
     751             :   {
     752       16919 :     b = negi(b);
     753       16919 :     z = u1; u1 = u2; u2 = negi(z);
     754             :   }
     755             :   /* Let q = (A,B,C). q o [u1,u2; v1,v2] = Q implies
     756             :    * [v1,v2] = (1/C) [(b-B)/2 u1 - a u2, c u1 - (b+B)/2 u2] */
     757     2159871 :   z = shifti(subii(b, gel(q,2)), -1);
     758     2159794 :   v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
     759     2159720 :   z = subii(z, b);
     760     2159807 :   v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
     761     2159619 :   *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2));
     762     2159912 :   Q = lg(q)==5 ? mkqfb(a,b,c,gel(q,4)): mkvec3(a,b,c);
     763     2159902 :   return gc_all(av, 2, &Q, U);
     764             : }
     765             : 
     766             : static GEN
     767     1057615 : setq_b0(ulong a, ulong c, GEN D)
     768     1057615 : { retmkqfb(utoipos(a), gen_0, utoipos(c), icopy(D)); }
     769             : /* assume |sb| = 1 */
     770             : static GEN
     771    75092828 : setq(ulong a, ulong b, ulong c, long sb, GEN D)
     772    75092828 : { retmkqfb(utoipos(a), sb==1? utoipos(b): utoineg(b), utoipos(c), icopy(D)); }
     773             : /* 0 < a, c < 2^BIL, b = 0 */
     774             : static GEN
     775      956997 : redimag_1_b0(ulong a, ulong c, GEN D)
     776      956997 : { return (a <= c)? setq_b0(a, c, D): setq_b0(c, a, D); }
     777             : 
     778             : /* 0 < a, c < 2^BIL: single word affair */
     779             : static GEN
     780    76242253 : redimag_1(pari_sp av, GEN a, GEN b, GEN c, GEN D)
     781             : {
     782             :   ulong ua, ub, uc;
     783             :   long sb;
     784             :   for(;;)
     785      140874 :   { /* at most twice */
     786    76242253 :     long lb = lgefint(b); /* <= 3 after first loop */
     787    76242253 :     if (lb == 2) return redimag_1_b0(a[2],c[2], D);
     788    75285256 :     if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
     789      139339 :     REDB(a,&b,&c);
     790      140142 :     if (uel(a,2) <= uel(c,2))
     791             :     { /* lg(b) <= 3 but may be too large for itos */
     792           0 :       long s = signe(b);
     793           0 :       set_avma(av);
     794           0 :       if (!s) return redimag_1_b0(a[2], c[2], D);
     795           0 :       if (a[2] == c[2]) s = 1;
     796           0 :       return setq(a[2], b[2], c[2], s, D);
     797             :     }
     798      140142 :     swap(a,c); b = negi(b);
     799             :   }
     800             :   /* b != 0 */
     801    75145917 :   set_avma(av);
     802    75156932 :   ua = a[2];
     803    75156932 :   ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
     804    75156932 :   uc = c[2];
     805    75156932 :   if (ua < ub)
     806    27366541 :     sREDB(ua, &sb, &uc);
     807    47790391 :   else if (ua == ub && sb < 0) sb = (long)ub;
     808   141489499 :   while(ua > uc)
     809             :   {
     810    66290246 :     lswap(ua,uc); sb = -sb;
     811    66290246 :     sREDB(ua, &sb, &uc);
     812             :   }
     813    75199253 :   if (!sb) return setq_b0(ua, uc, D);
     814             :   else
     815             :   {
     816    75098634 :     long s = 1;
     817    75098634 :     if (sb < 0)
     818             :     {
     819    28332231 :       sb = -sb;
     820    28332231 :       if (ua != uc) s = -1;
     821             :     }
     822    75098634 :     return setq(ua, sb, uc, s, D);
     823             :   }
     824             : }
     825             : 
     826             : static GEN
     827    76047599 : redimag_av(pari_sp av, GEN q)
     828             : {
     829    76047599 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3), D = gel(q,4);
     830    76047599 :   long cmp, lc = lgefint(c);
     831             : 
     832    76047599 :   if (lgefint(a) == 3 && lc == 3) return redimag_1(av, a, b, c, D);
     833      906428 :   cmp = abscmpii(a, b);
     834      914075 :   if (cmp < 0)
     835      434190 :     REDB(a,&b,&c);
     836      479885 :   else if (cmp == 0 && signe(b) < 0)
     837          30 :     b = negi(b);
     838             :   for(;;)
     839             :   {
     840     3108567 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     841     2918664 :     lc = lgefint(a); /* lg(future c): we swap a & c next */
     842     2918664 :     if (lc == 3) return redimag_1(av, a, b, c, D);
     843     2194492 :     swap(a,c); b = negi(b); /* apply rho */
     844     2194492 :     REDB(a,&b,&c);
     845             :   }
     846      189903 :   if (cmp == 0 && signe(b) < 0) b = negi(b);
     847      189903 :   return gerepilecopy(av, mkqfb(a, b, c, D));
     848             : }
     849             : static GEN
     850    16709859 : redimag(GEN q) { return redimag_av(avma, q); }
     851             : 
     852             : static GEN
     853           7 : rhoimag(GEN x)
     854             : {
     855           7 :   pari_sp av = avma;
     856           7 :   GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
     857           7 :   int fl = abscmpii(a, c);
     858           7 :   if (fl <= 0)
     859             :   {
     860           7 :     int fg = abscmpii(a, b);
     861           7 :     if (fg >= 0)
     862             :     {
     863           7 :       x = gcopy(x);
     864           7 :       if ((!fl || !fg) && signe(gel(x,2)) < 0) setsigne(gel(x,2), 1);
     865           7 :       return x;
     866             :     }
     867             :   }
     868           0 :   swap(a,c); b = negi(b);
     869           0 :   REDB(a, &b, &c);
     870           0 :   return gerepilecopy(av, mkqfb(a,b,c, qfb_disc(x)));
     871             : }
     872             : 
     873             : /* qfr3 / qfr5 */
     874             : 
     875             : /* t_QFB are unusable: D, sqrtD, isqrtD are recomputed all the time and the
     876             :  * logarithmic Shanks's distance is costly and hard to control.
     877             :  * qfr3 / qfr5 routines take a container of t_INTs (e.g a t_VEC) as argument,
     878             :  * at least 3 (resp. 5) components [it is a feature that they do not check the
     879             :  * precise type or length of the input]. They return a vector of length 3
     880             :  * (resp. 5). A qfr3 [a,b,c] contains the form coeffs, in a qfr5 [a,b,c, e,d]
     881             :  * the t_INT e is a binary exponent, d a t_REAL, coding the distance in
     882             :  * multiplicative form: the true distance is obtained from qfr5_dist.
     883             :  * All other qfr routines are obsolete (inefficient) wrappers */
     884             : 
     885             : /* static functions are not stack-clean. Unless mentionned otherwise, public
     886             :  * functions are. */
     887             : 
     888             : #define EMAX 22
     889             : static void
     890    12022731 : fix_expo(GEN x)
     891             : {
     892    12022731 :   if (expo(gel(x,5)) >= (1L << EMAX)) {
     893           0 :     gel(x,4) = addiu(gel(x,4), 1);
     894           0 :     shiftr_inplace(gel(x,5), - (1L << EMAX));
     895             :   }
     896    12022731 : }
     897             : 
     898             : /* (1/2) log (d * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
     899             : GEN
     900      184254 : qfr5_dist(GEN e, GEN d, long prec)
     901             : {
     902      184254 :   GEN t = logr_abs(d);
     903      184254 :   if (signe(e)) {
     904           0 :     GEN u = mulir(e, mplog2(prec));
     905           0 :     shiftr_inplace(u, EMAX); t = addrr(t, u);
     906             :   }
     907      184254 :   shiftr_inplace(t, -1); return t;
     908             : }
     909             : 
     910             : static void
     911    17290498 : rho_get_BC(GEN *B, GEN *C, GEN b, GEN c, struct qfr_data *S)
     912             : {
     913             :   GEN t, u;
     914    17290498 :   u = shifti(c,1);
     915    17290498 :   t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: c;
     916    17290498 :   u = remii(addii_sign(t,1, b,signe(b)), u);
     917    17290498 :   *B = addii_sign(t, 1, u, -signe(u)); /* |t| - (|t|+b) % |2c| */
     918    17290498 :   if (*B == gen_0)
     919        3406 :   { u = shifti(S->D, -2); setsigne(u, -1); }
     920             :   else
     921    17287092 :     u = shifti(addii_sign(sqri(*B),1, S->D,-1), -2);
     922    17290498 :   *C = diviiexact(u, c); /* = (B^2-D)/4c */
     923    17290498 : }
     924             : /* Not stack-clean */
     925             : GEN
     926     6994478 : qfr3_rho(GEN x, struct qfr_data *S)
     927             : {
     928     6994478 :   GEN B, C, b = gel(x,2), c = gel(x,3);
     929     6994478 :   rho_get_BC(&B, &C, b, c, S);
     930     6994478 :   return mkvec3(c,B,C);
     931             : }
     932             : /* Not stack-clean */
     933             : GEN
     934    10296020 : qfr5_rho(GEN x, struct qfr_data *S)
     935             : {
     936    10296020 :   GEN B, C, y, b = gel(x,2), c = gel(x,3);
     937    10296020 :   long sb = signe(b);
     938             : 
     939    10296020 :   rho_get_BC(&B, &C, b, c, S);
     940    10296020 :   y = mkvec5(c,B,C, gel(x,4), gel(x,5));
     941    10296020 :   if (sb) {
     942    10292324 :     GEN t = subii(sqri(b), S->D);
     943    10292324 :     if (sb < 0)
     944     2299122 :       t = divir(t, sqrr(subir(b,S->sqrtD)));
     945             :     else
     946     7993202 :       t = divri(sqrr(addir(b,S->sqrtD)), t);
     947             :     /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
     948    10292324 :     gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
     949             :   }
     950    10296020 :   return y;
     951             : }
     952             : 
     953             : /* Not stack-clean */
     954             : GEN
     955      217070 : qfr_to_qfr5(GEN x, long prec)
     956      217070 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
     957             : 
     958             : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
     959             : GEN
     960         448 : qfr5_to_qfr(GEN x, GEN D, GEN d0)
     961             : {
     962         448 :   if (d0)
     963             :   {
     964         126 :     GEN n = gel(x,4), d = absr(gel(x,5));
     965         126 :     if (signe(n))
     966             :     {
     967           0 :       n = addis(shifti(n, EMAX), expo(d));
     968           0 :       setexpo(d, 0); d = logr_abs(d);
     969           0 :       if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
     970           0 :       shiftr_inplace(d, -1);
     971           0 :       d0 = addrr(d0, d);
     972             :     }
     973         126 :     else if (!gequal1(d)) /* avoid loss of precision */
     974             :     {
     975          91 :       d = logr_abs(d);
     976          91 :       shiftr_inplace(d, -1);
     977          91 :       d0 = addrr(d0, d);
     978             :     }
     979             :   }
     980         448 :   x = qfr3_to_qfr(x, D);
     981         448 :   return d0 ? mkvec2(x,d0): x;
     982             : }
     983             : 
     984             : /* Not stack-clean */
     985             : GEN
     986       31899 : qfr3_to_qfr(GEN x, GEN d) { retmkqfb(gel(x,1), gel(x,2), gel(x,3), d); }
     987             : 
     988             : static int
     989    20823927 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
     990             : {
     991    20823927 :   if (signe(b) > 0 && abscmpii(b, isqrtD) <= 0)
     992             :   {
     993     5505658 :     GEN t = addii_sign(isqrtD,1, shifti(a,1),-1);
     994     5505658 :     long l = abscmpii(b, t); /* compare |b| and |floor(sqrt(D)) - |2a|| */
     995     5505658 :     if (l > 0 || (l == 0 && signe(t) < 0)) return 1;
     996             :   }
     997    16726977 :   return 0;
     998             : }
     999             : 
    1000             : INLINE int
    1001    17275760 : qfr_isreduced(GEN x, GEN isqrtD)
    1002             : {
    1003    17275760 :   return ab_isreduced(gel(x,1),gel(x,2),isqrtD);
    1004             : }
    1005             : 
    1006             : /* Not stack-clean */
    1007             : GEN
    1008     1947414 : qfr5_red(GEN x, struct qfr_data *S)
    1009             : {
    1010     1947414 :   pari_sp av = avma;
    1011    10247839 :   while (!qfr_isreduced(x, S->isqrtD))
    1012             :   {
    1013     8300425 :     x = qfr5_rho(x, S);
    1014     8300425 :     if (gc_needed(av,2))
    1015             :     {
    1016           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr5_red");
    1017           0 :       x = gerepilecopy(av, x);
    1018             :     }
    1019             :   }
    1020     1947414 :   return x;
    1021             : }
    1022             : /* Not stack-clean */
    1023             : GEN
    1024     1169914 : qfr3_red(GEN x, struct qfr_data *S)
    1025             : {
    1026     1169914 :   pari_sp av = avma;
    1027     7027921 :   while (!qfr_isreduced(x, S->isqrtD))
    1028             :   {
    1029     5858007 :     x = qfr3_rho(x, S);
    1030     5858007 :     if (gc_needed(av,2))
    1031             :     {
    1032           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"qfr3_red");
    1033           0 :       x = gerepilecopy(av, x);
    1034             :     }
    1035             :   }
    1036     1169914 :   return x;
    1037             : }
    1038             : 
    1039             : void
    1040        2184 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
    1041             : {
    1042        2184 :   S->D = D;
    1043        2184 :   S->sqrtD = sqrtr(itor(S->D,prec));
    1044        2184 :   S->isqrtD = truncr(S->sqrtD);
    1045        2184 : }
    1046             : 
    1047             : static GEN
    1048         126 : qfr5_init(GEN x, GEN d, struct qfr_data *S)
    1049             : {
    1050         126 :   long prec = realprec(d), l = -expo(d);
    1051         126 :   if (l < BITS_IN_LONG) l = BITS_IN_LONG;
    1052         126 :   prec = maxss(prec, nbits2prec(l));
    1053         126 :   S->D = qfb_disc(x);
    1054         126 :   x = qfr_to_qfr5(x,prec);
    1055         126 :   if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
    1056           0 :   else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
    1057             : 
    1058         126 :   if (!S->isqrtD)
    1059             :   {
    1060         112 :     pari_sp av=avma;
    1061             :     long e;
    1062         112 :     S->isqrtD = gcvtoi(S->sqrtD,&e);
    1063         112 :     if (e>-2) { set_avma(av); S->isqrtD = sqrti(S->D); }
    1064             :   }
    1065          14 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
    1066         126 :   return x;
    1067             : }
    1068             : static GEN
    1069         364 : qfr3_init(GEN x, struct qfr_data *S)
    1070             : {
    1071         364 :   S->D = qfb_disc(x);
    1072         364 :   if (!S->isqrtD) S->isqrtD = sqrti(S->D);
    1073         280 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
    1074         364 :   return x;
    1075             : }
    1076             : 
    1077             : #define qf_NOD  2
    1078             : #define qf_STEP 1
    1079             : 
    1080             : static GEN
    1081         392 : redreal_i(GEN x, long flag, GEN isqrtD, GEN sqrtD)
    1082             : {
    1083             :   struct qfr_data S;
    1084         392 :   GEN d = NULL, y;
    1085         392 :   if (typ(x)==t_VEC) { d = gel(x,2); x = gel(x,1); } else flag |= qf_NOD;
    1086         392 :   S.sqrtD = sqrtD;
    1087         392 :   S.isqrtD = isqrtD;
    1088         392 :   y = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, d, &S);
    1089         392 :   switch(flag) {
    1090          49 :     case 0:              y = qfr5_red(y,&S); break;
    1091         315 :     case qf_NOD:         y = qfr3_red(y,&S); break;
    1092          21 :     case qf_STEP:        y = qfr5_rho(y,&S); break;
    1093           7 :     case qf_STEP|qf_NOD: y = qfr3_rho(y,&S); break;
    1094           0 :     default: pari_err_FLAG("qfbred");
    1095             :   }
    1096         392 :   return qfr5_to_qfr(y, qfb_disc(x), d);
    1097             : }
    1098             : static GEN
    1099          42 : redreal(GEN x) { return redreal_i(x,0,NULL,NULL); }
    1100             : 
    1101             : GEN
    1102       53993 : qfbred0(GEN x, long flag, GEN isqrtD, GEN sqrtD)
    1103             : {
    1104             :   pari_sp av;
    1105       53993 :   GEN q = check_qfbext("qfbred",x);
    1106       53993 :   if (qfb_is_qfi(q)) return (flag & qf_STEP)? rhoimag(x): redimag(x);
    1107         350 :   if (typ(x)==t_QFB) flag |= qf_NOD;
    1108          42 :   else               flag &= ~qf_NOD;
    1109         350 :   av = avma;
    1110         350 :   return gerepilecopy(av, redreal_i(x,flag,isqrtD,sqrtD));
    1111             : }
    1112             : /* t_QFB */
    1113             : GEN
    1114    16656223 : qfbred_i(GEN x) { return qfb_is_qfi(x)? redimag(x): redreal(x); }
    1115             : GEN
    1116       52201 : qfbred(GEN x) { return qfbred0(x, 0, NULL, NULL); }
    1117             : 
    1118             : /* Not stack-clean */
    1119             : GEN
    1120     1730407 : qfr5_compraw(GEN x, GEN y)
    1121             : {
    1122     1730407 :   GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
    1123     1730407 :   if (x == y)
    1124             :   {
    1125       34097 :     gel(z,4) = shifti(gel(x,4),1);
    1126       34097 :     gel(z,5) = sqrr(gel(x,5));
    1127             :   }
    1128             :   else
    1129             :   {
    1130     1696310 :     gel(z,4) = addii(gel(x,4),gel(y,4));
    1131     1696310 :     gel(z,5) = mulrr(gel(x,5),gel(y,5));
    1132             :   }
    1133     1730407 :   fix_expo(z); return z;
    1134             : }
    1135             : GEN
    1136     1730393 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
    1137     1730393 : { return qfr5_red(qfr5_compraw(x, y), S); }
    1138             : /* Not stack-clean */
    1139             : GEN
    1140     1006919 : qfr3_compraw(GEN x, GEN y)
    1141             : {
    1142     1006919 :   GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
    1143     1006919 :   return z;
    1144             : }
    1145             : GEN
    1146     1006919 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
    1147     1006919 : { return qfr3_red(qfr3_compraw(x,y), S); }
    1148             : 
    1149             : /* m > 0. Not stack-clean */
    1150             : static GEN
    1151           7 : qfr5_powraw(GEN x, long m)
    1152             : {
    1153           7 :   GEN y = NULL;
    1154          14 :   for (; m; m >>= 1)
    1155             :   {
    1156          14 :     if (m&1) y = y? qfr5_compraw(y,x): x;
    1157          14 :     if (m == 1) break;
    1158           7 :     x = qfr5_compraw(x,x);
    1159             :   }
    1160           7 :   return y;
    1161             : }
    1162             : 
    1163             : /* return x^n. Not stack-clean */
    1164             : GEN
    1165          21 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
    1166             : {
    1167          21 :   GEN y = NULL;
    1168          21 :   long i, m, s = signe(n);
    1169          21 :   if (!s) return qfr5_1(S, lg(gel(x,5)));
    1170          21 :   if (s < 0) x = qfb_inv(x);
    1171          42 :   for (i=lgefint(n)-1; i>1; i--)
    1172             :   {
    1173          21 :     m = n[i];
    1174          56 :     for (; m; m>>=1)
    1175             :     {
    1176          56 :       if (m&1) y = y? qfr5_comp(y,x,S): x;
    1177          56 :       if (m == 1 && i == 2) break;
    1178          35 :       x = qfr5_comp(x,x,S);
    1179             :     }
    1180             :   }
    1181          21 :   return y;
    1182             : }
    1183             : /* m > 0; return x^m. Not stack-clean */
    1184             : static GEN
    1185           0 : qfr3_powraw(GEN x, long m)
    1186             : {
    1187           0 :   GEN y = NULL;
    1188           0 :   for (; m; m>>=1)
    1189             :   {
    1190           0 :     if (m&1) y = y? qfr3_compraw(y,x): x;
    1191           0 :     if (m == 1) break;
    1192           0 :     x = qfr3_compraw(x,x);
    1193             :   }
    1194           0 :   return y;
    1195             : }
    1196             : /* return x^n. Not stack-clean */
    1197             : GEN
    1198        4529 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
    1199             : {
    1200        4529 :   GEN y = NULL;
    1201        4529 :   long i, m, s = signe(n);
    1202        4529 :   if (!s) return qfr3_1(S);
    1203        4529 :   if (s < 0) x = qfb_inv(x);
    1204        9074 :   for (i=lgefint(n)-1; i>1; i--)
    1205             :   {
    1206        4545 :     m = n[i];
    1207        5137 :     for (; m; m>>=1)
    1208             :     {
    1209        5117 :       if (m&1) y = y? qfr3_comp(y,x,S): x;
    1210        5117 :       if (m == 1 && i == 2) break;
    1211         592 :       x = qfr3_comp(x,x,S);
    1212             :     }
    1213             :   }
    1214        4529 :   return y;
    1215             : }
    1216             : 
    1217             : static GEN
    1218           7 : qfrinvraw(GEN x)
    1219             : {
    1220           7 :   if (typ(x) == t_VEC) retmkvec2(qfbinv(gel(x,1)), negr(gel(x,2)));
    1221           7 :  return qfbinv(x);
    1222             : }
    1223             : static GEN
    1224          14 : qfrpowraw(GEN x, long n)
    1225             : {
    1226          14 :   struct qfr_data S = { NULL, NULL, NULL };
    1227          14 :   pari_sp av = avma;
    1228          14 :   if (n==1) return gcopy(x);
    1229          14 :   if (n==-1) return qfrinvraw(x);
    1230           7 :   if (typ(x)==t_QFB)
    1231             :   {
    1232           0 :     GEN D = qfb_disc(x);
    1233           0 :     if (!n) return qfr_1(x);
    1234           0 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1235           0 :     x = qfr3_powraw(x, n);
    1236           0 :     x = qfr3_to_qfr(x, D);
    1237             :   }
    1238             :   else
    1239             :   {
    1240           7 :     GEN d0 = gel(x,2);
    1241           7 :     x = gel(x,1);
    1242           7 :     if (!n) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1243           7 :     if (n < 0) { x = qfb_inv(x); n = -n; }
    1244           7 :     x = qfr5_init(x, d0, &S);
    1245           7 :     if (labs(n) != 1) x = qfr5_powraw(x, n);
    1246           7 :     x = qfr5_to_qfr(x, S.D, mulrs(d0,n));
    1247             :   }
    1248           7 :   return gerepilecopy(av, x);
    1249             : }
    1250             : static GEN
    1251         119 : qfrpow(GEN x, GEN n)
    1252             : {
    1253         119 :   struct qfr_data S = { NULL, NULL, NULL };
    1254         119 :   long s = signe(n);
    1255         119 :   pari_sp av = avma;
    1256         119 :   if (typ(x)==t_QFB)
    1257             :   {
    1258          56 :     if (!s) return qfr_1(x);
    1259          42 :     if (s < 0) x = qfb_inv(x);
    1260          42 :     x = qfr3_init(x, &S);
    1261          42 :     x = is_pm1(n)? qfr3_red(x, &S): qfr3_pow(x, n, &S);
    1262          42 :     x = qfr3_to_qfr(x, S.D);
    1263             :   }
    1264             :   else
    1265             :   {
    1266          63 :     GEN d0 = gel(x,2);
    1267          63 :     x = gel(x,1);
    1268          63 :     if (!s) retmkvec2(qfr_1(x), real_0(precision(d0)));
    1269          49 :     if (s < 0) x = qfb_inv(x);
    1270          49 :     x = qfr5_init(x, d0, &S);
    1271          49 :     x = is_pm1(n)? qfr5_red(x, &S): qfr5_pow(x, n, &S);
    1272          49 :     x = qfr5_to_qfr(x, S.D, mulri(d0,n));
    1273             :   }
    1274          91 :   return gerepilecopy(av, x);
    1275             : }
    1276             : GEN
    1277          21 : qfbpowraw(GEN x, long n)
    1278             : {
    1279          21 :   GEN q = check_qfbext("qfbpowraw",x);
    1280          21 :   return qfb_is_qfi(q)? qfipowraw(x,n): qfrpowraw(x,n);
    1281             : }
    1282             : /* same discriminant, no distance, no checks */
    1283             : GEN
    1284     9210069 : qfbpow_i(GEN x, GEN n) { return qfb_is_qfi(x)? qfipow(x,n): qfrpow(x,n); }
    1285             : GEN
    1286     1250362 : qfbpow(GEN x, GEN n)
    1287             : {
    1288     1250362 :   GEN q = check_qfbext("qfbpow",x);
    1289     1250361 :   return qfb_is_qfi(q)? qfipow(x,n): qfrpow(x,n);
    1290             : }
    1291             : GEN
    1292     1151174 : qfbpows(GEN x, long n)
    1293             : {
    1294     1151174 :   long N[] = { evaltyp(t_INT) | _evallg(3), 0, 0};
    1295     1151174 :   affsi(n, N); return qfbpow(x, N);
    1296             : }
    1297             : 
    1298             : /* Prime forms attached to prime ideals of degree 1 */
    1299             : 
    1300             : /* assume x != 0 a t_INT, p > 0
    1301             :  * Return a t_QFB, but discriminant sign is not checked: can be used for
    1302             :  * real forms as well */
    1303             : GEN
    1304    11775505 : primeform_u(GEN x, ulong p)
    1305             : {
    1306    11775505 :   GEN c, y = cgetg(5, t_QFB);
    1307    11775519 :   pari_sp av = avma;
    1308             :   ulong b;
    1309             :   long s;
    1310             : 
    1311    11775519 :   s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
    1312             :   /* 2 or 3 mod 4 */
    1313    11775503 :   if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
    1314    11775925 :   if (p == 2) {
    1315     3580388 :     switch(s) {
    1316      588811 :       case 0: b = 0; break;
    1317     2640924 :       case 1: b = 1; break;
    1318      350654 :       case 4: b = 2; break;
    1319           0 :       default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1320           0 :                b = 0; /* -Wall */
    1321             :     }
    1322     3580389 :     c = shifti(subsi(s,x), -3);
    1323             :   } else {
    1324     8195537 :     b = Fl_sqrt(umodiu(x,p), p);
    1325     8198292 :     if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1326             :     /* mod(b) != mod2(x) ? */
    1327     8198529 :     if ((b ^ s) & 1) b = p - b;
    1328     8198529 :     c = diviuexact(shifti(subii(sqru(b), x), -2), p);
    1329             :   }
    1330    11760382 :   gel(y,3) = gerepileuptoint(av, c);
    1331    11775891 :   gel(y,4) = icopy(x);
    1332    11775857 :   gel(y,2) = utoi(b);
    1333    11776091 :   gel(y,1) = utoipos(p); return y;
    1334             : }
    1335             : 
    1336             : /* special case: p = 1 return unit form */
    1337             : GEN
    1338       70042 : primeform(GEN x, GEN p)
    1339             : {
    1340       70042 :   const char *f = "primeform";
    1341             :   pari_sp av;
    1342       70042 :   long s, sx = signe(x), sp = signe(p);
    1343             :   GEN y, b, absp;
    1344             : 
    1345       70042 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
    1346       70042 :   if (typ(p) != t_INT) pari_err_TYPE(f,p);
    1347       70042 :   if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
    1348       70042 :   if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
    1349       70042 :   if (lgefint(p) == 3)
    1350             :   {
    1351       70028 :     ulong pp = p[2];
    1352       70028 :     if (pp == 1) {
    1353       13433 :       if (sx < 0) {
    1354             :         long r;
    1355        6601 :         if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1356        6601 :         r = mod4(x);
    1357        6601 :         if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
    1358        6601 :         return qfi_1_by_disc(x);
    1359             :       }
    1360        6832 :       y = qfr_1_by_disc(x);
    1361        6832 :       if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
    1362        6832 :       return y;
    1363             :     }
    1364       56595 :     y = primeform_u(x, pp);
    1365       56588 :     if (sx < 0) {
    1366       28889 :       if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1367       28889 :       return y;
    1368             :     }
    1369       27699 :     if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
    1370       27699 :     return gcopy( qfr3_to_qfr(y, x) );
    1371             :   }
    1372          14 :   s = mod8(x);
    1373          14 :   if (sx < 0)
    1374             :   {
    1375           7 :     if (sp < 0) pari_err_IMPL("negative definite t_QFB");
    1376           7 :     if (s) s = 8-s;
    1377             :   }
    1378          14 :   y = cgetg(5, t_QFB);
    1379             :   /* 2 or 3 mod 4 */
    1380          14 :   if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
    1381          14 :   absp = absi_shallow(p); av = avma;
    1382          14 :   b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
    1383          14 :   s &= 1; /* s = x mod 2 */
    1384             :   /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
    1385          14 :   if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(absp,b));
    1386             : 
    1387          14 :   av = avma;
    1388          14 :   gel(y,3) = gerepileuptoint(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
    1389          14 :   gel(y,4) = icopy(x);
    1390          14 :   gel(y,2) = b;
    1391          14 :   gel(y,1) = icopy(p);
    1392          14 :   return y;
    1393             : }
    1394             : 
    1395             : static GEN
    1396     2843836 : normforms(GEN D, GEN fa)
    1397             : {
    1398             :   long i, j, k, lB, aN, sa;
    1399             :   GEN a, L, V, B, N, N2;
    1400     2843836 :   int D_odd = mpodd(D);
    1401     2843836 :   a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
    1402     2843835 :   sa = signe(a);
    1403     2843835 :   if (sa==0 || (signe(D)<0 && sa<0)) return NULL;
    1404     1299046 :   V = D_odd? Zn_quad_roots(fa, gen_1, shifti(subsi(1, D), -2))
    1405     2774829 :            : Zn_quad_roots(fa, gen_0, negi(shifti(D, -2)));
    1406     2774824 :   if (!V) return NULL;
    1407      638935 :   N = gel(V,1); B = gel(V,2); lB = lg(B);
    1408      638935 :   N2 = shifti(N,1);
    1409      638911 :   aN = itou(diviiexact(a, N)); /* |a|/N */
    1410      638905 :   L = cgetg((lB-1)*aN+1, t_VEC);
    1411     2738595 :   for (k = 1, i = 1; i < lB; i++)
    1412             :   {
    1413     2099675 :     GEN b = shifti(gel(B,i), 1), c, C;
    1414     2099566 :     if (D_odd) b = addiu(b , 1);
    1415     2099566 :     c = diviiexact(shifti(subii(sqri(b), D), -2), a);
    1416     2099504 :     for (j = 0;; b = addii(b, N2))
    1417             :     {
    1418     2467585 :       gel(L, k++) = mkqfb(a, b, c, D);
    1419     2467718 :       if (++j == aN) break;
    1420      368054 :       C = addii(b, N); if (aN > 1) C = diviuexact(C, aN);
    1421      368081 :       c = sa > 0? addii(c, C): subii(c, C);
    1422             :     }
    1423             :   }
    1424      638920 :   return L;
    1425             : }
    1426             : 
    1427             : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
    1428             : static GEN
    1429      412752 : SL2_div_mul_e1(GEN N, GEN M)
    1430             : {
    1431      412752 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1432      412752 :   GEN p = subii(mulii(gcoeff(N,1,1), d), mulii(gcoeff(N,1,2), b));
    1433      412576 :   GEN q = subii(mulii(gcoeff(N,2,1), d), mulii(gcoeff(N,2,2), b));
    1434      412625 :   return mkvec2(p,q);
    1435             : }
    1436             : 
    1437             : /* Let M and N in SL2(Z), return (N*[1,0;0,-1]*M^-1)[,1] */
    1438             : static GEN
    1439           0 : SL2_swap_div_mul_e1(GEN N, GEN M)
    1440             : {
    1441           0 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1442           0 :   GEN p = addii(mulii(gcoeff(N,1,1), d), mulii(gcoeff(N,1,2), b));
    1443           0 :   GEN q = addii(mulii(gcoeff(N,2,1), d), mulii(gcoeff(N,2,2), b));
    1444           0 :   return mkvec2(p,q);
    1445             : }
    1446             : 
    1447             : /* Test equality modulo GL2 of two reduced forms */
    1448             : static int
    1449     1690154 : GL2_qfb_equal(GEN a, GEN b)
    1450             : {
    1451     1690154 :   return equalii(gel(a,1),gel(b,1))
    1452      254999 :    && absequalii(gel(a,2),gel(b,2))
    1453     1945117 :    &&    equalii(gel(a,3),gel(b,3));
    1454             : }
    1455             : 
    1456             : static GEN
    1457     1690122 : qfisolve_normform(GEN Q, GEN P)
    1458             : {
    1459     1690122 :   GEN a = gel(Q,1), N = gel(Q,2);
    1460     1690122 :   GEN M, b = redimagsl2(P, &M);
    1461     1690168 :   if (!GL2_qfb_equal(a,b) || signe(gel(a,2))!=signe(gel(b,2)))
    1462     1509088 :     return NULL;
    1463      181017 :   return SL2_div_mul_e1(N,M);
    1464             : }
    1465             : 
    1466             : static GEN
    1467           0 : qfbsolve_cornacchia(GEN c, GEN p, int swap)
    1468             : {
    1469           0 :   pari_sp av = avma;
    1470             :   GEN M, N;
    1471           0 :   if (kronecker(negi(c), p) < 0 || !cornacchia(c, p, &M,&N))
    1472           0 :   { set_avma(av); return gen_0; }
    1473           0 :   return gerepilecopy(av, swap? mkvec2(N,M): mkvec2(M,N));
    1474             : }
    1475             : 
    1476             : GEN
    1477           0 : qfisolvep(GEN Q, GEN p)
    1478             : {
    1479             :   GEN M, N, x,y, a,b,c, d;
    1480           0 :   pari_sp av = avma;
    1481           0 :   if (!signe(gel(Q,2)))
    1482             :   {
    1483           0 :     a = gel(Q,1);
    1484           0 :     c = gel(Q,3); /* if principal form, use faster cornacchia */
    1485           0 :     if (equali1(a)) return qfbsolve_cornacchia(c, p, 0);
    1486           0 :     if (equali1(c)) return qfbsolve_cornacchia(a, p, 1);
    1487             :   }
    1488           0 :   d = qfb_disc(Q); if (kronecker(d,p) < 0) return gen_0;
    1489           0 :   a = redimagsl2(Q, &N);
    1490           0 :   if (equali1(gel(a,1))) /* principal form */
    1491             :   {
    1492             :     long r;
    1493           0 :     if (!signe(gel(a,2)))
    1494             :     {
    1495           0 :       a = qfbsolve_cornacchia(gel(a,3), p, 0);
    1496           0 :       if (a == gen_0) { set_avma(av); return gen_0; }
    1497           0 :       a = ZM_ZC_mul(N, a);
    1498           0 :       a[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1499           0 :       return gerepileupto(av, a);
    1500             :     }
    1501             :     /* x^2 + xy + ((1-d)/4)y^2 = p <==> (2x + y)^2 - d y^2 = 4p */
    1502           0 :     if (!cornacchia2(negi(d), p, &x, &y)) { set_avma(av); return gen_0; }
    1503           0 :     x = divis_rem(subii(x,y), 2, &r); if (r) { set_avma(av); return gen_0; }
    1504           0 :     a = ZM_ZC_mul(N, mkvec2(x,y));
    1505           0 :     a[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1506           0 :     return gerepileupto(av, a);
    1507             :   }
    1508           0 :   b = redimagsl2(primeform(d, p), &M);
    1509           0 :   if (!GL2_qfb_equal(a,b)) { set_avma(av); return gen_0; }
    1510           0 :   if (signe(gel(a,2))==signe(gel(b,2)))
    1511           0 :     x = SL2_div_mul_e1(N,M);
    1512             :   else
    1513           0 :     x = SL2_swap_div_mul_e1(N,M);
    1514           0 :   return gerepilecopy(av, x);
    1515             : }
    1516             : 
    1517             : static GEN
    1518    10810674 : redrealsl2step(GEN A, GEN rd)
    1519             : {
    1520    10810674 :   GEN N, V = gel(A,1), M = gel(A,2);
    1521    10810674 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
    1522    10810674 :   GEN C = absi_shallow(c);
    1523    10810674 :   GEN t = addii(b, gmax_shallow(rd, C));
    1524    10810674 :   GEN r, q = truedvmdii(t, shifti(C,1), &r);
    1525    10810674 :   b = subii(t, addii(r,b));
    1526    10810674 :   a = c;
    1527    10810674 :   c = truedivii(subii(sqri(b), d), shifti(c,2));
    1528    10810674 :   if (signe(a) < 0) togglesign(q);
    1529    21621348 :   N = mkmat2(gel(M,2),
    1530    10810674 :              mkcol2(subii(mulii(q, gcoeff(M, 1, 2)), gcoeff(M, 1, 1)),
    1531    10810674 :                     subii(mulii(q, gcoeff(M, 2, 2)), gcoeff(M, 2, 1))));
    1532    10810674 :   return mkvec2(mkqfb(a,b,c,gel(V,4)), N);
    1533             : }
    1534             : 
    1535             : static GEN
    1536      979622 : redrealsl2(GEN V, GEN rd)
    1537             : {
    1538      979622 :   pari_sp ltop = avma;
    1539             :   GEN M, u1, u2, v1, v2;
    1540      979622 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3), d = qfb_disc(V);
    1541      979622 :   u1 = v2 = gen_1; v1 = u2 = gen_0;
    1542     3548167 :   while (!ab_isreduced(a,b,rd))
    1543             :   {
    1544     2568545 :     GEN C = mpabs_shallow(c);
    1545     2568545 :     GEN t = addii(b, gmax_shallow(rd,C));
    1546     2568545 :     GEN r, q = truedvmdii(t, shifti(C,1), &r);
    1547     2568545 :     b = subii(t, addii(r,b));
    1548     2568545 :     a = c;
    1549     2568545 :     c = truedivii(subii(sqri(b), d), shifti(c,2));
    1550     2568545 :     if (signe(a) < 0) togglesign(q);
    1551     2568545 :     r = u1; u1 = v1; v1 = subii(mulii(q, v1), r);
    1552     2568545 :     r = u2; u2 = v2; v2 = subii(mulii(q, v2), r);
    1553     2568545 :     if (gc_needed(ltop, 1))
    1554             :     {
    1555           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"redrealsl2");
    1556           0 :       gerepileall(ltop, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
    1557             :     }
    1558             :   }
    1559      979622 :   M = mkmat2(mkcol2(u1,u2), mkcol2(v1,v2));
    1560      979622 :   return gerepilecopy(ltop, mkvec2(lg(V)==5?mkqfb(a,b,c,gel(V,4))
    1561           0 :                                            :mkvec3(a,b,c), M));
    1562             : }
    1563             : 
    1564             : GEN
    1565      646376 : qfbredsl2(GEN q, GEN isD)
    1566             : {
    1567             :   pari_sp av;
    1568      646376 :   if (typ(q) != t_QFB) pari_err_TYPE("qfbredsl2",q);
    1569      646376 :   if (qfb_is_qfi(q))
    1570             :   {
    1571      436873 :     GEN v = cgetg(3,t_VEC);
    1572      436869 :     if (isD) pari_err_TYPE("qfbredsl2", isD);
    1573      436869 :     gel(v,1) = redimagsl2(q, &gel(v,2)); return v;
    1574             :   }
    1575      209503 :   av = avma;
    1576      209503 :   if (!isD) isD = sqrti(qfb_disc(q));
    1577      208068 :   else if (typ(isD) != t_INT) pari_err_TYPE("qfbredsl2",isD);
    1578      209496 :   return gerepileupto(av, redrealsl2(q, isD));
    1579             : }
    1580             : 
    1581             : static GEN
    1582      770126 : qfrsolve_normform(GEN N, GEN Ps, GEN rd)
    1583             : {
    1584      770126 :   pari_sp av = avma, btop;
    1585      770126 :   GEN M = N, P = redrealsl2(Ps, rd), Q = P;
    1586             : 
    1587      770126 :   btop = avma;
    1588             :   for(;;)
    1589             :   {
    1590     5840681 :     if (qfb_equal(gel(M,1), gel(P,1)))
    1591      154084 :       return gerepilecopy(av, SL2_div_mul_e1(gel(M,2),gel(P,2)));
    1592     5686597 :     if (qfb_equal(gel(N,1), gel(Q,1)))
    1593       77658 :       return gerepilecopy(av, SL2_div_mul_e1(gel(N,2),gel(Q,2)));
    1594     5608939 :     M = redrealsl2step(M, rd);
    1595     5608939 :     if (qfb_equal(gel(M,1), gel(N,1))) return gc_NULL(av);
    1596     5201735 :     Q = redrealsl2step(Q, rd);
    1597     5201735 :     if (qfb_equal(gel(P,1), gel(Q,1))) return gc_NULL(av);
    1598     5070555 :     if (gc_needed(btop, 1)) gerepileall(btop, 2, &M, &Q);
    1599             :   }
    1600             : }
    1601             : 
    1602             : GEN
    1603           0 : qfrsolvep(GEN Q, GEN p)
    1604             : {
    1605           0 :   pari_sp av = avma;
    1606           0 :   GEN N, x, rd, d = qfb_disc(Q);
    1607             : 
    1608           0 :   if (kronecker(d, p) < 0) return gc_const(av, gen_0);
    1609           0 :   rd = sqrti(d);
    1610           0 :   N = redrealsl2(Q, rd);
    1611           0 :   x = qfrsolve_normform(N, primeform(d, p), rd);
    1612           0 :   return x? gerepileupto(av, x): gc_const(av, gen_0);
    1613             : }
    1614             : 
    1615             : static GEN
    1616     2460247 : qfsolve_normform(GEN Q, GEN f, GEN rd)
    1617     2460247 : { return rd? qfrsolve_normform(Q, f, rd): qfisolve_normform(Q, f); }
    1618             : static GEN
    1619     2843839 : qfbsolve_primitive_i(GEN Q, GEN rd, GEN *Qr, GEN fa, long all)
    1620             : {
    1621     2843839 :   GEN x, W, F = normforms(qfb_disc(Q), fa);
    1622             :   long i, j, l;
    1623     2843815 :   if (!F) return NULL;
    1624      638920 :   if (!*Qr) *Qr = qfbredsl2(Q, rd);
    1625      638942 :   l = lg(F); W = all? cgetg(l, t_VEC): NULL;
    1626     3091929 :   for (j = i = 1; i < l; i++)
    1627     2460247 :     if ((x = qfsolve_normform(*Qr, gel(F,i), rd)))
    1628             :     {
    1629      412677 :       if (!all) return x;
    1630      405530 :       gel(W,j++) = x;
    1631             :     }
    1632      631682 :   if (j == 1) return NULL;
    1633      168744 :   setlg(W,j); return W;
    1634             : }
    1635             : 
    1636             : static GEN
    1637     2838553 : qfb_initrd(GEN Q) { GEN d = qfb_disc(Q); return signe(d) > 0? sqrti(d): NULL; }
    1638             : static GEN
    1639     2828341 : qfbsolve_primitive(GEN Q, GEN fa, long all)
    1640             : {
    1641     2828341 :   GEN x, Qr = NULL, rdQ = qfb_initrd(Q);
    1642     2828340 :   x = qfbsolve_primitive_i(Q, rdQ, &Qr, fa, all);
    1643     2828304 :   if (!x) return cgetg(1, t_VEC);
    1644      174765 :   return x;
    1645             : }
    1646             : 
    1647             : /* f / g^2 */
    1648             : static GEN
    1649        5285 : famat_divsqr(GEN f, GEN g)
    1650        5285 : { return famat_reduce(famat_div_shallow(f, famat_pows_shallow(g,2))); }
    1651             : static GEN
    1652       10213 : qfbsolve_all(GEN Q, GEN n, long all)
    1653             : {
    1654       10213 :   GEN W, Qr = NULL, fa = factorint(n, 0), rdQ = qfb_initrd(Q);
    1655       10213 :   GEN D = divisors_factored(mkmat2(gel(fa,1), gshift(gel(fa,2),-1)));
    1656       10213 :   long i, j, l = lg(D);
    1657       10213 :   W = all? cgetg(l, t_VEC): NULL;
    1658       25109 :   for (i = j = 1; i < l; i++)
    1659             :   {
    1660       15498 :     GEN w, d = gel(D,i), FA = i == 1? fa: famat_divsqr(fa, gel(d,2));
    1661       15498 :     if ((w = qfbsolve_primitive_i(Q, rdQ, &Qr, FA, all)))
    1662             :     {
    1663        1204 :       if (i != 1) w = RgV_Rg_mul(w, gel(d,1));
    1664        1204 :       if (!all) return w;
    1665         602 :       gel(W,j++) = w;
    1666             :     }
    1667             :   }
    1668        9611 :   if (j == 1) return cgetg(1, t_VEC);
    1669         511 :   setlg(W,j); return shallowconcat1(W);
    1670             : }
    1671             : 
    1672             : GEN
    1673     2838562 : qfbsolve(GEN Q, GEN n, long fl)
    1674             : {
    1675     2838562 :   pari_sp av = avma;
    1676     2838562 :   if (typ(Q) != t_QFB) pari_err_TYPE("qfbsolve",Q);
    1677     2838562 :   if (fl < 0 || fl > 3) pari_err_FLAG("qfbsolve");
    1678     5666858 :   return gerepilecopy(av, (fl & 2)? qfbsolve_all(Q, n, fl & 1)
    1679     2828341 :                                   : qfbsolve_primitive(Q, n, fl & 1));
    1680             : }
    1681             : 
    1682             : /* 1 if there exists x,y such that x^2 + dy^2 = p, 0 otherwise;
    1683             :  * Assume d > 0 and p is prime */
    1684             : long
    1685        7483 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
    1686             : {
    1687        7483 :   pari_sp av = avma;
    1688             :   GEN b, c, r;
    1689             : 
    1690        7483 :   *px = *py = gen_0;
    1691        7483 :   b = subii(p, d);
    1692        7483 :   if (signe(b) < 0) return gc_long(av,0);
    1693        7483 :   if (signe(b) == 0) { *py = gen_1; return gc_long(av,1); }
    1694        7476 :   b = Fp_sqrt(b, p); /* sqrt(-d) */
    1695        7476 :   if (!b) return gc_long(av,0);
    1696        3745 :   b = gmael(halfgcdii(p, b), 2, 2);
    1697        3745 :   c = dvmdii(subii(p, sqri(b)), d, &r);
    1698        3745 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    1699         231 :   set_avma(av);
    1700         231 :   *px = icopy(b);
    1701         231 :   *py = icopy(c); return 1;
    1702             : }
    1703             : 
    1704             : static GEN
    1705     1111217 : lastqi(GEN Q)
    1706             : {
    1707     1111217 :   GEN s = gcoeff(Q,1,1), q = gcoeff(Q,1,2), p = absi_shallow(gcoeff(Q,2,2));
    1708     1111216 :   if (!signe(q)) return gen_0;
    1709     1111209 :   if (!signe(s)) return p;
    1710     1106330 :   if (is_pm1(q)) return subiu(p,1);
    1711     1106330 :   return divii(p, absi_shallow(q));
    1712             : }
    1713             : 
    1714             : static long
    1715     1111236 : cornacchia2_i(long av, GEN d, GEN p, GEN b, GEN px4, GEN *px, GEN *py)
    1716             : {
    1717             :   GEN M, Q, V, c, r, b2;
    1718     1111236 :   if (!signe(b)) { /* d = p,2p,3p,4p */
    1719           7 :     set_avma(av);
    1720           7 :     if (absequalii(d, px4)){ *py = gen_1; return 1; }
    1721           7 :     if (absequalii(d, p))  { *py = gen_2; return 1; }
    1722           0 :     return 0;
    1723             :   }
    1724     1111229 :   if (mod2(b) != mod2(d)) b = subii(p,b);
    1725     1111187 :   M = halfgcdii(shifti(p,1), b); Q = gel(M,1); V = gel(M, 2);
    1726     1111223 :   b = addii(mulii(gel(V,1), lastqi(Q)), gel(V,2));
    1727     1111074 :   b2 = sqri(b);
    1728     1111079 :   if (cmpii(b2,px4) > 0)
    1729             :   {
    1730     1104374 :     b = gel(V,1); b2 = sqri(b);
    1731     1104335 :     if (cmpii(b2,px4) > 0) { b = gel(V,2); b2 = sqri(b); }
    1732             :   }
    1733     1111120 :   c = dvmdii(subii(px4, b2), d, &r);
    1734     1111111 :   if (r != gen_0 || !Z_issquareall(c, &c)) return gc_long(av,0);
    1735     1081392 :   set_avma(av);
    1736     1081388 :   *px = icopy(b);
    1737     1081392 :   *py = icopy(c); return 1;
    1738             : }
    1739             : 
    1740             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p, 0 otherwise;
    1741             :  * Assume d > 0 is congruent to 0 or 3 mod 4 and p is prime */
    1742             : long
    1743     1070597 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
    1744             : {
    1745     1070597 :   pari_sp av = avma;
    1746     1070597 :   GEN b, p4 = shifti(p,2);
    1747             : 
    1748     1070510 :   *px = *py = gen_0;
    1749     1070510 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    1750     1070595 :   if (absequaliu(p, 2))
    1751             :   {
    1752           7 :     set_avma(av);
    1753           7 :     switch (itou_or_0(d)) {
    1754           0 :       case 4: *px = gen_2; break;
    1755           0 :       case 7: *px = gen_1; break;
    1756           7 :       default: return 0;
    1757             :     }
    1758           0 :     *py = gen_1; return 1;
    1759             :   }
    1760     1070591 :   b = Fp_sqrt(negi(d), p);
    1761     1070642 :   if (!b) return gc_long(av,0);
    1762     1070558 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    1763             : }
    1764             : 
    1765             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
    1766             : long
    1767       40677 : cornacchia2_sqrt(GEN d, GEN p, GEN b, GEN *px, GEN *py)
    1768             : {
    1769       40677 :   pari_sp av = avma;
    1770       40677 :   GEN p4 = shifti(p,2);
    1771       40677 :   *px = *py = gen_0;
    1772       40677 :   if (abscmpii(p4, d) < 0) return gc_long(av,0);
    1773       40677 :   return cornacchia2_i(av, d, p, b, p4, px, py);
    1774             : }
    1775             : 
    1776             : GEN
    1777        7630 : qfbcornacchia(GEN d, GEN p)
    1778             : {
    1779        7630 :   pari_sp av = avma;
    1780             :   GEN x, y;
    1781        7630 :   if (typ(d) != t_INT || signe(d) <= 0) pari_err_TYPE("qfbcornacchia", d);
    1782        7630 :   if (typ(p) != t_INT || cmpiu(p, 2) < 0) pari_err_TYPE("qfbcornacchia", p);
    1783        7630 :   if (mod4(p)? cornacchia(d, p, &x, &y): cornacchia2(d, shifti(p, -2), &x, &y))
    1784         287 :     return gerepilecopy(av, mkvec2(x, y));
    1785        7343 :   set_avma(av); return cgetg(1, t_VEC);
    1786             : }

Generated by: LCOV version 1.13