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

Generated by: LCOV version 1.11