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

Generated by: LCOV version 1.13