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 - kernel/none - mp.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16741-1378b1c) Lines: 1108 1148 96.5 %
Date: 2014-08-17 Functions: 69 69 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 924 1110 83.2 %

           Branch data     Line data    Source code
       1                 :            : #line 2 "../src/kernel/none/mp.c"
       2                 :            : /* Copyright (C) 2000-2003 The PARI group.
       3                 :            : 
       4                 :            : This file is part of the PARI/GP package.
       5                 :            : 
       6                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       7                 :            : terms of the GNU General Public License as published by the Free Software
       8                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       9                 :            : ANY WARRANTY WHATSOEVER.
      10                 :            : 
      11                 :            : Check the License for details. You should have received a copy of it, along
      12                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      13                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14                 :            : 
      15                 :            : /***********************************************************************/
      16                 :            : /**                                                                   **/
      17                 :            : /**                       MULTIPRECISION KERNEL                       **/
      18                 :            : /**                                                                   **/
      19                 :            : /***********************************************************************/
      20                 :            : #include "pari.h"
      21                 :            : #include "paripriv.h"
      22                 :            : #include "../src/kernel/none/tune-gen.h"
      23                 :            : 
      24                 :        498 : void pari_kernel_init(void) { }
      25                 :        498 : void pari_kernel_close(void) { }
      26                 :            : 
      27                 :            : /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
      28                 :            :  * GENs but pairs (long *a, long na) representing a list of digits (in basis
      29                 :            :  * BITS_IN_LONG) : a[0], ..., a[na-1]. [ In ordre to facilitate splitting: no
      30                 :            :  * need to reintroduce codewords ] */
      31                 :            : 
      32                 :            : #define LIMBS(x)  ((x)+2)
      33                 :            : #define NLIMBS(x) (lgefint(x)-2)
      34                 :            : 
      35                 :            : /* Normalize a non-negative integer */
      36                 :            : GEN
      37                 :  149581212 : int_normalize(GEN x, long known_zero_words)
      38                 :            : {
      39                 :  149581212 :   long i, lx = lgefint(x);
      40                 :            :   GEN x0;
      41         [ -  + ]:  149581212 :   if (lx == 2) { x[1] = evalsigne(0) | evallgefint(2); return x; }
      42 [ +  + ][ +  + ]:  149581212 :   if (!known_zero_words && x[2]) return x;
      43         [ +  + ]:  627865503 :   for (i = 2+known_zero_words; i < lx; i++)
      44         [ +  + ]:  580677777 :     if (x[i]) break;
      45                 :   96806337 :   x0 = x; i -= 2; x += i;
      46         [ +  + ]:   96806337 :   if (x0 == (GEN)avma) avma = (pari_sp)x;
      47                 :   47429658 :   else stackdummy((pari_sp)(x0+i), (pari_sp)x0);
      48                 :   96806337 :   lx -= i;
      49                 :   96806337 :   x[0] = evaltyp(t_INT) | evallg(lx);
      50         [ +  + ]:   96806337 :   if (lx == 2) x[1] = evalsigne(0) | evallgefint(lx);
      51                 :   49618611 :   else         x[1] = evalsigne(1) | evallgefint(lx);
      52                 :  149581212 :   return x;
      53                 :            : }
      54                 :            : 
      55                 :            : /***********************************************************************/
      56                 :            : /**                                                                   **/
      57                 :            : /**                      ADDITION / SUBTRACTION                       **/
      58                 :            : /**                                                                   **/
      59                 :            : /***********************************************************************/
      60                 :            : 
      61                 :            : GEN
      62                 :    1098567 : setloop(GEN a)
      63                 :            : {
      64                 :    1098567 :   pari_sp av = avma;
      65                 :    1098567 :   (void)cgetg(lgefint(a) + 3, t_VECSMALL);
      66                 :    1098567 :   return icopy_avma(a, av); /* two cells of extra space before a */
      67                 :            : }
      68                 :            : 
      69                 :            : /* we had a = setloop(?), then some incloops. Reset a to b */
      70                 :            : GEN
      71                 :     109113 : resetloop(GEN a, GEN b) {
      72                 :     109113 :   long lb = lgefint(b);
      73                 :     109113 :   a += lgefint(a) - lb;
      74                 :     109113 :   a[0] = evaltyp(t_INT) | evallg(lb);
      75                 :     109113 :   affii(b, a); return a;
      76                 :            : }
      77                 :            : 
      78                 :            : /* assume a > 0, initialized by setloop. Do a++ */
      79                 :            : static GEN
      80                 :   10402101 : incpos(GEN a)
      81                 :            : {
      82                 :   10402101 :   long i, l = lgefint(a);
      83         [ +  + ]:   10402104 :   for (i=l-1; i>1; i--)
      84         [ +  + ]:   10402101 :     if (++a[i]) return a;
      85                 :          3 :   l++; a--; /* use extra cell */
      86                 :          3 :   a[0]=evaltyp(t_INT) | _evallg(l);
      87                 :          3 :   a[1]=evalsigne(1) | evallgefint(l);
      88                 :   10402101 :   a[2]=1; return a;
      89                 :            : }
      90                 :            : 
      91                 :            : /* assume a < 0, initialized by setloop. Do a++ */
      92                 :            : static GEN
      93                 :       4056 : incneg(GEN a)
      94                 :            : {
      95                 :       4056 :   long i, l = lgefint(a)-1;
      96         [ +  + ]:       4056 :   if (a[l]--)
      97                 :            :   {
      98 [ +  + ][ +  + ]:       4053 :     if (l == 2 && !a[2])
      99                 :            :     {
     100                 :        315 :       a++; /* save one cell */
     101                 :        315 :       a[0] = evaltyp(t_INT) | _evallg(2);
     102                 :        315 :       a[1] = evalsigne(0) | evallgefint(2);
     103                 :            :     }
     104                 :       4053 :     return a;
     105                 :            :   }
     106                 :          3 :   for (i = l-1;; i--) /* finishes since a[2] != 0 */
     107         [ +  - ]:          3 :     if (a[i]--) break;
     108         [ +  - ]:          3 :   if (!a[2])
     109                 :            :   {
     110                 :          3 :     a++; /* save one cell */
     111                 :          3 :     a[0] = evaltyp(t_INT) | _evallg(l);
     112                 :          3 :     a[1] = evalsigne(-1) | evallgefint(l);
     113                 :            :   }
     114                 :       4056 :   return a;
     115                 :            : }
     116                 :            : 
     117                 :            : /* assume a initialized by setloop. Do a++ */
     118                 :            : GEN
     119                 :   10412730 : incloop(GEN a)
     120                 :            : {
     121      [ +  +  + ]:   10412730 :   switch(signe(a))
     122                 :            :   {
     123                 :       6573 :     case 0: a--; /* use extra cell */
     124                 :       6573 :       a[0]=evaltyp(t_INT) | _evallg(3);
     125                 :       6573 :       a[1]=evalsigne(1) | evallgefint(3);
     126                 :       6573 :       a[2]=1; return a;
     127                 :       4056 :     case -1: return incneg(a);
     128                 :   10412730 :     default: return incpos(a);
     129                 :            :   }
     130                 :            : }
     131                 :            : 
     132                 :            : INLINE GEN
     133                 :  719356827 : adduispec(ulong s, GEN x, long nx)
     134                 :            : {
     135                 :  719356827 :   GEN xd, zd = (GEN)avma;
     136                 :            :   long lz;
     137                 :            : 
     138         [ +  + ]:  719356827 :   if (nx == 1) return adduu(s, uel(x,0));
     139                 :  152134959 :   lz = nx+3; (void)new_chunk(lz);
     140                 :  152134959 :   xd = x + nx;
     141                 :  152134959 :   *--zd = (ulong)*--xd + s;
     142         [ +  + ]:  152134959 :   if ((ulong)*zd < s)
     143                 :            :     for(;;)
     144                 :            :     {
     145         [ +  + ]:   11239230 :       if (xd == x) { *--zd = 1; break; } /* enlarge z */
     146                 :    7622196 :       *--zd = ((ulong)*--xd) + 1;
     147         [ +  + ]:    7622196 :       if (*zd) { lz--; break; }
     148                 :   11239230 :     }
     149                 :  146557977 :   else lz--;
     150         [ +  + ]:  441822360 :   while (xd > x) *--zd = *--xd;
     151                 :  152134959 :   *--zd = evalsigne(1) | evallgefint(lz);
     152                 :  152134959 :   *--zd = evaltyp(t_INT) | evallg(lz);
     153                 :  719356827 :   avma=(pari_sp)zd; return zd;
     154                 :            : }
     155                 :            : 
     156                 :            : GEN
     157                 :  183910515 : adduispec_offset(ulong s, GEN x, long offset, long nx)
     158                 :            : {
     159                 :  183910515 :   GEN xd = x+lgefint(x)-nx-offset;
     160 [ +  + ][ +  + ]:  221549322 :   while (nx && *xd==0) {xd++; nx--;}
     161         [ +  + ]:  183910515 :   if (!nx) return utoi(s);
     162                 :  183910515 :   return adduispec(s,xd,nx);
     163                 :            : }
     164                 :            : 
     165                 :            : static GEN
     166                 : 1109436027 : addiispec(GEN x, GEN y, long nx, long ny)
     167                 :            : {
     168                 :            :   GEN xd, yd, zd;
     169                 : 1109436027 :   long lz, i = -2;
     170                 :            :   LOCAL_OVERFLOW;
     171                 :            : 
     172         [ +  + ]: 1109436027 :   if (nx < ny) swapspec(x,y, nx,ny);
     173         [ +  + ]: 1109436027 :   if (ny == 1) return adduispec(*y,x,nx);
     174                 :  581302665 :   zd = (GEN)avma;
     175                 :  581302665 :   lz = nx+3; (void)new_chunk(lz);
     176                 :  581302665 :   xd = x + nx;
     177                 :  581302665 :   yd = y + ny;
     178                 :  581302665 :   zd[-1] = addll(xd[-1], yd[-1]);
     179                 :            : #ifdef addllx8
     180         [ +  + ]:  318441625 :   for (  ; i-8 > -ny; i-=8)
     181                 :  124674070 :     addllx8(xd+i, yd+i, zd+i, overflow);
     182                 :            : #endif
     183         [ +  + ]: 4918661659 :   for (  ; i >= -ny; i--) zd[i] = addllx(xd[i], yd[i]);
              [ +  +  + ]
     184         [ +  + ]:  581302665 :   if (overflow)
     185                 :            :     for(;;)
     186                 :            :     {
     187         [ +  + ]:   78068856 :       if (i < -nx) { zd[i] = 1; i--; break; } /* enlarge z */
     188                 :   59481903 :       zd[i] = uel(xd,i) + 1;
     189         [ +  + ]:   59481903 :       if (zd[i]) { i--; lz--; break; }
     190                 :   24774774 :       i--;
     191                 :   78068856 :     }
     192                 :  528008583 :   else lz--;
     193         [ +  + ]: 2476736634 :   for (; i >= -nx; i--) zd[i] = xd[i];
     194                 :  581302665 :   zd += i+1;
     195                 :  581302665 :   *--zd = evalsigne(1) | evallgefint(lz);
     196                 :  581302665 :   *--zd = evaltyp(t_INT) | evallg(lz);
     197                 : 1109436027 :   avma=(pari_sp)zd; return zd;
     198                 :            : }
     199                 :            : 
     200                 :            : /* assume x >= s */
     201                 :            : INLINE GEN
     202                 :  704436831 : subiuspec(GEN x, ulong s, long nx)
     203                 :            : {
     204                 :  704436831 :   GEN xd, zd = (GEN)avma;
     205                 :            :   long lz;
     206                 :            :   LOCAL_OVERFLOW;
     207                 :            : 
     208         [ +  + ]:  704436831 :   if (nx == 1) return utoi(x[0] - s);
     209                 :            : 
     210                 :   49985754 :   lz = nx+2; (void)new_chunk(lz);
     211                 :   49985754 :   xd = x + nx;
     212                 :   49985754 :   *--zd = subll(*--xd, s);
     213         [ +  + ]:   49985754 :   if (overflow)
     214                 :            :     for(;;)
     215                 :            :     {
     216                 :   24529017 :       *--zd = ((ulong)*--xd) - 1;
     217         [ +  + ]:   24529017 :       if (*xd) break;
     218                 :    1754175 :     }
     219         [ +  + ]:   49985754 :   if (xd == x)
     220         [ +  + ]:   35187930 :     while (*zd == 0) { zd++; lz--; } /* shorten z */
     221                 :            :   else
     222         [ +  + ]:  163881891 :     do  *--zd = *--xd; while (xd > x);
     223                 :   49985754 :   *--zd = evalsigne(1) | evallgefint(lz);
     224                 :   49985754 :   *--zd = evaltyp(t_INT) | evallg(lz);
     225                 :  704436831 :   avma=(pari_sp)zd; return zd;
     226                 :            : }
     227                 :            : 
     228                 :            : /* assume x > y */
     229                 :            : static GEN
     230                 : 1036302180 : subiispec(GEN x, GEN y, long nx, long ny)
     231                 :            : {
     232                 :            :   GEN xd,yd,zd;
     233                 : 1036302180 :   long lz, i = -2;
     234                 :            :   LOCAL_OVERFLOW;
     235                 :            : 
     236         [ +  + ]: 1036302180 :   if (ny==1) return subiuspec(x,*y,nx);
     237                 :  342180903 :   zd = (GEN)avma;
     238                 :  342180903 :   lz = nx+2; (void)new_chunk(lz);
     239                 :  342180903 :   xd = x + nx;
     240                 :  342180903 :   yd = y + ny;
     241                 :  342180903 :   zd[-1] = subll(xd[-1], yd[-1]);
     242                 :            : #ifdef subllx8
     243         [ +  + ]:  203058757 :   for (  ; i-8 > -ny; i-=8)
     244                 :   88998456 :     subllx8(xd+i, yd+i, zd+i, overflow);
     245                 :            : #endif
     246         [ +  + ]: 2827542444 :   for (  ; i >= -ny; i--) zd[i] = subllx(xd[i], yd[i]);
              [ +  +  + ]
     247         [ +  + ]:  342180903 :   if (overflow)
     248                 :            :     for(;;)
     249                 :            :     {
     250                 :   58761693 :       zd[i] = uel(xd,i) - 1;
     251         [ +  + ]:   58761693 :       if (xd[i--]) break;
     252                 :   28892913 :     }
     253         [ +  + ]:  342180903 :   if (i>=-nx)
     254         [ +  + ]:  227503691 :     for (; i >= -nx; i--) zd[i] = xd[i];
     255                 :            :   else
     256         [ +  + ]:  378011682 :     while (zd[i+1] == 0) { i++; lz--; } /* shorten z */
     257                 :  342180903 :   zd += i+1;
     258                 :  342180903 :   *--zd = evalsigne(1) | evallgefint(lz);
     259                 :  342180903 :   *--zd = evaltyp(t_INT) | evallg(lz);
     260                 : 1036302180 :   avma=(pari_sp)zd; return zd;
     261                 :            : }
     262                 :            : 
     263                 :            : static void
     264                 :   73849439 : roundr_up_ip(GEN x, long l)
     265                 :            : {
     266                 :   73849439 :   long i = l;
     267                 :            :   for(;;)
     268                 :            :   {
     269         [ +  + ]:   73937015 :     if (++uel(x,--i)) break;
     270         [ +  + ]:     110268 :     if (i == 2) { x[2] = (long)HIGHBIT; shiftr_inplace(x, 1); break; }
     271                 :      87576 :   }
     272                 :   73849439 : }
     273                 :            : 
     274                 :            : void
     275                 :   76691925 : affir(GEN x, GEN y)
     276                 :            : {
     277                 :   76691925 :   const long s = signe(x), ly = lg(y);
     278                 :            :   long lx, sh, i;
     279                 :            : 
     280         [ +  + ]:   76691925 :   if (!s)
     281                 :            :   {
     282                 :     651291 :     y[1] = evalexpo(-bit_accuracy(ly));
     283                 :     651291 :     return;
     284                 :            :   }
     285                 :            : 
     286 [ +  + ][ +  + ]:   76040634 :   lx = lgefint(x); sh = bfffo(x[2]);
         [ +  + ][ +  + ]
     287                 :   76040634 :   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
     288         [ +  + ]:   76040634 :   if (sh) {
     289         [ +  + ]:   75363570 :     if (lx <= ly)
     290                 :            :     {
     291         [ +  + ]:  153578850 :       for (i=lx; i<ly; i++) y[i]=0;
     292                 :   54368064 :       shift_left(y,x,2,lx-1, 0,sh);
     293                 :   54368064 :       return;
     294                 :            :     }
     295                 :   20995506 :     shift_left(y,x,2,ly-1, x[ly],sh);
     296                 :            :     /* lx > ly: round properly */
     297         [ +  + ]:   20995506 :     if ((x[ly]<<sh) & HIGHBIT) roundr_up_ip(y, ly);
     298                 :            :   }
     299                 :            :   else {
     300         [ +  + ]:     677064 :     if (lx <= ly)
     301                 :            :     {
     302         [ +  + ]:    1223091 :       for (i=2; i<lx; i++) y[i]=x[i];
     303         [ +  + ]:    1212183 :       for (   ; i<ly; i++) y[i]=0;
     304                 :     464010 :       return;
     305                 :            :     }
     306         [ +  + ]:     617508 :     for (i=2; i<ly; i++) y[i]=x[i];
     307                 :            :     /* lx > ly: round properly */
     308         [ +  + ]:   76691925 :     if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
     309                 :            :   }
     310                 :            : }
     311                 :            : 
     312                 :            : INLINE GEN
     313                 :  213946032 : shiftispec(GEN x, long nx, long n)
     314                 :            : {
     315                 :            :   long ny, i, m;
     316                 :            :   GEN y, yd;
     317         [ +  + ]:  213946032 :   if (!n)  return icopyspec(x, nx);
     318                 :            : 
     319         [ +  + ]:  203676012 :   if (n > 0)
     320                 :            :   {
     321                 :  147166335 :     GEN z = (GEN)avma;
     322                 :  147166335 :     long d = dvmdsBIL(n, &m);
     323                 :            : 
     324                 :  147166335 :     ny = nx+d; y = new_chunk(ny + 2); yd = y + 2;
     325         [ +  + ]:  303899658 :     for ( ; d; d--) *--z = 0;
     326 [ +  + ][ +  + ]:  164985687 :     if (!m) for (i=0; i<nx; i++) yd[i]=x[i];
     327                 :            :     else
     328                 :            :     {
     329                 :  146422776 :       register const ulong sh = BITS_IN_LONG - m;
     330                 :  146422776 :       shift_left(yd,x, 0,nx-1, 0,m);
     331                 :  146422776 :       i = uel(x,0) >> sh;
     332                 :            :       /* Extend y on the left? */
     333         [ +  + ]:  146422776 :       if (i) { ny++; y = new_chunk(1); y[2] = i; }
     334                 :            :     }
     335                 :            :   }
     336                 :            :   else
     337                 :            :   {
     338                 :   56509677 :     ny = nx - dvmdsBIL(-n, &m);
     339         [ +  + ]:   56509677 :     if (ny<1) return gen_0;
     340                 :   56362005 :     y = new_chunk(ny + 2); yd = y + 2;
     341         [ +  + ]:   56362005 :     if (m) {
     342                 :   51712404 :       shift_right(yd,x, 0,ny, 0,m);
     343         [ +  + ]:   51712404 :       if (yd[0] == 0)
     344                 :            :       {
     345         [ +  + ]:    4198977 :         if (ny==1) { avma = (pari_sp)(y+3); return gen_0; }
     346                 :    3687852 :         ny--; avma = (pari_sp)(++y);
     347                 :            :       }
     348                 :            :     } else {
     349         [ +  + ]:  193586676 :       for (i=0; i<ny; i++) yd[i]=x[i];
     350                 :            :     }
     351                 :            :   }
     352                 :  203017215 :   y[1] = evalsigne(1)|evallgefint(ny + 2);
     353                 :  213946032 :   y[0] = evaltyp(t_INT)|evallg(ny + 2); return y;
     354                 :            : }
     355                 :            : 
     356                 :            : GEN
     357                 :    7778949 : mantissa2nr(GEN x, long n)
     358                 :            : { /*This is a kludge since x is not an integer*/
     359                 :    7778949 :   long s = signe(x);
     360                 :            :   GEN y;
     361                 :            : 
     362         [ +  + ]:    7778949 :   if(s == 0) return gen_0;
     363                 :    7778079 :   y = shiftispec(x + 2, lg(x) - 2, n);
     364         [ +  - ]:    7778079 :   if (signe(y)) setsigne(y, s);
     365                 :    7778949 :   return y;
     366                 :            : }
     367                 :            : 
     368                 :            : GEN
     369                 :     338802 : truncr(GEN x)
     370                 :            : {
     371                 :            :   long d,e,i,s,m;
     372                 :            :   GEN y;
     373                 :            : 
     374 [ +  + ][ +  + ]:     338802 :   if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
     375                 :     156438 :   d = nbits2prec(e+1); m = remsBIL(e);
     376         [ -  + ]:     156438 :   if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
     377                 :            : 
     378                 :     156438 :   y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
     379         [ +  + ]:     156438 :   if (++m == BITS_IN_LONG)
     380         [ +  + ]:          6 :     for (i=2; i<d; i++) y[i]=x[i];
     381                 :            :   else
     382                 :     156435 :     shift_right(y,x, 2,d,0, BITS_IN_LONG - m);
     383                 :     338802 :   return y;
     384                 :            : }
     385                 :            : 
     386                 :            : /* integral part */
     387                 :            : GEN
     388                 :     348228 : floorr(GEN x)
     389                 :            : {
     390                 :            :   long d,e,i,lx,m;
     391                 :            :   GEN y;
     392                 :            : 
     393         [ +  + ]:     348228 :   if (signe(x) >= 0) return truncr(x);
     394         [ +  + ]:     120621 :   if ((e=expo(x)) < 0) return gen_m1;
     395                 :      47262 :   d = nbits2prec(e+1); m = remsBIL(e);
     396         [ -  + ]:      47262 :   lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
     397                 :      47262 :   y = new_chunk(d);
     398         [ -  + ]:      47262 :   if (++m == BITS_IN_LONG)
     399                 :            :   {
     400         [ #  # ]:          0 :     for (i=2; i<d; i++) y[i]=x[i];
     401 [ #  # ][ #  # ]:          0 :     i=d; while (i<lx && !x[i]) i++;
     402         [ #  # ]:          0 :     if (i==lx) goto END;
     403                 :            :   }
     404                 :            :   else
     405                 :            :   {
     406                 :      47262 :     shift_right(y,x, 2,d,0, BITS_IN_LONG - m);
     407         [ +  + ]:      47262 :     if (x[d-1]<<m == 0)
     408                 :            :     {
     409 [ +  + ][ +  + ]:      65949 :       i=d; while (i<lx && !x[i]) i++;
     410         [ +  + ]:      15576 :       if (i==lx) goto END;
     411                 :            :     }
     412                 :            :   }
     413                 :            :   /* set y:=y+1 */
     414 [ +  - ][ +  - ]:      36429 :   for (i=d-1; i>=2; i--) { y[i]++; if (y[i]) goto END; }
     415                 :          0 :   y=new_chunk(1); y[2]=1; d++;
     416                 :            : END:
     417                 :      47262 :   y[1] = evalsigne(-1) | evallgefint(d);
     418                 :     348228 :   y[0] = evaltyp(t_INT) | evallg(d); return y;
     419                 :            : }
     420                 :            : 
     421                 :            : INLINE int
     422                 : 1197922557 : cmpiispec(GEN x, GEN y, long lx, long ly)
     423                 :            : {
     424                 :            :   long i;
     425         [ +  + ]: 1197922557 :   if (lx < ly) return -1;
     426         [ +  + ]: 1143931560 :   if (lx > ly) return  1;
     427 [ +  + ][ +  + ]: 1184909871 :   i = 0; while (i<lx && x[i]==y[i]) i++;
     428         [ +  + ]: 1091844495 :   if (i==lx) return 0;
     429         [ +  + ]: 1197922557 :   return (uel(x,i) > uel(y,i))? 1: -1;
     430                 :            : }
     431                 :            : 
     432                 :            : INLINE int
     433                 :   17662065 : equaliispec(GEN x, GEN y, long lx, long ly)
     434                 :            : {
     435                 :            :   long i;
     436         [ +  + ]:   17662065 :   if (lx != ly) return 0;
     437 [ +  + ][ +  + ]:   28541286 :   i = ly-1; while (i>=0 && x[i]==y[i]) i--;
     438                 :   17662065 :   return i < 0;
     439                 :            : }
     440                 :            : 
     441                 :            : /***********************************************************************/
     442                 :            : /**                                                                   **/
     443                 :            : /**                          MULTIPLICATION                           **/
     444                 :            : /**                                                                   **/
     445                 :            : /***********************************************************************/
     446                 :            : /* assume ny > 0 */
     447                 :            : INLINE GEN
     448                 :  888114567 : muluispec(ulong x, GEN y, long ny)
     449                 :            : {
     450                 :  888114567 :   GEN yd, z = (GEN)avma;
     451                 :  888114567 :   long lz = ny+3;
     452                 :            :   LOCAL_HIREMAINDER;
     453                 :            : 
     454                 :  888114567 :   (void)new_chunk(lz);
     455                 :  888114567 :   yd = y + ny; *--z = mulll(x, *--yd);
     456      [ +  +  + ]: 1608297096 :   while (yd > y) *--z = addmul(x,*--yd);
     457         [ +  + ]:  888114567 :   if (hiremainder) *--z = hiremainder; else lz--;
     458                 :  888114567 :   *--z = evalsigne(1) | evallgefint(lz);
     459                 :  888114567 :   *--z = evaltyp(t_INT) | evallg(lz);
     460                 :  888114567 :   avma=(pari_sp)z; return z;
     461                 :            : }
     462                 :            : 
     463                 :            : /* a + b*|Y| */
     464                 :            : GEN
     465                 :     542115 : addumului(ulong a, ulong b, GEN Y)
     466                 :            : {
     467                 :            :   GEN yd,y,z;
     468                 :            :   long ny,lz;
     469                 :            :   LOCAL_HIREMAINDER;
     470                 :            :   LOCAL_OVERFLOW;
     471                 :            : 
     472         [ +  + ]:     542115 :   if (!signe(Y)) return utoi(a);
     473                 :            : 
     474                 :     541623 :   y = LIMBS(Y); z = (GEN)avma;
     475                 :     541623 :   ny = NLIMBS(Y);
     476                 :     541623 :   lz = ny+3;
     477                 :            : 
     478                 :     541623 :   (void)new_chunk(lz);
     479                 :     541623 :   yd = y + ny; *--z = addll(a, mulll(b, *--yd));
     480         [ +  + ]:     541623 :   if (overflow) hiremainder++; /* can't overflow */
     481      [ +  +  + ]:    1153965 :   while (yd > y) *--z = addmul(b,*--yd);
     482         [ +  + ]:     541623 :   if (hiremainder) *--z = hiremainder; else lz--;
     483                 :     541623 :   *--z = evalsigne(1) | evallgefint(lz);
     484                 :     541623 :   *--z = evaltyp(t_INT) | evallg(lz);
     485                 :     542115 :   avma=(pari_sp)z; return z;
     486                 :            : }
     487                 :            : 
     488                 :            : /***********************************************************************/
     489                 :            : /**                                                                   **/
     490                 :            : /**                          DIVISION                                 **/
     491                 :            : /**                                                                   **/
     492                 :            : /***********************************************************************/
     493                 :            : 
     494                 :            : ulong
     495                 :  548013381 : umodiu(GEN y, ulong x)
     496                 :            : {
     497                 :  548013381 :   long sy=signe(y),ly,i;
     498                 :            :   ulong xi;
     499                 :            :   LOCAL_HIREMAINDER;
     500                 :            : 
     501         [ -  + ]:  548013381 :   if (!x) pari_err_INV("umodiu",gen_0);
     502         [ +  + ]:  548013381 :   if (!sy) return 0;
     503                 :  503277912 :   ly = lgefint(y);
     504         [ +  + ]:  503277912 :   if (x <= uel(y,2))
     505                 :            :   {
     506                 :  119932935 :     hiremainder=0;
     507         [ +  + ]:  119932935 :     if (ly==3)
     508                 :            :     {
     509                 :  110772267 :       hiremainder=uel(y,2)%x;
     510         [ +  + ]:  110772267 :       if (!hiremainder) return 0;
     511         [ +  + ]:   93673584 :       return (sy > 0)? hiremainder: x - hiremainder;
     512                 :            :     }
     513                 :            :   }
     514                 :            :   else
     515                 :            :   {
     516 [ +  + ][ +  + ]:  383344977 :     if (ly==3) return (sy > 0)? uel(y,2): x - uel(y,2);
     517                 :   98944428 :     hiremainder=uel(y,2); ly--; y++;
     518                 :            :   }
     519                 :  108105096 :   xi = get_Fl_red(x);
     520 [ +  + ][ +  + ]:  354912456 :   for (i=2; i<ly; i++) (void)divll_pre(y[i],x,xi);
            [ +  + ][ + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     521         [ +  + ]:  108105096 :   if (!hiremainder) return 0;
     522         [ +  + ]:  548013381 :   return (sy > 0)? hiremainder: x - hiremainder;
     523                 :            : }
     524                 :            : 
     525                 :            : /* return |y| \/ x */
     526                 :            : GEN
     527                 :   76420104 : diviu_rem(GEN y, ulong x, ulong *rem)
     528                 :            : {
     529                 :            :   long ly,i;
     530                 :            :   GEN z;
     531                 :            :   ulong xi;
     532                 :            :   LOCAL_HIREMAINDER;
     533                 :            : 
     534         [ -  + ]:   76420104 :   if (!x) pari_err_INV("diviu_rem",gen_0);
     535         [ +  + ]:   76420104 :   if (!signe(y)) { *rem = 0; return gen_0; }
     536                 :            : 
     537                 :   76023789 :   ly = lgefint(y);
     538         [ +  + ]:   76023789 :   if (x <= uel(y,2))
     539                 :            :   {
     540                 :   63790440 :     hiremainder=0;
     541         [ +  + ]:   63790440 :     if (ly==3)
     542                 :            :     {
     543                 :   16023069 :       z = cgetipos(3);
     544 [ +  - ][ #  # ]:   16023069 :       z[2] = divll(uel(y,2),x);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     545                 :   16023069 :       *rem = hiremainder; return z;
     546                 :            :     }
     547                 :            :   }
     548                 :            :   else
     549                 :            :   {
     550         [ +  + ]:   12233349 :     if (ly==3) { *rem = uel(y,2); return gen_0; }
     551                 :   10085574 :     hiremainder = uel(y,2); ly--; y++;
     552                 :            :   }
     553                 :   57852945 :   xi = get_Fl_red(x);
     554                 :   57852945 :   z = cgetipos(ly);
     555 [ +  + ][ +  + ]:  380977722 :   for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x,xi);
            [ +  + ][ + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     556                 :   76420104 :   *rem = hiremainder; return z;
     557                 :            : }
     558                 :            : 
     559                 :            : GEN
     560                 :   27271746 : divis_rem(GEN y, long x, long *rem)
     561                 :            : {
     562                 :   27271746 :   long sy=signe(y),ly,s,i;
     563                 :            :   GEN z;
     564                 :            :   ulong xi;
     565                 :            :   LOCAL_HIREMAINDER;
     566                 :            : 
     567         [ -  + ]:   27271746 :   if (!x) pari_err_INV("divis_rem",gen_0);
     568         [ +  + ]:   27271746 :   if (!sy) { *rem=0; return gen_0; }
     569         [ +  + ]:   22483737 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     570                 :            : 
     571                 :   22483737 :   ly = lgefint(y);
     572         [ +  + ]:   22483737 :   if ((ulong)x <= uel(y,2))
     573                 :            :   {
     574                 :   18370935 :     hiremainder=0;
     575         [ +  + ]:   18370935 :     if (ly==3)
     576                 :            :     {
     577                 :   18275169 :       z = cgeti(3); z[1] = evallgefint(3) | evalsigne(s);
     578 [ +  - ][ #  # ]:   18275169 :       z[2] = divll(uel(y,2),x);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     579 [ +  + ][ +  + ]:   18275169 :       if (sy<0) hiremainder = - ((long)hiremainder);
     580                 :   18275169 :       *rem = (long)hiremainder; return z;
     581                 :            :     }
     582                 :            :   }
     583                 :            :   else
     584                 :            :   {
     585         [ +  + ]:    4112802 :     if (ly==3) { *rem = itos(y); return gen_0; }
     586                 :      41373 :     hiremainder = uel(y,2); ly--; y++;
     587                 :            :   }
     588                 :     137139 :   xi = get_Fl_red(x);
     589                 :     137139 :   z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
     590 [ +  + ][ +  + ]:     510258 :   for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x,xi);
            [ +  + ][ + ]
         [ +  + ][ +  - ]
         [ +  + ][ +  + ]
     591         [ +  + ]:     137139 :   if (sy<0) hiremainder = - ((long)hiremainder);
     592                 :   27271746 :   *rem = (long)hiremainder; return z;
     593                 :            : }
     594                 :            : 
     595                 :            : GEN
     596                 :      18168 : divis(GEN y, long x)
     597                 :            : {
     598                 :      18168 :   long sy=signe(y),ly,s,i;
     599                 :            :   ulong xi;
     600                 :            :   GEN z;
     601                 :            :   LOCAL_HIREMAINDER;
     602                 :            : 
     603         [ -  + ]:      18168 :   if (!x) pari_err_INV("divis",gen_0);
     604         [ -  + ]:      18168 :   if (!sy) return gen_0;
     605         [ -  + ]:      18168 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     606                 :            : 
     607                 :      18168 :   ly = lgefint(y);
     608         [ +  + ]:      18168 :   if ((ulong)x <= uel(y,2))
     609                 :            :   {
     610                 :      15015 :     hiremainder=0;
     611         [ +  + ]:      15015 :     if (ly==3)
     612                 :            :     {
     613                 :      14874 :       z = cgeti(3); z[1] = evallgefint(3) | evalsigne(s);
     614 [ +  - ][ #  # ]:      14874 :       z[2] = divll(y[2],x);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     615                 :      14874 :       return z;
     616                 :            :     }
     617                 :            :   }
     618                 :            :   else
     619                 :            :   {
     620         [ -  + ]:       3153 :     if (ly==3) return gen_0;
     621                 :       3153 :     hiremainder=y[2]; ly--; y++;
     622                 :            :   }
     623                 :       3294 :   xi = get_Fl_red(x);
     624                 :       3294 :   z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
     625 [ +  + ][ +  + ]:       9471 :   for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x, xi);
            [ -  + ][ + ]
         [ +  + ][ +  - ]
         [ +  + ][ +  + ]
     626                 :      18168 :   return z;
     627                 :            : }
     628                 :            : 
     629                 :            : GEN
     630                 :   58650198 : divrr(GEN x, GEN y)
     631                 :            : {
     632                 :   58650198 :   long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
     633                 :            :   ulong y0,y1;
     634                 :            :   GEN r, r1;
     635                 :            : 
     636         [ -  + ]:   58650198 :   if (!x) pari_err_INV("divrr",y);
     637                 :   58650198 :   e = expo(x) - expo(y);
     638         [ +  + ]:   58650198 :   if (!sx) return real_0_bit(e);
     639         [ +  + ]:   57367182 :   if (sy<0) sx = -sx;
     640                 :            : 
     641                 :   57367182 :   lx=lg(x); ly=lg(y);
     642         [ +  + ]:   57367182 :   if (ly==3)
     643                 :            :   {
     644         [ +  + ]:   46882803 :     ulong k = x[2], l = (lx>3)? x[3]: 0;
     645                 :            :     LOCAL_HIREMAINDER;
     646         [ +  + ]:   46882803 :     if (k < uel(y,2)) e--;
     647                 :            :     else
     648                 :            :     {
     649         [ +  + ]:   21723717 :       l >>= 1; if (k&1) l |= HIGHBIT;
     650                 :   21723717 :       k >>= 1;
     651                 :            :     }
     652 [ -  + ][ -  + ]:   46882803 :     hiremainder = k; k = divll(l,y[2]);
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
                 [ +  + ]
     653 [ +  + ][ +  + ]:   46882803 :     if (hiremainder & HIGHBIT)
     654                 :            :     {
     655                 :   12968940 :       k++;
     656         [ -  + ]:   12968940 :       if (!k) { k = HIGHBIT; e++; }
     657                 :            :     }
     658                 :   46882803 :     r = cgetr(3);
     659                 :   46882803 :     r[1] = evalsigne(sx) | evalexpo(e);
     660                 :   46882803 :     r[2] = k; return r;
     661                 :            :   }
     662                 :            : 
     663                 :   10484379 :   lr = minss(lx,ly); r = new_chunk(lr);
     664                 :   10484379 :   r1 = r-1;
     665         [ +  + ]:   81149655 :   r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
     666         [ +  + ]:   10484379 :   r1[lr] = (lx>ly)? x[lr]: 0;
     667                 :   10484379 :   y0 = y[2]; y1 = y[3];
     668         [ +  + ]:   91634034 :   for (i=0; i<lr-1; i++)
     669                 :            :   { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
     670                 :            :     ulong k, qp;
     671                 :            :     LOCAL_HIREMAINDER;
     672                 :            :     LOCAL_OVERFLOW;
     673                 :            : 
     674         [ +  + ]:   81149655 :     if (uel(r1,1) == y0)
     675                 :            :     {
     676                 :     246117 :       qp = ULONG_MAX; k = addll(y0,r1[2]);
     677                 :            :     }
     678                 :            :     else
     679                 :            :     {
     680         [ -  + ]:   80903538 :       if (uel(r1,1) > y0) /* can't happen if i=0 */
     681                 :            :       {
     682                 :          0 :         GEN y1 = y+1;
     683                 :          0 :         j = lr-i; r1[j] = subll(r1[j],y1[j]);
     684 [ #  # ][ #  # ]:          0 :         for (j--; j>0; j--) r1[j] = subllx(r1[j],y1[j]);
     685 [ #  # ][ #  # ]:          0 :         j=i; do uel(r,--j)++; while (j && !uel(r,j));
     686                 :            :       }
     687                 :   80903538 :       hiremainder = r1[1]; overflow = 0;
     688 [ +  + ][ -  + ]:   80903538 :       qp = divll(r1[2],y0); k = hiremainder;
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
                 [ +  + ]
     689                 :            :     }
     690                 :   81149655 :     j = lr-i+1;
     691         [ +  + ]:   81149655 :     if (!overflow)
     692                 :            :     {
     693                 :            :       long k3, k4;
     694                 :   80986503 :       k3 = mulll(qp,y1);
     695 [ +  + ][ +  + ]:   80986503 :       if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
     696                 :   10474239 :         k4 = subll(hiremainder,k);
     697                 :            :       else
     698                 :            :       {
     699                 :   70512264 :         k3 = subll(k3, r1[3]);
     700                 :   70512264 :         k4 = subllx(hiremainder,k);
     701                 :            :       }
     702 [ +  + ][ +  + ]:   98786454 :       while (!overflow && k4) { qp--; k3 = subll(k3,y1); k4 = subllx(k4,y0); }
         [ +  + ][ +  + ]
     703                 :            :     }
     704         [ +  + ]:   81149655 :     if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
     705 [ +  + ][ +  + ]:  578356305 :     for (j--; j>1; j--)
     706                 :            :     {
     707                 :  497206650 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     708                 :  497206650 :       hiremainder += overflow;
     709                 :            :     }
     710         [ +  + ]:   81149655 :     if (uel(r1,1) != hiremainder)
     711                 :            :     {
     712         [ +  - ]:     150165 :       if (uel(r1,1) < hiremainder)
     713                 :            :       {
     714                 :     150165 :         qp--;
     715                 :     150165 :         j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
     716 [ +  + ][ +  + ]:     985431 :         for (j--; j>1; j--) r1[j] = addllx(r1[j], y[j]);
     717                 :            :       }
     718                 :            :       else
     719                 :            :       {
     720                 :          0 :         r1[1] -= hiremainder;
     721         [ #  # ]:          0 :         while (r1[1])
     722                 :            :         {
     723 [ #  # ][ #  # ]:          0 :           qp++; if (!qp) { j=i; do ((ulong*)r)[--j]++; while (j && !r[j]); }
                 [ #  # ]
     724                 :          0 :           j = lr-i-(lr-i>=ly); r1[j] = subll(r1[j],y[j]);
     725 [ #  # ][ #  # ]:          0 :           for (j--; j>1; j--) r1[j] = subllx(r1[j],y[j]);
     726                 :          0 :           r1[1] -= overflow;
     727                 :            :         }
     728                 :            :       }
     729                 :            :     }
     730                 :   81149655 :     *++r1 = qp;
     731                 :            :   }
     732                 :            :   /* i = lr-1 */
     733                 :            :   /* round correctly */
     734         [ +  + ]:   10484379 :   if (uel(r1,1) > (y0>>1))
     735                 :            :   {
     736 [ +  + ][ +  + ]:    5098718 :     j=i; do uel(r,--j)++; while (j && !r[j]);
     737                 :            :   }
     738         [ +  + ]:   81149655 :   r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
     739         [ +  + ]:   10484379 :   if (r[0] == 0) e--;
     740         [ +  + ]:    3855912 :   else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
     741                 :            :   else { /* possible only when rounding up to 0x2 0x0 ... */
     742                 :          3 :     r[2] = (long)HIGHBIT; e++;
     743                 :            :   }
     744                 :   10484379 :   r[0] = evaltyp(t_REAL)|evallg(lr);
     745                 :   10484379 :   r[1] = evalsigne(sx) | evalexpo(e);
     746                 :   58650198 :   return r;
     747                 :            : }
     748                 :            : 
     749                 :            : GEN
     750                 :   10771665 : divri(GEN x, GEN y)
     751                 :            : {
     752                 :   10771665 :   long lx, s = signe(y);
     753                 :            :   pari_sp av;
     754                 :            :   GEN z;
     755                 :            : 
     756         [ -  + ]:   10771665 :   if (!s) pari_err_INV("divri",y);
     757         [ +  + ]:   10771665 :   if (!signe(x)) return real_0_bit(expo(x) - expi(y));
     758         [ +  + ]:   10764327 :   if (!is_bigint(y)) {
     759                 :    9089445 :     GEN z = divru(x, y[2]);
     760         [ +  + ]:    9089445 :     if (s < 0) togglesign(z);
     761                 :    9089445 :     return z;
     762                 :            :   }
     763                 :    1674882 :   lx = lg(x); z = cgetr(lx); av = avma;
     764                 :    1674882 :   affrr(divrr(x, itor(y, lx+1)), z);
     765                 :   10771665 :   avma = av; return z;
     766                 :            : }
     767                 :            : 
     768                 :            : /* Integer division x / y: such that sign(r) = sign(x)
     769                 :            :  *   if z = ONLY_REM return remainder, otherwise return quotient
     770                 :            :  *   if z != NULL set *z to remainder
     771                 :            :  *   *z is the last object on stack (and thus can be disposed of with cgiv
     772                 :            :  *   instead of gerepile)
     773                 :            :  * If *z is zero, we put gen_0 here and no copy.
     774                 :            :  * space needed: lx + ly */
     775                 :            : GEN
     776                 :  373396857 : dvmdii(GEN x, GEN y, GEN *z)
     777                 :            : {
     778                 :  373396857 :   long sx = signe(x), sy = signe(y);
     779                 :  373396857 :   long lx, ly = lgefint(y), lz, i, j, sh, lq, lr;
     780                 :            :   pari_sp av;
     781                 :            :   ulong y0,y0i,y1, *xd,*rd,*qd;
     782                 :            :   GEN q, r, r1;
     783                 :            : 
     784         [ +  + ]:  373396857 :   if (!sx)
     785                 :            :   {
     786         [ +  + ]:   10440207 :     if (ly < 3) pari_err_INV("dvmdii",gen_0);
     787 [ +  + ][ +  + ]:   10440192 :     if (!z || z == ONLY_REM) return gen_0;
     788                 :    4431042 :     *z=gen_0; return gen_0;
     789                 :            :   }
     790         [ +  + ]:  362956650 :   if (ly <= 3)
     791                 :            :   {
     792                 :            :     ulong rem;
     793         [ -  + ]:  228196038 :     if (ly < 3) pari_err_INV("dvmdii",gen_0);
     794         [ +  + ]:  228196038 :     if (z == ONLY_REM)
     795                 :            :     {
     796                 :  208672011 :       rem = umodiu(x,uel(y,2));
     797         [ +  + ]:  208672011 :       if (!rem) return gen_0;
     798         [ +  + ]:  193770705 :       return (sx < 0)? utoineg(uel(y,2) - rem): utoipos(rem);
     799                 :            :     }
     800                 :   19524027 :     q = diviu_rem(x, uel(y,2), &rem);
     801         [ +  + ]:   19524027 :     if (sx != sy) togglesign(q);
     802         [ +  + ]:   19524027 :     if (!z) return q;
     803         [ +  + ]:   19142892 :     if (!rem) *z = gen_0;
     804         [ +  + ]:    8215584 :     else *z = sx < 0? utoineg(rem): utoipos(rem);
     805                 :  228196038 :     return q;
     806                 :            :   }
     807                 :  134760612 :   lx=lgefint(x);
     808                 :  134760612 :   lz=lx-ly;
     809         [ +  + ]:  134760612 :   if (lz <= 0)
     810                 :            :   {
     811         [ +  + ]:   54527700 :     if (lz == 0)
     812                 :            :     {
     813         [ +  + ]:   45829083 :       for (i=2; i<lx; i++)
     814         [ +  + ]:   45693288 :         if (x[i] != y[i])
     815                 :            :         {
     816         [ +  + ]:   44742303 :           if (uel(x,i) > uel(y,i)) goto DIVIDE;
     817                 :   11561661 :           goto TRIVIAL;
     818                 :            :         }
     819         [ +  + ]:     135795 :       if (z == ONLY_REM) return gen_0;
     820         [ +  + ]:       5526 :       if (z) *z = gen_0;
     821         [ +  + ]:       5526 :       if (sx < 0) sy = -sy;
     822                 :       5526 :       return stoi(sy);
     823                 :            :     }
     824                 :            : TRIVIAL:
     825         [ +  + ]:   21211263 :     if (z == ONLY_REM) return icopy(x);
     826         [ +  + ]:     772971 :     if (z) *z = icopy(x);
     827                 :     772971 :     return gen_0;
     828                 :            :   }
     829                 :            : DIVIDE: /* quotient is non-zero */
     830         [ +  + ]:  113413554 :   av=avma; if (sx<0) sy = -sy;
     831 [ +  + ][ +  + ]:  113413554 :   r1 = new_chunk(lx); sh = bfffo(y[2]);
         [ +  + ][ +  + ]
     832 [ +  + ][ +  + ]:  113413554 :   if (sh)
     833                 :            :   { /* normalize so that highbit(y) = 1 (shift left x and y by sh bits)*/
     834                 :  106391448 :     register const ulong m = BITS_IN_LONG - sh;
     835                 :  106391448 :     r = new_chunk(ly);
     836                 :  106391448 :     shift_left(r, y,2,ly-1, 0,sh); y = r;
     837                 :  106391448 :     shift_left(r1,x,2,lx-1, 0,sh);
     838                 :  106391448 :     r1[1] = uel(x,2) >> m;
     839                 :            :   }
     840                 :            :   else
     841                 :            :   {
     842         [ +  + ]:   49723797 :     r1[1] = 0; for (j=2; j<lx; j++) r1[j] = x[j];
     843                 :            :   }
     844                 :  113413554 :   x = r1;
     845                 :  113413554 :   y0 = y[2]; y0i = get_Fl_red(y0);
     846                 :  113413554 :   y1 = y[3];
     847         [ +  + ]:  492847452 :   for (i=0; i<=lz; i++)
     848                 :            :   { /* r1 = x + i */
     849                 :            :     ulong k, qp;
     850                 :            :     LOCAL_HIREMAINDER;
     851                 :            :     LOCAL_OVERFLOW;
     852                 :            : 
     853         [ +  + ]:  379433898 :     if (uel(r1,1) == y0)
     854                 :            :     {
     855                 :      17880 :       qp = ULONG_MAX; k = addll(y0,r1[2]);
     856                 :            :     }
     857                 :            :     else
     858                 :            :     {
     859                 :  379416018 :       hiremainder = r1[1]; overflow = 0;
     860 [ +  + ][ +  - ]:  379416018 :       qp = divll_pre(r1[2],y0,y0i); k = hiremainder;
         [ +  - ][ +  - ]
                 [ -  + ]
     861                 :            :     }
     862         [ +  + ]:  379433898 :     if (!overflow)
     863                 :            :     {
     864                 :  379418685 :       long k3 = subll(mulll(qp,y1), r1[3]);
     865                 :  379418685 :       long k4 = subllx(hiremainder,k);
     866 [ +  + ][ +  + ]:  449361397 :       while (!overflow && k4) { qp--; k3 = subll(k3,y1); k4 = subllx(k4,y0); }
         [ +  + ][ +  + ]
     867                 :            :     }
     868                 :  379433898 :     hiremainder = 0; j = ly;
     869 [ +  + ][ +  + ]: 4437697824 :     for (j--; j>1; j--)
     870                 :            :     {
     871                 : 4058263926 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     872                 : 4058263926 :       hiremainder += overflow;
     873                 :            :     }
     874         [ +  + ]:  379433898 :     if (uel(r1,1) < hiremainder)
     875                 :            :     {
     876                 :     104866 :       qp--;
     877                 :     104866 :       j = ly-1; r1[j] = addll(r1[j],y[j]);
     878 [ +  + ][ +  + ]:     879201 :       for (j--; j>1; j--) r1[j] = addllx(r1[j],y[j]);
     879                 :            :     }
     880                 :  379433898 :     *++r1 = qp;
     881                 :            :   }
     882                 :            : 
     883                 :  113413554 :   lq = lz+2;
     884         [ +  + ]:  113413554 :   if (!z)
     885                 :            :   {
     886                 :     410967 :     qd = (ulong*)av;
     887                 :     410967 :     xd = (ulong*)(x + lq);
     888         [ +  + ]:     410967 :     if (x[1]) { lz++; lq++; }
     889         [ +  + ]:     821967 :     while (lz--) *--qd = *--xd;
     890                 :     410967 :     *--qd = evalsigne(sy) | evallgefint(lq);
     891                 :     410967 :     *--qd = evaltyp(t_INT) | evallg(lq);
     892                 :     410967 :     avma = (pari_sp)qd; return (GEN)qd;
     893                 :            :   }
     894                 :            : 
     895 [ +  + ][ +  + ]:  123826818 :   j=lq; while (j<lx && !x[j]) j++;
     896                 :  113002587 :   lz = lx-j;
     897         [ +  + ]:  113002587 :   if (z == ONLY_REM)
     898                 :            :   {
     899         [ +  + ]:   86654487 :     if (lz==0) { avma = av; return gen_0; }
     900                 :   85252938 :     rd = (ulong*)av; lr = lz+2;
     901                 :   85252938 :     xd = (ulong*)(x + lx);
     902 [ +  + ][ +  + ]:  104488641 :     if (!sh) while (lz--) *--rd = *--xd;
     903                 :            :     else
     904                 :            :     { /* shift remainder right by sh bits */
     905                 :   78260727 :       const ulong shl = BITS_IN_LONG - sh;
     906                 :            :       ulong l;
     907                 :   78260727 :       xd--;
     908         [ +  + ]:  291561036 :       while (--lz) /* fill r[3..] */
     909                 :            :       {
     910                 :  213300309 :         l = *xd >> sh;
     911                 :  213300309 :         *--rd = l | (*--xd << shl);
     912                 :            :       }
     913                 :   78260727 :       l = *xd >> sh;
     914         [ +  + ]:   78260727 :       if (l) *--rd = l; else lr--;
     915                 :            :     }
     916                 :   85252938 :     *--rd = evalsigne(sx) | evallgefint(lr);
     917                 :   85252938 :     *--rd = evaltyp(t_INT) | evallg(lr);
     918                 :   85252938 :     avma = (pari_sp)rd; return (GEN)rd;
     919                 :            :   }
     920                 :            : 
     921                 :   26348100 :   lr = lz+2;
     922                 :   26348100 :   rd = NULL; /* gcc -Wall */
     923         [ +  + ]:   26348100 :   if (lz)
     924                 :            :   { /* non zero remainder: initialize rd */
     925                 :   25540143 :     xd = (ulong*)(x + lx);
     926         [ +  + ]:   25540143 :     if (!sh)
     927                 :            :     {
     928                 :       8439 :       rd = (ulong*)avma; (void)new_chunk(lr);
     929         [ +  + ]:      68262 :       while (lz--) *--rd = *--xd;
     930                 :            :     }
     931                 :            :     else
     932                 :            :     { /* shift remainder right by sh bits */
     933                 :   25531704 :       const ulong shl = BITS_IN_LONG - sh;
     934                 :            :       ulong l;
     935                 :   25531704 :       rd = (ulong*)x; /* overwrite shifted y */
     936                 :   25531704 :       xd--;
     937         [ +  + ]:   89385912 :       while (--lz)
     938                 :            :       {
     939                 :   63854208 :         l = *xd >> sh;
     940                 :   63854208 :         *--rd = l | (*--xd << shl);
     941                 :            :       }
     942                 :   25531704 :       l = *xd >> sh;
     943         [ +  + ]:   25531704 :       if (l) *--rd = l; else lr--;
     944                 :            :     }
     945                 :   25540143 :     *--rd = evalsigne(sx) | evallgefint(lr);
     946                 :   25540143 :     *--rd = evaltyp(t_INT) | evallg(lr);
     947                 :   25540143 :     rd += lr;
     948                 :            :   }
     949                 :   26348100 :   qd = (ulong*)av;
     950                 :   26348100 :   xd = (ulong*)(x + lq);
     951         [ +  + ]:   26348100 :   if (x[1]) lq++;
     952         [ +  + ]:   82328700 :   j = lq-2; while (j--) *--qd = *--xd;
     953                 :   26348100 :   *--qd = evalsigne(sy) | evallgefint(lq);
     954                 :   26348100 :   *--qd = evaltyp(t_INT) | evallg(lq);
     955                 :   26348100 :   q = (GEN)qd;
     956         [ +  + ]:   26348100 :   if (lr==2) *z = gen_0;
     957                 :            :   else
     958                 :            :   { /* rd has been properly initialized: we had lz > 0 */
     959         [ +  + ]:  148425107 :     while (lr--) *--qd = *--rd;
     960                 :   25540143 :     *z = (GEN)qd;
     961                 :            :   }
     962                 :  373396842 :   avma = (pari_sp)qd; return q;
     963                 :            : }
     964                 :            : 
     965                 :            : /* Montgomery reduction.
     966                 :            :  * N has k words, assume T >= 0 has less than 2k.
     967                 :            :  * Return res := T / B^k mod N, where B = 2^BIL
     968                 :            :  * such that 0 <= res < T/B^k + N  and  res has less than k words */
     969                 :            : GEN
     970                 :    6207447 : red_montgomery(GEN T, GEN N, ulong inv)
     971                 :            : {
     972                 :            :   pari_sp av;
     973                 :            :   GEN Te, Td, Ne, Nd, scratch;
     974                 :    6207447 :   ulong i, j, m, t, d, k = NLIMBS(N);
     975                 :            :   int carry;
     976                 :            :   LOCAL_HIREMAINDER;
     977                 :            :   LOCAL_OVERFLOW;
     978                 :            : 
     979         [ -  + ]:    6207447 :   if (k == 0) return gen_0;
     980                 :    6207447 :   d = NLIMBS(T); /* <= 2*k */
     981         [ +  + ]:    6207447 :   if (d == 0) return gen_0;
     982                 :            : #ifdef DEBUG
     983                 :            :   if (d > 2*k) pari_err_BUG("red_montgomery");
     984                 :            : #endif
     985         [ -  + ]:    6207177 :   if (k == 1)
     986                 :            :   { /* as below, special cased for efficiency */
     987                 :          0 :     ulong n = uel(N,2);
     988         [ #  # ]:          0 :     if (d == 1) {
     989                 :          0 :       hiremainder = uel(T,2);
     990                 :          0 :       m = hiremainder * inv;
     991                 :          0 :       (void)addmul(m, n); /* t + m*n = 0 */
     992                 :          0 :       return utoi(hiremainder);
     993                 :            :     } else { /* d = 2 */
     994                 :          0 :       hiremainder = uel(T,3);
     995                 :          0 :       m = hiremainder * inv;
     996                 :          0 :       (void)addmul(m, n); /* t + m*n = 0 */
     997                 :          0 :       t = addll(hiremainder, uel(T,2));
     998 [ #  # ][ #  # ]:          0 :       if (overflow) t -= n; /* t > n doesn't fit in 1 word */
     999                 :          0 :       return utoi(t);
    1000                 :            :     }
    1001                 :            :   }
    1002                 :            :   /* assume k >= 2 */
    1003                 :    6207177 :   av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
    1004                 :            : 
    1005                 :            :   /* copy T to scratch space (pad with zeroes to 2k words) */
    1006                 :    6207177 :   Td = (GEN)av;
    1007                 :    6207177 :   Te = T + (d+2);
    1008         [ +  + ]:   52592814 :   for (i=0; i < d     ; i++) *--Td = *--Te;
    1009         [ +  + ]:   15271746 :   for (   ; i < (k<<1); i++) *--Td = 0;
    1010                 :            : 
    1011                 :    6207177 :   Te = (GEN)av; /* 1 beyond end of current T mantissa (in scratch) */
    1012                 :    6207177 :   Ne = N + k+2; /* 1 beyond end of N mantissa */
    1013                 :            : 
    1014                 :    6207177 :   carry = 0;
    1015         [ +  + ]:   33932280 :   for (i=0; i<k; i++) /* set T := T/B nod N, k times */
    1016                 :            :   {
    1017                 :   27725103 :     Td = Te; /* one beyond end of (new) T mantissa */
    1018                 :   27725103 :     Nd = Ne;
    1019                 :   27725103 :     hiremainder = *--Td;
    1020                 :   27725103 :     m = hiremainder * inv; /* solve T + m N = O(B) */
    1021                 :            : 
    1022                 :            :     /* set T := (T + mN) / B */
    1023                 :   27725103 :     Te = Td;
    1024                 :   27725103 :     (void)addmul(m, *--Nd); /* = 0 */
    1025 [ +  + ][ +  + ]:  283452585 :     for (j=1; j<k; j++)
    1026                 :            :     {
    1027                 :  255727482 :       t = addll(addmul(m, *--Nd), *--Td);
    1028                 :  255727482 :       *Td = t;
    1029                 :  255727482 :       hiremainder += overflow;
    1030                 :            :     }
    1031                 :   27725103 :     t = addll(hiremainder, *--Td); *Td = t + carry;
    1032 [ +  + ][ +  + ]:   27725103 :     carry = (overflow || (carry && *Td == 0));
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
    1033                 :            :   }
    1034         [ +  + ]:    6207177 :   if (carry)
    1035                 :            :   { /* Td > N overflows (k+1 words), set Td := Td - N */
    1036                 :      87213 :     Td = Te;
    1037                 :      87213 :     Nd = Ne;
    1038                 :      87213 :     t = subll(*--Td, *--Nd); *Td = t;
    1039 [ +  + ][ +  + ]:    1251240 :     while (Td > scratch) { t = subllx(*--Td, *--Nd); *Td = t; }
    1040                 :            :   }
    1041                 :            : 
    1042                 :            :   /* copy result */
    1043                 :    6207177 :   Td = (GEN)av;
    1044 [ +  + ][ +  - ]:    9798351 :   while (*scratch == 0 && Te > scratch) scratch++; /* strip leading 0s */
    1045         [ +  + ]:   30341106 :   while (Te > scratch) *--Td = *--Te;
    1046         [ -  + ]:    6207177 :   k = (GEN)av - Td; if (!k) { avma = av; return gen_0; }
    1047                 :    6207177 :   k += 2;
    1048                 :    6207177 :   *--Td = evalsigne(1) | evallgefint(k);
    1049                 :    6207177 :   *--Td = evaltyp(t_INT) | evallg(k);
    1050                 :            : #ifdef DEBUG
    1051                 :            : {
    1052                 :            :   long l = NLIMBS(N), s = BITS_IN_LONG*l;
    1053                 :            :   GEN R = int2n(s);
    1054                 :            :   GEN res = remii(mulii(T, Fp_inv(R, N)), N);
    1055                 :            :   if (k > lgefint(N)
    1056                 :            :     || !equalii(remii(Td,N),res)
    1057                 :            :     || cmpii(Td, addii(shifti(T, -s), N)) >= 0) pari_err_BUG("red_montgomery");
    1058                 :            : }
    1059                 :            : #endif
    1060                 :    6207447 :   avma = (pari_sp)Td; return Td;
    1061                 :            : }
    1062                 :            : 
    1063                 :            : /* EXACT INTEGER DIVISION */
    1064                 :            : 
    1065                 :            : /* assume xy>0, the division is exact and y is odd. Destroy x */
    1066                 :            : static GEN
    1067                 :   10681032 : diviuexact_i(GEN x, ulong y)
    1068                 :            : {
    1069                 :            :   long i, lz, lx;
    1070                 :            :   ulong q, yinv;
    1071                 :            :   GEN z, z0, x0, x0min;
    1072                 :            : 
    1073         [ +  + ]:   10681032 :   if (y == 1) return icopy(x);
    1074                 :    9707304 :   lx = lgefint(x);
    1075         [ +  + ]:    9707304 :   if (lx == 3) return utoipos(uel(x,2) / y);
    1076                 :    9325338 :   yinv = invmod2BIL(y);
    1077         [ +  + ]:    9325338 :   lz = (y <= uel(x,2)) ? lx : lx-1;
    1078                 :    9325338 :   z = new_chunk(lz);
    1079                 :    9325338 :   z0 = z + lz;
    1080                 :    9325338 :   x0 = x + lx; x0min = x + lx-lz+2;
    1081                 :            : 
    1082         [ +  + ]:   23000097 :   while (x0 > x0min)
    1083                 :            :   {
    1084                 :   13674759 :     *--z0 = q = yinv*((ulong)*--x0); /* i-th quotient */
    1085         [ +  + ]:   13674759 :     if (!q) continue;
    1086                 :            :     /* x := x - q * y */
    1087                 :            :     { /* update neither lowest word (could set it to 0) nor highest ones */
    1088                 :   13593543 :       register GEN x1 = x0 - 1;
    1089                 :            :       LOCAL_HIREMAINDER;
    1090                 :   13593543 :       (void)mulll(q,y);
    1091 [ +  + ][ +  + ]:   13593543 :       if (hiremainder)
    1092                 :            :       {
    1093         [ +  + ]:   11484033 :         if ((ulong)*x1 < hiremainder)
    1094                 :            :         {
    1095                 :      33810 :           *x1 -= hiremainder;
    1096         [ +  + ]:      36072 :           do (*--x1)--; while ((ulong)*x1 == ULONG_MAX);
    1097                 :            :         }
    1098                 :            :         else
    1099                 :   11450223 :           *x1 -= hiremainder;
    1100                 :            :       }
    1101                 :            :     }
    1102                 :            :   }
    1103         [ -  + ]:    9325338 :   i=2; while(!z[i]) i++;
    1104                 :    9325338 :   z += i-2; lz -= i-2;
    1105                 :    9325338 :   z[0] = evaltyp(t_INT)|evallg(lz);
    1106                 :    9325338 :   z[1] = evalsigne(1)|evallg(lz);
    1107                 :   10681032 :   avma = (pari_sp)z; return z;
    1108                 :            : }
    1109                 :            : 
    1110                 :            : /* assume y != 0 and the division is exact */
    1111                 :            : GEN
    1112                 :    3517665 : diviuexact(GEN x, ulong y)
    1113                 :            : {
    1114                 :            :   pari_sp av;
    1115                 :    3517665 :   long lx, vy, s = signe(x);
    1116                 :            :   GEN z;
    1117                 :            : 
    1118         [ +  + ]:    3517665 :   if (!s) return gen_0;
    1119         [ +  + ]:    3320568 :   if (y == 1) return icopy(x);
    1120                 :    3316821 :   lx = lgefint(x);
    1121         [ +  + ]:    3316821 :   if (lx == 3) {
    1122                 :    3169824 :     ulong q = uel(x,2) / y;
    1123         [ +  + ]:    3169824 :     return (s > 0)? utoipos(q): utoineg(q);
    1124                 :            :   }
    1125                 :     146997 :   av = avma; (void)new_chunk(lx); vy = vals(y);
    1126         [ +  + ]:     146997 :   if (vy) {
    1127                 :      92319 :     y >>= vy;
    1128         [ +  + ]:      92319 :     if (y == 1) { avma = av; return shifti(x, -vy); }
    1129                 :      70416 :     x = shifti(x, -vy);
    1130         [ -  + ]:      70416 :     if (lx == 3) {
    1131                 :          0 :       ulong q = uel(x,2) / y;
    1132                 :          0 :       avma = av;
    1133         [ #  # ]:          0 :       return (s > 0)? utoipos(q): utoineg(q);
    1134                 :            :     }
    1135                 :      54678 :   } else x = icopy(x);
    1136                 :     125094 :   avma = av;
    1137                 :     125094 :   z = diviuexact_i(x, y);
    1138                 :    3517665 :   setsigne(z, s); return z;
    1139                 :            : }
    1140                 :            : 
    1141                 :            : /* Find z such that x=y*z, knowing that y | x (unchecked)
    1142                 :            :  * Method: y0 z0 = x0 mod B = 2^BITS_IN_LONG ==> z0 = 1/y0 mod B.
    1143                 :            :  *    Set x := (x - z0 y) / B, updating only relevant words, and repeat */
    1144                 :            : GEN
    1145                 :  163832727 : diviiexact(GEN x, GEN y)
    1146                 :            : {
    1147                 :  163832727 :   long lx, ly, lz, vy, i, ii, sx = signe(x), sy = signe(y);
    1148                 :            :   pari_sp av;
    1149                 :            :   ulong y0inv,q;
    1150                 :            :   GEN z;
    1151                 :            : 
    1152         [ -  + ]:  163832727 :   if (!sy) pari_err_INV("diviiexact",gen_0);
    1153         [ +  + ]:  163832727 :   if (!sx) return gen_0;
    1154                 :   64527690 :   lx = lgefint(x);
    1155         [ +  + ]:   64527690 :   if (lx == 3) {
    1156                 :   49711818 :     q = uel(x,2) / uel(y,2);
    1157         [ +  + ]:   49711818 :     return (sx+sy) ? utoipos(q): utoineg(q);
    1158                 :            :   }
    1159                 :   14815872 :   vy = vali(y); av = avma;
    1160                 :   14815872 :   (void)new_chunk(lx); /* enough room for z */
    1161         [ +  + ]:   14815872 :   if (vy)
    1162                 :            :   { /* make y odd */
    1163                 :    6849438 :     y = shifti(y,-vy);
    1164                 :    6849438 :     x = shifti(x,-vy); lx = lgefint(x);
    1165                 :            :   }
    1166                 :    7966434 :   else x = icopy(x); /* necessary because we destroy x */
    1167                 :   14815872 :   avma = av; /* will erase our x,y when exiting */
    1168                 :            :   /* now y is odd */
    1169                 :   14815872 :   ly = lgefint(y);
    1170         [ +  + ]:   14815872 :   if (ly == 3)
    1171                 :            :   {
    1172                 :   10555938 :     x = diviuexact_i(x,uel(y,2)); /* x != 0 */
    1173         [ +  + ]:   10555938 :     setsigne(x, (sx+sy)? 1: -1); return x;
    1174                 :            :   }
    1175                 :    4259934 :   y0inv = invmod2BIL(y[ly-1]);
    1176 [ +  + ][ +  + ]:    7657695 :   i=2; while (i<ly && y[i]==x[i]) i++;
    1177 [ +  + ][ +  + ]:    4259934 :   lz = (i==ly || uel(y,i) < uel(x,i)) ? lx-ly+3 : lx-ly+2;
    1178                 :    4259934 :   z = new_chunk(lz);
    1179                 :            : 
    1180                 :    4259934 :   y += ly - 1; /* now y[-i] = i-th word of y */
    1181         [ +  + ]:   16301106 :   for (ii=lx-1,i=lz-1; i>=2; i--,ii--)
    1182                 :            :   {
    1183                 :            :     long limj;
    1184                 :            :     LOCAL_HIREMAINDER;
    1185                 :            :     LOCAL_OVERFLOW;
    1186                 :            : 
    1187                 :   12041172 :     z[i] = q = y0inv*uel(x,ii); /* i-th quotient */
    1188         [ +  + ]:   12041172 :     if (!q) continue;
    1189                 :            : 
    1190                 :            :     /* x := x - q * y */
    1191                 :   12041130 :     (void)mulll(q,y[0]); limj = maxss(lx - lz, ii+3-ly);
    1192                 :            :     { /* update neither lowest word (could set it to 0) nor highest ones */
    1193                 :   12041130 :       register GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj;
    1194 [ +  + ][ +  + ]:   59328993 :       for (; x0 >= xlim; x0--, y0--)
    1195                 :            :       {
    1196                 :   47287863 :         *x0 = subll(*x0, addmul(q,*y0));
    1197                 :   47287863 :         hiremainder += overflow;
    1198                 :            :       }
    1199 [ +  + ][ +  + ]:   12041130 :       if (hiremainder && limj != lx - lz)
    1200                 :            :       {
    1201         [ +  + ]:    5445825 :         if ((ulong)*x0 < hiremainder)
    1202                 :            :         {
    1203                 :      41400 :           *x0 -= hiremainder;
    1204         [ -  + ]:      41400 :           do (*--x0)--; while ((ulong)*x0 == ULONG_MAX);
    1205                 :            :         }
    1206                 :            :         else
    1207                 :    5404425 :           *x0 -= hiremainder;
    1208                 :            :       }
    1209                 :            :     }
    1210                 :            :   }
    1211         [ -  + ]:    4259934 :   i=2; while(!z[i]) i++;
    1212                 :    4259934 :   z += i-2; lz -= (i-2);
    1213                 :    4259934 :   z[0] = evaltyp(t_INT)|evallg(lz);
    1214         [ +  + ]:    4259934 :   z[1] = evalsigne((sx+sy)? 1: -1) | evallg(lz);
    1215                 :  163832727 :   avma = (pari_sp)z; return z;
    1216                 :            : }
    1217                 :            : 
    1218                 :            : /* assume yz != and yz | x */
    1219                 :            : GEN
    1220                 :      39897 : diviuuexact(GEN x, ulong y, ulong z)
    1221                 :            : {
    1222                 :            :   long tmp[4];
    1223                 :            :   ulong t;
    1224                 :            :   LOCAL_HIREMAINDER;
    1225                 :      39897 :   t = mulll(y, z);
    1226 [ +  - ][ +  - ]:      39897 :   if (!hiremainder) return diviuexact(x, t);
    1227                 :          0 :   tmp[0] = evaltyp(t_INT)|_evallg(4);
    1228                 :          0 :   tmp[1] = evalsigne(1)|evallgefint(4);
    1229                 :          0 :   tmp[2] = hiremainder;
    1230                 :          0 :   tmp[3] = t;
    1231                 :      39897 :   return diviiexact(x, tmp);
    1232                 :            : }
    1233                 :            : 
    1234                 :            : /********************************************************************/
    1235                 :            : /**                                                                **/
    1236                 :            : /**               INTEGER MULTIPLICATION (BASECASE)                **/
    1237                 :            : /**                                                                **/
    1238                 :            : /********************************************************************/
    1239                 :            : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
    1240                 :            : INLINE GEN
    1241                 :  928193850 : muliispec_basecase(GEN x, GEN y, long nx, long ny)
    1242                 :            : {
    1243                 :            :   GEN z2e,z2d,yd,xd,ye,zd;
    1244                 :            :   long p1,lz;
    1245                 :            :   LOCAL_HIREMAINDER;
    1246                 :            : 
    1247         [ +  + ]:  928193850 :   if (ny == 1) return muluispec((ulong)*y, x, nx);
    1248         [ +  + ]:  368718522 :   if (ny == 0) return gen_0;
    1249                 :  368185980 :   zd = (GEN)avma; lz = nx+ny+2;
    1250                 :  368185980 :   (void)new_chunk(lz);
    1251                 :  368185980 :   xd = x + nx;
    1252                 :  368185980 :   yd = y + ny;
    1253                 :  368185980 :   ye = yd; p1 = *--xd;
    1254                 :            : 
    1255                 :  368185980 :   *--zd = mulll(p1, *--yd); z2e = zd;
    1256      [ +  +  + ]: 1789268541 :   while (yd > y) *--zd = addmul(p1, *--yd);
    1257                 :  368185980 :   *--zd = hiremainder;
    1258                 :            : 
    1259         [ +  + ]: 1978078167 :   while (xd > x)
    1260                 :            :   {
    1261                 :            :     LOCAL_OVERFLOW;
    1262                 : 1609892187 :     yd = ye; p1 = *--xd;
    1263                 :            : 
    1264                 : 1609892187 :     z2d = --z2e;
    1265                 : 1609892187 :     *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
    1266 [ +  + ][ +  + ]:15218154825 :     while (yd > y)
    1267                 :            :     {
    1268                 :13608262638 :       hiremainder += overflow;
    1269                 :13608262638 :       *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
    1270                 :            :     }
    1271                 : 1609892187 :     *--zd = hiremainder + overflow;
    1272                 :            :   }
    1273         [ +  + ]:  368185980 :   if (*zd == 0) { zd++; lz--; } /* normalize */
    1274                 :  368185980 :   *--zd = evalsigne(1) | evallgefint(lz);
    1275                 :  368185980 :   *--zd = evaltyp(t_INT) | evallg(lz);
    1276                 :  928193850 :   avma=(pari_sp)zd; return zd;
    1277                 :            : }
    1278                 :            : 
    1279                 :            : INLINE GEN
    1280                 :  452824479 : sqrispec_basecase(GEN x, long nx)
    1281                 :            : {
    1282                 :            :   GEN z2e,z2d,yd,xd,zd,x0,z0;
    1283                 :            :   long p1,lz;
    1284                 :            :   LOCAL_HIREMAINDER;
    1285                 :            :   LOCAL_OVERFLOW;
    1286                 :            : 
    1287         [ +  + ]:  452824479 :   if (nx == 1) return sqru((ulong)*x);
    1288         [ +  + ]:  401057367 :   if (nx == 0) return gen_0;
    1289                 :   40397355 :   zd = (GEN)avma; lz = (nx+1) << 1;
    1290                 :   40397355 :   z0 = new_chunk(lz);
    1291         [ -  + ]:   40397355 :   if (nx == 1)
    1292                 :            :   {
    1293                 :          0 :     *--zd = mulll(*x, *x);
    1294                 :          0 :     *--zd = hiremainder; goto END;
    1295                 :            :   }
    1296                 :   40397355 :   xd = x + nx;
    1297                 :            : 
    1298                 :            :   /* compute double products --> zd */
    1299                 :   40397355 :   p1 = *--xd; yd = xd; --zd;
    1300                 :   40397355 :   *--zd = mulll(p1, *--yd); z2e = zd;
    1301 [ +  + ][ +  + ]:  201813951 :   while (yd > x) *--zd = addmul(p1, *--yd);
    1302                 :   40397355 :   *--zd = hiremainder;
    1303                 :            : 
    1304                 :   40397355 :   x0 = x+1;
    1305         [ +  + ]:  201813951 :   while (xd > x0)
    1306                 :            :   {
    1307                 :            :     LOCAL_OVERFLOW;
    1308                 :  161416596 :     p1 = *--xd; yd = xd;
    1309                 :            : 
    1310                 :  161416596 :     z2e -= 2; z2d = z2e;
    1311                 :  161416596 :     *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
    1312 [ +  + ][ +  + ]: 1378800474 :     while (yd > x)
    1313                 :            :     {
    1314                 : 1217383878 :       hiremainder += overflow;
    1315                 : 1217383878 :       *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
    1316                 :            :     }
    1317                 :  161416596 :     *--zd = hiremainder + overflow;
    1318                 :            :   }
    1319                 :            :   /* multiply zd by 2 (put result in zd - 1) */
    1320                 :   40397355 :   zd[-1] = ((*zd & HIGHBIT) != 0);
    1321                 :   40397355 :   shift_left(zd, zd, 0, (nx<<1)-3, 0, 1);
    1322                 :            : 
    1323                 :            :   /* add the squares */
    1324                 :   40397355 :   xd = x + nx; zd = z0 + lz;
    1325                 :   40397355 :   p1 = *--xd;
    1326                 :   40397355 :   zd--; *zd = mulll(p1,p1);
    1327                 :   40397355 :   zd--; *zd = addll(hiremainder, *zd);
    1328         [ +  + ]:  242211306 :   while (xd > x)
    1329                 :            :   {
    1330                 :  201813951 :     p1 = *--xd;
    1331                 :  201813951 :     zd--; *zd = addll(mulll(p1,p1)+ overflow, *zd);
    1332                 :  201813951 :     zd--; *zd = addll(hiremainder + overflow, *zd);
    1333                 :            :   }
    1334                 :            : 
    1335                 :            : END:
    1336         [ +  + ]:   40397355 :   if (*zd == 0) { zd++; lz--; } /* normalize */
    1337                 :   40397355 :   *--zd = evalsigne(1) | evallgefint(lz);
    1338                 :   40397355 :   *--zd = evaltyp(t_INT) | evallg(lz);
    1339                 :  452824479 :   avma=(pari_sp)zd; return zd;
    1340                 :            : }
    1341                 :            : 
    1342                 :            : /********************************************************************/
    1343                 :            : /**                                                                **/
    1344                 :            : /**               INTEGER MULTIPLICATION (FFT)                     **/
    1345                 :            : /**                                                                **/
    1346                 :            : /********************************************************************/
    1347                 :            : 
    1348                 :            : /*
    1349                 :            :  Compute parameters for FFT:
    1350                 :            :    len: result length
    1351                 :            :    k: FFT depth.
    1352                 :            :    n: number of blocks (2^k)
    1353                 :            :    bs: block size
    1354                 :            :    mod: Modulus is M=2^(BIL*mod)+1
    1355                 :            :    ord: order of 2 in Z/MZ.
    1356                 :            :  We must have:
    1357                 :            :    bs*n >= l
    1358                 :            :    2^(BIL*mod) > nb*2^(2*BIL*bs)
    1359                 :            :    2^k | 2*BIL*mod
    1360                 :            : */
    1361                 :            : static void
    1362                 :       2532 : mulliifft_params(long len, long *k, long *mod, long *bs, long *n, ulong *ord)
    1363                 :            : {
    1364                 :            :   long r;
    1365                 :       2532 :   *k = expu((3*len)>>2)-3;
    1366                 :            :   do {
    1367                 :       2532 :     (*k)--;
    1368                 :       2532 :     r  = *k-(TWOPOTBITS_IN_LONG+2);
    1369                 :       2532 :     *n = 1L<<*k;
    1370                 :       2532 :     *bs = (len+*n-1)>>*k;
    1371                 :       2532 :     *mod= 2**bs+1;
    1372         [ +  + ]:       2532 :     if (r>0)
    1373                 :        153 :       *mod=((*mod+(1L<<r)-1)>>r)<<r;
    1374         [ -  + ]:       2532 :   } while(*mod>=3**bs);
    1375                 :       2532 :   *ord= 4**mod*BITS_IN_LONG;
    1376                 :       2532 : }
    1377                 :            : 
    1378                 :            : /* Zf_: arithmetic in ring Z/MZ where M= 2^(BITS_IN_LONG*mod)+1
    1379                 :            :  * for some mod.
    1380                 :            :  * Do not garbage collect.
    1381                 :            :  */
    1382                 :            : 
    1383                 :            : static GEN
    1384                 :    5954496 : Zf_add(GEN a, GEN b, GEN M)
    1385                 :            : {
    1386                 :    5954496 :   GEN y, z = addii(a,b);
    1387                 :    5954496 :   long mod = lgefint(M)-3;
    1388                 :    5954496 :   long l = NLIMBS(z);
    1389         [ +  + ]:    5954496 :   if (l<=mod) return z;
    1390                 :    1795956 :   y = subis(z, 1);
    1391         [ -  + ]:    1795956 :   if (NLIMBS(y)<=mod) return z;
    1392                 :    5954496 :   return int_normalize(y,1);
    1393                 :            : }
    1394                 :            : 
    1395                 :            : static GEN
    1396                 :    6093231 : Zf_sub(GEN a, GEN b, GEN M)
    1397                 :            : {
    1398                 :    6093231 :   GEN z = subii(a,b);
    1399         [ +  + ]:    6093231 :   return signe(z)>=0? z: addii(M,z);
    1400                 :            : }
    1401                 :            : 
    1402                 :            : /* destroy z */
    1403                 :            : static GEN
    1404                 :   10683852 : Zf_red_destroy(GEN z, GEN M)
    1405                 :            : {
    1406                 :   10683852 :   long mod = lgefint(M)-3;
    1407                 :   10683852 :   long l = NLIMBS(z);
    1408                 :            :   GEN y;
    1409         [ +  + ]:   10683852 :   if (l<=mod) return z;
    1410                 :    4537806 :   y = shifti(z, -mod*BITS_IN_LONG);
    1411                 :    4537806 :   z = int_normalize(z, NLIMBS(y));
    1412                 :    4537806 :   y = Zf_red_destroy(y, M);
    1413                 :    4537806 :   z = subii(z, y);
    1414         [ +  + ]:    4537806 :   if (signe(z)<0) z = addii(z, M);
    1415                 :   10683852 :   return z;
    1416                 :            : }
    1417                 :            : 
    1418                 :            : INLINE GEN
    1419                 :    5615742 : Zf_shift(GEN a, ulong s, GEN M) { return Zf_red_destroy(shifti(a, s), M); }
    1420                 :            : 
    1421                 :            : /*
    1422                 :            :  Multiply by sqrt(2)^s
    1423                 :            :  We use the formula sqrt(2)=z_8*(1-z_4)) && z_8=2^(ord/16) [2^(ord/4)+1]
    1424                 :            : */
    1425                 :            : 
    1426                 :            : static GEN
    1427                 :    5954496 : Zf_mulsqrt2(GEN a, ulong s, ulong ord, GEN M)
    1428                 :            : {
    1429                 :    5954496 :   ulong hord = ord>>1;
    1430         [ +  + ]:    5954496 :   if (!signe(a)) return gen_0;
    1431         [ +  + ]:    4807968 :   if (odd(s)) /* Multiply by 2^(s/2) */
    1432                 :            :   {
    1433                 :     138735 :     GEN az8  = Zf_shift(a,   ord>>4, M);
    1434                 :     138735 :     GEN az83 = Zf_shift(az8, ord>>3, M);
    1435                 :     138735 :     a = Zf_sub(az8, az83, M);
    1436                 :     138735 :     s--;
    1437                 :            :   }
    1438         [ +  + ]:    4807968 :   if (s < hord)
    1439                 :    3530406 :     return Zf_shift(a, s>>1, M);
    1440                 :            :   else
    1441                 :    5954496 :     return subii(M,Zf_shift(a, (s-hord)>>1, M));
    1442                 :            : }
    1443                 :            : 
    1444                 :            : INLINE GEN
    1445                 :      83712 : Zf_sqr(GEN a, GEN M) { return Zf_red_destroy(sqri(a), M); }
    1446                 :            : 
    1447                 :            : INLINE GEN
    1448                 :     446592 : Zf_mul(GEN a, GEN b, GEN M) { return Zf_red_destroy(mulii(a,b), M); }
    1449                 :            : 
    1450                 :            : /* In place, bit reversing FFT */
    1451                 :            : static void
    1452                 :     972147 : muliifft_dit(ulong o, ulong ord, GEN M, GEN FFT, long d, long step)
    1453                 :            : {
    1454                 :     972147 :   pari_sp av = avma;
    1455                 :            :   long i;
    1456                 :     972147 :   ulong j, no = (o<<1)%ord;
    1457                 :     972147 :   long hstep=step>>1;
    1458         [ +  + ]:    4826355 :   for (i = d+1, j = 0; i <= d+hstep; ++i, j =(j+o)%ord)
    1459                 :            :   {
    1460                 :    3854208 :     GEN a = Zf_add(gel(FFT,i), gel(FFT,i+hstep), M);
    1461                 :    3854208 :     GEN b = Zf_mulsqrt2(Zf_sub(gel(FFT,i), gel(FFT,i+hstep), M), j, ord, M);
    1462                 :    3854208 :     affii(a,gel(FFT,i));
    1463                 :    3854208 :     affii(b,gel(FFT,i+hstep));
    1464                 :    3854208 :     avma = av;
    1465                 :            :   }
    1466         [ +  + ]:     972147 :   if (hstep>1)
    1467                 :            :   {
    1468                 :     483699 :     muliifft_dit(no, ord, M, FFT, d, hstep);
    1469                 :     483699 :     muliifft_dit(no, ord, M, FFT, d+hstep, hstep);
    1470                 :            :   }
    1471                 :     972147 : }
    1472                 :            : 
    1473                 :            : /* In place, bit reversed FFT, inverse of muliifft_dit */
    1474                 :            : static void
    1475                 :     527772 : muliifft_dis(ulong o, ulong ord, GEN M, GEN FFT, long d, long step)
    1476                 :            : {
    1477                 :     527772 :   pari_sp av = avma;
    1478                 :            :   long i;
    1479                 :     527772 :   ulong j, no = (o<<1)%ord;
    1480                 :     527772 :   long hstep=step>>1;
    1481         [ +  + ]:     527772 :   if (hstep>1)
    1482                 :            :   {
    1483                 :     262620 :     muliifft_dis(no, ord, M, FFT, d, hstep);
    1484                 :     262620 :     muliifft_dis(no, ord, M, FFT, d+hstep, hstep);
    1485                 :            :   }
    1486         [ +  + ]:    2628060 :   for (i = d+1, j = 0; i <= d+hstep; ++i, j =(j+o)%ord)
    1487                 :            :   {
    1488                 :    2100288 :     GEN z = Zf_mulsqrt2(gel(FFT,i+hstep), j, ord, M);
    1489                 :    2100288 :     GEN a = Zf_add(gel(FFT,i), z, M);
    1490                 :    2100288 :     GEN b = Zf_sub(gel(FFT,i), z, M);
    1491                 :    2100288 :     affii(a,gel(FFT,i));
    1492                 :    2100288 :     affii(b,gel(FFT,i+hstep));
    1493                 :    2100288 :     avma = av;
    1494                 :            :   }
    1495                 :     527772 : }
    1496                 :            : 
    1497                 :            : static GEN
    1498                 :       4749 : muliifft_spliti(GEN a, long na, long bs, long n, long mod)
    1499                 :            : {
    1500                 :       4749 :   GEN ap = a+na-1;
    1501                 :       4749 :   GEN c  = cgetg(n+1, t_VEC);
    1502                 :            :   long i,j;
    1503         [ +  + ]:     981645 :   for(i=1;i<=n;i++)
    1504                 :            :   {
    1505                 :     976896 :     GEN z = cgeti(mod+3);
    1506         [ +  + ]:     976896 :     if (na)
    1507                 :            :     {
    1508                 :     481404 :       long m = minss(bs, na), v=0;
    1509                 :     481404 :       GEN zp, aa=ap-m+1;
    1510 [ +  + ][ +  + ]:    8552346 :       while (!*aa && v<m) {aa++; v++;}
    1511                 :     481404 :       zp = z+m-v+1;
    1512         [ +  + ]:    6781143 :       for (j=v; j < m; j++)
    1513                 :    6299739 :         *zp-- = *ap--;
    1514                 :     481404 :       ap -= v; na -= m;
    1515         [ +  + ]:     481404 :       z[1] = evalsigne(m!=v) | evallgefint(m-v+2);
    1516                 :            :     } else
    1517                 :     495492 :       z[1] = evalsigne(0) | evallgefint(2);
    1518                 :     976896 :     gel(c, i) = z;
    1519                 :            :   }
    1520                 :       4749 :   return c;
    1521                 :            : }
    1522                 :            : 
    1523                 :            : static GEN
    1524                 :       2532 : muliifft_unspliti(GEN V, long bs, long len)
    1525                 :            : {
    1526                 :       2532 :   long s, i, j, l = lg(V);
    1527                 :       2532 :   GEN a = cgeti(len);
    1528                 :       2532 :   a[1] = evalsigne(1)|evallgefint(len);
    1529         [ +  + ]:   15650619 :   for(i=2;i<len;i++)
    1530                 :   15648087 :     a[i] = 0;
    1531         [ +  + ]:     532836 :   for(i=1, s=0; i<l; i++, s+=bs)
    1532                 :            :   {
    1533                 :     530304 :     GEN  u = gel(V,i);
    1534         [ +  + ]:     530304 :     if (signe(u))
    1535                 :            :     {
    1536                 :     393462 :       GEN ap = int_W(a,s);
    1537                 :     393462 :       GEN up = int_LSW(u);
    1538                 :     393462 :       long lu = NLIMBS(u);
    1539                 :            :       LOCAL_OVERFLOW;
    1540                 :     393462 :       *ap = addll(*ap, *up--); ap--;
    1541 [ +  + ][ +  + ]:   20010216 :       for (j=1; j<lu; j++)
    1542                 :   19616754 :        { *ap = addllx(*ap, *up--); ap--; }
    1543         [ +  + ]:     393624 :       while (overflow)
    1544                 :        162 :        { *ap = addll(*ap, 1); ap--; }
    1545                 :            :     }
    1546                 :            :   }
    1547                 :       2532 :   return int_normalize(a,0);
    1548                 :            : }
    1549                 :            : 
    1550                 :            : static GEN
    1551                 :        315 : sqrispec_fft(GEN a, long na)
    1552                 :            : {
    1553                 :        315 :   pari_sp av, ltop = avma;
    1554                 :        315 :   long len = 2*na;
    1555                 :            :   long k, mod, bs, n;
    1556                 :            :   GEN  FFT, M;
    1557                 :            :   long i;
    1558                 :            :   ulong o, ord;
    1559                 :            : 
    1560                 :        315 :   mulliifft_params(len,&k,&mod,&bs,&n,&ord);
    1561                 :        315 :   o = ord>>k;
    1562                 :        315 :   M = int2n(mod*BITS_IN_LONG);
    1563                 :        315 :   M[2+mod] = 1;
    1564                 :        315 :   FFT = muliifft_spliti(a, na, bs, n, mod);
    1565                 :        315 :   muliifft_dit(o, ord, M, FFT, 0, n);
    1566                 :        315 :   av = avma;
    1567         [ +  + ]:      84027 :   for(i=1; i<=n; i++)
    1568                 :            :   {
    1569                 :      83712 :     affii(Zf_sqr(gel(FFT,i), M), gel(FFT,i));
    1570                 :      83712 :     avma=av;
    1571                 :            :   }
    1572                 :        315 :   muliifft_dis(ord-o, ord, M, FFT, 0, n);
    1573         [ +  + ]:      84027 :   for(i=1; i<=n; i++)
    1574                 :            :   {
    1575                 :      83712 :     affii(Zf_shift(gel(FFT,i), (ord>>1)-k, M), gel(FFT,i));
    1576                 :      83712 :     avma=av;
    1577                 :            :   }
    1578                 :        315 :   return gerepileuptoint(ltop, muliifft_unspliti(FFT,bs,2+len));
    1579                 :            : }
    1580                 :            : 
    1581                 :            : static GEN
    1582                 :       2217 : muliispec_fft(GEN a, GEN b, long na, long nb)
    1583                 :            : {
    1584                 :       2217 :   pari_sp av, av2, ltop = avma;
    1585                 :       2217 :   long len = na+nb;
    1586                 :            :   long k, mod, bs, n;
    1587                 :            :   GEN FFT, FFTb, M;
    1588                 :            :   long i;
    1589                 :            :   ulong o, ord;
    1590                 :            : 
    1591                 :       2217 :   mulliifft_params(len,&k,&mod,&bs,&n,&ord);
    1592                 :       2217 :   o = ord>>k;
    1593                 :       2217 :   M = int2n(mod*BITS_IN_LONG);
    1594                 :       2217 :   M[2+mod] = 1;
    1595                 :       2217 :   FFT = muliifft_spliti(a, na, bs, n, mod);
    1596                 :       2217 :   av=avma;
    1597                 :       2217 :   muliifft_dit(o, ord, M, FFT, 0, n);
    1598                 :       2217 :   FFTb = muliifft_spliti(b, nb, bs, n, mod);
    1599                 :       2217 :   av2 = avma;
    1600                 :       2217 :   muliifft_dit(o, ord, M, FFTb, 0, n);
    1601         [ +  + ]:     448809 :   for(i=1; i<=n; i++)
    1602                 :            :   {
    1603                 :     446592 :     affii(Zf_mul(gel(FFT,i), gel(FFTb,i), M), gel(FFT,i));
    1604                 :     446592 :     avma=av2;
    1605                 :            :   }
    1606                 :       2217 :   avma=av;
    1607                 :       2217 :   muliifft_dis(ord-o, ord, M, FFT, 0, n);
    1608         [ +  + ]:     448809 :   for(i=1; i<=n; i++)
    1609                 :            :   {
    1610                 :     446592 :     affii(Zf_shift(gel(FFT,i),(ord>>1)-k,M), gel(FFT,i));
    1611                 :     446592 :     avma=av;
    1612                 :            :   }
    1613                 :       2217 :   return gerepileuptoint(ltop, muliifft_unspliti(FFT,bs,2+len));
    1614                 :            : }
    1615                 :            : 
    1616                 :            : /********************************************************************/
    1617                 :            : /**                                                                **/
    1618                 :            : /**               INTEGER MULTIPLICATION (KARATSUBA)               **/
    1619                 :            : /**                                                                **/
    1620                 :            : /********************************************************************/
    1621                 :            : 
    1622                 :            : /* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */
    1623                 :            : static GEN
    1624                 :   66968466 : addshiftw(GEN x, GEN y, long d)
    1625                 :            : {
    1626                 :   66968466 :   GEN z,z0,y0,yd, zd = (GEN)avma;
    1627                 :   66968466 :   long a,lz,ly = lgefint(y);
    1628                 :            : 
    1629                 :   66968466 :   z0 = new_chunk(d);
    1630                 :   66968466 :   a = ly-2; yd = y+ly;
    1631         [ +  + ]:   66968466 :   if (a >= d)
    1632                 :            :   {
    1633         [ +  + ]: 1291530930 :     y0 = yd-d; while (yd > y0) *--zd = *--yd; /* copy last d words of y */
    1634                 :   65140542 :     a -= d;
    1635         [ +  + ]:   65140542 :     if (a)
    1636                 :   40700604 :       z = addiispec(LIMBS(x), LIMBS(y), NLIMBS(x), a);
    1637                 :            :     else
    1638                 :   24439938 :       z = icopy(x);
    1639                 :            :   }
    1640                 :            :   else
    1641                 :            :   {
    1642         [ +  + ]:    7816140 :     y0 = yd-a; while (yd > y0) *--zd = *--yd; /* copy last a words of y */
    1643         [ +  + ]:   33466038 :     while (zd > z0) *--zd = 0;    /* complete with 0s */
    1644                 :    1827924 :     z = icopy(x);
    1645                 :            :   }
    1646                 :   66968466 :   lz = lgefint(z)+d;
    1647                 :   66968466 :   z[1] = evalsigne(1) | evallgefint(lz);
    1648                 :   66968466 :   z[0] = evaltyp(t_INT) | evallg(lz); return z;
    1649                 :            : }
    1650                 :            : 
    1651                 :            : /* Fast product (Karatsuba) of integers. a and b are "special" GENs
    1652                 :            :  * c,c0,c1,c2 are genuine GENs.
    1653                 :            :  */
    1654                 :            : GEN
    1655                 :  948490566 : muliispec(GEN a, GEN b, long na, long nb)
    1656                 :            : {
    1657                 :            :   GEN a0,c,c0;
    1658                 :            :   long n0, n0a, i;
    1659                 :            :   pari_sp av;
    1660                 :            : 
    1661         [ +  + ]:  948490566 :   if (na < nb) swapspec(a,b, na,nb);
    1662         [ +  + ]:  948490566 :   if (nb < MULII_KARATSUBA_LIMIT) return muliispec_basecase(a,b,na,nb);
    1663         [ +  + ]:   20296716 :   if (nb >= MULII_FFT_LIMIT)      return muliispec_fft(a,b,na,nb);
    1664                 :   20294499 :   i=(na>>1); n0=na-i; na=i;
    1665                 :   20294499 :   av=avma; a0=a+na; n0a=n0;
    1666 [ +  + ][ +  + ]:   58794378 :   while (n0a && !*a0) { a0++; n0a--; }
    1667                 :            : 
    1668 [ +  + ][ +  + ]:   20294499 :   if (n0a && nb > n0)
    1669                 :   19247088 :   { /* nb <= na <= n0 */
    1670                 :            :     GEN b0,c1,c2;
    1671                 :            :     long n0b;
    1672                 :            : 
    1673                 :   19247088 :     nb -= n0;
    1674                 :   19247088 :     c = muliispec(a,b,na,nb);
    1675                 :   19247088 :     b0 = b+nb; n0b = n0;
    1676 [ +  + ][ +  + ]:   57939237 :     while (n0b && !*b0) { b0++; n0b--; }
    1677         [ +  + ]:   19247088 :     if (n0b)
    1678                 :            :     {
    1679                 :   18667380 :       c0 = muliispec(a0,b0, n0a,n0b);
    1680                 :            : 
    1681                 :   18667380 :       c2 = addiispec(a0,a, n0a,na);
    1682                 :   18667380 :       c1 = addiispec(b0,b, n0b,nb);
    1683                 :   18667380 :       c1 = muliispec(LIMBS(c1),LIMBS(c2), NLIMBS(c1),NLIMBS(c2));
    1684                 :   18667380 :       c2 = addiispec(LIMBS(c0),LIMBS(c),  NLIMBS(c0),NLIMBS(c));
    1685                 :            : 
    1686                 :   18667380 :       c1 = subiispec(LIMBS(c1),LIMBS(c2), NLIMBS(c1),NLIMBS(c2));
    1687                 :            :     }
    1688                 :            :     else
    1689                 :            :     {
    1690                 :     579708 :       c0 = gen_0;
    1691                 :     579708 :       c1 = muliispec(a0,b, n0a,nb);
    1692                 :            :     }
    1693                 :   19247088 :     c = addshiftw(c,c1, n0);
    1694                 :            :   }
    1695                 :            :   else
    1696                 :            :   {
    1697                 :    1047411 :     c = muliispec(a,b,na,nb);
    1698                 :    1047411 :     c0 = muliispec(a0,b,n0a,nb);
    1699                 :            :   }
    1700                 :  948490566 :   return gerepileuptoint(av, addshiftw(c,c0, n0));
    1701                 :            : }
    1702                 :            : GEN
    1703                 :      56553 : muluui(ulong x, ulong y, GEN z)
    1704                 :            : {
    1705                 :      56553 :   long t, s = signe(z);
    1706                 :            :   GEN r;
    1707                 :            :   LOCAL_HIREMAINDER;
    1708                 :            : 
    1709 [ +  - ][ +  - ]:      56553 :   if (!x || !y || !signe(z)) return gen_0;
                 [ +  + ]
    1710                 :      56256 :   t = mulll(x,y);
    1711 [ +  - ][ +  - ]:      56256 :   if (!hiremainder)
    1712                 :      56256 :     r = muluispec(t, z+2, lgefint(z)-2);
    1713                 :            :   else
    1714                 :            :   {
    1715                 :            :     long tmp[2];
    1716                 :          0 :     tmp[0] = hiremainder;
    1717                 :          0 :     tmp[1] = t;
    1718                 :          0 :     r = muliispec(z+2,tmp,lgefint(z)-2,2);
    1719                 :            :   }
    1720                 :      56553 :   setsigne(r,s); return r;
    1721                 :            : }
    1722                 :            : 
    1723                 :            : #define sqrispec_mirror sqrispec
    1724                 :            : #define muliispec_mirror muliispec
    1725                 :            : 
    1726                 :            : /* x % (2^n), assuming n >= 0 */
    1727                 :            : GEN
    1728                 :    5384187 : remi2n(GEN x, long n)
    1729                 :            : {
    1730                 :    5384187 :   long hi,l,k,lx,ly, sx = signe(x);
    1731                 :            :   GEN z, xd, zd;
    1732                 :            : 
    1733 [ +  + ][ -  + ]:    5384187 :   if (!sx || !n) return gen_0;
    1734                 :            : 
    1735                 :    5374215 :   k = dvmdsBIL(n, &l);
    1736                 :    5374215 :   lx = lgefint(x);
    1737         [ +  + ]:    5374215 :   if (lx < k+3) return icopy(x);
    1738                 :            : 
    1739                 :    5348514 :   xd = x + (lx-k-1);
    1740                 :            :   /* x = |_|...|#|1|...|k| : copy the last l bits of # and the last k words
    1741                 :            :    *            ^--- initial xd  */
    1742                 :    5348514 :   hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
    1743         [ +  + ]:    5348514 :   if (!hi)
    1744                 :            :   { /* strip leading zeroes from result */
    1745 [ +  + ][ +  + ]:     142719 :     xd++; while (k && !*xd) { k--; xd++; }
    1746         [ +  + ]:     139917 :     if (!k) return gen_0;
    1747                 :      38052 :     ly = k+2; xd--;
    1748                 :            :   }
    1749                 :            :   else
    1750                 :    5208597 :     ly = k+3;
    1751                 :            : 
    1752                 :    5246649 :   zd = z = cgeti(ly);
    1753                 :    5246649 :   *++zd = evalsigne(sx) | evallgefint(ly);
    1754         [ +  + ]:    5246649 :   if (hi) *++zd = hi;
    1755         [ +  + ]:    9321912 :   for ( ;k; k--) *++zd = *++xd;
    1756                 :    5384187 :   return z;
    1757                 :            : }
    1758                 :            : 
    1759                 :            : GEN
    1760                 :  454145103 : sqrispec(GEN a, long na)
    1761                 :            : {
    1762                 :            :   GEN a0,c;
    1763                 :            :   long n0, n0a, i;
    1764                 :            :   pari_sp av;
    1765                 :            : 
    1766         [ +  + ]:  454145103 :   if (na < SQRI_KARATSUBA_LIMIT) return sqrispec_basecase(a,na);
    1767         [ +  + ]:    1320624 :   if (na >= SQRI_FFT_LIMIT) return sqrispec_fft(a,na);
    1768                 :    1320309 :   i=(na>>1); n0=na-i; na=i;
    1769                 :    1320309 :   av=avma; a0=a+na; n0a=n0;
    1770 [ +  + ][ +  + ]:    2622906 :   while (n0a && !*a0) { a0++; n0a--; }
    1771                 :    1320309 :   c = sqrispec(a,na);
    1772         [ +  + ]:    1320309 :   if (n0a)
    1773                 :            :   {
    1774                 :    1308126 :     GEN t, c1, c0 = sqrispec(a0,n0a);
    1775                 :            : #if 0
    1776                 :            :     c1 = shifti(muliispec(a0,a, n0a,na),1);
    1777                 :            : #else /* faster */
    1778                 :    1308126 :     t = addiispec(a0,a,n0a,na);
    1779                 :    1308126 :     t = sqrispec(LIMBS(t),NLIMBS(t));
    1780                 :    1308126 :     c1= addiispec(LIMBS(c0),LIMBS(c), NLIMBS(c0), NLIMBS(c));
    1781                 :    1308126 :     c1= subiispec(LIMBS(t),LIMBS(c1), NLIMBS(t), NLIMBS(c1));
    1782                 :            : #endif
    1783                 :    1308126 :     c = addshiftw(c,c1, n0);
    1784                 :    1308126 :     c = addshiftw(c,c0, n0);
    1785                 :            :   }
    1786                 :            :   else
    1787                 :      12183 :     c = addshiftw(c,gen_0,n0<<1);
    1788                 :  454145103 :   return gerepileuptoint(av, c);
    1789                 :            : }
    1790                 :            : 
    1791                 :            : /********************************************************************/
    1792                 :            : /**                                                                **/
    1793                 :            : /**                    KARATSUBA SQUARE ROOT                       **/
    1794                 :            : /**      adapted from Paul Zimmermann's implementation of          **/
    1795                 :            : /**      his algorithm in GMP (mpn_sqrtrem)                        **/
    1796                 :            : /**                                                                **/
    1797                 :            : /********************************************************************/
    1798                 :            : 
    1799                 :            : /* Square roots table */
    1800                 :            : static const unsigned char approx_tab[192] = {
    1801                 :            :   128,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,
    1802                 :            :   143,144,144,145,146,147,148,149,150,150,151,152,153,154,155,155,
    1803                 :            :   156,157,158,159,160,160,161,162,163,163,164,165,166,167,167,168,
    1804                 :            :   169,170,170,171,172,173,173,174,175,176,176,177,178,178,179,180,
    1805                 :            :   181,181,182,183,183,184,185,185,186,187,187,188,189,189,190,191,
    1806                 :            :   192,192,193,193,194,195,195,196,197,197,198,199,199,200,201,201,
    1807                 :            :   202,203,203,204,204,205,206,206,207,208,208,209,209,210,211,211,
    1808                 :            :   212,212,213,214,214,215,215,216,217,217,218,218,219,219,220,221,
    1809                 :            :   221,222,222,223,224,224,225,225,226,226,227,227,228,229,229,230,
    1810                 :            :   230,231,231,232,232,233,234,234,235,235,236,236,237,237,238,238,
    1811                 :            :   239,240,240,241,241,242,242,243,243,244,244,245,245,246,246,247,
    1812                 :            :   247,248,248,249,249,250,250,251,251,252,252,253,253,254,254,255
    1813                 :            : };
    1814                 :            : 
    1815                 :            : /* N[0], assume N[0] >= 2^(BIL-2).
    1816                 :            :  * Return r,s such that s^2 + r = N, 0 <= r <= 2s */
    1817                 :            : static void
    1818                 :   10356480 : p_sqrtu1(ulong *N, ulong *ps, ulong *pr)
    1819                 :            : {
    1820                 :   10356480 :   ulong prec, r, s, q, u, n0 = N[0];
    1821                 :            : 
    1822                 :   10356480 :   q = n0 >> (BITS_IN_LONG - 8);
    1823                 :            :   /* 2^6 = 64 <= q < 256 = 2^8 */
    1824                 :   10356480 :   s = approx_tab[q - 64];                                /* 128 <= s < 255 */
    1825                 :   10356480 :   r = (n0 >> (BITS_IN_LONG - 16)) - s * s;                /* r <= 2*s */
    1826         [ +  + ]:   10356480 :   if (r > (s << 1)) { r -= (s << 1) | 1; s++; }
    1827                 :            : 
    1828                 :            :   /* 8-bit approximation from the high 8-bits of N[0] */
    1829                 :   10356480 :   prec = 8;
    1830                 :   10356480 :   n0 <<= 2 * prec;
    1831         [ +  + ]:   31069440 :   while (2 * prec < BITS_IN_LONG)
    1832                 :            :   { /* invariant: s has prec bits, and r <= 2*s */
    1833                 :   20712960 :     r = (r << prec) + (n0 >> (BITS_IN_LONG - prec));
    1834                 :   20712960 :     n0 <<= prec;
    1835                 :   20712960 :     u = 2 * s;
    1836                 :   20712960 :     q = r / u; u = r - q * u;
    1837                 :   20712960 :     s = (s << prec) + q;
    1838                 :   20712960 :     u = (u << prec) + (n0 >> (BITS_IN_LONG - prec));
    1839                 :   20712960 :     q = q * q;
    1840                 :   20712960 :     r = u - q;
    1841         [ +  + ]:   20712960 :     if (u < q) { s--; r += (s << 1) | 1; }
    1842                 :   20712960 :     n0 <<= prec;
    1843                 :   20712960 :     prec = 2 * prec;
    1844                 :            :   }
    1845                 :   10356480 :   *ps = s;
    1846                 :   10356480 :   *pr = r;
    1847                 :   10356480 : }
    1848                 :            : 
    1849                 :            : /* N[0..1], assume N[0] >= 2^(BIL-2).
    1850                 :            :  * Return 1 if remainder overflows, 0 otherwise */
    1851                 :            : static int
    1852                 :    9313734 : p_sqrtu2(ulong *N, ulong *ps, ulong *pr)
    1853                 :            : {
    1854                 :    9313734 :   ulong cc, qhl, r, s, q, u, n1 = N[1];
    1855                 :            :   LOCAL_OVERFLOW;
    1856                 :            : 
    1857                 :    9313734 :   p_sqrtu1(N, &s, &r); /* r <= 2s */
    1858         [ +  + ]:   13952583 :   qhl = 0; while (r >= s) { qhl++; r -= s; }
    1859                 :            :   /* now r < s < 2^(BIL/2) */
    1860                 :    9313734 :   r = (r << BITS_IN_HALFULONG) | (n1 >> BITS_IN_HALFULONG);
    1861                 :    9313734 :   u = s << 1;
    1862                 :    9313734 :   q = r / u; u = r - q * u;
    1863                 :    9313734 :   q += (qhl & 1) << (BITS_IN_HALFULONG - 1);
    1864                 :    9313734 :   qhl >>= 1;
    1865                 :            :   /* (initial r)<<(BIL/2) + n1>>(BIL/2) = (qhl<<(BIL/2) + q) * 2s + u */
    1866                 :    9313734 :   s = ((s + qhl) << BITS_IN_HALFULONG) + q;
    1867                 :    9313734 :   cc = u >> BITS_IN_HALFULONG;
    1868                 :    9313734 :   r = (u << BITS_IN_HALFULONG) | (n1 & LOWMASK);
    1869                 :    9313734 :   r = subll(r, q * q);
    1870                 :    9313734 :   cc -= overflow + qhl;
    1871                 :            :   /* now subtract 2*q*2^(BIL/2) + 2^BIL if qhl is set */
    1872 [ +  + ][ +  + ]:    9313734 :   if ((long)cc < 0)
    1873                 :            :   {
    1874         [ +  + ]:    2293773 :     if (s) {
    1875                 :    2286732 :       r = addll(r, s);
    1876                 :    2286732 :       cc += overflow;
    1877                 :    2286732 :       s--;
    1878                 :            :     } else {
    1879                 :       7041 :       cc++;
    1880                 :       7041 :       s = ~0UL;
    1881                 :            :     }
    1882                 :    2293773 :     r = addll(r, s);
    1883                 :    2293773 :     cc += overflow;
    1884                 :            :   }
    1885                 :    9313734 :   *ps = s;
    1886                 :    9313734 :   *pr = r; return cc;
    1887                 :            : }
    1888                 :            : 
    1889                 :            : static void
    1890                 :    9303204 : xmpn_zero(GEN x, long n)
    1891                 :            : {
    1892         [ +  + ]:   68545833 :   while (--n >= 0) x[n]=0;
    1893                 :    9303204 : }
    1894                 :            : static void
    1895                 :  113353290 : xmpn_copy(GEN z, GEN x, long n)
    1896                 :            : {
    1897                 :  113353290 :   long k = n;
    1898         [ +  + ]:  431721386 :   while (--k >= 0) z[k] = x[k];
    1899                 :  113353290 : }
    1900                 :            : /* a[0..la-1] * 2^(lb BIL) | b[0..lb-1] */
    1901                 :            : static GEN
    1902                 :   49596888 : catii(GEN a, long la, GEN b, long lb)
    1903                 :            : {
    1904                 :   49596888 :   long l = la + lb + 2;
    1905                 :   49596888 :   GEN z = cgetipos(l);
    1906                 :   49596888 :   xmpn_copy(LIMBS(z), a, la);
    1907                 :   49596888 :   xmpn_copy(LIMBS(z) + la, b, lb);
    1908                 :   49596888 :   return int_normalize(z, 0);
    1909                 :            : }
    1910                 :            : 
    1911                 :            : /* sqrt n[0..1], assume n normalized */
    1912                 :            : static GEN
    1913                 :    9310737 : sqrtispec2(GEN n, GEN *pr)
    1914                 :            : {
    1915                 :            :   ulong s, r;
    1916                 :    9310737 :   int hi = p_sqrtu2((ulong*)n, &s, &r);
    1917                 :    9310737 :   GEN S = utoi(s);
    1918         [ +  + ]:    9310737 :   *pr = hi? uutoi(1,r): utoi(r);
    1919                 :    9310737 :   return S;
    1920                 :            : }
    1921                 :            : 
    1922                 :            : /* sqrt n[0], _dont_ assume n normalized */
    1923                 :            : static GEN
    1924                 :    1042746 : sqrtispec1_sh(GEN n, GEN *pr)
    1925                 :            : {
    1926                 :            :   GEN S;
    1927                 :    1042746 :   ulong r, s, u0 = uel(n,0);
    1928 [ +  + ][ +  + ]:    1042746 :   int sh = bfffo(u0) & ~1UL;
         [ +  + ][ +  + ]
    1929 [ +  + ][ +  + ]:    1042746 :   if (sh) u0 <<= sh;
    1930                 :    1042746 :   p_sqrtu1(&u0, &s, &r);
    1931                 :            :   /* s^2 + r = u0, s < 2^(BIL/2). Rescale back:
    1932                 :            :    * 2^(2k) n = S^2 + R
    1933                 :            :    * so 2^(2k) n = (S - s0)^2 + (2*S*s0 - s0^2 + R), s0 = S mod 2^k. */
    1934         [ +  + ]:    1042746 :   if (sh) {
    1935                 :    1042704 :     int k = sh >> 1;
    1936                 :    1042704 :     ulong s0 = s & ((1L<<k) - 1);
    1937                 :    1042704 :     r += s * (s0<<1);
    1938                 :    1042704 :     s >>= k;
    1939                 :    1042704 :     r >>= sh;
    1940                 :            :   }
    1941                 :    1042746 :   S = utoi(s);
    1942         [ -  + ]:    1042746 :   if (pr) *pr = utoi(r);
    1943                 :    1042746 :   return S;
    1944                 :            : }
    1945                 :            : 
    1946                 :            : /* sqrt n[0..1], _dont_ assume n normalized */
    1947                 :            : static GEN
    1948                 :       2997 : sqrtispec2_sh(GEN n, GEN *pr)
    1949                 :            : {
    1950                 :            :   GEN S;
    1951                 :       2997 :   ulong U[2], r, s, u0 = uel(n,0), u1 = uel(n,1);
    1952 [ +  + ][ +  + ]:       2997 :   int hi, sh = bfffo(u0) & ~1UL;
         [ +  + ][ +  + ]
    1953 [ +  + ][ +  + ]:       2997 :   if (sh) {
    1954                 :       2919 :     u0 = (u0 << sh) | (u1 >> (BITS_IN_LONG-sh));
    1955                 :       2919 :     u1 <<= sh;
    1956                 :            :   }
    1957                 :       2997 :   U[0] = u0;
    1958                 :       2997 :   U[1] = u1; hi = p_sqrtu2(U, &s, &r);
    1959                 :            :   /* s^2 + R = u0|u1. Rescale back:
    1960                 :            :    * 2^(2k) n = S^2 + R
    1961                 :            :    * so 2^(2k) n = (S - s0)^2 + (2*S*s0 - s0^2 + R), s0 = S mod 2^k. */
    1962         [ +  + ]:       2997 :   if (sh) {
    1963                 :       2919 :     int k = sh >> 1;
    1964                 :       2919 :     ulong s0 = s & ((1L<<k) - 1);
    1965                 :            :     LOCAL_HIREMAINDER;
    1966                 :            :     LOCAL_OVERFLOW;
    1967                 :       2919 :     r = addll(r, mulll(s, (s0<<1)));
    1968 [ +  + ][ +  + ]:       2919 :     if (overflow) hiremainder++;
    1969                 :       2919 :     hiremainder += hi; /* + 0 or 1 */
    1970                 :       2919 :     s >>= k;
    1971                 :       2919 :     r = (r>>sh) | (hiremainder << (BITS_IN_LONG-sh));
    1972                 :       2919 :     hi = (hiremainder & (1L<<sh));
    1973                 :            :   }
    1974                 :       2997 :   S = utoi(s);
    1975 [ +  + ][ -  + ]:       2997 :   if (pr) *pr = hi? uutoi(1,r): utoi(r);
    1976                 :       2997 :   return S;
    1977                 :            : }
    1978                 :            : 
    1979                 :            : /* Let N = N[0..2n-1]. Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S
    1980                 :            :  * Assume N normalized */
    1981                 :            : static GEN
    1982                 :   34109181 : sqrtispec(GEN N, long n, GEN *r)
    1983                 :            : {
    1984                 :            :   GEN S, R, q, z, u;
    1985                 :            :   long l, h;
    1986                 :            : 
    1987         [ +  + ]:   34109181 :   if (n == 1) return sqrtispec2(N, r);
    1988                 :   24798444 :   l = n >> 1;
    1989                 :   24798444 :   h = n - l; /* N = a3(h) | a2(h) | a1(l) | a0(l words) */
    1990                 :   24798444 :   S = sqrtispec(N, h, &R); /* S^2 + R = a3|a2 */
    1991                 :            : 
    1992                 :   24798444 :   z = catii(LIMBS(R), NLIMBS(R), N + 2*h, l); /* = R | a1(l) */
    1993                 :   24798444 :   q = dvmdii(z, shifti(S,1), &u);
    1994                 :   24798444 :   z = catii(LIMBS(u), NLIMBS(u), N + n + h, l); /* = u | a0(l) */
    1995                 :            : 
    1996                 :   24798444 :   S = addshiftw(S, q, l);
    1997                 :   24798444 :   R = subii(z, sqri(q));
    1998         [ +  + ]:   24798444 :   if (signe(R) < 0)
    1999                 :            :   {
    2000                 :    3435018 :     GEN S2 = shifti(S,1);
    2001                 :    3435018 :     R = addis(subiispec(LIMBS(S2),LIMBS(R), NLIMBS(S2),NLIMBS(R)), -1);
    2002                 :    3435018 :     S = addis(S, -1);
    2003                 :            :   }
    2004                 :   34109181 :   *r = R; return S;
    2005                 :            : }
    2006                 :            : 
    2007                 :            : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
    2008                 :            :  * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
    2009                 :            :  * remainder is 0. R = NULL is allowed. */
    2010                 :            : GEN
    2011                 :    1053276 : sqrtremi(GEN N, GEN *r)
    2012                 :            : {
    2013                 :            :   pari_sp av;
    2014                 :    1053276 :   GEN S, R, n = N+2;
    2015                 :    1053276 :   long k, l2, ln = NLIMBS(N);
    2016                 :            :   int sh;
    2017                 :            : 
    2018         [ +  + ]:    1053276 :   if (ln <= 2)
    2019                 :            :   {
    2020         [ +  + ]:    1045743 :     if (ln == 2) return sqrtispec2_sh(n, r);
    2021         [ +  - ]:    1042746 :     if (ln == 1) return sqrtispec1_sh(n, r);
    2022         [ #  # ]:          0 :     if (r) *r = gen_0;
    2023                 :          0 :     return gen_0;
    2024                 :            :   }
    2025                 :       7533 :   av = avma;
    2026 [ +  + ][ +  + ]:       7533 :   sh = bfffo(n[0]) >> 1;
         [ +  + ][ +  + ]
    2027                 :       7533 :   l2 = (ln + 1) >> 1;
    2028 [ +  + ][ +  + ]:      15003 :   if (sh || (ln & 1)) { /* normalize n, so that n[0] >= 2^BIL / 4 */
         [ +  + ][ +  + ]
    2029                 :       7470 :     GEN s0, t = new_chunk(ln + 1);
    2030                 :       7470 :     t[ln] = 0;
    2031         [ +  + ]:       7470 :     if (sh)
    2032                 :       7317 :       shift_left(t, n, 0,ln-1, 0, sh << 1);
    2033                 :            :     else
    2034                 :        153 :       xmpn_copy(t, n, ln);
    2035                 :       7470 :     S = sqrtispec(t, l2, &R); /* t normalized, 2 * l2 words */
    2036                 :            :     /* Rescale back:
    2037                 :            :      * 2^(2k) n = S^2 + R, k = sh + (ln & 1)*BIL/2
    2038                 :            :      * so 2^(2k) n = (S - s0)^2 + (2*S*s0 - s0^2 + R), s0 = S mod 2^k. */
    2039                 :       7470 :     k = sh + (ln & 1) * (BITS_IN_LONG/2);
    2040                 :       7470 :     s0 = remi2n(S, k);
    2041                 :       7470 :     R = addii(shifti(R,-1), mulii(s0, S));
    2042                 :       7470 :     R = shifti(R, 1 - (k<<1));
    2043                 :       7470 :     S = shifti(S, -k);
    2044                 :            :   }
    2045                 :            :   else
    2046                 :         63 :     S = sqrtispec(n, l2, &R);
    2047                 :            : 
    2048         [ +  + ]:       7533 :   if (!r) { avma = (pari_sp)S; return gerepileuptoint(av, S); }
    2049                 :    1053276 :   gerepileall(av, 2, &S, &R); *r = R; return S;
    2050                 :            : }
    2051                 :            : 
    2052                 :            : /* compute sqrt(|a|), assuming a != 0 */
    2053                 :            : 
    2054                 :            : #if 1
    2055                 :            : GEN
    2056                 :    9303204 : sqrtr_abs(GEN x)
    2057                 :            : {
    2058                 :    9303204 :   long l = lg(x) - 2, e = expo(x), er = e>>1;
    2059                 :    9303204 :   GEN b, c, res = cgetr(2 + l);
    2060                 :    9303204 :   res[1] = evalsigne(1) | evalexpo(er);
    2061         [ +  + ]:    9303204 :   if (e&1) {
    2062                 :    4856157 :     b = new_chunk(l << 1);
    2063                 :    4856157 :     xmpn_copy(b, x+2, l);
    2064                 :    4856157 :     xmpn_zero(b + l,l);
    2065                 :    4856157 :     b = sqrtispec(b, l, &c);
    2066                 :    4856157 :     xmpn_copy(res+2, b+2, l);
    2067         [ +  + ]:    4856157 :     if (cmpii(c, b) > 0) roundr_up_ip(res, l+2);
    2068                 :            :   } else {
    2069                 :            :     ulong u;
    2070                 :    4447047 :     b = new_chunk(2 + (l << 1));
    2071                 :    4447047 :     shift_left(b+1, x+2, 0,l-1, 0, BITS_IN_LONG-1);
    2072                 :    4447047 :     b[0] = uel(x,2)>>1;
    2073                 :    4447047 :     xmpn_zero(b + l+1,l+1);
    2074                 :    4447047 :     b = sqrtispec(b, l+1, &c);
    2075                 :    4447047 :     xmpn_copy(res+2, b+2, l);
    2076                 :    4447047 :     u = uel(b,l+2);
    2077 [ +  + ][ +  + ]:    4447047 :     if ( u&HIGHBIT || (u == ~HIGHBIT && cmpii(c,b) > 0))
                 [ +  - ]
    2078                 :    2186638 :       roundr_up_ip(res, l+2);
    2079                 :            :   }
    2080                 :    9303204 :   avma = (pari_sp)res; return res;
    2081                 :            : }
    2082                 :            : 
    2083                 :            : #else /* use t_REAL: currently much slower (quadratic division) */
    2084                 :            : 
    2085                 :            : #ifdef LONG_IS_64BIT
    2086                 :            : /* 64 bits of b = sqrt(a[0] * 2^64 + a[1])  [ up to 1ulp ] */
    2087                 :            : static ulong
    2088                 :            : sqrtu2(ulong *a)
    2089                 :            : {
    2090                 :            :   ulong c, b = dblmantissa( sqrt((double)a[0]) );
    2091                 :            :   LOCAL_HIREMAINDER;
    2092                 :            :   LOCAL_OVERFLOW;
    2093                 :            : 
    2094                 :            :   /* > 32 correct bits, 1 Newton iteration to reach 64 */
    2095                 :            :   if (b <= a[0]) return HIGHBIT | (a[0] >> 1);
    2096                 :            :   hiremainder = a[0]; c = divll(a[1], b);
    2097                 :            :   return (addll(c, b) >> 1) | HIGHBIT;
    2098                 :            : }
    2099                 :            : /* 64 bits of sqrt(a[0] * 2^63) */
    2100                 :            : static ulong
    2101                 :            : sqrtu2_1(ulong *a)
    2102                 :            : {
    2103                 :            :   ulong t[2];
    2104                 :            :   t[0] = (a[0] >> 1);
    2105                 :            :   t[1] = (a[0] << (BITS_IN_LONG-1)) | (a[1] >> 1);
    2106                 :            :   return sqrtu2(t);
    2107                 :            : }
    2108                 :            : #else
    2109                 :            : /* 32 bits of sqrt(a[0] * 2^32) */
    2110                 :            : static ulong
    2111                 :            : sqrtu2(ulong *a)   { return dblmantissa( sqrt((double)a[0]) ); }
    2112                 :            : /* 32 bits of sqrt(a[0] * 2^31) */
    2113                 :            : static ulong
    2114                 :            : sqrtu2_1(ulong *a) { return dblmantissa( sqrt(2. * a[0]) ); }
    2115                 :            : #endif
    2116                 :            : 
    2117                 :            : GEN
    2118                 :            : sqrtr_abs(GEN x)
    2119                 :            : {
    2120                 :            :   long l1, i, l = lg(x), ex = expo(x);
    2121                 :            :   GEN a, t, y = cgetr(l);
    2122                 :            :   pari_sp av, av0 = avma;
    2123                 :            : 
    2124                 :            :   a = rtor(x, l+1);
    2125                 :            :   t = cgetr(l+1);
    2126                 :            :   if (ex & 1) { /* odd exponent */
    2127                 :            :     a[1] = evalsigne(1) | _evalexpo(1);
    2128                 :            :     t[2] = (long)sqrtu2((ulong*)a + 2);
    2129                 :            :   } else { /* even exponent */
    2130                 :            :     a[1] = evalsigne(1) | _evalexpo(0);
    2131                 :            :     t[2] = (long)sqrtu2_1((ulong*)a + 2);
    2132                 :            :   }
    2133                 :            :   t[1] = evalsigne(1) | _evalexpo(0);
    2134                 :            :   for (i = 3; i <= l; i++) t[i] = 0;
    2135                 :            : 
    2136                 :            :   /* |x| = 2^(ex/2) a, t ~ sqrt(a) */
    2137                 :            :   l--; l1 = 1; av = avma;
    2138                 :            :   while (l1 < l) { /* let t := (t + a/t)/2 */
    2139                 :            :     l1 <<= 1; if (l1 > l) l1 = l;
    2140                 :            :     setlg(a, l1 + 2);
    2141                 :            :     setlg(t, l1 + 2);
    2142                 :            :     affrr(addrr(t, divrr(a,t)), t); shiftr_inplace(t, -1);
    2143                 :            :     avma = av;
    2144                 :            :   }
    2145                 :            :   affrr(t,y); shiftr_inplace(y, (ex>>1));
    2146                 :            :   avma = av0; return y;
    2147                 :            : }
    2148                 :            : 
    2149                 :            : #endif
    2150                 :            : 
    2151                 :            : /*******************************************************************
    2152                 :            :  *                                                                 *
    2153                 :            :  *                           Base Conversion                       *
    2154                 :            :  *                                                                 *
    2155                 :            :  *******************************************************************/
    2156                 :            : 
    2157                 :            : static void
    2158                 :    1223598 : convi_dac(GEN x, ulong l, ulong *res)
    2159                 :            : {
    2160                 :    1223598 :   pari_sp ltop=avma;
    2161                 :            :   ulong m;
    2162                 :            :   GEN x1,x2;
    2163         [ +  + ]:    1773060 :   if (l==1) { *res=itou(x); return; }
    2164                 :     549462 :   m=l>>1;
    2165                 :     549462 :   x1=dvmdii(x,powuu(1000000000UL,m),&x2);
    2166                 :     549462 :   convi_dac(x1,l-m,res+m);
    2167                 :     549462 :   convi_dac(x2,m,res);
    2168                 :     549462 :   avma=ltop;
    2169                 :            : }
    2170                 :            : 
    2171                 :            : /* convert integer --> base 10^9 [not memory clean] */
    2172                 :            : ulong *
    2173                 :     261669 : convi(GEN x, long *l)
    2174                 :            : {
    2175                 :     261669 :   long lz, lx = lgefint(x);
    2176                 :            :   ulong *z;
    2177 [ +  + ][ +  + ]:     261669 :   if (lx == 3 && uel(x,2) < 1000000000UL) {
    2178                 :     136995 :     z = (ulong*)new_chunk(1);
    2179                 :     136995 :     *z = x[2];
    2180                 :     136995 :     *l = 1; return z+1;
    2181                 :            :   }
    2182                 :     124674 :   lz = 1 + (long)bit_accuracy_mul(lx, LOG10_2/9);
    2183                 :     124674 :   z = (ulong*)new_chunk(lz);
    2184                 :     124674 :   convi_dac(x,(ulong)lz,z);
    2185         [ +  + ]:     268014 :   while (z[lz-1]==0) lz--;
    2186                 :     261669 :   *l=lz; return z+lz;
    2187                 :            : }
    2188                 :            : 

Generated by: LCOV version 1.9