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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - Qfb.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17097-9391e68) Lines: 835 904 92.4 %
Date: 2014-11-21 Functions: 91 97 93.8 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 422 567 74.4 %

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

Generated by: LCOV version 1.9