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 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.8.0 lcov report (development 19357-d770f77) Lines: 833 905 92.0 %
Date: 2016-08-27 06:11:27 Functions: 91 97 93.8 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : #include "pari.h"
      14             : #include "paripriv.h"
      15             : /*******************************************************************/
      16             : /*                                                                 */
      17             : /*         QUADRATIC POLYNOMIAL ASSOCIATED TO A DISCRIMINANT       */
      18             : /*                                                                 */
      19             : /*******************************************************************/
      20             : 
      21             : void
      22     1827292 : check_quaddisc(GEN x, long *s, long *r, const char *f)
      23             : {
      24     1827292 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
      25     1827292 :   *s = signe(x);
      26     1827292 :   if (Z_issquare(x)) pari_err_DOMAIN(f,"issquare(disc)","=", gen_1,x);
      27     1827285 :   *r = mod4(x); if (*s < 0 && *r) *r = 4 - *r;
      28     1827285 :   if (*r > 1) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
      29     1827271 : }
      30             : void
      31        6608 : check_quaddisc_real(GEN x, long *r, const char *f)
      32             : {
      33        6608 :   long sx; check_quaddisc(x, &sx, r, f);
      34        6608 :   if (sx < 0) pari_err_DOMAIN(f, "disc","<",gen_0,x);
      35        6608 : }
      36             : void
      37         609 : check_quaddisc_imag(GEN x, long *r, const char *f)
      38             : {
      39         609 :   long sx; check_quaddisc(x, &sx, r, f);
      40         602 :   if (sx > 0) pari_err_DOMAIN(f, "disc",">",gen_0,x);
      41         602 : }
      42             : 
      43             : /* X^2 + b X + c is the canonical quadratic t_POL of discriminant D.
      44             :  * Dodd is non-zero iff D is odd */
      45             : static void
      46     2211811 : quadpoly_bc(GEN D, long Dodd, GEN *b, GEN *c)
      47             : {
      48     2211811 :   if (Dodd)
      49             :   {
      50     2207968 :     pari_sp av = avma;
      51     2207968 :     *b = gen_m1;
      52     2207968 :     *c = gerepileuptoint(av, shifti(subui(1,D), -2));
      53             :   }
      54             :   else
      55             :   {
      56        3843 :     *b = gen_0;
      57        3843 :     *c = shifti(D,-2); togglesign(*c);
      58             :   }
      59     2211811 : }
      60             : /* X^2 - X - (D-1)/4 or X^2 - D/4 */
      61             : GEN
      62        1323 : quadpoly(GEN D)
      63             : {
      64             :   long Dmod4, s;
      65        1323 :   GEN b, c, y = cgetg(5,t_POL);
      66        1323 :   check_quaddisc(D, &s, &Dmod4, "quadpoly");
      67        1316 :   y[1] = evalsigne(1) | evalvarn(0);
      68        1316 :   quadpoly_bc(D, Dmod4, &b,&c);
      69        1316 :   gel(y,2) = c;
      70        1316 :   gel(y,3) = b;
      71        1316 :   gel(y,4) = gen_1; return y;
      72             : }
      73             : 
      74             : GEN
      75         441 : quadpoly0(GEN x, long v)
      76             : {
      77         441 :   GEN T = quadpoly(x);
      78         434 :   if (v > 0) setvarn(T, v);
      79         434 :   return T;
      80             : }
      81             : 
      82             : GEN
      83         259 : quadgen(GEN x)
      84         259 : { retmkquad(quadpoly(x), gen_0, gen_1); }
      85             : 
      86             : /***********************************************************************/
      87             : /**                                                                   **/
      88             : /**                      BINARY QUADRATIC FORMS                       **/
      89             : /**                                                                   **/
      90             : /***********************************************************************/
      91             : GEN
      92      105476 : qfi(GEN x, GEN y, GEN z)
      93             : {
      94      105476 :   if (signe(x) < 0) pari_err_IMPL("negative definite t_QFI");
      95      105476 :   retmkqfi(icopy(x),icopy(y),icopy(z));
      96             : }
      97             : GEN
      98       34244 : qfr(GEN x, GEN y, GEN z, GEN d)
      99             : {
     100       34244 :   if (typ(d) != t_REAL) pari_err_TYPE("qfr",d);
     101       34244 :   retmkqfr(icopy(x),icopy(y),icopy(z),rcopy(d));
     102             : }
     103             : 
     104             : GEN
     105       91840 : Qfb0(GEN x, GEN y, GEN z, GEN d, long prec)
     106             : {
     107       91840 :   pari_sp av = avma;
     108             :   GEN D;
     109             :   long s, r;
     110       91840 :   if (typ(x)!=t_INT) pari_err_TYPE("Qfb",x);
     111       91840 :   if (typ(y)!=t_INT) pari_err_TYPE("Qfb",y);
     112       91840 :   if (typ(z)!=t_INT) pari_err_TYPE("Qfb",z);
     113       91840 :   D = qfb_disc3(x,y,z);
     114       91840 :   check_quaddisc(D, &s, &r, "Qfb");
     115       91833 :   avma = av;
     116       91833 :   if (s < 0) return qfi(x, y, z);
     117             : 
     118       34244 :   d = d? gtofp(d,prec): real_0(prec);
     119       34244 :   return qfr(x,y,z,d);
     120             : }
     121             : 
     122             : /* composition */
     123             : static void
     124    65977490 : qfb_sqr(GEN z, GEN x)
     125             : {
     126             :   GEN c, d1, x2, v1, v2, c3, m, p1, r;
     127             : 
     128    65977490 :   d1 = bezout(gel(x,2),gel(x,1),&x2, NULL); /* usually 1 */
     129    65977490 :   c = gel(x,3);
     130    65977490 :   m = mulii(c,x2);
     131    65977490 :   if (is_pm1(d1))
     132    60418799 :     v1 = v2 = gel(x,1);
     133             :   else
     134             :   {
     135     5558691 :     v1 = diviiexact(gel(x,1),d1);
     136     5558691 :     v2 = mulii(v1, gcdii(d1,c)); /* = v1 iff x primitive */
     137     5558691 :     c = mulii(c, d1);
     138             :   }
     139    65977490 :   togglesign(m);
     140    65977490 :   r = modii(m,v2);
     141    65977490 :   p1 = mulii(r, v1);
     142    65977490 :   c3 = addii(c, mulii(r,addii(gel(x,2),p1)));
     143    65977490 :   gel(z,1) = mulii(v1,v2);
     144    65977490 :   gel(z,2) = addii(gel(x,2), shifti(p1,1));
     145    65977490 :   gel(z,3) = diviiexact(c3,v2);
     146    65977490 : }
     147             : /* z <- x * y */
     148             : static void
     149    45003763 : qfb_comp(GEN z, GEN x, GEN y)
     150             : {
     151             :   GEN n, c, d, y1, v1, v2, c3, m, p1, r;
     152             : 
     153    90007526 :   if (x == y) { qfb_sqr(z,x); return; }
     154    42678727 :   n = shifti(subii(gel(y,2),gel(x,2)), -1);
     155    42678727 :   v1 = gel(x,1);
     156    42678727 :   v2 = gel(y,1);
     157    42678727 :   c  = gel(y,3);
     158    42678727 :   d = bezout(v2,v1,&y1,NULL);
     159    42678727 :   if (is_pm1(d))
     160    19928543 :     m = mulii(y1,n);
     161             :   else
     162             :   {
     163    22750184 :     GEN s = subii(gel(y,2), n);
     164    22750184 :     GEN x2, y2, d1 = bezout(s,d,&x2,&y2); /* x2 s + y2 (x1 v1 + y1 v2) = d1 */
     165    22750184 :     if (!is_pm1(d1))
     166             :     {
     167    11687947 :       v1 = diviiexact(v1,d1);
     168    11687947 :       v2 = diviiexact(v2,d1); /* gcd = 1 iff x or y primitive */
     169    11687947 :       v1 = mulii(v1, gcdii(c,gcdii(gel(x,3),gcdii(d1,n))));
     170    11687947 :       c = mulii(c, d1);
     171             :     }
     172    22750184 :     m = addii(mulii(mulii(y1,y2),n), mulii(gel(y,3),x2));
     173             :   }
     174    42678727 :   togglesign(m);
     175    42678727 :   r = modii(m, v1);
     176    42678727 :   p1 = mulii(r, v2);
     177    42678727 :   c3 = addii(c, mulii(r,addii(gel(y,2),p1)));
     178    42678727 :   gel(z,1) = mulii(v1,v2);
     179    42678727 :   gel(z,2) = addii(gel(y,2), shifti(p1,1));
     180    42678727 :   gel(z,3) = diviiexact(c3,v1);
     181             : }
     182             : 
     183             : static GEN redimag_av(pari_sp av, GEN q);
     184             : static GEN
     185    42266416 : qficomp0(GEN x, GEN y, int raw)
     186             : {
     187    42266416 :   pari_sp av = avma;
     188    42266416 :   GEN z = cgetg(4,t_QFI);
     189    42266416 :   qfb_comp(z, x,y);
     190    42266416 :   if (raw) return gerepilecopy(av,z);
     191    42264687 :   return redimag_av(av, z);
     192             : }
     193             : static GEN
     194         427 : qfrcomp0(GEN x, GEN y, int raw)
     195             : {
     196         427 :   pari_sp av = avma;
     197         427 :   GEN z = cgetg(5,t_QFR);
     198         427 :   qfb_comp(z, x,y); gel(z,4) = addrr(gel(x,4),gel(y,4));
     199         427 :   if (raw) return gerepilecopy(av,z);
     200           7 :   return gerepileupto(av, redreal(z));
     201             : }
     202             : GEN
     203           7 : qfrcomp(GEN x, GEN y) { return qfrcomp0(x,y,0); }
     204             : GEN
     205         420 : qfrcompraw(GEN x, GEN y) { return qfrcomp0(x,y,1); }
     206             : GEN
     207    42264687 : qficomp(GEN x, GEN y) { return qficomp0(x,y,0); }
     208             : GEN
     209        1729 : qficompraw(GEN x, GEN y) { return qficomp0(x,y,1); }
     210             : GEN
     211        2135 : qfbcompraw(GEN x, GEN y)
     212             : {
     213        2135 :   long tx = typ(x);
     214        2135 :   if (typ(y) != tx) pari_err_TYPE2("*",x,y);
     215        2135 :   switch(tx) {
     216        1722 :     case t_QFI: return qficompraw(x,y);
     217         413 :     case t_QFR: return qfrcompraw(x,y);
     218             :   }
     219           0 :   pari_err_TYPE("composition",x);
     220           0 :   return NULL; /* not reached */
     221             : }
     222             : 
     223             : static GEN
     224    63652433 : qfisqr0(GEN x, long raw)
     225             : {
     226    63652433 :   pari_sp av = avma;
     227    63652433 :   GEN z = cgetg(4,t_QFI);
     228             : 
     229    63652433 :   if (typ(x)!=t_QFI) pari_err_TYPE("composition",x);
     230    63652433 :   qfb_sqr(z,x);
     231    63652433 :   if (raw) return gerepilecopy(av,z);
     232    63652426 :   return redimag_av(av, z);
     233             : }
     234             : static GEN
     235          21 : qfrsqr0(GEN x, long raw)
     236             : {
     237          21 :   pari_sp av = avma;
     238          21 :   GEN z = cgetg(5,t_QFR);
     239             : 
     240          21 :   if (typ(x)!=t_QFR) pari_err_TYPE("composition",x);
     241          21 :   qfb_sqr(z,x); gel(z,4) = shiftr(gel(x,4),1);
     242          21 :   if (raw) return gerepilecopy(av,z);
     243          14 :   return gerepileupto(av, redreal(z));
     244             : }
     245             : GEN
     246          14 : qfrsqr(GEN x) { return qfrsqr0(x,0); }
     247             : GEN
     248           7 : qfrsqrraw(GEN x) { return qfrsqr0(x,1); }
     249             : GEN
     250    63652426 : qfisqr(GEN x) { return qfisqr0(x,0); }
     251             : GEN
     252           7 : qfisqrraw(GEN x) { return qfisqr0(x,1); }
     253             : 
     254             : static GEN
     255        6580 : qfr_1_by_disc(GEN D, long prec)
     256             : {
     257        6580 :   GEN y = cgetg(5,t_QFR), isqrtD;
     258        6580 :   pari_sp av = avma;
     259             :   long r;
     260             : 
     261        6580 :   check_quaddisc_real(D, &r, "qfr_1_by_disc");
     262        6580 :   gel(y,1) = gen_1; isqrtD = sqrti(D);
     263        6580 :   if ((r & 1) != mod2(isqrtD)) /* we know isqrtD > 0 */
     264        3129 :     isqrtD = gerepileuptoint(av, addsi(-1,isqrtD));
     265        6580 :   gel(y,2) = isqrtD; av = avma;
     266        6580 :   gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(isqrtD), D),-2));
     267        6580 :   gel(y,4) = real_0(prec); return y;
     268             : }
     269             : GEN
     270          14 : qfr_1(GEN x)
     271             : {
     272          14 :   if (typ(x) != t_QFR) pari_err_TYPE("qfr_1",x);
     273          14 :   return qfr_1_by_disc(qfb_disc(x), precision(gel(x,4)));
     274             : }
     275             : 
     276             : static void
     277           0 : qfr_1_fill(GEN y, struct qfr_data *S)
     278             : {
     279           0 :   pari_sp av = avma;
     280           0 :   GEN y2 = S->isqrtD;
     281           0 :   gel(y,1) = gen_1;
     282           0 :   if (mod2(S->D) != mod2(y2)) y2 = addsi(-1,y2);
     283           0 :   gel(y,2) = y2; av = avma;
     284           0 :   gel(y,3) = gerepileuptoint(av, shifti(subii(sqri(y2), S->D),-2));
     285           0 : }
     286             : static GEN
     287           0 : qfr5_1(struct qfr_data *S, long prec)
     288             : {
     289           0 :   GEN y = cgetg(6, t_VEC);
     290           0 :   qfr_1_fill(y, S);
     291           0 :   gel(y,4) = gen_0;
     292           0 :   gel(y,5) = real_1(prec); return y;
     293             : }
     294             : static GEN
     295           0 : qfr3_1(struct qfr_data *S)
     296             : {
     297           0 :   GEN y = cgetg(4, t_VEC);
     298           0 :   qfr_1_fill(y, S); return y;
     299             : }
     300             : 
     301             : /* Assume D < 0 is the discriminant of a t_QFI */
     302             : static GEN
     303     2210495 : qfi_1_by_disc(GEN D)
     304             : {
     305     2210495 :   GEN b,c, y = cgetg(4,t_QFI);
     306     2210495 :   quadpoly_bc(D, mod2(D), &b,&c);
     307     2210495 :   if (b == gen_m1) b = gen_1;
     308     2210495 :   gel(y,1) = gen_1;
     309     2210495 :   gel(y,2) = b;
     310     2210495 :   gel(y,3) = c; return y;
     311             : }
     312             : GEN
     313     2203852 : qfi_1(GEN x)
     314             : {
     315     2203852 :   if (typ(x) != t_QFI) pari_err_TYPE("qfi_1",x);
     316     2203852 :   return qfi_1_by_disc(qfb_disc(x));
     317             : }
     318             : 
     319             : static GEN
     320           7 : invraw(GEN x)
     321             : {
     322           7 :   GEN y = gcopy(x);
     323           7 :   if (typ(y) == t_QFR) togglesign(gel(y,4));
     324           7 :   togglesign(gel(y,2)); return y;
     325             : }
     326             : GEN
     327          14 : qfrpowraw(GEN x, long n)
     328             : {
     329          14 :   pari_sp av = avma;
     330             :   long m;
     331             :   GEN y;
     332             : 
     333          14 :   if (typ(x) != t_QFR) pari_err_TYPE("qfrpowraw",x);
     334          14 :   if (!n) return qfr_1(x);
     335          14 :   if (n== 1) return gcopy(x);
     336          14 :   if (n==-1) return invraw(x);
     337             : 
     338           7 :   y = NULL; m = labs(n);
     339          14 :   for (; m>1; m>>=1)
     340             :   {
     341           7 :     if (m&1) y = y? qfrcompraw(y,x): x;
     342           7 :     x = qfrsqrraw(x);
     343             :   }
     344           7 :   y = y? qfrcompraw(y,x): x;
     345           7 :   if (n < 0) y = invraw(y);
     346           7 :   return gerepileupto(av,y);
     347             : }
     348             : GEN
     349           7 : qfipowraw(GEN x, long n)
     350             : {
     351           7 :   pari_sp av = avma;
     352             :   long m;
     353             :   GEN y;
     354             : 
     355           7 :   if (typ(x) != t_QFI) pari_err_TYPE("qfipow",x);
     356           7 :   if (!n) return qfi_1(x);
     357           7 :   if (n== 1) return gcopy(x);
     358           7 :   if (n==-1) return invraw(x);
     359             : 
     360           7 :   y = NULL; m = labs(n);
     361          14 :   for (; m>1; m>>=1)
     362             :   {
     363           7 :     if (m&1) y = y? qficompraw(y,x): x;
     364           7 :     x = qfisqrraw(x);
     365             :   }
     366           7 :   y = y? qficompraw(y,x): x;
     367           7 :   if (n < 0) y = invraw(y);
     368           7 :   return gerepileupto(av,y);
     369             : }
     370             : 
     371             : GEN
     372          21 : qfbpowraw(GEN x, long n)
     373          21 : { return (typ(x)==t_QFI)? qfipowraw(x,n): qfrpowraw(x,n); }
     374             : 
     375             : static long
     376      322777 : parteucl(GEN L, GEN *d, GEN *v3, GEN *v, GEN *v2)
     377             : {
     378             :   long z;
     379      322777 :   *v = gen_0; *v2 = gen_1;
     380     3534545 :   for (z=0; abscmpii(*v3,L) > 0; z++)
     381             :   {
     382     3211768 :     GEN t3, t2 = subii(*v, mulii(truedvmdii(*d,*v3,&t3),*v2));
     383     3211768 :     *v = *v2; *d = *v3; *v2 = t2; *v3 = t3;
     384             :   }
     385      322777 :   return z;
     386             : }
     387             : 
     388             : /* composition: Shanks' NUCOMP & NUDUPL */
     389             : /* L = floor((|d|/4)^(1/4)) */
     390             : GEN
     391      319956 : nucomp(GEN x, GEN y, GEN L)
     392             : {
     393      319956 :   pari_sp av = avma;
     394             :   long z;
     395             :   GEN a, a1, a2, b2, b, d, d1, g, n, p1, q1, q2, s, u, u1, v, v2, v3, Q;
     396             : 
     397      319956 :   if (x==y) return nudupl(x,L);
     398      319935 :   if (typ(x) != t_QFI) pari_err_TYPE("nucomp",x);
     399      319935 :   if (typ(y) != t_QFI) pari_err_TYPE("nucomp",y);
     400             : 
     401      319935 :   if (abscmpii(gel(x,1),gel(y,1)) < 0) swap(x, y);
     402      319935 :   s = shifti(addii(gel(x,2),gel(y,2)), -1);
     403      319935 :   n = subii(gel(y,2), s);
     404      319935 :   a1 = gel(x,1);
     405      319935 :   a2 = gel(y,1); d = bezout(a2,a1,&u,&v);
     406      319935 :   if (is_pm1(d)) { a = negi(mulii(u,n)); d1 = d; }
     407             :   else
     408      120141 :     if (remii(s,d) == gen_0) /* d | s */
     409             :     {
     410       60109 :       a = negi(mulii(u,n)); d1 = d;
     411       60109 :       a1 = diviiexact(a1, d1);
     412       60109 :       a2 = diviiexact(a2, d1);
     413       60109 :       s = diviiexact(s, d1);
     414             :     }
     415             :     else
     416             :     {
     417             :       GEN p2, l;
     418       60032 :       d1 = bezout(s,d,&u1,NULL);
     419       60032 :       if (!is_pm1(d1))
     420             :       {
     421          28 :         a1 = diviiexact(a1,d1);
     422          28 :         a2 = diviiexact(a2,d1);
     423          28 :         s = diviiexact(s,d1);
     424          28 :         d = diviiexact(d,d1);
     425             :       }
     426       60032 :       p1 = remii(gel(x,3),d);
     427       60032 :       p2 = remii(gel(y,3),d);
     428       60032 :       l = modii(mulii(negi(u1), addii(mulii(u,p1),mulii(v,p2))), d);
     429       60032 :       a = subii(mulii(l,diviiexact(a1,d)), mulii(u,diviiexact(n,d)));
     430             :     }
     431      319935 :   a = modii(a,a1); p1 = subii(a,a1); if (abscmpii(a,p1) > 0) a = p1;
     432      319935 :   d = a1; v3 = a; z = parteucl(L, &d,&v3, &v,&v2);
     433      319935 :   Q = cgetg(4,t_QFI);
     434      319935 :   if (!z) {
     435       37317 :     g = diviiexact(addii(mulii(v3,s),gel(y,3)), d);
     436       37317 :     b = a2;
     437       37317 :     b2 = gel(y,2);
     438       37317 :     v2 = d1;
     439       37317 :     gel(Q,1) = mulii(d,b);
     440             :   } else {
     441             :     GEN e, q3, q4;
     442      282618 :     if (z&1) { v3 = negi(v3); v2 = negi(v2); }
     443      282618 :     b = diviiexact(addii(mulii(a2,d), mulii(n,v)), a1);
     444      282618 :     e = diviiexact(addii(mulii(s,d),mulii(gel(y,3),v)), a1);
     445      282618 :     q3 = mulii(e,v2);
     446      282618 :     q4 = subii(q3,s);
     447      282618 :     b2 = addii(q3,q4);
     448      282618 :     g = diviiexact(q4,v);
     449      282618 :     if (!is_pm1(d1)) { v2 = mulii(d1,v2); v = mulii(d1,v); b2 = mulii(d1,b2); }
     450      282618 :     gel(Q,1) = addii(mulii(d,b), mulii(e,v));
     451             :   }
     452      319935 :   q1 = mulii(b, v3);
     453      319935 :   q2 = addii(q1,n);
     454      319935 :   gel(Q,2) = addii(b2, z? addii(q1,q2): shifti(q1, 1));
     455      319935 :   gel(Q,3) = addii(mulii(v3,diviiexact(q2,d)), mulii(g,v2));
     456      319935 :   return redimag_av(av, Q);
     457             : }
     458             : 
     459             : GEN
     460        2842 : nudupl(GEN x, GEN L)
     461             : {
     462        2842 :   pari_sp av = avma;
     463             :   long z;
     464             :   GEN u, v, d, d1, p1, a, b, c, a2, b2, c2, Q, v2, v3, g;
     465             : 
     466        2842 :   if (typ(x) != t_QFI) pari_err_TYPE("nudupl",x);
     467        2842 :   a = gel(x,1);
     468        2842 :   b = gel(x,2);
     469        2842 :   d1 = bezout(b,a, &u,NULL);
     470        2842 :   if (!is_pm1(d1))
     471             :   {
     472        1141 :     a = diviiexact(a, d1);
     473        1141 :     b = diviiexact(b, d1);
     474             :   }
     475        2842 :   c = modii(negi(mulii(u,gel(x,3))), a);
     476        2842 :   p1 = subii(c,a); if (abscmpii(c,p1) > 0) c = p1;
     477        2842 :   d = a; v3 = c; z = parteucl(L, &d,&v3, &v,&v2);
     478        2842 :   a2 = sqri(d);
     479        2842 :   c2 = sqri(v3);
     480        2842 :   Q = cgetg(4,t_QFI);
     481        2842 :   if (!z) {
     482         273 :     g = diviiexact(addii(mulii(v3,b),gel(x,3)), d);
     483         273 :     b2 = gel(x,2);
     484         273 :     v2 = d1;
     485         273 :     gel(Q,1) = a2;
     486             :   } else {
     487             :     GEN e;
     488        2569 :     if (z&1) { v = negi(v); d = negi(d); }
     489        2569 :     e = diviiexact(addii(mulii(gel(x,3),v), mulii(b,d)), a);
     490        2569 :     g = diviiexact(subii(mulii(e,v2), b), v);
     491        2569 :     b2 = addii(mulii(e,v2), mulii(v,g));
     492        2569 :     if (!is_pm1(d1)) { b2 = mulii(d1,b2); v = mulii(d1,v); v2 = mulii(d1,v2); }
     493        2569 :     gel(Q,1) = addii(a2, mulii(e,v));
     494             :   }
     495        2842 :   gel(Q,2) = addii(b2, subii(sqri(addii(d,v3)), addii(a2,c2)));
     496        2842 :   gel(Q,3) = addii(c2, mulii(g,v2));
     497        2842 :   return redimag_av(av, Q);
     498             : }
     499             : 
     500             : static GEN
     501        1029 : mul_nucomp(void *l, GEN x, GEN y) { return nucomp(x, y, (GEN)l); }
     502             : static GEN
     503        2821 : mul_nudupl(void *l, GEN x) { return nudupl(x, (GEN)l); }
     504             : GEN
     505         168 : nupow(GEN x, GEN n, GEN L)
     506             : {
     507             :   pari_sp av;
     508             :   GEN y, D;
     509             : 
     510         168 :   if (typ(n) != t_INT) pari_err_TYPE("nupow",n);
     511         168 :   if (typ(x) != t_QFI) pari_err_TYPE("nupow",x);
     512         168 :   if (gequal1(n)) return gcopy(x);
     513         168 :   av = avma;
     514         168 :   D = qfb_disc(x);
     515         168 :   y = qfi_1_by_disc(D);
     516         168 :   if (!signe(n)) return y;
     517         154 :   if (!L) L = sqrtnint(absi(D), 4);
     518         154 :   y = gen_pow(x, n, (void*)L, &mul_nudupl, &mul_nucomp);
     519         154 :   if (signe(n) < 0
     520          14 :   && !absequalii(gel(y,1),gel(y,2))
     521          14 :   && !absequalii(gel(y,1),gel(y,3))) togglesign(gel(y,2));
     522         154 :   return gerepileupto(av, y);
     523             : }
     524             : 
     525             : /* Reduction */
     526             : 
     527             : /* assume a > 0. Write b = q*2a + r, with -a < r <= a */
     528             : static GEN
     529     6118140 : dvmdii_round(GEN b, GEN a, GEN *r)
     530             : {
     531     6118140 :   GEN a2 = shifti(a, 1), q = dvmdii(b, a2, r);
     532     6118140 :   if (signe(b) >= 0) {
     533     2747939 :     if (abscmpii(*r, a) > 0) { q = addis(q,  1); *r = subii(*r, a2); }
     534             :   } else { /* r <= 0 */
     535     3370201 :     if (abscmpii(*r, a) >= 0){ q = addis(q, -1); *r = addii(*r, a2); }
     536             :   }
     537     6118140 :   return q;
     538             : }
     539             : /* Assume 0 < a <= LONG_MAX. Ensure no overflow */
     540             : static long
     541   181344339 : dvmdsu_round(long b, ulong a, long *r)
     542             : {
     543   181344339 :   ulong a2 = a << 1, q, ub, ur;
     544   181344339 :   if (b >= 0) {
     545   112525408 :     ub = b;
     546   112525408 :     q = ub / a2;
     547   112525408 :     ur = ub % a2;
     548   112525408 :     if (ur > a) { ur -= a; q++; *r = (long)ur; *r -= (long)a; }
     549    39763181 :     else *r = (long)ur;
     550   112525408 :     return (long)q;
     551             :   } else { /* r <= 0 */
     552    68818931 :     ub = (ulong)-b; /* |b| */
     553    68818931 :     q = ub / a2;
     554    68818931 :     ur = ub % a2;
     555    68818931 :     if (ur >= a) { ur -= a; q++; *r = (long)ur; *r = (long)a - *r; }
     556    37411314 :     else *r = -(long)ur;
     557    68818931 :     return -(long)q;
     558             :   }
     559             : }
     560             : /* reduce b mod 2*a. Update b,c */
     561             : static void
     562     2696491 : REDB(GEN a, GEN *b, GEN *c)
     563             : {
     564     2696491 :   GEN r, q = dvmdii_round(*b, a, &r);
     565     5392982 :   if (!signe(q)) return;
     566     2632075 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     567     2632075 :   *b = r;
     568             : }
     569             : /* Assume a > 0. Reduce b mod 2*a. Update b,c */
     570             : static void
     571   181344339 : sREDB(ulong a, long *b, ulong *c)
     572             : {
     573             :   long r, q;
     574             :   ulong uz;
     575   197772334 :   if (a > LONG_MAX) return; /* b already reduced */
     576   181344339 :   q = dvmdsu_round(*b, a, &r);
     577   181344339 :   if (q == 0) return;
     578             :   /* Final (a,r,c2) satisfies |r| <= |b| hence c2 <= c, c2 = c - q*z,
     579             :    * where z = (b+r) / 2, representable as long, has the same sign as q. */
     580   164916344 :   if (*b < 0)
     581             :   { /* uz = -z >= 0, q < 0 */
     582    60451316 :     if (r >= 0) /* different signs=>no overflow, exact division */
     583    31581152 :       uz = (ulong)-((*b + r)>>1);
     584             :     else
     585             :     {
     586    28870164 :       ulong ub = (ulong)-*b, ur = (ulong)-r;
     587    28870164 :       uz = (ub + ur) >> 1;
     588             :     }
     589    60451316 :     *c -= (-q) * uz; /* c -= qz */
     590             :   }
     591             :   else
     592             :   { /* uz = z >= 0, q > 0 */
     593   104465028 :     if (r <= 0)
     594    72820612 :       uz = (*b + r)>>1;
     595             :     else
     596             :     {
     597    31644416 :       ulong ub = (ulong)*b, ur = (ulong)r;
     598    31644416 :       uz = ((ub + ur) >> 1);
     599             :     }
     600   104465028 :     *c -= q * uz; /* c -= qz */
     601             :   }
     602   164916344 :   *b = r;
     603             : }
     604             : static void
     605     3421649 : REDBU(GEN a, GEN *b, GEN *c, GEN u1, GEN *u2)
     606             : { /* REDB(a,b,c) */
     607     3421649 :   GEN r, q = dvmdii_round(*b, a, &r);
     608     3421649 :   *c = subii(*c, mulii(q, shifti(addii(*b, r),-1)));
     609     3421649 :   *b = r;
     610     3421649 :   *u2 = subii(*u2, mulii(q, u1));
     611     3421649 : }
     612             : 
     613             : /* q t_QFI, return reduced representative and set base change U in Sl2(Z) */
     614             : GEN
     615     1967301 : redimagsl2(GEN q, GEN *U)
     616             : {
     617     1967301 :   GEN Q = cgetg(4, t_QFI);
     618     1967301 :   pari_sp av = avma, av2;
     619     1967301 :   GEN z, u1,u2,v1,v2, a = gel(q,1), b = gel(q,2), c = gel(q,3);
     620             :   long cmp;
     621             :   /* upper bound for size of final (a,b,c) */
     622     1967301 :   (void)new_chunk(2*(lgefint(a) + lgefint(b) + lgefint(c) + 3));
     623     1967301 :   av2 = avma;
     624     1967301 :   u1 = gen_1; u2 = gen_0;
     625     1967301 :   cmp = abscmpii(a, b);
     626     1967301 :   if (cmp < 0)
     627      216293 :     REDBU(a,&b,&c, u1,&u2);
     628     1751008 :   else if (cmp == 0 && signe(b) < 0)
     629             :   { /* b = -a */
     630           0 :     b = negi(b);
     631           0 :     u2 = gen_1;
     632             :   }
     633             :   for(;;)
     634             :   {
     635     5172657 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     636     3205356 :     swap(a,c); b = negi(b);
     637     3205356 :     z = u1; u1 = u2; u2 = negi(z);
     638     3205356 :     REDBU(a,&b,&c, u1,&u2);
     639     3205356 :     if (gc_needed(av, 1)) {
     640           0 :       if (DEBUGMEM>1) pari_warn(warnmem, "redimagsl2");
     641           0 :       gerepileall(av2, 5, &a,&b,&c, &u1,&u2);
     642             :     }
     643     3205356 :   }
     644     1967301 :   if (cmp == 0 && signe(b) < 0)
     645             :   {
     646        7168 :     b = negi(b);
     647        7168 :     z = u1; u1 = u2; u2 = negi(z);
     648             :   }
     649     1967301 :   avma = av;
     650     1967301 :   a = icopy(a); gel(Q,1) = a;
     651     1967301 :   b = icopy(b); gel(Q,2) = b;
     652     1967301 :   c = icopy(c); gel(Q,3) = c;
     653     1967301 :   u1 = icopy(u1);
     654     1967301 :   u2 = icopy(u2); av = avma;
     655             : 
     656             :   /* Let q = (A,B,C). q o [u1,u2; v1,v2] = Q implies
     657             :    * [v1,v2] = (1/C) [(b-B)/2 u1 - a u2, c u1 - (b+B)/2 u2] */
     658     1967301 :   z = shifti(subii(b, gel(q,2)), -1);
     659     1967301 :   v1 = subii(mulii(z, u1), mulii(a, u2)); v1 = diviiexact(v1, gel(q,3));
     660     1967301 :   z = subii(z, b);
     661     1967301 :   v2 = addii(mulii(z, u2), mulii(c, u1)); v2 = diviiexact(v2, gel(q,3));
     662     1967301 :   avma = av;
     663     1967301 :   v1 = icopy(v1);
     664     1967301 :   v2 = icopy(v2);
     665     1967301 :   *U = mkmat2(mkcol2(u1,v1), mkcol2(u2,v2)); return Q;
     666             : }
     667             : 
     668             : static GEN
     669      309027 : setq_b0(ulong a, ulong c)
     670      309027 : { retmkqfi( utoipos(a), gen_0, utoipos(c) ); }
     671             : /* assume |sb| = 1 */
     672             : static GEN
     673   141903020 : setq(ulong a, ulong b, ulong c, long sb)
     674   141903020 : { retmkqfi( utoipos(a), sb == 1? utoipos(b): utoineg(b), utoipos(c) ); }
     675             : /* 0 < a, c < 2^BIL, b = 0 */
     676             : static GEN
     677       77107 : redimag_1_b0(ulong a, ulong c)
     678       77107 : { return (a <= c)? setq_b0(a, c): setq_b0(c, a); }
     679             : 
     680             : /* 0 < a, c < 2^BIL: single word affair */
     681             : static GEN
     682   142349935 : redimag_1(pari_sp av, GEN a, GEN b, GEN c)
     683             : {
     684             :   ulong ua, ub, uc;
     685             :   long sb;
     686             :   for(;;)
     687             :   { /* at most twice */
     688   142349935 :     long lb = lgefint(b); /* <= 3 after first loop */
     689   142349935 :     if (lb == 2) return redimag_1_b0(a[2],c[2]);
     690   142272828 :     if (lb == 3 && uel(b,2) <= (ulong)LONG_MAX) break;
     691      137888 :     REDB(a,&b,&c);
     692      137888 :     if (uel(a,2) <= uel(c,2))
     693             :     { /* lg(b) <= 3 but may be too large for itos */
     694           0 :       long s = signe(b);
     695           0 :       avma = av;
     696           0 :       if (!s) return redimag_1_b0(a[2], c[2]);
     697           0 :       if (a[2] == c[2]) s = 1;
     698           0 :       return setq(a[2], b[2], c[2], s);
     699             :     }
     700      137888 :     swap(a,c); b = negi(b);
     701      137888 :   }
     702             :   /* b != 0 */
     703   142134940 :   avma = av;
     704   142134940 :   ua = a[2];
     705   142134940 :   ub = sb = b[2]; if (signe(b) < 0) sb = -sb;
     706   142134940 :   uc = c[2];
     707   142134940 :   if (ua < ub)
     708    48513489 :     sREDB(ua, &sb, &uc);
     709    93621451 :   else if (ua == ub && sb < 0) sb = (long)ub;
     710   417100730 :   while(ua > uc)
     711             :   {
     712   132830850 :     lswap(ua,uc); sb = -sb;
     713   132830850 :     sREDB(ua, &sb, &uc);
     714             :   }
     715   142134940 :   if (!sb) return setq_b0(ua, uc);
     716             :   else
     717             :   {
     718   141903020 :     long s = 1;
     719   141903020 :     if (sb < 0)
     720             :     {
     721    64396621 :       sb = -sb;
     722    64396621 :       if (ua != uc) s = -1;
     723             :     }
     724   141903020 :     return setq(ua, sb, uc, s);
     725             :   }
     726             : }
     727             : 
     728             : static GEN
     729   142386023 : redimag_av(pari_sp av, GEN q)
     730             : {
     731   142386023 :   GEN a = gel(q,1), b = gel(q,2), c = gel(q,3);
     732   142386023 :   long cmp, lc = lgefint(c);
     733             : 
     734   142386023 :   if (lgefint(a) == 3 && lc == 3) return redimag_1(av, a, b, c);
     735      892885 :   cmp = abscmpii(a, b);
     736      892885 :   if (cmp < 0)
     737      421986 :     REDB(a,&b,&c);
     738      470899 :   else if (cmp == 0 && signe(b) < 0)
     739          17 :     b = negi(b);
     740             :   for(;;)
     741             :   {
     742     3029502 :     cmp = abscmpii(a, c); if (cmp <= 0) break;
     743     2855526 :     lc = lgefint(a); /* lg(future c): we swap a & c next */
     744     2855526 :     if (lc == 3) return redimag_1(av, a, b, c);
     745     2136617 :     swap(a,c); b = negi(b); /* apply rho */
     746     2136617 :     REDB(a,&b,&c);
     747     2136617 :   }
     748      173976 :   if (cmp == 0 && signe(b) < 0) b = negi(b);
     749             :   /* size of reduced Qfb(a,b,c) <= 3 lg(c) + 4 <= 4 lg(c) */
     750      173976 :   (void)new_chunk(lc<<2);
     751      173976 :   a = icopy(a); b = icopy(b); c = icopy(c);
     752      173976 :   avma = av;
     753      173976 :   retmkqfi(icopy(a), icopy(b), icopy(c));
     754             : }
     755             : GEN
     756    36146133 : redimag(GEN q) { return redimag_av(avma, q); }
     757             : 
     758             : static GEN
     759           7 : rhoimag(GEN x)
     760             : {
     761           7 :   GEN a = gel(x,1), b = gel(x,2), c = gel(x,3);
     762           7 :   int fl = abscmpii(a, c);
     763           7 :   if (fl <= 0) {
     764           7 :     int fg = abscmpii(a, b);
     765           7 :     if (fg >= 0) {
     766           7 :       x = qfi(a,b,c);
     767           7 :       if ((!fl || !fg) && signe(gel(x,2)) < 0) setsigne(gel(x,2), 1);
     768           7 :       return x;
     769             :     }
     770             :   }
     771           0 :   x = cgetg(4, t_QFI);
     772           0 :   (void)new_chunk(lgefint(a) + lgefint(b) + lgefint(c) + 3);
     773           0 :   swap(a,c); b = negi(b);
     774           0 :   REDB(a, &b, &c); avma = (pari_sp)x;
     775           0 :   gel(x,1) = icopy(a);
     776           0 :   gel(x,2) = icopy(b);
     777           0 :   gel(x,3) = icopy(c); return x;
     778             : }
     779             : 
     780             : /* qfr3 / qfr5 */
     781             : 
     782             : /* t_QFR are unusable: D, sqrtD, isqrtD are recomputed all the time and the
     783             :  * logarithmic Shanks's distance is costly and hard to control.
     784             :  * qfr3 / qfr5 routines take a container of t_INTs (e.g a t_VEC) as argument,
     785             :  * at least 3 (resp. 5) components [it is a feature that they do not check the
     786             :  * precise type or length of the input]. They return a vector of length 3
     787             :  * (resp. 5). A qfr3 [a,b,c] contains the form coeffs, in a qfr5 [a,b,c, e,d]
     788             :  * the t_INT e is a binary exponent, d a t_REAL, coding the distance in
     789             :  * multiplicative form: the true distance is obtained from qfr5_dist.
     790             :  * All other qfr routines are obsolete (inefficient) wrappers */
     791             : 
     792             : /* static functions are not stack-clean. Unless mentionned otherwise, public
     793             :  * functions are. */
     794             : 
     795             : #define EMAX 22
     796             : static void
     797    12022346 : fix_expo(GEN x)
     798             : {
     799    12022346 :   if (expo(gel(x,5)) >= (1L << EMAX)) {
     800           0 :     gel(x,4) = addsi(1, gel(x,4));
     801           0 :     shiftr_inplace(gel(x, 5), - (1L << EMAX));
     802             :   }
     803    12022346 : }
     804             : 
     805             : /* (1/2) log (d * 2^{e * 2^EMAX}). Not stack clean if e != 0 */
     806             : GEN
     807      184254 : qfr5_dist(GEN e, GEN d, long prec)
     808             : {
     809      184254 :   GEN t = logr_abs(d);
     810      184254 :   if (signe(e)) {
     811           0 :     GEN u = mulir(e, mplog2(prec));
     812           0 :     shiftr_inplace(u, EMAX); t = addrr(t, u);
     813             :   }
     814      184254 :   shiftr_inplace(t, -1); return t;
     815             : }
     816             : 
     817             : static void
     818    17288489 : rho_get_BC(GEN *B, GEN *C, GEN b, GEN c, struct qfr_data *S)
     819             : {
     820             :   GEN t, u;
     821    17288489 :   u = shifti(c,1);
     822    17288489 :   t = (abscmpii(S->isqrtD,c) >= 0)? S->isqrtD: c;
     823    17288489 :   u = remii(addii_sign(t,1, b,signe(b)), u);
     824    17288489 :   *B = addii_sign(t, 1, u, -signe(u)); /* |t| - (|t|+b) % |2c| */
     825    17288489 :   if (*B == gen_0)
     826        3343 :   { u = shifti(S->D, -2); setsigne(u, -1); }
     827             :   else
     828    17285146 :     u = shifti(addii_sign(sqri(*B),1, S->D,-1), -2);
     829    17288489 :   *C = diviiexact(u, c); /* = (B^2-D)/4c */
     830    17288489 : }
     831             : /* Not stack-clean */
     832             : GEN
     833     6992791 : qfr3_rho(GEN x, struct qfr_data *S)
     834             : {
     835     6992791 :   GEN B, C, b = gel(x,2), c = gel(x,3);
     836     6992791 :   rho_get_BC(&B, &C, b, c, S);
     837     6992791 :   return mkvec3(c,B,C);
     838             : }
     839             : /* Not stack-clean */
     840             : GEN
     841    10295698 : qfr5_rho(GEN x, struct qfr_data *S)
     842             : {
     843    10295698 :   GEN B, C, y, b = gel(x,2), c = gel(x,3);
     844    10295698 :   long sb = signe(b);
     845             : 
     846    10295698 :   rho_get_BC(&B, &C, b, c, S);
     847    10295698 :   y = mkvec5(c,B,C, gel(x,4), gel(x,5));
     848    10295698 :   if (sb) {
     849    10291995 :     GEN t = subii(sqri(b), S->D);
     850    10291995 :     if (sb < 0)
     851     2299003 :       t = divir(t, sqrr(subir(b,S->sqrtD)));
     852             :     else
     853     7992992 :       t = divri(sqrr(addir(b,S->sqrtD)), t);
     854             :     /* t = (b + sqrt(D)) / (b - sqrt(D)), evaluated stably */
     855    10291995 :     gel(y,5) = mulrr(t, gel(y,5)); fix_expo(y);
     856             :   }
     857    10295698 :   return y;
     858             : }
     859             : 
     860             : /* Not stack-clean */
     861             : GEN
     862      217021 : qfr_to_qfr5(GEN x, long prec)
     863      217021 : { return mkvec5(gel(x,1),gel(x,2),gel(x,3),gen_0,real_1(prec)); }
     864             : 
     865             : /* d0 = initial distance, x = [a,b,c, expo(d), d], d = exp(2*distance) */
     866             : GEN
     867          84 : qfr5_to_qfr(GEN x, GEN d0)
     868             : {
     869             :   GEN y;
     870          84 :   if (lg(x) ==  6)
     871             :   {
     872          70 :     GEN n = gel(x,4), d = absr(gel(x,5));
     873          70 :     if (signe(n))
     874             :     {
     875           0 :       n = addis(shifti(n, EMAX), expo(d));
     876           0 :       setexpo(d, 0); d = logr_abs(d);
     877           0 :       if (signe(n)) d = addrr(d, mulir(n, mplog2(lg(d0))));
     878           0 :       shiftr_inplace(d, -1);
     879           0 :       d0 = addrr(d0, d);
     880             :     }
     881          70 :     else if (!gequal1(d)) /* avoid loss of precision */
     882             :     {
     883          35 :       d = logr_abs(d);
     884          35 :       shiftr_inplace(d, -1);
     885          35 :       d0 = addrr(d0, d);
     886             :     }
     887             :   }
     888          84 :   y = cgetg(5, t_QFR);
     889          84 :   gel(y,1) = gel(x,1);
     890          84 :   gel(y,2) = gel(x,2);
     891          84 :   gel(y,3) = gel(x,3);
     892          84 :   gel(y,4) = d0; return y;
     893             : }
     894             : 
     895             : /* Not stack-clean */
     896             : GEN
     897      819882 : qfr3_to_qfr(GEN x, GEN d)
     898             : {
     899      819882 :   GEN z = cgetg(5, t_QFR);
     900      819882 :   gel(z,1) = gel(x,1);
     901      819882 :   gel(z,2) = gel(x,2);
     902      819882 :   gel(z,3) = gel(x,3);
     903      819882 :   gel(z,4) = d; return z;
     904             : }
     905             : 
     906             : static int
     907    24899124 : ab_isreduced(GEN a, GEN b, GEN isqrtD)
     908             : {
     909    24899124 :   if (signe(b) > 0 && abscmpii(b, isqrtD) <= 0)
     910             :   {
     911     7108637 :     GEN t = addii_sign(isqrtD,1, shifti(a,1),-1);
     912     7108637 :     long l = abscmpii(b, t); /* compare |b| and |floor(sqrt(D)) - |2a|| */
     913     7108637 :     if (l > 0 || (l == 0 && signe(t) < 0)) return 1;
     914             :   }
     915    19416706 :   return 0;
     916             : }
     917             : 
     918             : INLINE int
     919    17273016 : qfr_isreduced(GEN x, GEN isqrtD)
     920             : {
     921    17273016 :   return ab_isreduced(gel(x,1),gel(x,2),isqrtD);
     922             : }
     923             : 
     924             : /* Not stack-clean */
     925             : GEN
     926     1947344 : qfr5_red(GEN x, struct qfr_data *S) {
     927     1947344 :   while (!qfr_isreduced(x,S->isqrtD)) x = qfr5_rho(x,S);
     928     1947344 :   return x;
     929             : }
     930             : /* Not stack-clean */
     931             : GEN
     932     1169242 : qfr3_red(GEN x, struct qfr_data *S) {
     933     1169242 :   while (!qfr_isreduced(x, S->isqrtD)) x = qfr3_rho(x, S);
     934     1169242 :   return x;
     935             : }
     936             : 
     937             : static void
     938          91 : get_disc(GEN x, struct qfr_data *S)
     939             : {
     940          91 :   if (!S->D) S->D = qfb_disc(x);
     941           0 :   else if (typ(S->D) != t_INT) pari_err_TYPE("qfr_init",S->D);
     942          91 :   if (!signe(S->D)) pari_err_DOMAIN("qfr_init", "disc", "=", gen_0,x);
     943          91 : }
     944             : 
     945             : void
     946        2184 : qfr_data_init(GEN D, long prec, struct qfr_data *S)
     947             : {
     948        2184 :   S->D = D;
     949        2184 :   S->sqrtD = sqrtr(itor(S->D,prec));
     950        2184 :   S->isqrtD = truncr(S->sqrtD);
     951        2184 : }
     952             : 
     953             : static GEN
     954          70 : qfr5_init(GEN x, struct qfr_data *S)
     955             : {
     956          70 :   GEN d = gel(x,4);
     957          70 :   long prec = realprec(d), l = -expo(d);
     958          70 :   if (l < BITS_IN_LONG) l = BITS_IN_LONG;
     959          70 :   prec = maxss(prec, nbits2prec(l));
     960          70 :   x = qfr_to_qfr5(x,prec);
     961             : 
     962          70 :   get_disc(x, S);
     963          70 :   if (!S->sqrtD) S->sqrtD = sqrtr(itor(S->D,prec));
     964           0 :   else if (typ(S->sqrtD) != t_REAL) pari_err_TYPE("qfr_init",S->sqrtD);
     965             : 
     966          70 :   if (!S->isqrtD)
     967             :   {
     968          70 :     pari_sp av=avma;
     969             :     long e;
     970          70 :     S->isqrtD = gcvtoi(S->sqrtD,&e);
     971          70 :     if (e>-2) { avma = av; S->isqrtD = sqrti(S->D); }
     972             :   }
     973           0 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
     974          70 :   return x;
     975             : }
     976             : static GEN
     977          21 : qfr3_init(GEN x, struct qfr_data *S)
     978             : {
     979          21 :   get_disc(x, S);
     980          21 :   if (!S->isqrtD) S->isqrtD = sqrti(S->D);
     981          14 :   else if (typ(S->isqrtD) != t_INT) pari_err_TYPE("qfr_init",S->isqrtD);
     982          21 :   return x;
     983             : }
     984             : 
     985             : #define qf_NOD  2
     986             : #define qf_STEP 1
     987             : 
     988             : static GEN
     989          77 : redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
     990             : {
     991          77 :   pari_sp av = avma;
     992             :   struct qfr_data S;
     993             :   GEN d;
     994          77 :   if (typ(x) != t_QFR) pari_err_TYPE("redreal",x);
     995          77 :   d = gel(x,4);
     996          77 :   S.D = D;
     997          77 :   S.sqrtD = sqrtD;
     998          77 :   S.isqrtD = isqrtD;
     999          77 :   x = (flag & qf_NOD)? qfr3_init(x, &S): qfr5_init(x, &S);
    1000          77 :   switch(flag) {
    1001          49 :     case 0:              x = qfr5_red(x,&S); break;
    1002           7 :     case qf_NOD:         x = qfr3_red(x,&S); break;
    1003          14 :     case qf_STEP:        x = qfr5_rho(x,&S); break;
    1004           7 :     case qf_STEP|qf_NOD: x = qfr3_rho(x,&S); break;
    1005           0 :     default: pari_err_FLAG("qfbred");
    1006             :   }
    1007          77 :   return gerepilecopy(av, qfr5_to_qfr(x,d));
    1008             : }
    1009             : GEN
    1010          42 : redreal(GEN x)
    1011          42 : { return redreal0(x,0,NULL,NULL,NULL); }
    1012             : GEN
    1013           0 : rhoreal(GEN x)
    1014           0 : { return redreal0(x,qf_STEP,NULL,NULL,NULL); }
    1015             : GEN
    1016           0 : redrealnod(GEN x, GEN isqrtD)
    1017           0 : { return redreal0(x,qf_NOD,NULL,isqrtD,NULL); }
    1018             : GEN
    1019           0 : rhorealnod(GEN x, GEN isqrtD)
    1020           0 : { return redreal0(x,qf_STEP|qf_NOD,NULL,isqrtD,NULL); }
    1021             : GEN
    1022          56 : qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
    1023             : {
    1024          56 :   if (typ(x) == t_QFI)
    1025          21 :     return (flag & qf_STEP)? rhoimag(x): redimag(x);
    1026          35 :   return redreal0(x,flag,D,isqrtD,sqrtD);
    1027             : }
    1028             : 
    1029             : GEN
    1030     1730351 : qfr5_comp(GEN x, GEN y, struct qfr_data *S)
    1031             : {
    1032     1730351 :   pari_sp av = avma;
    1033     1730351 :   GEN z = cgetg(6,t_VEC); qfb_comp(z,x,y);
    1034     1730351 :   if (x == y)
    1035             :   {
    1036       34062 :     gel(z,4) = shifti(gel(x,4),1);
    1037       34062 :     gel(z,5) = sqrr(gel(x,5));
    1038             :   }
    1039             :   else
    1040             :   {
    1041     1696289 :     gel(z,4) = addii(gel(x,4),gel(y,4));
    1042     1696289 :     gel(z,5) = mulrr(gel(x,5),gel(y,5));
    1043             :   }
    1044     1730351 :   fix_expo(z); z = qfr5_red(z,S);
    1045     1730351 :   return gerepilecopy(av,z);
    1046             : }
    1047             : /* Not stack-clean */
    1048             : GEN
    1049     1006569 : qfr3_comp(GEN x, GEN y, struct qfr_data *S)
    1050             : {
    1051     1006569 :   GEN z = cgetg(4,t_VEC); qfb_comp(z,x,y);
    1052     1006569 :   return qfr3_red(z, S);
    1053             : }
    1054             : 
    1055             : /* valid for t_QFR, qfr3, qfr5 */
    1056             : static GEN
    1057        3577 : qfr_inv(GEN x) {
    1058        3577 :   GEN z = shallowcopy(x);
    1059        3577 :   gel(z,2) = negi(gel(z,2));
    1060        3577 :   return z;
    1061             : }
    1062             : 
    1063             : /* return x^n. Not stack-clean */
    1064             : GEN
    1065           7 : qfr5_pow(GEN x, GEN n, struct qfr_data *S)
    1066             : {
    1067           7 :   GEN y = NULL;
    1068           7 :   long i, m, s = signe(n);
    1069           7 :   if (!s) return qfr5_1(S, lg(gel(x,5)));
    1070          22 :   for (i=lgefint(n)-1; i>1; i--)
    1071             :   {
    1072          15 :     m = n[i];
    1073          22 :     for (; m; m>>=1)
    1074             :     {
    1075          14 :       if (m&1) y = y? qfr5_comp(y,x,S): x;
    1076          14 :       if (m == 1 && i == 2) break;
    1077           7 :       x = qfr5_comp(x,x,S);
    1078             :     }
    1079             :   }
    1080           7 :   return y;
    1081             : }
    1082             : /* return x^n. Not stack-clean */
    1083             : GEN
    1084        4396 : qfr3_pow(GEN x, GEN n, struct qfr_data *S)
    1085             : {
    1086        4396 :   GEN y = NULL;
    1087        4396 :   long i, m, s = signe(n);
    1088        4396 :   if (!s) return qfr3_1(S);
    1089        4396 :   if (s < 0) x = qfr_inv(x);
    1090        8800 :   for (i=lgefint(n)-1; i>1; i--)
    1091             :   {
    1092        4404 :     m = n[i];
    1093        4863 :     for (; m; m>>=1)
    1094             :     {
    1095        4851 :       if (m&1) y = y? qfr3_comp(y,x,S): x;
    1096        4851 :       if (m == 1 && i == 2) break;
    1097         459 :       x = qfr3_comp(x,x,S);
    1098             :     }
    1099             :   }
    1100        4396 :   return y;
    1101             : }
    1102             : GEN
    1103          14 : qfrpow(GEN x, GEN n)
    1104             : {
    1105          14 :   struct qfr_data S = { NULL, NULL, NULL };
    1106          14 :   long s = signe(n);
    1107          14 :   pari_sp av = avma;
    1108             :   GEN d0;
    1109             : 
    1110          14 :   if (!s) return qfr_1(x);
    1111          14 :   if (is_pm1(n)) return s > 0? redreal(x): ginv(x);
    1112          14 :   if (s < 0) x = qfr_inv(x);
    1113          14 :   d0 = gel(x,4);
    1114          14 :   if (!signe(d0)) {
    1115           7 :     x = qfr3_init(x, &S);
    1116           7 :     x = qfr3_pow(x, n, &S);
    1117           7 :     x = qfr3_to_qfr(x, d0);
    1118             :   } else {
    1119           7 :     x = qfr5_init(x, &S);
    1120           7 :     x = qfr5_pow(qfr_to_qfr5(x, lg(S.sqrtD)), n, &S);
    1121           7 :     x = qfr5_to_qfr(x, mulri(d0,n));
    1122             :   }
    1123          14 :   return gerepilecopy(av, x);
    1124             : }
    1125             : 
    1126             : /* Prime forms attached to prime ideals of degree 1 */
    1127             : 
    1128             : /* assume x != 0 a t_INT, p > 0
    1129             :  * Return a t_QFI, but discriminant sign is not checked: can be used for
    1130             :  * real forms as well */
    1131             : GEN
    1132    39409794 : primeform_u(GEN x, ulong p)
    1133             : {
    1134    39409794 :   GEN c, y = cgetg(4, t_QFI);
    1135    39409794 :   pari_sp av = avma;
    1136             :   ulong b;
    1137             :   long s;
    1138             : 
    1139    39409794 :   s = mod8(x); if (signe(x) < 0 && s) s = 8-s;
    1140             :   /* 2 or 3 mod 4 */
    1141    39409794 :   if (s & 2) pari_err_DOMAIN("primeform", "disc % 4", ">",gen_1, x);
    1142    39409787 :   if (p == 2) {
    1143     3334247 :     switch(s) {
    1144       25967 :       case 0: b = 0; break;
    1145     3296342 :       case 1: b = 1; break;
    1146       11938 :       case 4: b = 2; break;
    1147           0 :       default: pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1148           0 :                b = 0; /* -Wall */
    1149             :     }
    1150     3334247 :     c = shifti(subsi(s,x), -3);
    1151             :   } else {
    1152    36075540 :     b = Fl_sqrt(umodiu(x,p), p);
    1153    36075540 :     if (b == ~0UL) pari_err_SQRTN("primeform", mkintmod(x,utoi(p)) );
    1154             :     /* mod(b) != mod2(x) ? */
    1155    36075540 :     if ((b ^ s) & 1) b = p - b;
    1156    36075540 :     c = diviuexact(shifti(subii(sqru(b), x), -2), p);
    1157             :   }
    1158    39409787 :   gel(y,3) = gerepileuptoint(av, c);
    1159    39409787 :   gel(y,2) = utoi(b);
    1160    39409787 :   gel(y,1) = utoipos(p); return y;
    1161             : }
    1162             : 
    1163             : /* special case: p = 1 return unit form */
    1164             : GEN
    1165     1758099 : primeform(GEN x, GEN p, long prec)
    1166             : {
    1167     1758099 :   const char *f = "primeform";
    1168             :   pari_sp av;
    1169     1758099 :   long s, sx = signe(x), sp = signe(p);
    1170             :   GEN y, b, absp;
    1171             : 
    1172     1758099 :   if (typ(x) != t_INT) pari_err_TYPE(f,x);
    1173     1758099 :   if (typ(p) != t_INT) pari_err_TYPE(f,p);
    1174     1758099 :   if (!sp) pari_err_DOMAIN(f,"p","=",gen_0,p);
    1175     1758099 :   if (!sx) pari_err_DOMAIN(f,"D","=",gen_0,x);
    1176     1758099 :   if (lgefint(p) == 3)
    1177             :   {
    1178     1758085 :     ulong pp = p[2];
    1179     1758085 :     if (pp == 1) {
    1180       13041 :       if (sx < 0) {
    1181             :         long r;
    1182        6475 :         if (sp < 0) pari_err_IMPL("negative definite t_QFI");
    1183        6475 :         r = mod4(x);
    1184        6475 :         if (r && r != 3) pari_err_DOMAIN(f,"disc % 4",">", gen_1,x);
    1185        6475 :         return qfi_1_by_disc(x);
    1186             :       }
    1187        6566 :       y = qfr_1_by_disc(x,prec);
    1188        6566 :       if (sp < 0) { gel(y,1) = gen_m1; togglesign(gel(y,3)); }
    1189        6566 :       return y;
    1190             :     }
    1191     1745044 :     y = primeform_u(x, pp);
    1192     1745037 :     if (sx < 0) {
    1193      928886 :       if (sp < 0) pari_err_IMPL("negative definite t_QFI");
    1194      928886 :       return y;
    1195             :     }
    1196      816151 :     if (sp < 0) { togglesign(gel(y,1)); togglesign(gel(y,3)); }
    1197      816151 :     return gcopy( qfr3_to_qfr(y, real_0(prec)) );
    1198             :   }
    1199          14 :   s = mod8(x);
    1200          14 :   if (sx < 0)
    1201             :   {
    1202           7 :     if (sp < 0) pari_err_IMPL("negative definite t_QFI");
    1203           7 :     if (s) s = 8-s;
    1204           7 :     y = cgetg(4, t_QFI);
    1205             :   }
    1206             :   else
    1207             :   {
    1208           7 :     y = cgetg(5, t_QFR);
    1209           7 :     gel(y,4) = real_0(prec);
    1210             :   }
    1211             :   /* 2 or 3 mod 4 */
    1212          14 :   if (s & 2) pari_err_DOMAIN(f, "disc % 4", ">",gen_1, x);
    1213          14 :   absp = absi(p); av = avma;
    1214          14 :   b = Fp_sqrt(x, absp); if (!b) pari_err_SQRTN(f, mkintmod(x,absp));
    1215          14 :   s &= 1; /* s = x mod 2 */
    1216             :   /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
    1217          14 :   if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(absp,b));
    1218             : 
    1219          14 :   av = avma;
    1220          14 :   gel(y,3) = gerepileuptoint(av, diviiexact(shifti(subii(sqri(b), x), -2), p));
    1221          14 :   gel(y,2) = b;
    1222          14 :   gel(y,1) = icopy(p);
    1223          14 :   return y;
    1224             : }
    1225             : 
    1226             : /* Let M and N in SL2(Z), return (N*M^-1)[,1] */
    1227             : static GEN
    1228       99575 : SL2_div_mul_e1(GEN N, GEN M)
    1229             : {
    1230       99575 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1231       99575 :   GEN p = subii(mulii(gcoeff(N,1,1), d), mulii(gcoeff(N,1,2), b));
    1232       99575 :   GEN q = subii(mulii(gcoeff(N,2,1), d), mulii(gcoeff(N,2,2), b));
    1233       99575 :   return mkvec2(p,q);
    1234             : }
    1235             : /* Let M and N in SL2(Z), return (N*[1,0;0,-1]*M^-1)[,1] */
    1236             : static GEN
    1237       21049 : SL2_swap_div_mul_e1(GEN N, GEN M)
    1238             : {
    1239       21049 :   GEN b = gcoeff(M,2,1), d = gcoeff(M,2,2);
    1240       21049 :   GEN p = addii(mulii(gcoeff(N,1,1), d), mulii(gcoeff(N,1,2), b));
    1241       21049 :   GEN q = addii(mulii(gcoeff(N,2,1), d), mulii(gcoeff(N,2,2), b));
    1242       21049 :   return mkvec2(p,q);
    1243             : }
    1244             : 
    1245             : /* Test equality modulo GL2 of two reduced forms */
    1246             : static int
    1247      899458 : GL2_qfb_equal(GEN a, GEN b)
    1248             : {
    1249     1798916 :   return equalii(gel(a,1),gel(b,1))
    1250       58590 :    && absequalii(gel(a,2),gel(b,2))
    1251      951874 :    &&    equalii(gel(a,3),gel(b,3));
    1252             : }
    1253             : 
    1254             : static GEN
    1255      193942 : qfbsolve_cornacchia(GEN c, GEN p, int swap)
    1256             : {
    1257      193942 :   pari_sp av = avma;
    1258             :   GEN M, N;
    1259      193942 :   if (kronecker(negi(c), p) < 0 || !cornacchia(c, p, &M,&N)) {
    1260      162946 :     avma = av; return gen_0;
    1261             :   }
    1262       30996 :   return gerepilecopy(av, swap? mkvec2(N,M): mkvec2(M,N));
    1263             : }
    1264             : 
    1265             : GEN
    1266     2259880 : qfisolvep(GEN Q, GEN p)
    1267             : {
    1268             :   GEN M, N, x,y, a,b,c, d;
    1269     2259880 :   pari_sp av = avma;
    1270     2259880 :   if (!signe(gel(Q,2)))
    1271             :   {
    1272      177142 :     a = gel(Q,1);
    1273      177142 :     c = gel(Q,3); /* if principal form, use faster cornacchia */
    1274      177142 :     if (equali1(a)) return qfbsolve_cornacchia(c, p, 0);
    1275       78932 :     if (equali1(c)) return qfbsolve_cornacchia(a, p, 1);
    1276             :   }
    1277     2085006 :   d = qfb_disc(Q); if (kronecker(d,p) < 0) return gen_0;
    1278     1035524 :   a = redimagsl2(Q, &N);
    1279     1035524 :   if (is_pm1(gel(a,1))) /* principal form */
    1280             :   {
    1281             :     long r;
    1282      136066 :     if (!signe(gel(a,2)))
    1283             :     {
    1284       19068 :       a = qfbsolve_cornacchia(gel(a,3), p, 0);
    1285       19068 :       if (a == gen_0) { avma = av; return gen_0; }
    1286        3178 :       a = ZM_ZC_mul(N, a);
    1287        3178 :       a[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1288        3178 :       return gerepileupto(av, a);
    1289             :     }
    1290             :     /* x^2 + xy + ((1-d)/4)y^2 = p <==> (2x + y)^2 - d y^2 = 4p */
    1291      116998 :     if (!cornacchia2(negi(d), p, &x, &y)) { avma = av; return gen_0; }
    1292        9926 :     x = divis_rem(subii(x,y), 2, &r); if (r) { avma = av; return gen_0; }
    1293        9926 :     a = ZM_ZC_mul(N, mkvec2(x,y));
    1294        9926 :     a[0] = evaltyp(t_VEC) | _evallg(3); /* transpose */
    1295        9926 :     return gerepileupto(av, a);
    1296             :   }
    1297      899458 :   b = redimagsl2(primeform(d, p, 0), &M);
    1298      899458 :   if (!GL2_qfb_equal(a,b)) { avma = av; return gen_0; }
    1299       52416 :   if (signe(gel(a,2))==signe(gel(b,2)))
    1300       31367 :     x = SL2_div_mul_e1(N,M);
    1301             :   else
    1302       21049 :     x = SL2_swap_div_mul_e1(N,M);
    1303       52416 :   return gerepilecopy(av, x);
    1304             : }
    1305             : 
    1306             : GEN
    1307     2550268 : redrealsl2step(GEN A, GEN d, GEN rd)
    1308             : {
    1309     2550268 :   pari_sp ltop = avma;
    1310     2550268 :   GEN N, V = gel(A,1), M = gel(A,2);
    1311     2550268 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3);
    1312     2550268 :   GEN C = mpabs(c);
    1313     2550268 :   GEN t = addii(b, gmax(rd, C));
    1314     2550268 :   GEN r, q = truedvmdii(t, shifti(C,1), &r);
    1315     2550268 :   b = subii(t, addii(r,b));
    1316     2550268 :   a = c;
    1317     2550268 :   c = truedivii(subii(sqri(b), d), shifti(c,2));
    1318     2550268 :   if (signe(a) < 0) togglesign(q);
    1319    10201072 :   N = mkmat2(gel(M,2),
    1320     5100536 :              mkcol2(subii(mulii(q, gcoeff(M, 1, 2)), gcoeff(M, 1, 1)),
    1321     5100536 :                     subii(mulii(q, gcoeff(M, 2, 2)), gcoeff(M, 2, 1))));
    1322     2550268 :   return gerepilecopy(ltop, mkvec2(mkvec3(a,b,c),N));
    1323             : }
    1324             : 
    1325             : GEN
    1326     2365832 : redrealsl2(GEN V, GEN d, GEN rd)
    1327             : {
    1328     2365832 :   pari_sp ltop = avma;
    1329             :   GEN M, u1, u2, v1, v2;
    1330     2365832 :   GEN a = gel(V,1), b = gel(V,2), c = gel(V,3);
    1331     2365832 :   u1 = v2 = gen_1; v1 = u2 = gen_0;
    1332     9991940 :   while (!ab_isreduced(a,b,rd))
    1333             :   {
    1334     5260276 :     GEN C = mpabs(c);
    1335     5260276 :     GEN t = addii(b, gmax(rd,C));
    1336     5260276 :     GEN r, q = truedvmdii(t, shifti(C,1), &r);
    1337     5260276 :     b = subii(t, addii(r,b));
    1338     5260276 :     a = c;
    1339     5260276 :     c = truedivii(subii(sqri(b), d), shifti(c,2));
    1340     5260276 :     if (signe(a) < 0) togglesign(q);
    1341     5260276 :     r = u1; u1 = v1; v1 = subii(mulii(q, v1), r);
    1342     5260276 :     r = u2; u2 = v2; v2 = subii(mulii(q, v2), r);
    1343     5260276 :     if (gc_needed(ltop, 1))
    1344             :     {
    1345           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"redrealsl2");
    1346           0 :       gerepileall(ltop, 7, &a,&b,&c,&u1,&u2,&v1,&v2);
    1347             :     }
    1348             :   }
    1349     2365832 :   M = mkmat2(mkcol2(u1,u2), mkcol2(v1,v2));
    1350     2365832 :   return gerepilecopy(ltop, mkvec2(mkvec3(a,b,c), M));
    1351             : }
    1352             : 
    1353             : GEN
    1354          35 : qfbredsl2(GEN q, GEN S)
    1355             : {
    1356             :   GEN v, D, isD;
    1357             :   pari_sp av;
    1358          35 :   switch(typ(q))
    1359             :   {
    1360             :     case t_QFI:
    1361           7 :       if (S) pari_err_TYPE("qfbredsl2",S);
    1362           7 :       v = cgetg(3,t_VEC);
    1363           7 :       gel(v,1) = redimagsl2(q, &gel(v,2));
    1364           7 :       return v;
    1365             :     case t_QFR:
    1366          28 :       av = avma;
    1367          28 :       if (S) {
    1368          21 :         if (typ(S) != t_VEC || lg(S) != 3) pari_err_TYPE("qfbredsl2",S);
    1369          14 :         D = gel(S,1);
    1370          14 :         isD = gel(S,2);
    1371          14 :         if (typ(D) != t_INT || signe(D) <= 0 || typ(isD) != t_INT)
    1372           7 :           pari_err_TYPE("qfbredsl2",S);
    1373             :       }
    1374             :       else
    1375             :       {
    1376           7 :         D = qfb_disc(q);
    1377           7 :         isD = sqrtint(D);
    1378             :       }
    1379          14 :       v = redrealsl2(q,D,isD);
    1380          14 :       gel(v,1) = qfr3_to_qfr(gel(v,1), real_0(precision(gel(q,4))));
    1381          14 :       return gerepilecopy(av, v);
    1382             : 
    1383             :     default:
    1384           0 :         pari_err_TYPE("qfbredsl2",q);
    1385           0 :         return NULL;
    1386             :   }
    1387             : }
    1388             : 
    1389             : GEN
    1390     1582882 : qfrsolvep(GEN Q, GEN p)
    1391             : {
    1392     1582882 :   pari_sp ltop = avma, btop;
    1393     1582882 :   GEN N, P, P1, P2, M, rd, d = qfb_disc(Q);
    1394     1582882 :   if (kronecker(d, p) < 0) { avma = ltop; return gen_0; }
    1395      788606 :   rd = sqrti(d);
    1396      788606 :   M = N = redrealsl2(Q, d,rd);
    1397      788606 :   P = primeform(d, p, DEFAULTPREC);
    1398      788606 :   P1 = redrealsl2(P, d,rd);
    1399      788606 :   togglesign( gel(P,2) );
    1400      788606 :   P2 = redrealsl2(P, d,rd);
    1401      788606 :   btop = avma;
    1402             :   for(;;)
    1403             :   {
    1404     2618476 :     if (ZV_equal(gel(M,1), gel(P1,1))) { N = gel(P1,2); break; }
    1405     2570575 :     if (ZV_equal(gel(M,1), gel(P2,1))) { N = gel(P2,2); break; }
    1406     2550268 :     M = redrealsl2step(M, d,rd);
    1407     2550268 :     if (ZV_equal(gel(M,1), gel(N,1))) { avma = ltop; return gen_0; }
    1408     1829870 :     if (gc_needed(btop, 1)) M = gerepileupto(btop, M);
    1409     1829870 :   }
    1410       68208 :   return gerepilecopy(ltop, SL2_div_mul_e1(gel(M,2),N));
    1411             : }
    1412             : 
    1413             : GEN
    1414     3842762 : qfbsolve(GEN Q,GEN n)
    1415             : {
    1416     3842762 :   if (typ(n)!=t_INT) pari_err_TYPE("qfbsolve",n);
    1417     3842762 :   switch(typ(Q))
    1418             :   {
    1419     2259880 :   case t_QFI: return qfisolvep(Q,n);
    1420     1582882 :   case t_QFR: return qfrsolvep(Q,n);
    1421             :   default:
    1422           0 :     pari_err_TYPE("qfbsolve",Q);
    1423           0 :     return NULL; /* NOT REACHED */
    1424             :   }
    1425             : }
    1426             : 
    1427             : /* 1 if there exists x,y such that x^2 + dy^2 = p [prime], 0 otherwise */
    1428             : long
    1429      116865 : cornacchia(GEN d, GEN p, GEN *px, GEN *py)
    1430             : {
    1431      116865 :   pari_sp av = avma, av2;
    1432             :   GEN a, b, c, L, r;
    1433             : 
    1434      116865 :   if (typ(d) != t_INT) pari_err_TYPE("cornacchia", d);
    1435      116865 :   if (typ(p) != t_INT) pari_err_TYPE("cornacchia", p);
    1436      116865 :   if (signe(d) <= 0) pari_err_DOMAIN("cornacchia", "d","<=",gen_0,d);
    1437      116865 :   *px = *py = gen_0;
    1438      116865 :   b = subii(p, d);
    1439      116865 :   if (signe(b) < 0) return 0;
    1440      115871 :   if (signe(b) == 0) { avma = av; *py = gen_1; return 1; }
    1441      115871 :   b = Fp_sqrt(b, p); /* sqrt(-d) */
    1442      115871 :   if (!b) { avma = av; return 0; }
    1443      115871 :   if (abscmpii(shifti(b,1), p) > 0) b = subii(b,p);
    1444      115871 :   a = p; L = sqrti(p);
    1445      115871 :   av2 = avma;
    1446      622244 :   while (abscmpii(b, L) > 0)
    1447             :   {
    1448      390502 :     r = remii(a, b); a = b; b = r;
    1449      390502 :     if (gc_needed(av2, 1)) {
    1450           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"cornacchia");
    1451           0 :       gerepileall(av2, 2, &a,&b);
    1452             :     }
    1453             :   }
    1454      115871 :   a = subii(p, sqri(b));
    1455      115871 :   c = dvmdii(a, d, &r);
    1456      115871 :   if (r != gen_0 || !Z_issquareall(c, &c)) { avma = av; return 0; }
    1457       31031 :   avma = av;
    1458       31031 :   *px = icopy(b);
    1459       31031 :   *py = icopy(c); return 1;
    1460             : }
    1461             : /* 1 if there exists x,y such that x^2 + dy^2 = 4p [p prime], 0 otherwise */
    1462             : long
    1463     1447439 : cornacchia2(GEN d, GEN p, GEN *px, GEN *py)
    1464             : {
    1465     1447439 :   pari_sp av = avma, av2;
    1466             :   GEN a, b, c, L, r, px4;
    1467             :   long k;
    1468             : 
    1469     1447439 :   if (typ(d) != t_INT) pari_err_TYPE("cornacchia2", d);
    1470     1447439 :   if (typ(p) != t_INT) pari_err_TYPE("cornacchia2", p);
    1471     1447439 :   if (signe(d) <= 0) pari_err_DOMAIN("cornacchia2", "d","<=",gen_0,d);
    1472     1447439 :   *px = *py = gen_0;
    1473     1447439 :   k = mod4(d);
    1474     1447439 :   if (k == 1 || k == 2) pari_err_DOMAIN("cornacchia2","-d mod 4", ">",gen_1,d);
    1475     1447439 :   px4 = shifti(p,2);
    1476     1447439 :   if (abscmpii(px4, d) < 0) { avma = av; return 0; }
    1477     1446277 :   if (absequaliu(p, 2))
    1478             :   {
    1479           0 :     avma = av;
    1480           0 :     switch (itou_or_0(d)) {
    1481           0 :       case 4: *px = gen_2; break;
    1482           0 :       case 7: *px = gen_1; break;
    1483           0 :       default: return 0;
    1484             :     }
    1485           0 :     *py = gen_1; return 1;
    1486             :   }
    1487     1446277 :   b = Fp_sqrt(negi(d), p);
    1488     1446277 :   if (!b) { avma = av; return 0; }
    1489     1446277 :   if (!signe(b)) { /* d = p,2p,3p,4p */
    1490          28 :     avma = av;
    1491          28 :     if (absequalii(d, px4)){ *py = gen_1; return 1; }
    1492          28 :     if (absequalii(d, p))  { *py = gen_2; return 1; }
    1493           0 :     return 0;
    1494             :   }
    1495     1446249 :   if (mod2(b) != (k & 1)) b = subii(p,b);
    1496     1446249 :   a = shifti(p,1); L = sqrti(px4);
    1497     1446249 :   av2 = avma;
    1498     9808666 :   while (cmpii(b, L) > 0)
    1499             :   {
    1500     6916168 :     r = remii(a, b); a = b; b = r;
    1501     6916168 :     if (gc_needed(av2, 1)) {
    1502           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"cornacchia");
    1503           0 :       gerepileall(av2, 2, &a,&b);
    1504             :     }
    1505             :   }
    1506     1446249 :   a = subii(px4, sqri(b));
    1507     1446249 :   c = dvmdii(a, d, &r);
    1508     1446249 :   if (r != gen_0 || !Z_issquareall(c, &c)) { avma = av; return 0; }
    1509     1340339 :   avma = av;
    1510     1340339 :   *px = icopy(b);
    1511     1340339 :   *py = icopy(c); return 1;
    1512             : }

Generated by: LCOV version 1.11