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 17953-c39f2e6) Lines: 1108 1148 96.5 %
Date: 2015-08-29 Functions: 69 69 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 925 1110 83.3 %

           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                 :        543 : void pari_kernel_init(void) { }
      25                 :        543 : 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                 :  200791800 : int_normalize(GEN x, long known_zero_words)
      38                 :            : {
      39                 :  200791800 :   long i, lx = lgefint(x);
      40                 :            :   GEN x0;
      41         [ -  + ]:  200791800 :   if (lx == 2) { x[1] = evalsigne(0) | evallgefint(2); return x; }
      42 [ +  + ][ +  + ]:  200791800 :   if (!known_zero_words && x[2]) return x;
      43         [ +  + ]: 1006392633 :   for (i = 2+known_zero_words; i < lx; i++)
      44         [ +  + ]:  952161207 :     if (x[i]) break;
      45                 :  145720746 :   x0 = x; i -= 2; x += i;
      46         [ +  + ]:  145720746 :   if (x0 == (GEN)avma) avma = (pari_sp)x;
      47                 :   74147088 :   else stackdummy((pari_sp)(x0+i), (pari_sp)x0);
      48                 :  145720746 :   lx -= i;
      49                 :  145720746 :   x[0] = evaltyp(t_INT) | evallg(lx);
      50         [ +  + ]:  145720746 :   if (lx == 2) x[1] = evalsigne(0) | evallgefint(lx);
      51                 :   91489320 :   else         x[1] = evalsigne(1) | evallgefint(lx);
      52                 :  200791800 :   return x;
      53                 :            : }
      54                 :            : 
      55                 :            : /***********************************************************************/
      56                 :            : /**                                                                   **/
      57                 :            : /**                      ADDITION / SUBTRACTION                       **/
      58                 :            : /**                                                                   **/
      59                 :            : /***********************************************************************/
      60                 :            : 
      61                 :            : GEN
      62                 :    1631913 : setloop(GEN a)
      63                 :            : {
      64                 :    1631913 :   pari_sp av = avma;
      65                 :    1631913 :   (void)cgetg(lgefint(a) + 3, t_VECSMALL);
      66                 :    1631913 :   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                 :     109206 : resetloop(GEN a, GEN b) {
      72                 :     109206 :   long lb = lgefint(b);
      73                 :     109206 :   a += lgefint(a) - lb;
      74                 :     109206 :   a[0] = evaltyp(t_INT) | evallg(lb);
      75                 :     109206 :   affii(b, a); return a;
      76                 :            : }
      77                 :            : 
      78                 :            : /* assume a > 0, initialized by setloop. Do a++ */
      79                 :            : static GEN
      80                 :   17959638 : incpos(GEN a)
      81                 :            : {
      82                 :   17959638 :   long i, l = lgefint(a);
      83         [ +  + ]:   17959641 :   for (i=l-1; i>1; i--)
      84         [ +  + ]:   17959638 :     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                 :   17959638 :   a[2]=1; return a;
      89                 :            : }
      90                 :            : 
      91                 :            : /* assume a < 0, initialized by setloop. Do a++ */
      92                 :            : static GEN
      93                 :       4638 : incneg(GEN a)
      94                 :            : {
      95                 :       4638 :   long i, l = lgefint(a)-1;
      96         [ +  + ]:       4638 :   if (a[l]--)
      97                 :            :   {
      98 [ +  + ][ +  + ]:       4635 :     if (l == 2 && !a[2])
      99                 :            :     {
     100                 :        399 :       a++; /* save one cell */
     101                 :        399 :       a[0] = evaltyp(t_INT) | _evallg(2);
     102                 :        399 :       a[1] = evalsigne(0) | evallgefint(2);
     103                 :            :     }
     104                 :       4635 :     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                 :       4638 :   return a;
     115                 :            : }
     116                 :            : 
     117                 :            : /* assume a initialized by setloop. Do a++ */
     118                 :            : GEN
     119                 :   17971050 : incloop(GEN a)
     120                 :            : {
     121      [ +  +  + ]:   17971050 :   switch(signe(a))
     122                 :            :   {
     123                 :       6774 :     case 0: a--; /* use extra cell */
     124                 :       6774 :       a[0]=evaltyp(t_INT) | _evallg(3);
     125                 :       6774 :       a[1]=evalsigne(1) | evallgefint(3);
     126                 :       6774 :       a[2]=1; return a;
     127                 :       4638 :     case -1: return incneg(a);
     128                 :   17971050 :     default: return incpos(a);
     129                 :            :   }
     130                 :            : }
     131                 :            : 
     132                 :            : INLINE GEN
     133                 :  885403038 : adduispec(ulong s, GEN x, long nx)
     134                 :            : {
     135                 :  885403038 :   GEN xd, zd = (GEN)avma;
     136                 :            :   long lz;
     137                 :            : 
     138         [ +  + ]:  885403038 :   if (nx == 1) return adduu(s, uel(x,0));
     139                 :  189369564 :   lz = nx+3; (void)new_chunk(lz);
     140                 :  189369564 :   xd = x + nx;
     141                 :  189369564 :   *--zd = (ulong)*--xd + s;
     142         [ +  + ]:  189369564 :   if ((ulong)*zd < s)
     143                 :            :     for(;;)
     144                 :            :     {
     145         [ +  + ]:   34018362 :       if (xd == x) { *--zd = 1; break; } /* enlarge z */
     146                 :   30317151 :       *--zd = ((ulong)*--xd) + 1;
     147         [ +  + ]:   30317151 :       if (*zd) { lz--; break; }
     148                 :   34018362 :     }
     149                 :  179571630 :   else lz--;
     150         [ +  + ]:  577251846 :   while (xd > x) *--zd = *--xd;
     151                 :  189369564 :   *--zd = evalsigne(1) | evallgefint(lz);
     152                 :  189369564 :   *--zd = evaltyp(t_INT) | evallg(lz);
     153                 :  885403038 :   avma=(pari_sp)zd; return zd;
     154                 :            : }
     155                 :            : 
     156                 :            : GEN
     157                 :  218225145 : adduispec_offset(ulong s, GEN x, long offset, long nx)
     158                 :            : {
     159                 :  218225145 :   GEN xd = x+lgefint(x)-nx-offset;
     160 [ +  + ][ +  + ]:  293665986 :   while (nx && *xd==0) {xd++; nx--;}
     161         [ +  + ]:  218225145 :   if (!nx) return utoi(s);
     162                 :  218225145 :   return adduispec(s,xd,nx);
     163                 :            : }
     164                 :            : 
     165                 :            : static GEN
     166                 : 1729826766 : addiispec(GEN x, GEN y, long nx, long ny)
     167                 :            : {
     168                 :            :   GEN xd, yd, zd;
     169                 : 1729826766 :   long lz, i = -2;
     170                 :            :   LOCAL_OVERFLOW;
     171                 :            : 
     172         [ +  + ]: 1729826766 :   if (nx < ny) swapspec(x,y, nx,ny);
     173         [ +  + ]: 1729826766 :   if (ny == 1) return adduispec(*y,x,nx);
     174                 : 1083984819 :   zd = (GEN)avma;
     175                 : 1083984819 :   lz = nx+3; (void)new_chunk(lz);
     176                 : 1083984819 :   xd = x + nx;
     177                 : 1083984819 :   yd = y + ny;
     178                 : 1083984819 :   zd[-1] = addll(xd[-1], yd[-1]);
     179                 :            : #ifdef addllx8
     180         [ +  + ]:  710455806 :   for (  ; i-8 > -ny; i-=8)
     181                 :  349127533 :     addllx8(xd+i, yd+i, zd+i, overflow);
     182                 :            : #endif
     183         [ +  + ]:11177585671 :   for (  ; i >= -ny; i--) zd[i] = addllx(xd[i], yd[i]);
              [ +  +  + ]
     184         [ +  + ]: 1083984819 :   if (overflow)
     185                 :            :     for(;;)
     186                 :            :     {
     187         [ +  + ]:  190683009 :       if (i < -nx) { zd[i] = 1; i--; break; } /* enlarge z */
     188                 :  136309293 :       zd[i] = uel(xd,i) + 1;
     189         [ +  + ]:  136309293 :       if (zd[i]) { i--; lz--; break; }
     190                 :   68657892 :       i--;
     191                 :  190683009 :     }
     192                 :  961959702 :   else lz--;
     193         [ +  + ]: 5239419399 :   for (; i >= -nx; i--) zd[i] = xd[i];
     194                 : 1083984819 :   zd += i+1;
     195                 : 1083984819 :   *--zd = evalsigne(1) | evallgefint(lz);
     196                 : 1083984819 :   *--zd = evaltyp(t_INT) | evallg(lz);
     197                 : 1729826766 :   avma=(pari_sp)zd; return zd;
     198                 :            : }
     199                 :            : 
     200                 :            : /* assume x >= s */
     201                 :            : INLINE GEN
     202                 :  852182314 : subiuspec(GEN x, ulong s, long nx)
     203                 :            : {
     204                 :  852182314 :   GEN xd, zd = (GEN)avma;
     205                 :            :   long lz;
     206                 :            :   LOCAL_OVERFLOW;
     207                 :            : 
     208         [ +  + ]:  852182314 :   if (nx == 1) return utoi(x[0] - s);
     209                 :            : 
     210                 :   87159640 :   lz = nx+2; (void)new_chunk(lz);
     211                 :   87159640 :   xd = x + nx;
     212                 :   87159640 :   *--zd = subll(*--xd, s);
     213         [ +  + ]:   87159640 :   if (overflow)
     214                 :            :     for(;;)
     215                 :            :     {
     216                 :   41825169 :       *--zd = ((ulong)*--xd) - 1;
     217         [ +  + ]:   41825169 :       if (*xd) break;
     218                 :    3086118 :     }
     219         [ +  + ]:   87159640 :   if (xd == x)
     220         [ +  + ]:   63738006 :     while (*zd == 0) { zd++; lz--; } /* shorten z */
     221                 :            :   else
     222         [ +  + ]:  894339764 :     do  *--zd = *--xd; while (xd > x);
     223                 :   87159640 :   *--zd = evalsigne(1) | evallgefint(lz);
     224                 :   87159640 :   *--zd = evaltyp(t_INT) | evallg(lz);
     225                 :  852182314 :   avma=(pari_sp)zd; return zd;
     226                 :            : }
     227                 :            : 
     228                 :            : /* assume x > y */
     229                 :            : static GEN
     230                 : 1584477744 : subiispec(GEN x, GEN y, long nx, long ny)
     231                 :            : {
     232                 :            :   GEN xd,yd,zd;
     233                 : 1584477744 :   long lz, i = -2;
     234                 :            :   LOCAL_OVERFLOW;
     235                 :            : 
     236         [ +  + ]: 1584477744 :   if (ny==1) return subiuspec(x,*y,nx);
     237                 :  754778339 :   zd = (GEN)avma;
     238                 :  754778339 :   lz = nx+2; (void)new_chunk(lz);
     239                 :  754778339 :   xd = x + nx;
     240                 :  754778339 :   yd = y + ny;
     241                 :  754778339 :   zd[-1] = subll(xd[-1], yd[-1]);
     242                 :            : #ifdef subllx8
     243         [ +  + ]:  575431187 :   for (  ; i-8 > -ny; i-=8)
     244                 :  323838407 :     subllx8(xd+i, yd+i, zd+i, overflow);
     245                 :            : #endif
     246         [ +  + ]: 8865665794 :   for (  ; i >= -ny; i--) zd[i] = subllx(xd[i], yd[i]);
              [ +  +  + ]
     247         [ +  + ]:  754778339 :   if (overflow)
     248                 :            :     for(;;)
     249                 :            :     {
     250                 :  189272081 :       zd[i] = uel(xd,i) - 1;
     251         [ +  + ]:  189272081 :       if (xd[i--]) break;
     252                 :  112755378 :     }
     253         [ +  + ]:  754778339 :   if (i>=-nx)
     254         [ +  + ]:  881570764 :     for (; i >= -nx; i--) zd[i] = xd[i];
     255                 :            :   else
     256         [ +  + ]:  856907978 :     while (zd[i+1] == 0) { i++; lz--; } /* shorten z */
     257                 :  754778339 :   zd += i+1;
     258                 :  754778339 :   *--zd = evalsigne(1) | evallgefint(lz);
     259                 :  754778339 :   *--zd = evaltyp(t_INT) | evallg(lz);
     260                 : 1584477744 :   avma=(pari_sp)zd; return zd;
     261                 :            : }
     262                 :            : 
     263                 :            : static void
     264                 :  113064383 : roundr_up_ip(GEN x, long l)
     265                 :            : {
     266                 :  113064383 :   long i = l;
     267                 :            :   for(;;)
     268                 :            :   {
     269         [ +  + ]:  113158823 :     if (++uel(x,--i)) break;
     270         [ +  + ]:     126396 :     if (i == 2) { x[2] = (long)HIGHBIT; shiftr_inplace(x, 1); break; }
     271                 :      94440 :   }
     272                 :  113064383 : }
     273                 :            : 
     274                 :            : void
     275                 :   96752457 : affir(GEN x, GEN y)
     276                 :            : {
     277                 :   96752457 :   const long s = signe(x), ly = lg(y);
     278                 :            :   long lx, sh, i;
     279                 :            : 
     280         [ +  + ]:   96752457 :   if (!s)
     281                 :            :   {
     282                 :    1021959 :     y[1] = evalexpo(-bit_accuracy(ly));
     283                 :    1021959 :     return;
     284                 :            :   }
     285                 :            : 
     286 [ +  + ][ +  + ]:   95730498 :   lx = lgefint(x); sh = bfffo(x[2]);
         [ +  + ][ +  + ]
     287                 :   95730498 :   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
     288         [ +  + ]:   95730498 :   if (sh) {
     289         [ +  + ]:   94896585 :     if (lx <= ly)
     290                 :            :     {
     291         [ +  + ]:  190692588 :       for (i=lx; i<ly; i++) y[i]=0;
     292                 :   64581756 :       shift_left(y,x,2,lx-1, 0,sh);
     293                 :   64581756 :       return;
     294                 :            :     }
     295                 :   30314829 :     shift_left(y,x,2,ly-1, x[ly],sh);
     296                 :            :     /* lx > ly: round properly */
     297         [ +  + ]:   30314829 :     if ((x[ly]<<sh) & HIGHBIT) roundr_up_ip(y, ly);
     298                 :            :   }
     299                 :            :   else {
     300         [ +  + ]:     833913 :     if (lx <= ly)
     301                 :            :     {
     302         [ +  + ]:    1476951 :       for (i=2; i<lx; i++) y[i]=x[i];
     303         [ +  + ]:    1362264 :       for (   ; i<ly; i++) y[i]=0;
     304                 :     515412 :       return;
     305                 :            :     }
     306         [ +  + ]:    1044888 :     for (i=2; i<ly; i++) y[i]=x[i];
     307                 :            :     /* lx > ly: round properly */
     308         [ +  + ]:   96752457 :     if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
     309                 :            :   }
     310                 :            : }
     311                 :            : 
     312                 :            : INLINE GEN
     313                 :  317609004 : shiftispec(GEN x, long nx, long n)
     314                 :            : {
     315                 :            :   long ny, i, m;
     316                 :            :   GEN y, yd;
     317         [ +  + ]:  317609004 :   if (!n)  return icopyspec(x, nx);
     318                 :            : 
     319         [ +  + ]:  295808223 :   if (n > 0)
     320                 :            :   {
     321                 :  191604195 :     GEN z = (GEN)avma;
     322                 :  191604195 :     long d = dvmdsBIL(n, &m);
     323                 :            : 
     324                 :  191604195 :     ny = nx+d; y = new_chunk(ny + 2); yd = y + 2;
     325         [ +  + ]: 1212217089 :     for ( ; d; d--) *--z = 0;
     326 [ +  + ][ +  + ]:  918785895 :     if (!m) for (i=0; i<nx; i++) yd[i]=x[i];
     327                 :            :     else
     328                 :            :     {
     329                 :  179550939 :       register const ulong sh = BITS_IN_LONG - m;
     330                 :  179550939 :       shift_left(yd,x, 0,nx-1, 0,m);
     331                 :  179550939 :       i = uel(x,0) >> sh;
     332                 :            :       /* Extend y on the left? */
     333         [ +  + ]:  179550939 :       if (i) { ny++; y = new_chunk(1); y[2] = i; }
     334                 :            :     }
     335                 :            :   }
     336                 :            :   else
     337                 :            :   {
     338                 :  104204028 :     ny = nx - dvmdsBIL(-n, &m);
     339         [ +  + ]:  104204028 :     if (ny<1) return gen_0;
     340                 :  104015472 :     y = new_chunk(ny + 2); yd = y + 2;
     341         [ +  + ]:  104015472 :     if (m) {
     342                 :   74480001 :       shift_right(yd,x, 0,ny, 0,m);
     343         [ +  + ]:   74480001 :       if (yd[0] == 0)
     344                 :            :       {
     345         [ +  + ]:    6763239 :         if (ny==1) { avma = (pari_sp)(y+3); return gen_0; }
     346                 :    6141765 :         ny--; avma = (pari_sp)(++y);
     347                 :            :       }
     348                 :            :     } else {
     349         [ +  + ]: 1216960917 :       for (i=0; i<ny; i++) yd[i]=x[i];
     350                 :            :     }
     351                 :            :   }
     352                 :  294998193 :   y[1] = evalsigne(1)|evallgefint(ny + 2);
     353                 :  317609004 :   y[0] = evaltyp(t_INT)|evallg(ny + 2); return y;
     354                 :            : }
     355                 :            : 
     356                 :            : GEN
     357                 :    8879316 : mantissa2nr(GEN x, long n)
     358                 :            : { /*This is a kludge since x is not an integer*/
     359                 :    8879316 :   long s = signe(x);
     360                 :            :   GEN y;
     361                 :            : 
     362         [ +  + ]:    8879316 :   if(s == 0) return gen_0;
     363                 :    8878443 :   y = shiftispec(x + 2, lg(x) - 2, n);
     364         [ +  - ]:    8878443 :   if (signe(y)) setsigne(y, s);
     365                 :    8879316 :   return y;
     366                 :            : }
     367                 :            : 
     368                 :            : GEN
     369                 :     493479 : truncr(GEN x)
     370                 :            : {
     371                 :            :   long d,e,i,s,m;
     372                 :            :   GEN y;
     373                 :            : 
     374 [ +  + ][ +  + ]:     493479 :   if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
     375                 :     382398 :   d = nbits2prec(e+1); m = remsBIL(e);
     376         [ -  + ]:     382398 :   if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
     377                 :            : 
     378                 :     382398 :   y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
     379         [ +  + ]:     382398 :   if (++m == BITS_IN_LONG)
     380         [ +  + ]:          6 :     for (i=2; i<d; i++) y[i]=x[i];
     381                 :            :   else
     382                 :     382395 :     shift_right(y,x, 2,d,0, BITS_IN_LONG - m);
     383                 :     493479 :   return y;
     384                 :            : }
     385                 :            : 
     386                 :            : /* integral part */
     387                 :            : GEN
     388                 :     268659 : floorr(GEN x)
     389                 :            : {
     390                 :            :   long d,e,i,lx,m;
     391                 :            :   GEN y;
     392                 :            : 
     393         [ +  + ]:     268659 :   if (signe(x) >= 0) return truncr(x);
     394         [ +  + ]:     109311 :   if ((e=expo(x)) < 0) return gen_m1;
     395                 :      46344 :   d = nbits2prec(e+1); m = remsBIL(e);
     396         [ -  + ]:      46344 :   lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
     397                 :      46344 :   y = new_chunk(d);
     398         [ -  + ]:      46344 :   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                 :      46344 :     shift_right(y,x, 2,d,0, BITS_IN_LONG - m);
     407         [ +  + ]:      46344 :     if (x[d-1]<<m == 0)
     408                 :            :     {
     409 [ +  + ][ +  + ]:      45108 :       i=d; while (i<lx && !x[i]) i++;
     410         [ +  + ]:      11694 :       if (i==lx) goto END;
     411                 :            :     }
     412                 :            :   }
     413                 :            :   /* set y:=y+1 */
     414 [ +  - ][ +  - ]:      40674 :   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                 :      46344 :   y[1] = evalsigne(-1) | evallgefint(d);
     418                 :     268659 :   y[0] = evaltyp(t_INT) | evallg(d); return y;
     419                 :            : }
     420                 :            : 
     421                 :            : INLINE int
     422                 : 2015071389 : cmpiispec(GEN x, GEN y, long lx, long ly)
     423                 :            : {
     424                 :            :   long i;
     425         [ +  + ]: 2015071389 :   if (lx < ly) return -1;
     426         [ +  + ]: 1928571792 :   if (lx > ly) return  1;
     427 [ +  + ][ +  + ]: 1943294737 :   i = 0; while (i<lx && x[i]==y[i]) i++;
     428         [ +  + ]: 1806506467 :   if (i==lx) return 0;
     429         [ +  + ]: 2015071389 :   return (uel(x,i) > uel(y,i))? 1: -1;
     430                 :            : }
     431                 :            : 
     432                 :            : INLINE int
     433                 :   23256534 : equaliispec(GEN x, GEN y, long lx, long ly)
     434                 :            : {
     435                 :            :   long i;
     436         [ +  + ]:   23256534 :   if (lx != ly) return 0;
     437 [ +  + ][ +  + ]:   38309883 :   i = ly-1; while (i>=0 && x[i]==y[i]) i--;
     438                 :   23256534 :   return i < 0;
     439                 :            : }
     440                 :            : 
     441                 :            : /***********************************************************************/
     442                 :            : /**                                                                   **/
     443                 :            : /**                          MULTIPLICATION                           **/
     444                 :            : /**                                                                   **/
     445                 :            : /***********************************************************************/
     446                 :            : /* assume ny > 0 */
     447                 :            : INLINE GEN
     448                 : 1180869222 : muluispec(ulong x, GEN y, long ny)
     449                 :            : {
     450                 : 1180869222 :   GEN yd, z = (GEN)avma;
     451                 : 1180869222 :   long lz = ny+3;
     452                 :            :   LOCAL_HIREMAINDER;
     453                 :            : 
     454                 : 1180869222 :   (void)new_chunk(lz);
     455                 : 1180869222 :   yd = y + ny; *--z = mulll(x, *--yd);
     456      [ +  +  + ]: 2304725676 :   while (yd > y) *--z = addmul(x,*--yd);
     457         [ +  + ]: 1180869222 :   if (hiremainder) *--z = hiremainder; else lz--;
     458                 : 1180869222 :   *--z = evalsigne(1) | evallgefint(lz);
     459                 : 1180869222 :   *--z = evaltyp(t_INT) | evallg(lz);
     460                 : 1180869222 :   avma=(pari_sp)z; return z;
     461                 :            : }
     462                 :            : 
     463                 :            : /* a + b*|Y| */
     464                 :            : GEN
     465                 :     543180 : 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         [ +  + ]:     543180 :   if (!signe(Y)) return utoi(a);
     473                 :            : 
     474                 :     542760 :   y = LIMBS(Y); z = (GEN)avma;
     475                 :     542760 :   ny = NLIMBS(Y);
     476                 :     542760 :   lz = ny+3;
     477                 :            : 
     478                 :     542760 :   (void)new_chunk(lz);
     479                 :     542760 :   yd = y + ny; *--z = addll(a, mulll(b, *--yd));
     480         [ +  + ]:     542760 :   if (overflow) hiremainder++; /* can't overflow */
     481      [ +  +  + ]:    1180626 :   while (yd > y) *--z = addmul(b,*--yd);
     482         [ +  + ]:     542760 :   if (hiremainder) *--z = hiremainder; else lz--;
     483                 :     542760 :   *--z = evalsigne(1) | evallgefint(lz);
     484                 :     542760 :   *--z = evaltyp(t_INT) | evallg(lz);
     485                 :     543180 :   avma=(pari_sp)z; return z;
     486                 :            : }
     487                 :            : 
     488                 :            : /***********************************************************************/
     489                 :            : /**                                                                   **/
     490                 :            : /**                          DIVISION                                 **/
     491                 :            : /**                                                                   **/
     492                 :            : /***********************************************************************/
     493                 :            : 
     494                 :            : ulong
     495                 :  779245074 : umodiu(GEN y, ulong x)
     496                 :            : {
     497                 :  779245074 :   long sy=signe(y),ly,i;
     498                 :            :   ulong xi;
     499                 :            :   LOCAL_HIREMAINDER;
     500                 :            : 
     501         [ -  + ]:  779245074 :   if (!x) pari_err_INV("umodiu",gen_0);
     502         [ +  + ]:  779245074 :   if (!sy) return 0;
     503                 :  648445794 :   ly = lgefint(y);
     504         [ +  + ]:  648445794 :   if (x <= uel(y,2))
     505                 :            :   {
     506                 :  180859005 :     hiremainder=0;
     507         [ +  + ]:  180859005 :     if (ly==3)
     508                 :            :     {
     509                 :  167209158 :       hiremainder=uel(y,2)%x;
     510         [ +  + ]:  167209158 :       if (!hiremainder) return 0;
     511         [ +  + ]:  146764155 :       return (sy > 0)? hiremainder: x - hiremainder;
     512                 :            :     }
     513                 :            :   }
     514                 :            :   else
     515                 :            :   {
     516 [ +  + ][ +  + ]:  467586789 :     if (ly==3) return (sy > 0)? uel(y,2): x - uel(y,2);
     517                 :  109808751 :     hiremainder=uel(y,2); ly--; y++;
     518                 :            :   }
     519                 :  123458598 :   xi = get_Fl_red(x);
     520 [ +  + ][ +  + ]:  398208309 :   for (i=2; i<ly; i++) (void)divll_pre(y[i],x,xi);
            [ +  + ][ + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     521         [ +  + ]:  123458598 :   if (!hiremainder) return 0;
     522         [ +  + ]:  779245074 :   return (sy > 0)? hiremainder: x - hiremainder;
     523                 :            : }
     524                 :            : 
     525                 :            : /* return |y| \/ x */
     526                 :            : GEN
     527                 :   96047394 : diviu_rem(GEN y, ulong x, ulong *rem)
     528                 :            : {
     529                 :            :   long ly,i;
     530                 :            :   GEN z;
     531                 :            :   ulong xi;
     532                 :            :   LOCAL_HIREMAINDER;
     533                 :            : 
     534         [ -  + ]:   96047394 :   if (!x) pari_err_INV("diviu_rem",gen_0);
     535         [ +  + ]:   96047394 :   if (!signe(y)) { *rem = 0; return gen_0; }
     536                 :            : 
     537                 :   95646849 :   ly = lgefint(y);
     538         [ +  + ]:   95646849 :   if (x <= uel(y,2))
     539                 :            :   {
     540                 :   80587299 :     hiremainder=0;
     541         [ +  + ]:   80587299 :     if (ly==3)
     542                 :            :     {
     543                 :   25168284 :       z = cgetipos(3);
     544 [ +  - ][ #  # ]:   25168284 :       z[2] = divll(uel(y,2),x);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     545                 :   25168284 :       *rem = hiremainder; return z;
     546                 :            :     }
     547                 :            :   }
     548                 :            :   else
     549                 :            :   {
     550         [ +  + ]:   15059550 :     if (ly==3) { *rem = uel(y,2); return gen_0; }
     551                 :   12238128 :     hiremainder = uel(y,2); ly--; y++;
     552                 :            :   }
     553                 :   67657143 :   xi = get_Fl_red(x);
     554                 :   67657143 :   z = cgetipos(ly);
     555 [ +  + ][ +  + ]:  417935706 :   for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x,xi);
            [ +  + ][ + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     556                 :   96047394 :   *rem = hiremainder; return z;
     557                 :            : }
     558                 :            : 
     559                 :            : GEN
     560                 :   29938350 : divis_rem(GEN y, long x, long *rem)
     561                 :            : {
     562                 :   29938350 :   long sy=signe(y),ly,s,i;
     563                 :            :   GEN z;
     564                 :            :   ulong xi;
     565                 :            :   LOCAL_HIREMAINDER;
     566                 :            : 
     567         [ -  + ]:   29938350 :   if (!x) pari_err_INV("divis_rem",gen_0);
     568         [ +  + ]:   29938350 :   if (!sy) { *rem=0; return gen_0; }
     569         [ +  + ]:   24877572 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     570                 :            : 
     571                 :   24877572 :   ly = lgefint(y);
     572         [ +  + ]:   24877572 :   if ((ulong)x <= uel(y,2))
     573                 :            :   {
     574                 :   20975988 :     hiremainder=0;
     575         [ +  + ]:   20975988 :     if (ly==3)
     576                 :            :     {
     577                 :   20637681 :       z = cgeti(3); z[1] = evallgefint(3) | evalsigne(s);
     578 [ +  - ][ #  # ]:   20637681 :       z[2] = divll(uel(y,2),x);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     579 [ +  + ][ +  + ]:   20637681 :       if (sy<0) hiremainder = - ((long)hiremainder);
     580                 :   20637681 :       *rem = (long)hiremainder; return z;
     581                 :            :     }
     582                 :            :   }
     583                 :            :   else
     584                 :            :   {
     585         [ +  + ]:    3901584 :     if (ly==3) { *rem = itos(y); return gen_0; }
     586                 :      67371 :     hiremainder = uel(y,2); ly--; y++;
     587                 :            :   }
     588                 :     405678 :   xi = get_Fl_red(x);
     589                 :     405678 :   z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
     590 [ +  + ][ +  + ]:    3509076 :   for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x,xi);
            [ +  + ][ + ]
         [ +  + ][ +  - ]
         [ +  + ][ +  + ]
     591         [ +  + ]:     405678 :   if (sy<0) hiremainder = - ((long)hiremainder);
     592                 :   29938350 :   *rem = (long)hiremainder; return z;
     593                 :            : }
     594                 :            : 
     595                 :            : GEN
     596                 :     195903 : divis(GEN y, long x)
     597                 :            : {
     598                 :     195903 :   long sy=signe(y),ly,s,i;
     599                 :            :   ulong xi;
     600                 :            :   GEN z;
     601                 :            :   LOCAL_HIREMAINDER;
     602                 :            : 
     603         [ -  + ]:     195903 :   if (!x) pari_err_INV("divis",gen_0);
     604         [ -  + ]:     195903 :   if (!sy) return gen_0;
     605         [ -  + ]:     195903 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     606                 :            : 
     607                 :     195903 :   ly = lgefint(y);
     608         [ +  + ]:     195903 :   if ((ulong)x <= uel(y,2))
     609                 :            :   {
     610                 :     178935 :     hiremainder=0;
     611         [ +  + ]:     178935 :     if (ly==3)
     612                 :            :     {
     613                 :      29145 :       z = cgeti(3); z[1] = evallgefint(3) | evalsigne(s);
     614 [ +  - ][ #  # ]:      29145 :       z[2] = divll(y[2],x);
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
     615                 :      29145 :       return z;
     616                 :            :     }
     617                 :            :   }
     618                 :            :   else
     619                 :            :   {
     620         [ -  + ]:      16968 :     if (ly==3) return gen_0;
     621                 :      16968 :     hiremainder=y[2]; ly--; y++;
     622                 :            :   }
     623                 :     166758 :   xi = get_Fl_red(x);
     624                 :     166758 :   z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
     625 [ +  + ][ +  + ]:    1888401 :   for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x, xi);
            [ -  + ][ + ]
         [ +  + ][ +  - ]
         [ +  + ][ +  + ]
     626                 :     195903 :   return z;
     627                 :            : }
     628                 :            : 
     629                 :            : GEN
     630                 :   72938664 : divrr(GEN x, GEN y)
     631                 :            : {
     632                 :   72938664 :   long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
     633                 :            :   ulong y0,y1;
     634                 :            :   GEN r, r1;
     635                 :            : 
     636         [ -  + ]:   72938664 :   if (!x) pari_err_INV("divrr",y);
     637                 :   72938664 :   e = expo(x) - expo(y);
     638         [ +  + ]:   72938664 :   if (!sx) return real_0_bit(e);
     639         [ +  + ]:   71671188 :   if (sy<0) sx = -sx;
     640                 :            : 
     641                 :   71671188 :   lx=lg(x); ly=lg(y);
     642         [ +  + ]:   71671188 :   if (ly==3)
     643                 :            :   {
     644         [ +  + ]:   58427856 :     ulong k = x[2], l = (lx>3)? x[3]: 0;
     645                 :            :     LOCAL_HIREMAINDER;
     646         [ +  + ]:   58427856 :     if (k < uel(y,2)) e--;
     647                 :            :     else
     648                 :            :     {
     649         [ +  + ]:   27095247 :       l >>= 1; if (k&1) l |= HIGHBIT;
     650                 :   27095247 :       k >>= 1;
     651                 :            :     }
     652 [ -  + ][ -  + ]:   58427856 :     hiremainder = k; k = divll(l,y[2]);
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
                 [ +  + ]
     653 [ +  + ][ +  + ]:   58427856 :     if (hiremainder & HIGHBIT)
     654                 :            :     {
     655                 :   16174497 :       k++;
     656         [ -  + ]:   16174497 :       if (!k) { k = HIGHBIT; e++; }
     657                 :            :     }
     658                 :   58427856 :     r = cgetr(3);
     659                 :   58427856 :     r[1] = evalsigne(sx) | evalexpo(e);
     660                 :   58427856 :     r[2] = k; return r;
     661                 :            :   }
     662                 :            : 
     663                 :   13243332 :   lr = minss(lx,ly); r = new_chunk(lr);
     664                 :   13243332 :   r1 = r-1;
     665         [ +  + ]:  105617259 :   r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
     666         [ +  + ]:   13243332 :   r1[lr] = (lx>ly)? x[lr]: 0;
     667                 :   13243332 :   y0 = y[2]; y1 = y[3];
     668         [ +  + ]:  118860591 :   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         [ +  + ]:  105617259 :     if (uel(r1,1) == y0)
     675                 :            :     {
     676                 :     326178 :       qp = ULONG_MAX; k = addll(y0,r1[2]);
     677                 :            :     }
     678                 :            :     else
     679                 :            :     {
     680         [ -  + ]:  105291081 :       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                 :  105291081 :       hiremainder = r1[1]; overflow = 0;
     688 [ +  + ][ -  + ]:  105291081 :       qp = divll(r1[2],y0); k = hiremainder;
         [ +  - ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
                 [ +  + ]
     689                 :            :     }
     690                 :  105617259 :     j = lr-i+1;
     691         [ +  + ]:  105617259 :     if (!overflow)
     692                 :            :     {
     693                 :            :       long k3, k4;
     694                 :  105361722 :       k3 = mulll(qp,y1);
     695 [ +  + ][ +  + ]:  105361722 :       if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
     696                 :   13234569 :         k4 = subll(hiremainder,k);
     697                 :            :       else
     698                 :            :       {
     699                 :   92127153 :         k3 = subll(k3, r1[3]);
     700                 :   92127153 :         k4 = subllx(hiremainder,k);
     701                 :            :       }
     702 [ +  + ][ +  + ]:  129261893 :       while (!overflow && k4) { qp--; k3 = subll(k3,y1); k4 = subllx(k4,y0); }
         [ +  + ][ +  + ]
     703                 :            :     }
     704         [ +  + ]:  105617259 :     if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
     705 [ +  + ][ +  + ]: 1085344467 :     for (j--; j>1; j--)
     706                 :            :     {
     707                 :  979727208 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     708                 :  979727208 :       hiremainder += overflow;
     709                 :            :     }
     710         [ +  + ]:  105617259 :     if (uel(r1,1) != hiremainder)
     711                 :            :     {
     712         [ +  - ]:     121551 :       if (uel(r1,1) < hiremainder)
     713                 :            :       {
     714                 :     121551 :         qp--;
     715                 :     121551 :         j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
     716 [ +  + ][ +  + ]:     850803 :         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                 :  105617259 :     *++r1 = qp;
     731                 :            :   }
     732                 :            :   /* i = lr-1 */
     733                 :            :   /* round correctly */
     734         [ +  + ]:   13243332 :   if (uel(r1,1) > (y0>>1))
     735                 :            :   {
     736 [ +  + ][ +  + ]:    6484472 :     j=i; do uel(r,--j)++; while (j && !r[j]);
     737                 :            :   }
     738         [ +  + ]:  105617259 :   r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
     739         [ +  + ]:   13243332 :   if (r[0] == 0) e--;
     740         [ +  + ]:    5339016 :   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                 :   13243332 :   r[0] = evaltyp(t_REAL)|evallg(lr);
     745                 :   13243332 :   r[1] = evalsigne(sx) | evalexpo(e);
     746                 :   72938664 :   return r;
     747                 :            : }
     748                 :            : 
     749                 :            : GEN
     750                 :   11851179 : divri(GEN x, GEN y)
     751                 :            : {
     752                 :   11851179 :   long lx, s = signe(y);
     753                 :            :   pari_sp av;
     754                 :            :   GEN z;
     755                 :            : 
     756         [ -  + ]:   11851179 :   if (!s) pari_err_INV("divri",y);
     757         [ +  + ]:   11851179 :   if (!signe(x)) return real_0_bit(expo(x) - expi(y));
     758         [ +  + ]:   11842863 :   if (!is_bigint(y)) {
     759                 :    9780486 :     GEN z = divru(x, y[2]);
     760         [ +  + ]:    9780486 :     if (s < 0) togglesign(z);
     761                 :    9780486 :     return z;
     762                 :            :   }
     763                 :    2062377 :   lx = lg(x); z = cgetr(lx); av = avma;
     764                 :    2062377 :   affrr(divrr(x, itor(y, lx+1)), z);
     765                 :   11851179 :   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                 :  469454454 : dvmdii(GEN x, GEN y, GEN *z)
     777                 :            : {
     778                 :  469454454 :   long sx = signe(x), sy = signe(y);
     779                 :  469454454 :   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         [ +  + ]:  469454454 :   if (!sx)
     785                 :            :   {
     786         [ +  + ]:   11663841 :     if (ly < 3) pari_err_INV("dvmdii",gen_0);
     787 [ +  + ][ +  + ]:   11663826 :     if (!z || z == ONLY_REM) return gen_0;
     788                 :    5288961 :     *z=gen_0; return gen_0;
     789                 :            :   }
     790         [ +  + ]:  457790613 :   if (ly <= 3)
     791                 :            :   {
     792                 :            :     ulong rem;
     793         [ -  + ]:  292950726 :     if (ly < 3) pari_err_INV("dvmdii",gen_0);
     794         [ +  + ]:  292950726 :     if (z == ONLY_REM)
     795                 :            :     {
     796                 :  262454661 :       rem = umodiu(x,uel(y,2));
     797         [ +  + ]:  262454661 :       if (!rem) return gen_0;
     798         [ +  + ]:  244468320 :       return (sx < 0)? utoineg(uel(y,2) - rem): utoipos(rem);
     799                 :            :     }
     800                 :   30496065 :     q = diviu_rem(x, uel(y,2), &rem);
     801         [ +  + ]:   30496065 :     if (sx != sy) togglesign(q);
     802         [ +  + ]:   30496065 :     if (!z) return q;
     803         [ +  + ]:   29127759 :     if (!rem) *z = gen_0;
     804         [ +  + ]:   12610395 :     else *z = sx < 0? utoineg(rem): utoipos(rem);
     805                 :  292950726 :     return q;
     806                 :            :   }
     807                 :  164839887 :   lx=lgefint(x);
     808                 :  164839887 :   lz=lx-ly;
     809         [ +  + ]:  164839887 :   if (lz <= 0)
     810                 :            :   {
     811         [ +  + ]:   69867201 :     if (lz == 0)
     812                 :            :     {
     813         [ +  + ]:   59287356 :       for (i=2; i<lx; i++)
     814         [ +  + ]:   59148711 :         if (x[i] != y[i])
     815                 :            :         {
     816         [ +  + ]:   57381474 :           if (uel(x,i) > uel(y,i)) goto DIVIDE;
     817                 :   14569953 :           goto TRIVIAL;
     818                 :            :         }
     819         [ +  + ]:     138645 :       if (z == ONLY_REM) return gen_0;
     820         [ +  + ]:       6012 :       if (z) *z = gen_0;
     821         [ +  + ]:       6012 :       if (sx < 0) sy = -sy;
     822                 :       6012 :       return stoi(sy);
     823                 :            :     }
     824                 :            : TRIVIAL:
     825         [ +  + ]:   26917035 :     if (z == ONLY_REM) return icopy(x);
     826         [ +  + ]:    2423043 :     if (z) *z = icopy(x);
     827                 :    2423043 :     return gen_0;
     828                 :            :   }
     829                 :            : DIVIDE: /* quotient is non-zero */
     830         [ +  + ]:  137784207 :   av=avma; if (sx<0) sy = -sy;
     831 [ +  + ][ +  + ]:  137784207 :   r1 = new_chunk(lx); sh = bfffo(y[2]);
         [ +  + ][ +  + ]
     832 [ +  + ][ +  + ]:  137784207 :   if (sh)
     833                 :            :   { /* normalize so that highbit(y) = 1 (shift left x and y by sh bits)*/
     834                 :  130390809 :     register const ulong m = BITS_IN_LONG - sh;
     835                 :  130390809 :     r = new_chunk(ly);
     836                 :  130390809 :     shift_left(r, y,2,ly-1, 0,sh); y = r;
     837                 :  130390809 :     shift_left(r1,x,2,lx-1, 0,sh);
     838                 :  130390809 :     r1[1] = uel(x,2) >> m;
     839                 :            :   }
     840                 :            :   else
     841                 :            :   {
     842         [ +  + ]:   52754289 :     r1[1] = 0; for (j=2; j<lx; j++) r1[j] = x[j];
     843                 :            :   }
     844                 :  137784207 :   x = r1;
     845                 :  137784207 :   y0 = y[2]; y0i = get_Fl_red(y0);
     846                 :  137784207 :   y1 = y[3];
     847         [ +  + ]:  592475934 :   for (i=0; i<=lz; i++)
     848                 :            :   { /* r1 = x + i */
     849                 :            :     ulong k, qp;
     850                 :            :     LOCAL_HIREMAINDER;
     851                 :            :     LOCAL_OVERFLOW;
     852                 :            : 
     853         [ +  + ]:  454691727 :     if (uel(r1,1) == y0)
     854                 :            :     {
     855                 :       4973 :       qp = ULONG_MAX; k = addll(y0,r1[2]);
     856                 :            :     }
     857                 :            :     else
     858                 :            :     {
     859                 :  454686754 :       hiremainder = r1[1]; overflow = 0;
     860 [ +  + ][ +  - ]:  454686754 :       qp = divll_pre(r1[2],y0,y0i); k = hiremainder;
         [ +  - ][ +  - ]
                 [ -  + ]
     861                 :            :     }
     862         [ +  + ]:  454691727 :     if (!overflow)
     863                 :            :     {
     864                 :  454687313 :       long k3 = subll(mulll(qp,y1), r1[3]);
     865                 :  454687313 :       long k4 = subllx(hiremainder,k);
     866 [ +  + ][ +  + ]:  540886234 :       while (!overflow && k4) { qp--; k3 = subll(k3,y1); k4 = subllx(k4,y0); }
         [ +  + ][ +  + ]
     867                 :            :     }
     868                 :  454691727 :     hiremainder = 0; j = ly;
     869 [ +  + ][ +  + ]: 5541574428 :     for (j--; j>1; j--)
     870                 :            :     {
     871                 : 5086882701 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     872                 : 5086882701 :       hiremainder += overflow;
     873                 :            :     }
     874         [ +  + ]:  454691727 :     if (uel(r1,1) < hiremainder)
     875                 :            :     {
     876                 :     138531 :       qp--;
     877                 :     138531 :       j = ly-1; r1[j] = addll(r1[j],y[j]);
     878 [ +  + ][ +  + ]:    1324058 :       for (j--; j>1; j--) r1[j] = addllx(r1[j],y[j]);
     879                 :            :     }
     880                 :  454691727 :     *++r1 = qp;
     881                 :            :   }
     882                 :            : 
     883                 :  137784207 :   lq = lz+2;
     884         [ +  + ]:  137784207 :   if (!z)
     885                 :            :   {
     886                 :     446988 :     qd = (ulong*)av;
     887                 :     446988 :     xd = (ulong*)(x + lq);
     888         [ +  + ]:     446988 :     if (x[1]) { lz++; lq++; }
     889         [ +  + ]:     908004 :     while (lz--) *--qd = *--xd;
     890                 :     446988 :     *--qd = evalsigne(sy) | evallgefint(lq);
     891                 :     446988 :     *--qd = evaltyp(t_INT) | evallg(lq);
     892                 :     446988 :     avma = (pari_sp)qd; return (GEN)qd;
     893                 :            :   }
     894                 :            : 
     895 [ +  + ][ +  + ]:  157473036 :   j=lq; while (j<lx && !x[j]) j++;
     896                 :  137337219 :   lz = lx-j;
     897         [ +  + ]:  137337219 :   if (z == ONLY_REM)
     898                 :            :   {
     899         [ +  + ]:  108360537 :     if (lz==0) { avma = av; return gen_0; }
     900                 :  106305480 :     rd = (ulong*)av; lr = lz+2;
     901                 :  106305480 :     xd = (ulong*)(x + lx);
     902 [ +  + ][ +  + ]:  126578328 :     if (!sh) while (lz--) *--rd = *--xd;
     903                 :            :     else
     904                 :            :     { /* shift remainder right by sh bits */
     905                 :   98975997 :       const ulong shl = BITS_IN_LONG - sh;
     906                 :            :       ulong l;
     907                 :   98975997 :       xd--;
     908         [ +  + ]:  358174737 :       while (--lz) /* fill r[3..] */
     909                 :            :       {
     910                 :  259198740 :         l = *xd >> sh;
     911                 :  259198740 :         *--rd = l | (*--xd << shl);
     912                 :            :       }
     913                 :   98975997 :       l = *xd >> sh;
     914         [ +  + ]:   98975997 :       if (l) *--rd = l; else lr--;
     915                 :            :     }
     916                 :  106305480 :     *--rd = evalsigne(sx) | evallgefint(lr);
     917                 :  106305480 :     *--rd = evaltyp(t_INT) | evallg(lr);
     918                 :  106305480 :     avma = (pari_sp)rd; return (GEN)rd;
     919                 :            :   }
     920                 :            : 
     921                 :   28976682 :   lr = lz+2;
     922                 :   28976682 :   rd = NULL; /* gcc -Wall */
     923         [ +  + ]:   28976682 :   if (lz)
     924                 :            :   { /* non zero remainder: initialize rd */
     925                 :   27334347 :     xd = (ulong*)(x + lx);
     926         [ +  + ]:   27334347 :     if (!sh)
     927                 :            :     {
     928                 :      17391 :       rd = (ulong*)avma; (void)new_chunk(lr);
     929         [ +  + ]:     158688 :       while (lz--) *--rd = *--xd;
     930                 :            :     }
     931                 :            :     else
     932                 :            :     { /* shift remainder right by sh bits */
     933                 :   27316956 :       const ulong shl = BITS_IN_LONG - sh;
     934                 :            :       ulong l;
     935                 :   27316956 :       rd = (ulong*)x; /* overwrite shifted y */
     936                 :   27316956 :       xd--;
     937         [ +  + ]:  105085101 :       while (--lz)
     938                 :            :       {
     939                 :   77768145 :         l = *xd >> sh;
     940                 :   77768145 :         *--rd = l | (*--xd << shl);
     941                 :            :       }
     942                 :   27316956 :       l = *xd >> sh;
     943         [ +  + ]:   27316956 :       if (l) *--rd = l; else lr--;
     944                 :            :     }
     945                 :   27334347 :     *--rd = evalsigne(sx) | evallgefint(lr);
     946                 :   27334347 :     *--rd = evaltyp(t_INT) | evallg(lr);
     947                 :   27334347 :     rd += lr;
     948                 :            :   }
     949                 :   28976682 :   qd = (ulong*)av;
     950                 :   28976682 :   xd = (ulong*)(x + lq);
     951         [ +  + ]:   28976682 :   if (x[1]) lq++;
     952         [ +  + ]:   94424622 :   j = lq-2; while (j--) *--qd = *--xd;
     953                 :   28976682 :   *--qd = evalsigne(sy) | evallgefint(lq);
     954                 :   28976682 :   *--qd = evaltyp(t_INT) | evallg(lq);
     955                 :   28976682 :   q = (GEN)qd;
     956         [ +  + ]:   28976682 :   if (lr==2) *z = gen_0;
     957                 :            :   else
     958                 :            :   { /* rd has been properly initialized: we had lz > 0 */
     959         [ +  + ]:  168389885 :     while (lr--) *--qd = *--rd;
     960                 :   27334347 :     *z = (GEN)qd;
     961                 :            :   }
     962                 :  469454439 :   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                 :    6389556 : red_montgomery(GEN T, GEN N, ulong inv)
     971                 :            : {
     972                 :            :   pari_sp av;
     973                 :            :   GEN Te, Td, Ne, Nd, scratch;
     974                 :    6389556 :   ulong i, j, m, t, d, k = NLIMBS(N);
     975                 :            :   int carry;
     976                 :            :   LOCAL_HIREMAINDER;
     977                 :            :   LOCAL_OVERFLOW;
     978                 :            : 
     979         [ -  + ]:    6389556 :   if (k == 0) return gen_0;
     980                 :    6389556 :   d = NLIMBS(T); /* <= 2*k */
     981         [ +  + ]:    6389556 :   if (d == 0) return gen_0;
     982                 :            : #ifdef DEBUG
     983                 :            :   if (d > 2*k) pari_err_BUG("red_montgomery");
     984                 :            : #endif
     985         [ -  + ]:    6389547 :   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                 :    6389547 :   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                 :    6389547 :   Td = (GEN)av;
    1007                 :    6389547 :   Te = T + (d+2);
    1008         [ +  + ]:   53456904 :   for (i=0; i < d     ; i++) *--Td = *--Te;
    1009         [ +  + ]:   15637104 :   for (   ; i < (k<<1); i++) *--Td = 0;
    1010                 :            : 
    1011                 :    6389547 :   Te = (GEN)av; /* 1 beyond end of current T mantissa (in scratch) */
    1012                 :    6389547 :   Ne = N + k+2; /* 1 beyond end of N mantissa */
    1013                 :            : 
    1014                 :    6389547 :   carry = 0;
    1015         [ +  + ]:   34547004 :   for (i=0; i<k; i++) /* set T := T/B nod N, k times */
    1016                 :            :   {
    1017                 :   28157457 :     Td = Te; /* one beyond end of (new) T mantissa */
    1018                 :   28157457 :     Nd = Ne;
    1019                 :   28157457 :     hiremainder = *--Td;
    1020                 :   28157457 :     m = hiremainder * inv; /* solve T + m N = O(B) */
    1021                 :            : 
    1022                 :            :     /* set T := (T + mN) / B */
    1023                 :   28157457 :     Te = Td;
    1024                 :   28157457 :     (void)addmul(m, *--Nd); /* = 0 */
    1025 [ +  + ][ +  + ]:  284758575 :     for (j=1; j<k; j++)
    1026                 :            :     {
    1027                 :  256601118 :       t = addll(addmul(m, *--Nd), *--Td);
    1028                 :  256601118 :       *Td = t;
    1029                 :  256601118 :       hiremainder += overflow;
    1030                 :            :     }
    1031                 :   28157457 :     t = addll(hiremainder, *--Td); *Td = t + carry;
    1032 [ +  + ][ +  + ]:   28157457 :     carry = (overflow || (carry && *Td == 0));
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
    1033                 :            :   }
    1034         [ +  + ]:    6389547 :   if (carry)
    1035                 :            :   { /* Td > N overflows (k+1 words), set Td := Td - N */
    1036                 :      87030 :     Td = Te;
    1037                 :      87030 :     Nd = Ne;
    1038                 :      87030 :     t = subll(*--Td, *--Nd); *Td = t;
    1039 [ +  + ][ +  + ]:    1250793 :     while (Td > scratch) { t = subllx(*--Td, *--Nd); *Td = t; }
    1040                 :            :   }
    1041                 :            : 
    1042                 :            :   /* copy result */
    1043                 :    6389547 :   Td = (GEN)av;
    1044 [ +  + ][ +  - ]:   10089870 :   while (*scratch == 0 && Te > scratch) scratch++; /* strip leading 0s */
    1045         [ +  + ]:   30846681 :   while (Te > scratch) *--Td = *--Te;
    1046         [ -  + ]:    6389547 :   k = (GEN)av - Td; if (!k) { avma = av; return gen_0; }
    1047                 :    6389547 :   k += 2;
    1048                 :    6389547 :   *--Td = evalsigne(1) | evallgefint(k);
    1049                 :    6389547 :   *--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                 :    6389556 :   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                 :   13995918 : diviuexact_i(GEN x, ulong y)
    1068                 :            : {
    1069                 :            :   long i, lz, lx;
    1070                 :            :   ulong q, yinv;
    1071                 :            :   GEN z, z0, x0, x0min;
    1072                 :            : 
    1073         [ +  + ]:   13995918 :   if (y == 1) return icopy(x);
    1074                 :   11757840 :   lx = lgefint(x);
    1075         [ +  + ]:   11757840 :   if (lx == 3) return utoipos(uel(x,2) / y);
    1076                 :   11266179 :   yinv = invmod2BIL(y);
    1077         [ +  + ]:   11266179 :   lz = (y <= uel(x,2)) ? lx : lx-1;
    1078                 :   11266179 :   z = new_chunk(lz);
    1079                 :   11266179 :   z0 = z + lz;
    1080                 :   11266179 :   x0 = x + lx; x0min = x + lx-lz+2;
    1081                 :            : 
    1082         [ +  + ]:   36840741 :   while (x0 > x0min)
    1083                 :            :   {
    1084                 :   25574562 :     *--z0 = q = yinv*((ulong)*--x0); /* i-th quotient */
    1085         [ +  + ]:   25574562 :     if (!q) continue;
    1086                 :            :     /* x := x - q * y */
    1087                 :            :     { /* update neither lowest word (could set it to 0) nor highest ones */
    1088                 :   23009064 :       register GEN x1 = x0 - 1;
    1089                 :            :       LOCAL_HIREMAINDER;
    1090                 :   23009064 :       (void)mulll(q,y);
    1091 [ +  + ][ +  + ]:   23009064 :       if (hiremainder)
    1092                 :            :       {
    1093         [ +  + ]:   19444887 :         if ((ulong)*x1 < hiremainder)
    1094                 :            :         {
    1095                 :      51678 :           *x1 -= hiremainder;
    1096         [ +  + ]:      53412 :           do (*--x1)--; while ((ulong)*x1 == ULONG_MAX);
    1097                 :            :         }
    1098                 :            :         else
    1099                 :   19393209 :           *x1 -= hiremainder;
    1100                 :            :       }
    1101                 :            :     }
    1102                 :            :   }
    1103         [ -  + ]:   11266179 :   i=2; while(!z[i]) i++;
    1104                 :   11266179 :   z += i-2; lz -= i-2;
    1105                 :   11266179 :   z[0] = evaltyp(t_INT)|evallg(lz);
    1106                 :   11266179 :   z[1] = evalsigne(1)|evallg(lz);
    1107                 :   13995918 :   avma = (pari_sp)z; return z;
    1108                 :            : }
    1109                 :            : 
    1110                 :            : /* assume y != 0 and the division is exact */
    1111                 :            : GEN
    1112                 :    7101501 : diviuexact(GEN x, ulong y)
    1113                 :            : {
    1114                 :            :   pari_sp av;
    1115                 :    7101501 :   long lx, vy, s = signe(x);
    1116                 :            :   GEN z;
    1117                 :            : 
    1118         [ +  + ]:    7101501 :   if (!s) return gen_0;
    1119         [ +  + ]:    6904299 :   if (y == 1) return icopy(x);
    1120                 :    6900519 :   lx = lgefint(x);
    1121         [ +  + ]:    6900519 :   if (lx == 3) {
    1122                 :    6558633 :     ulong q = uel(x,2) / y;
    1123         [ +  + ]:    6558633 :     return (s > 0)? utoipos(q): utoineg(q);
    1124                 :            :   }
    1125                 :     341886 :   av = avma; (void)new_chunk(lx); vy = vals(y);
    1126         [ +  + ]:     341886 :   if (vy) {
    1127                 :     231273 :     y >>= vy;
    1128         [ +  + ]:     231273 :     if (y == 1) { avma = av; return shifti(x, -vy); }
    1129                 :     204954 :     x = shifti(x, -vy);
    1130         [ -  + ]:     204954 :     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                 :     110613 :   } else x = icopy(x);
    1136                 :     315567 :   avma = av;
    1137                 :     315567 :   z = diviuexact_i(x, y);
    1138                 :    7101501 :   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                 :  277152057 : diviiexact(GEN x, GEN y)
    1146                 :            : {
    1147                 :  277152057 :   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         [ -  + ]:  277152057 :   if (!sy) pari_err_INV("diviiexact",gen_0);
    1153         [ +  + ]:  277152057 :   if (!sx) return gen_0;
    1154                 :   93707358 :   lx = lgefint(x);
    1155         [ +  + ]:   93707358 :   if (lx == 3) {
    1156                 :   74453370 :     q = uel(x,2) / uel(y,2);
    1157         [ +  + ]:   74453370 :     return (sx+sy) ? utoipos(q): utoineg(q);
    1158                 :            :   }
    1159                 :   19253988 :   vy = vali(y); av = avma;
    1160                 :   19253988 :   (void)new_chunk(lx); /* enough room for z */
    1161         [ +  + ]:   19253988 :   if (vy)
    1162                 :            :   { /* make y odd */
    1163                 :    9508764 :     y = shifti(y,-vy);
    1164                 :    9508764 :     x = shifti(x,-vy); lx = lgefint(x);
    1165                 :            :   }
    1166                 :    9745224 :   else x = icopy(x); /* necessary because we destroy x */
    1167                 :   19253988 :   avma = av; /* will erase our x,y when exiting */
    1168                 :            :   /* now y is odd */
    1169                 :   19253988 :   ly = lgefint(y);
    1170         [ +  + ]:   19253988 :   if (ly == 3)
    1171                 :            :   {
    1172                 :   13680351 :     x = diviuexact_i(x,uel(y,2)); /* x != 0 */
    1173         [ +  + ]:   13680351 :     setsigne(x, (sx+sy)? 1: -1); return x;
    1174                 :            :   }
    1175                 :    5573637 :   y0inv = invmod2BIL(y[ly-1]);
    1176 [ +  + ][ +  + ]:   10130820 :   i=2; while (i<ly && y[i]==x[i]) i++;
    1177 [ +  + ][ +  + ]:    5573637 :   lz = (i==ly || uel(y,i) < uel(x,i)) ? lx-ly+3 : lx-ly+2;
    1178                 :    5573637 :   z = new_chunk(lz);
    1179                 :            : 
    1180                 :    5573637 :   y += ly - 1; /* now y[-i] = i-th word of y */
    1181         [ +  + ]:   26358246 :   for (ii=lx-1,i=lz-1; i>=2; i--,ii--)
    1182                 :            :   {
    1183                 :            :     long limj;
    1184                 :            :     LOCAL_HIREMAINDER;
    1185                 :            :     LOCAL_OVERFLOW;
    1186                 :            : 
    1187                 :   20784609 :     z[i] = q = y0inv*uel(x,ii); /* i-th quotient */
    1188         [ +  + ]:   20784609 :     if (!q) continue;
    1189                 :            : 
    1190                 :            :     /* x := x - q * y */
    1191                 :   18686748 :     (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                 :   18686748 :       register GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj;
    1194 [ +  + ][ +  + ]:   87651444 :       for (; x0 >= xlim; x0--, y0--)
    1195                 :            :       {
    1196                 :   68964696 :         *x0 = subll(*x0, addmul(q,*y0));
    1197                 :   68964696 :         hiremainder += overflow;
    1198                 :            :       }
    1199 [ +  + ][ +  + ]:   18686748 :       if (hiremainder && limj != lx - lz)
    1200                 :            :       {
    1201         [ +  + ]:    9916650 :         if ((ulong)*x0 < hiremainder)
    1202                 :            :         {
    1203                 :      79686 :           *x0 -= hiremainder;
    1204         [ -  + ]:      79686 :           do (*--x0)--; while ((ulong)*x0 == ULONG_MAX);
    1205                 :            :         }
    1206                 :            :         else
    1207                 :    9836964 :           *x0 -= hiremainder;
    1208                 :            :       }
    1209                 :            :     }
    1210                 :            :   }
    1211         [ -  + ]:    5573637 :   i=2; while(!z[i]) i++;
    1212                 :    5573637 :   z += i-2; lz -= (i-2);
    1213                 :    5573637 :   z[0] = evaltyp(t_INT)|evallg(lz);
    1214         [ +  + ]:    5573637 :   z[1] = evalsigne((sx+sy)? 1: -1) | evallg(lz);
    1215                 :  277152057 :   avma = (pari_sp)z; return z;
    1216                 :            : }
    1217                 :            : 
    1218                 :            : /* assume yz != and yz | x */
    1219                 :            : GEN
    1220                 :     172902 : diviuuexact(GEN x, ulong y, ulong z)
    1221                 :            : {
    1222                 :            :   long tmp[4];
    1223                 :            :   ulong t;
    1224                 :            :   LOCAL_HIREMAINDER;
    1225                 :     172902 :   t = mulll(y, z);
    1226 [ +  - ][ +  - ]:     172902 :   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                 :     172902 :   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                 : 1201441893 : 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         [ +  + ]: 1201441893 :   if (ny == 1) return muluispec((ulong)*y, x, nx);
    1248         [ +  + ]:  469271715 :   if (ny == 0) return gen_0;
    1249                 :  467983605 :   zd = (GEN)avma; lz = nx+ny+2;
    1250                 :  467983605 :   (void)new_chunk(lz);
    1251                 :  467983605 :   xd = x + nx;
    1252                 :  467983605 :   yd = y + ny;
    1253                 :  467983605 :   ye = yd; p1 = *--xd;
    1254                 :            : 
    1255                 :  467983605 :   *--zd = mulll(p1, *--yd); z2e = zd;
    1256      [ +  +  + ]: 2808067839 :   while (yd > y) *--zd = addmul(p1, *--yd);
    1257                 :  467983605 :   *--zd = hiremainder;
    1258                 :            : 
    1259         [ +  + ]: 3157634301 :   while (xd > x)
    1260                 :            :   {
    1261                 :            :     LOCAL_OVERFLOW;
    1262                 : 2689650696 :     yd = ye; p1 = *--xd;
    1263                 :            : 
    1264                 : 2689650696 :     z2d = --z2e;
    1265                 : 2689650696 :     *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
    1266 [ +  + ][ +  + ]:30324878592 :     while (yd > y)
    1267                 :            :     {
    1268                 :27635227896 :       hiremainder += overflow;
    1269                 :27635227896 :       *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
    1270                 :            :     }
    1271                 : 2689650696 :     *--zd = hiremainder + overflow;
    1272                 :            :   }
    1273         [ +  + ]:  467983605 :   if (*zd == 0) { zd++; lz--; } /* normalize */
    1274                 :  467983605 :   *--zd = evalsigne(1) | evallgefint(lz);
    1275                 :  467983605 :   *--zd = evaltyp(t_INT) | evallg(lz);
    1276                 : 1201441893 :   avma=(pari_sp)zd; return zd;
    1277                 :            : }
    1278                 :            : 
    1279                 :            : INLINE GEN
    1280                 :  464520192 : 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         [ +  + ]:  464520192 :   if (nx == 1) return sqru((ulong)*x);
    1288         [ +  + ]:  411249003 :   if (nx == 0) return gen_0;
    1289                 :   57999048 :   zd = (GEN)avma; lz = (nx+1) << 1;
    1290                 :   57999048 :   z0 = new_chunk(lz);
    1291         [ -  + ]:   57999048 :   if (nx == 1)
    1292                 :            :   {
    1293                 :          0 :     *--zd = mulll(*x, *x);
    1294                 :          0 :     *--zd = hiremainder; goto END;
    1295                 :            :   }
    1296                 :   57999048 :   xd = x + nx;
    1297                 :            : 
    1298                 :            :   /* compute double products --> zd */
    1299                 :   57999048 :   p1 = *--xd; yd = xd; --zd;
    1300                 :   57999048 :   *--zd = mulll(p1, *--yd); z2e = zd;
    1301 [ +  + ][ +  + ]:  338888604 :   while (yd > x) *--zd = addmul(p1, *--yd);
    1302                 :   57999048 :   *--zd = hiremainder;
    1303                 :            : 
    1304                 :   57999048 :   x0 = x+1;
    1305         [ +  + ]:  338888604 :   while (xd > x0)
    1306                 :            :   {
    1307                 :            :     LOCAL_OVERFLOW;
    1308                 :  280889556 :     p1 = *--xd; yd = xd;
    1309                 :            : 
    1310                 :  280889556 :     z2e -= 2; z2d = z2e;
    1311                 :  280889556 :     *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
    1312 [ +  + ][ +  + ]: 2430986244 :     while (yd > x)
    1313                 :            :     {
    1314                 : 2150096688 :       hiremainder += overflow;
    1315                 : 2150096688 :       *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
    1316                 :            :     }
    1317                 :  280889556 :     *--zd = hiremainder + overflow;
    1318                 :            :   }
    1319                 :            :   /* multiply zd by 2 (put result in zd - 1) */
    1320                 :   57999048 :   zd[-1] = ((*zd & HIGHBIT) != 0);
    1321                 :   57999048 :   shift_left(zd, zd, 0, (nx<<1)-3, 0, 1);
    1322                 :            : 
    1323                 :            :   /* add the squares */
    1324                 :   57999048 :   xd = x + nx; zd = z0 + lz;
    1325                 :   57999048 :   p1 = *--xd;
    1326                 :   57999048 :   zd--; *zd = mulll(p1,p1);
    1327                 :   57999048 :   zd--; *zd = addll(hiremainder, *zd);
    1328         [ +  + ]:  396887652 :   while (xd > x)
    1329                 :            :   {
    1330                 :  338888604 :     p1 = *--xd;
    1331                 :  338888604 :     zd--; *zd = addll(mulll(p1,p1)+ overflow, *zd);
    1332                 :  338888604 :     zd--; *zd = addll(hiremainder + overflow, *zd);
    1333                 :            :   }
    1334                 :            : 
    1335                 :            : END:
    1336         [ +  + ]:   57999048 :   if (*zd == 0) { zd++; lz--; } /* normalize */
    1337                 :   57999048 :   *--zd = evalsigne(1) | evallgefint(lz);
    1338                 :   57999048 :   *--zd = evaltyp(t_INT) | evallg(lz);
    1339                 :  464520192 :   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                 :       3735 : mulliifft_params(long len, long *k, long *mod, long *bs, long *n, ulong *ord)
    1363                 :            : {
    1364                 :            :   long r;
    1365                 :       3735 :   *k = expu((3*len)>>2)-3;
    1366                 :            :   do {
    1367                 :       3735 :     (*k)--;
    1368                 :       3735 :     r  = *k-(TWOPOTBITS_IN_LONG+2);
    1369                 :       3735 :     *n = 1L<<*k;
    1370                 :       3735 :     *bs = (len+*n-1)>>*k;
    1371                 :       3735 :     *mod= 2**bs+1;
    1372         [ +  + ]:       3735 :     if (r>0)
    1373                 :        702 :       *mod=((*mod+(1L<<r)-1)>>r)<<r;
    1374         [ -  + ]:       3735 :   } while(*mod>=3**bs);
    1375                 :       3735 :   *ord= 4**mod*BITS_IN_LONG;
    1376                 :       3735 : }
    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                 :   32912448 : Zf_add(GEN a, GEN b, GEN M)
    1385                 :            : {
    1386                 :   32912448 :   GEN y, z = addii(a,b);
    1387                 :   32912448 :   long mod = lgefint(M)-3;
    1388                 :   32912448 :   long l = NLIMBS(z);
    1389         [ +  + ]:   32912448 :   if (l<=mod) return z;
    1390                 :   12920508 :   y = subis(z, 1);
    1391         [ -  + ]:   12920508 :   if (NLIMBS(y)<=mod) return z;
    1392                 :   32912448 :   return int_normalize(y,1);
    1393                 :            : }
    1394                 :            : 
    1395                 :            : static GEN
    1396                 :   33179391 : Zf_sub(GEN a, GEN b, GEN M)
    1397                 :            : {
    1398                 :   33179391 :   GEN z = subii(a,b);
    1399         [ +  + ]:   33179391 :   return signe(z)>=0? z: addii(M,z);
    1400                 :            : }
    1401                 :            : 
    1402                 :            : /* destroy z */
    1403                 :            : static GEN
    1404                 :   64713447 : Zf_red_destroy(GEN z, GEN M)
    1405                 :            : {
    1406                 :   64713447 :   long mod = lgefint(M)-3;
    1407                 :   64713447 :   long l = NLIMBS(z);
    1408                 :            :   GEN y;
    1409         [ +  + ]:   64713447 :   if (l<=mod) return z;
    1410                 :   29391486 :   y = shifti(z, -mod*BITS_IN_LONG);
    1411                 :   29391486 :   z = int_normalize(z, NLIMBS(y));
    1412                 :   29391486 :   y = Zf_red_destroy(y, M);
    1413                 :   29391486 :   z = subii(z, y);
    1414         [ +  + ]:   29391486 :   if (signe(z)<0) z = addii(z, M);
    1415                 :   64713447 :   return z;
    1416                 :            : }
    1417                 :            : 
    1418                 :            : INLINE GEN
    1419                 :   33170409 : 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                 :   32912448 : Zf_mulsqrt2(GEN a, ulong s, ulong ord, GEN M)
    1428                 :            : {
    1429                 :   32912448 :   ulong hord = ord>>1;
    1430         [ +  + ]:   32912448 :   if (!signe(a)) return gen_0;
    1431         [ +  + ]:   30484971 :   if (odd(s)) /* Multiply by 2^(s/2) */
    1432                 :            :   {
    1433                 :     266943 :     GEN az8  = Zf_shift(a,   ord>>4, M);
    1434                 :     266943 :     GEN az83 = Zf_shift(az8, ord>>3, M);
    1435                 :     266943 :     a = Zf_sub(az8, az83, M);
    1436                 :     266943 :     s--;
    1437                 :            :   }
    1438         [ +  + ]:   30484971 :   if (s < hord)
    1439                 :   21272694 :     return Zf_shift(a, s>>1, M);
    1440                 :            :   else
    1441                 :   32912448 :     return subii(M,Zf_shift(a, (s-hord)>>1, M));
    1442                 :            : }
    1443                 :            : 
    1444                 :            : INLINE GEN
    1445                 :     373248 : Zf_sqr(GEN a, GEN M) { return Zf_red_destroy(sqri(a), M); }
    1446                 :            : 
    1447                 :            : INLINE GEN
    1448                 :    1778304 : 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                 :    3922878 : muliifft_dit(ulong o, ulong ord, GEN M, GEN FFT, long d, long step)
    1453                 :            : {
    1454                 :    3922878 :   pari_sp av = avma;
    1455                 :            :   long i;
    1456                 :    3922878 :   ulong j, no = (o<<1)%ord;
    1457                 :    3922878 :   long hstep=step>>1;
    1458         [ +  + ]:   25174590 :   for (i = d+1, j = 0; i <= d+hstep; ++i, j =(j+o)%ord)
    1459                 :            :   {
    1460                 :   21251712 :     GEN a = Zf_add(gel(FFT,i), gel(FFT,i+hstep), M);
    1461                 :   21251712 :     GEN b = Zf_mulsqrt2(Zf_sub(gel(FFT,i), gel(FFT,i+hstep), M), j, ord, M);
    1462                 :   21251712 :     affii(a,gel(FFT,i));
    1463                 :   21251712 :     affii(b,gel(FFT,i+hstep));
    1464                 :   21251712 :     avma = av;
    1465                 :            :   }
    1466         [ +  + ]:    3922878 :   if (hstep>1)
    1467                 :            :   {
    1468                 :    1957950 :     muliifft_dit(no, ord, M, FFT, d, hstep);
    1469                 :    1957950 :     muliifft_dit(no, ord, M, FFT, d+hstep, hstep);
    1470                 :            :   }
    1471                 :    3922878 : }
    1472                 :            : 
    1473                 :            : /* In place, bit reversed FFT, inverse of muliifft_dit */
    1474                 :            : static void
    1475                 :    2147817 : muliifft_dis(ulong o, ulong ord, GEN M, GEN FFT, long d, long step)
    1476                 :            : {
    1477                 :    2147817 :   pari_sp av = avma;
    1478                 :            :   long i;
    1479                 :    2147817 :   ulong j, no = (o<<1)%ord;
    1480                 :    2147817 :   long hstep=step>>1;
    1481         [ +  + ]:    2147817 :   if (hstep>1)
    1482                 :            :   {
    1483                 :    1072041 :     muliifft_dis(no, ord, M, FFT, d, hstep);
    1484                 :    1072041 :     muliifft_dis(no, ord, M, FFT, d+hstep, hstep);
    1485                 :            :   }
    1486         [ +  + ]:   13808553 :   for (i = d+1, j = 0; i <= d+hstep; ++i, j =(j+o)%ord)
    1487                 :            :   {
    1488                 :   11660736 :     GEN z = Zf_mulsqrt2(gel(FFT,i+hstep), j, ord, M);
    1489                 :   11660736 :     GEN a = Zf_add(gel(FFT,i), z, M);
    1490                 :   11660736 :     GEN b = Zf_sub(gel(FFT,i), z, M);
    1491                 :   11660736 :     affii(a,gel(FFT,i));
    1492                 :   11660736 :     affii(b,gel(FFT,i+hstep));
    1493                 :   11660736 :     avma = av;
    1494                 :            :   }
    1495                 :    2147817 : }
    1496                 :            : 
    1497                 :            : static GEN
    1498                 :       6978 : muliifft_spliti(GEN a, long na, long bs, long n, long mod)
    1499                 :            : {
    1500                 :       6978 :   GEN ap = a+na-1;
    1501                 :       6978 :   GEN c  = cgetg(n+1, t_VEC);
    1502                 :            :   long i,j;
    1503         [ +  + ]:    3936834 :   for(i=1;i<=n;i++)
    1504                 :            :   {
    1505                 :    3929856 :     GEN z = cgeti(mod+3);
    1506         [ +  + ]:    3929856 :     if (na)
    1507                 :            :     {
    1508                 :    1934364 :       long m = minss(bs, na), v=0;
    1509                 :    1934364 :       GEN zp, aa=ap-m+1;
    1510 [ +  + ][ +  + ]:   26731164 :       while (!*aa && v<m) {aa++; v++;}
    1511                 :    1934364 :       zp = z+m-v+1;
    1512         [ +  + ]:   32466684 :       for (j=v; j < m; j++)
    1513                 :   30532320 :         *zp-- = *ap--;
    1514                 :    1934364 :       ap -= v; na -= m;
    1515         [ +  + ]:    1934364 :       z[1] = evalsigne(m!=v) | evallgefint(m-v+2);
    1516                 :            :     } else
    1517                 :    1995492 :       z[1] = evalsigne(0) | evallgefint(2);
    1518                 :    3929856 :     gel(c, i) = z;
    1519                 :            :   }
    1520                 :       6978 :   return c;
    1521                 :            : }
    1522                 :            : 
    1523                 :            : static GEN
    1524                 :       3735 : muliifft_unspliti(GEN V, long bs, long len)
    1525                 :            : {
    1526                 :       3735 :   long s, i, j, l = lg(V);
    1527                 :       3735 :   GEN a = cgeti(len);
    1528                 :       3735 :   a[1] = evalsigne(1)|evallgefint(len);
    1529         [ +  + ]:   60422082 :   for(i=2;i<len;i++)
    1530                 :   60418347 :     a[i] = 0;
    1531         [ +  + ]:    2155287 :   for(i=1, s=0; i<l; i++, s+=bs)
    1532                 :            :   {
    1533                 :    2151552 :     GEN  u = gel(V,i);
    1534         [ +  + ]:    2151552 :     if (signe(u))
    1535                 :            :     {
    1536                 :    1969188 :       GEN ap = int_W(a,s);
    1537                 :    1969188 :       GEN up = int_LSW(u);
    1538                 :    1969188 :       long lu = NLIMBS(u);
    1539                 :            :       LOCAL_OVERFLOW;
    1540                 :    1969188 :       *ap = addll(*ap, *up--); ap--;
    1541 [ +  + ][ +  + ]:  109601475 :       for (j=1; j<lu; j++)
    1542                 :  107632287 :        { *ap = addllx(*ap, *up--); ap--; }
    1543         [ +  + ]:    1969488 :       while (overflow)
    1544                 :        300 :        { *ap = addll(*ap, 1); ap--; }
    1545                 :            :     }
    1546                 :            :   }
    1547                 :       3735 :   return int_normalize(a,0);
    1548                 :            : }
    1549                 :            : 
    1550                 :            : static GEN
    1551                 :        492 : sqrispec_fft(GEN a, long na)
    1552                 :            : {
    1553                 :        492 :   pari_sp av, ltop = avma;
    1554                 :        492 :   long len = 2*na;
    1555                 :            :   long k, mod, bs, n;
    1556                 :            :   GEN  FFT, M;
    1557                 :            :   long i;
    1558                 :            :   ulong o, ord;
    1559                 :            : 
    1560                 :        492 :   mulliifft_params(len,&k,&mod,&bs,&n,&ord);
    1561                 :        492 :   o = ord>>k;
    1562                 :        492 :   M = int2n(mod*BITS_IN_LONG);
    1563                 :        492 :   M[2+mod] = 1;
    1564                 :        492 :   FFT = muliifft_spliti(a, na, bs, n, mod);
    1565                 :        492 :   muliifft_dit(o, ord, M, FFT, 0, n);
    1566                 :        492 :   av = avma;
    1567         [ +  + ]:     373740 :   for(i=1; i<=n; i++)
    1568                 :            :   {
    1569                 :     373248 :     affii(Zf_sqr(gel(FFT,i), M), gel(FFT,i));
    1570                 :     373248 :     avma=av;
    1571                 :            :   }
    1572                 :        492 :   muliifft_dis(ord-o, ord, M, FFT, 0, n);
    1573         [ +  + ]:     373740 :   for(i=1; i<=n; i++)
    1574                 :            :   {
    1575                 :     373248 :     affii(Zf_shift(gel(FFT,i), (ord>>1)-k, M), gel(FFT,i));
    1576                 :     373248 :     avma=av;
    1577                 :            :   }
    1578                 :        492 :   return gerepileuptoint(ltop, muliifft_unspliti(FFT,bs,2+len));
    1579                 :            : }
    1580                 :            : 
    1581                 :            : static GEN
    1582                 :       3243 : muliispec_fft(GEN a, GEN b, long na, long nb)
    1583                 :            : {
    1584                 :       3243 :   pari_sp av, av2, ltop = avma;
    1585                 :       3243 :   long len = na+nb;
    1586                 :            :   long k, mod, bs, n;
    1587                 :            :   GEN FFT, FFTb, M;
    1588                 :            :   long i;
    1589                 :            :   ulong o, ord;
    1590                 :            : 
    1591                 :       3243 :   mulliifft_params(len,&k,&mod,&bs,&n,&ord);
    1592                 :       3243 :   o = ord>>k;
    1593                 :       3243 :   M = int2n(mod*BITS_IN_LONG);
    1594                 :       3243 :   M[2+mod] = 1;
    1595                 :       3243 :   FFT = muliifft_spliti(a, na, bs, n, mod);
    1596                 :       3243 :   av=avma;
    1597                 :       3243 :   muliifft_dit(o, ord, M, FFT, 0, n);
    1598                 :       3243 :   FFTb = muliifft_spliti(b, nb, bs, n, mod);
    1599                 :       3243 :   av2 = avma;
    1600                 :       3243 :   muliifft_dit(o, ord, M, FFTb, 0, n);
    1601         [ +  + ]:    1781547 :   for(i=1; i<=n; i++)
    1602                 :            :   {
    1603                 :    1778304 :     affii(Zf_mul(gel(FFT,i), gel(FFTb,i), M), gel(FFT,i));
    1604                 :    1778304 :     avma=av2;
    1605                 :            :   }
    1606                 :       3243 :   avma=av;
    1607                 :       3243 :   muliifft_dis(ord-o, ord, M, FFT, 0, n);
    1608         [ +  + ]:    1781547 :   for(i=1; i<=n; i++)
    1609                 :            :   {
    1610                 :    1778304 :     affii(Zf_shift(gel(FFT,i),(ord>>1)-k,M), gel(FFT,i));
    1611                 :    1778304 :     avma=av;
    1612                 :            :   }
    1613                 :       3243 :   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                 :  124769085 : addshiftw(GEN x, GEN y, long d)
    1625                 :            : {
    1626                 :  124769085 :   GEN z,z0,y0,yd, zd = (GEN)avma;
    1627                 :  124769085 :   long a,lz,ly = lgefint(y);
    1628                 :            : 
    1629                 :  124769085 :   z0 = new_chunk(d);
    1630                 :  124769085 :   a = ly-2; yd = y+ly;
    1631         [ +  + ]:  124769085 :   if (a >= d)
    1632                 :            :   {
    1633         [ +  + ]: 2872837329 :     y0 = yd-d; while (yd > y0) *--zd = *--yd; /* copy last d words of y */
    1634                 :  121392096 :     a -= d;
    1635         [ +  + ]:  121392096 :     if (a)
    1636                 :   96422628 :       z = addiispec(LIMBS(x), LIMBS(y), NLIMBS(x), a);
    1637                 :            :     else
    1638                 :   24969468 :       z = icopy(x);
    1639                 :            :   }
    1640                 :            :   else
    1641                 :            :   {
    1642         [ +  + ]:   11707467 :     y0 = yd-a; while (yd > y0) *--zd = *--yd; /* copy last a words of y */
    1643         [ +  + ]:   74002869 :     while (zd > z0) *--zd = 0;    /* complete with 0s */
    1644                 :    3376989 :     z = icopy(x);
    1645                 :            :   }
    1646                 :  124769085 :   lz = lgefint(z)+d;
    1647                 :  124769085 :   z[1] = evalsigne(1) | evallgefint(lz);
    1648                 :  124769085 :   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                 : 1249793925 : 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         [ +  + ]: 1249793925 :   if (na < nb) swapspec(a,b, na,nb);
    1662         [ +  + ]: 1249793925 :   if (nb < MULII_KARATSUBA_LIMIT) return muliispec_basecase(a,b,na,nb);
    1663         [ +  + ]:   48352032 :   if (nb >= MULII_FFT_LIMIT)      return muliispec_fft(a,b,na,nb);
    1664                 :   48348789 :   i=(na>>1); n0=na-i; na=i;
    1665                 :   48348789 :   av=avma; a0=a+na; n0a=n0;
    1666 [ +  + ][ +  + ]:  129885033 :   while (n0a && !*a0) { a0++; n0a--; }
    1667                 :            : 
    1668 [ +  + ][ +  + ]:   48348789 :   if (n0a && nb > n0)
    1669                 :   45832614 :   { /* nb <= na <= n0 */
    1670                 :            :     GEN b0,c1,c2;
    1671                 :            :     long n0b;
    1672                 :            : 
    1673                 :   45832614 :     nb -= n0;
    1674                 :   45832614 :     c = muliispec(a,b,na,nb);
    1675                 :   45832614 :     b0 = b+nb; n0b = n0;
    1676 [ +  + ][ +  + ]:  111317766 :     while (n0b && !*b0) { b0++; n0b--; }
    1677         [ +  + ]:   45832614 :     if (n0b)
    1678                 :            :     {
    1679                 :   44652369 :       c0 = muliispec(a0,b0, n0a,n0b);
    1680                 :            : 
    1681                 :   44652369 :       c2 = addiispec(a0,a, n0a,na);
    1682                 :   44652369 :       c1 = addiispec(b0,b, n0b,nb);
    1683                 :   44652369 :       c1 = muliispec(LIMBS(c1),LIMBS(c2), NLIMBS(c1),NLIMBS(c2));
    1684                 :   44652369 :       c2 = addiispec(LIMBS(c0),LIMBS(c),  NLIMBS(c0),NLIMBS(c));
    1685                 :            : 
    1686                 :   44652369 :       c1 = subiispec(LIMBS(c1),LIMBS(c2), NLIMBS(c1),NLIMBS(c2));
    1687                 :            :     }
    1688                 :            :     else
    1689                 :            :     {
    1690                 :    1180245 :       c0 = gen_0;
    1691                 :    1180245 :       c1 = muliispec(a0,b, n0a,nb);
    1692                 :            :     }
    1693                 :   45832614 :     c = addshiftw(c,c1, n0);
    1694                 :            :   }
    1695                 :            :   else
    1696                 :            :   {
    1697                 :    2516175 :     c = muliispec(a,b,na,nb);
    1698                 :    2516175 :     c0 = muliispec(a0,b,n0a,nb);
    1699                 :            :   }
    1700                 : 1249793925 :   return gerepileuptoint(av, addshiftw(c,c0, n0));
    1701                 :            : }
    1702                 :            : GEN
    1703                 :     189558 : muluui(ulong x, ulong y, GEN z)
    1704                 :            : {
    1705                 :     189558 :   long t, s = signe(z);
    1706                 :            :   GEN r;
    1707                 :            :   LOCAL_HIREMAINDER;
    1708                 :            : 
    1709 [ +  - ][ +  - ]:     189558 :   if (!x || !y || !signe(z)) return gen_0;
                 [ +  + ]
    1710                 :     189261 :   t = mulll(x,y);
    1711 [ +  - ][ +  - ]:     189261 :   if (!hiremainder)
    1712                 :     189261 :     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                 :     189558 :   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                 :    6402387 : remi2n(GEN x, long n)
    1729                 :            : {
    1730                 :    6402387 :   long hi,l,k,lx,ly, sx = signe(x);
    1731                 :            :   GEN z, xd, zd;
    1732                 :            : 
    1733 [ +  + ][ -  + ]:    6402387 :   if (!sx || !n) return gen_0;
    1734                 :            : 
    1735                 :    6384243 :   k = dvmdsBIL(n, &l);
    1736                 :    6384243 :   lx = lgefint(x);
    1737         [ +  + ]:    6384243 :   if (lx < k+3) return icopy(x);
    1738                 :            : 
    1739                 :    6353439 :   xd = x + (lx-k-1);
    1740                 :            :   /* x = |_|...|#|1|...|k| : copy the last l bits of # and the last k words
    1741                 :            :    *            ^--- initial xd  */
    1742                 :    6353439 :   hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
    1743         [ +  + ]:    6353439 :   if (!hi)
    1744                 :            :   { /* strip leading zeroes from result */
    1745 [ +  + ][ +  + ]:     249387 :     xd++; while (k && !*xd) { k--; xd++; }
    1746         [ +  + ]:     246552 :     if (!k) return gen_0;
    1747                 :      49359 :     ly = k+2; xd--;
    1748                 :            :   }
    1749                 :            :   else
    1750                 :    6106887 :     ly = k+3;
    1751                 :            : 
    1752                 :    6156246 :   zd = z = cgeti(ly);
    1753                 :    6156246 :   *++zd = evalsigne(sx) | evallgefint(ly);
    1754         [ +  + ]:    6156246 :   if (hi) *++zd = hi;
    1755         [ +  + ]:   11140452 :   for ( ;k; k--) *++zd = *++xd;
    1756                 :    6402387 :   return z;
    1757                 :            : }
    1758                 :            : 
    1759                 :            : GEN
    1760                 :  467147943 : sqrispec(GEN a, long na)
    1761                 :            : {
    1762                 :            :   GEN a0,c;
    1763                 :            :   long n0, n0a, i;
    1764                 :            :   pari_sp av;
    1765                 :            : 
    1766         [ +  + ]:  467147943 :   if (na < SQRI_KARATSUBA_LIMIT) return sqrispec_basecase(a,na);
    1767         [ +  + ]:    2627751 :   if (na >= SQRI_FFT_LIMIT) return sqrispec_fft(a,na);
    1768                 :    2627259 :   i=(na>>1); n0=na-i; na=i;
    1769                 :    2627259 :   av=avma; a0=a+na; n0a=n0;
    1770 [ +  + ][ +  + ]:    4696131 :   while (n0a && !*a0) { a0++; n0a--; }
    1771                 :    2627259 :   c = sqrispec(a,na);
    1772         [ +  + ]:    2627259 :   if (n0a)
    1773                 :            :   {
    1774                 :    2595921 :     GEN t, c1, c0 = sqrispec(a0,n0a);
    1775                 :            : #if 0
    1776                 :            :     c1 = shifti(muliispec(a0,a, n0a,na),1);
    1777                 :            : #else /* faster */
    1778                 :    2595921 :     t = addiispec(a0,a,n0a,na);
    1779                 :    2595921 :     t = sqrispec(LIMBS(t),NLIMBS(t));
    1780                 :    2595921 :     c1= addiispec(LIMBS(c0),LIMBS(c), NLIMBS(c0), NLIMBS(c));
    1781                 :    2595921 :     c1= subiispec(LIMBS(t),LIMBS(c1), NLIMBS(t), NLIMBS(c1));
    1782                 :            : #endif
    1783                 :    2595921 :     c = addshiftw(c,c1, n0);
    1784                 :    2595921 :     c = addshiftw(c,c0, n0);
    1785                 :            :   }
    1786                 :            :   else
    1787                 :      31338 :     c = addshiftw(c,gen_0,n0<<1);
    1788                 :  467147943 :   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                 :   11723781 : p_sqrtu1(ulong *N, ulong *ps, ulong *pr)
    1819                 :            : {
    1820                 :   11723781 :   ulong prec, r, s, q, u, n0 = N[0];
    1821                 :            : 
    1822                 :   11723781 :   q = n0 >> (BITS_IN_LONG - 8);
    1823                 :            :   /* 2^6 = 64 <= q < 256 = 2^8 */
    1824                 :   11723781 :   s = approx_tab[q - 64];                                /* 128 <= s < 255 */
    1825                 :   11723781 :   r = (n0 >> (BITS_IN_LONG - 16)) - s * s;                /* r <= 2*s */
    1826         [ +  + ]:   11723781 :   if (r > (s << 1)) { r -= (s << 1) | 1; s++; }
    1827                 :            : 
    1828                 :            :   /* 8-bit approximation from the high 8-bits of N[0] */
    1829                 :   11723781 :   prec = 8;
    1830                 :   11723781 :   n0 <<= 2 * prec;
    1831         [ +  + ]:   35171343 :   while (2 * prec < BITS_IN_LONG)
    1832                 :            :   { /* invariant: s has prec bits, and r <= 2*s */
    1833                 :   23447562 :     r = (r << prec) + (n0 >> (BITS_IN_LONG - prec));
    1834                 :   23447562 :     n0 <<= prec;
    1835                 :   23447562 :     u = 2 * s;
    1836                 :   23447562 :     q = r / u; u = r - q * u;
    1837                 :   23447562 :     s = (s << prec) + q;
    1838                 :   23447562 :     u = (u << prec) + (n0 >> (BITS_IN_LONG - prec));
    1839                 :   23447562 :     q = q * q;
    1840                 :   23447562 :     r = u - q;
    1841         [ +  + ]:   23447562 :     if (u < q) { s--; r += (s << 1) | 1; }
    1842                 :   23447562 :     n0 <<= prec;
    1843                 :   23447562 :     prec = 2 * prec;
    1844                 :            :   }
    1845                 :   11723781 :   *ps = s;
    1846                 :   11723781 :   *pr = r;
    1847                 :   11723781 : }
    1848                 :            : 
    1849                 :            : /* N[0..1], assume N[0] >= 2^(BIL-2).
    1850                 :            :  * Return 1 if remainder overflows, 0 otherwise */
    1851                 :            : static int
    1852                 :   10498788 : p_sqrtu2(ulong *N, ulong *ps, ulong *pr)
    1853                 :            : {
    1854                 :   10498788 :   ulong cc, qhl, r, s, q, u, n1 = N[1];
    1855                 :            :   LOCAL_OVERFLOW;
    1856                 :            : 
    1857                 :   10498788 :   p_sqrtu1(N, &s, &r); /* r <= 2s */
    1858         [ +  + ]:   15801066 :   qhl = 0; while (r >= s) { qhl++; r -= s; }
    1859                 :            :   /* now r < s < 2^(BIL/2) */
    1860                 :   10498788 :   r = (r << BITS_IN_HALFULONG) | (n1 >> BITS_IN_HALFULONG);
    1861                 :   10498788 :   u = s << 1;
    1862                 :   10498788 :   q = r / u; u = r - q * u;
    1863                 :   10498788 :   q += (qhl & 1) << (BITS_IN_HALFULONG - 1);
    1864                 :   10498788 :   qhl >>= 1;
    1865                 :            :   /* (initial r)<<(BIL/2) + n1>>(BIL/2) = (qhl<<(BIL/2) + q) * 2s + u */
    1866                 :   10498788 :   s = ((s + qhl) << BITS_IN_HALFULONG) + q;
    1867                 :   10498788 :   cc = u >> BITS_IN_HALFULONG;
    1868                 :   10498788 :   r = (u << BITS_IN_HALFULONG) | (n1 & LOWMASK);
    1869                 :   10498788 :   r = subll(r, q * q);
    1870                 :   10498788 :   cc -= overflow + qhl;
    1871                 :            :   /* now subtract 2*q*2^(BIL/2) + 2^BIL if qhl is set */
    1872 [ +  + ][ +  + ]:   10498788 :   if ((long)cc < 0)
    1873                 :            :   {
    1874         [ +  + ]:    2736033 :     if (s) {
    1875                 :    2732061 :       r = addll(r, s);
    1876                 :    2732061 :       cc += overflow;
    1877                 :    2732061 :       s--;
    1878                 :            :     } else {
    1879                 :       3972 :       cc++;
    1880                 :       3972 :       s = ~0UL;
    1881                 :            :     }
    1882                 :    2736033 :     r = addll(r, s);
    1883                 :    2736033 :     cc += overflow;
    1884                 :            :   }
    1885                 :   10498788 :   *ps = s;
    1886                 :   10498788 :   *pr = r; return cc;
    1887                 :            : }
    1888                 :            : 
    1889                 :            : static void
    1890                 :   10487391 : xmpn_zero(GEN x, long n)
    1891                 :            : {
    1892         [ +  + ]:   76537752 :   while (--n >= 0) x[n]=0;
    1893                 :   10487391 : }
    1894                 :            : static void
    1895                 :  117103035 : xmpn_copy(GEN z, GEN x, long n)
    1896                 :            : {
    1897                 :  117103035 :   long k = n;
    1898         [ +  + ]:  463168673 :   while (--k >= 0) z[k] = x[k];
    1899                 :  117103035 : }
    1900                 :            : /* a[0..la-1] * 2^(lb BIL) | b[0..lb-1] */
    1901                 :            : static GEN
    1902                 :   50729004 : catii(GEN a, long la, GEN b, long lb)
    1903                 :            : {
    1904                 :   50729004 :   long l = la + lb + 2;
    1905                 :   50729004 :   GEN z = cgetipos(l);
    1906                 :   50729004 :   xmpn_copy(LIMBS(z), a, la);
    1907                 :   50729004 :   xmpn_copy(LIMBS(z) + la, b, lb);
    1908                 :   50729004 :   return int_normalize(z, 0);
    1909                 :            : }
    1910                 :            : 
    1911                 :            : /* sqrt n[0..1], assume n normalized */
    1912                 :            : static GEN
    1913                 :   10495065 : sqrtispec2(GEN n, GEN *pr)
    1914                 :            : {
    1915                 :            :   ulong s, r;
    1916                 :   10495065 :   int hi = p_sqrtu2((ulong*)n, &s, &r);
    1917                 :   10495065 :   GEN S = utoi(s);
    1918         [ +  + ]:   10495065 :   *pr = hi? uutoi(1,r): utoi(r);
    1919                 :   10495065 :   return S;
    1920                 :            : }
    1921                 :            : 
    1922                 :            : /* sqrt n[0], _dont_ assume n normalized */
    1923                 :            : static GEN
    1924                 :    1224993 : sqrtispec1_sh(GEN n, GEN *pr)
    1925                 :            : {
    1926                 :            :   GEN S;
    1927                 :    1224993 :   ulong r, s, u0 = uel(n,0);
    1928 [ +  + ][ +  + ]:    1224993 :   int sh = bfffo(u0) & ~1UL;
         [ +  + ][ +  + ]
    1929 [ +  + ][ +  + ]:    1224993 :   if (sh) u0 <<= sh;
    1930                 :    1224993 :   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         [ +  + ]:    1224993 :   if (sh) {
    1935                 :    1224933 :     int k = sh >> 1;
    1936                 :    1224933 :     ulong s0 = s & ((1L<<k) - 1);
    1937                 :    1224933 :     r += s * (s0<<1);
    1938                 :    1224933 :     s >>= k;
    1939                 :    1224933 :     r >>= sh;
    1940                 :            :   }
    1941                 :    1224993 :   S = utoi(s);
    1942         [ +  + ]:    1224993 :   if (pr) *pr = utoi(r);
    1943                 :    1224993 :   return S;
    1944                 :            : }
    1945                 :            : 
    1946                 :            : /* sqrt n[0..1], _dont_ assume n normalized */
    1947                 :            : static GEN
    1948                 :       3723 : sqrtispec2_sh(GEN n, GEN *pr)
    1949                 :            : {
    1950                 :            :   GEN S;
    1951                 :       3723 :   ulong U[2], r, s, u0 = uel(n,0), u1 = uel(n,1);
    1952 [ +  + ][ +  + ]:       3723 :   int hi, sh = bfffo(u0) & ~1UL;
         [ +  + ][ +  + ]
    1953 [ +  + ][ +  + ]:       3723 :   if (sh) {
    1954                 :       3639 :     u0 = (u0 << sh) | (u1 >> (BITS_IN_LONG-sh));
    1955                 :       3639 :     u1 <<= sh;
    1956                 :            :   }
    1957                 :       3723 :   U[0] = u0;
    1958                 :       3723 :   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         [ +  + ]:       3723 :   if (sh) {
    1963                 :       3639 :     int k = sh >> 1;
    1964                 :       3639 :     ulong s0 = s & ((1L<<k) - 1);
    1965                 :            :     LOCAL_HIREMAINDER;
    1966                 :            :     LOCAL_OVERFLOW;
    1967                 :       3639 :     r = addll(r, mulll(s, (s0<<1)));
    1968 [ +  + ][ +  + ]:       3639 :     if (overflow) hiremainder++;
    1969                 :       3639 :     hiremainder += hi; /* + 0 or 1 */
    1970                 :       3639 :     s >>= k;
    1971                 :       3639 :     r = (r>>sh) | (hiremainder << (BITS_IN_LONG-sh));
    1972                 :       3639 :     hi = (hiremainder & (1L<<sh));
    1973                 :            :   }
    1974                 :       3723 :   S = utoi(s);
    1975 [ +  + ][ -  + ]:       3723 :   if (pr) *pr = hi? uutoi(1,r): utoi(r);
    1976                 :       3723 :   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                 :   35859567 : sqrtispec(GEN N, long n, GEN *r)
    1983                 :            : {
    1984                 :            :   GEN S, R, q, z, u;
    1985                 :            :   long l, h;
    1986                 :            : 
    1987         [ +  + ]:   35859567 :   if (n == 1) return sqrtispec2(N, r);
    1988                 :   25364502 :   l = n >> 1;
    1989                 :   25364502 :   h = n - l; /* N = a3(h) | a2(h) | a1(l) | a0(l words) */
    1990                 :   25364502 :   S = sqrtispec(N, h, &R); /* S^2 + R = a3|a2 */
    1991                 :            : 
    1992                 :   25364502 :   z = catii(LIMBS(R), NLIMBS(R), N + 2*h, l); /* = R | a1(l) */
    1993                 :   25364502 :   q = dvmdii(z, shifti(S,1), &u);
    1994                 :   25364502 :   z = catii(LIMBS(u), NLIMBS(u), N + n + h, l); /* = u | a0(l) */
    1995                 :            : 
    1996                 :   25364502 :   S = addshiftw(S, q, l);
    1997                 :   25364502 :   R = subii(z, sqri(q));
    1998         [ +  + ]:   25364502 :   if (signe(R) < 0)
    1999                 :            :   {
    2000                 :    4022721 :     GEN S2 = shifti(S,1);
    2001                 :    4022721 :     R = addis(subiispec(LIMBS(S2),LIMBS(R), NLIMBS(S2),NLIMBS(R)), -1);
    2002                 :    4022721 :     S = addis(S, -1);
    2003                 :            :   }
    2004                 :   35859567 :   *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                 :    1236390 : sqrtremi(GEN N, GEN *r)
    2012                 :            : {
    2013                 :            :   pari_sp av;
    2014                 :    1236390 :   GEN S, R, n = N+2;
    2015                 :    1236390 :   long k, l2, ln = NLIMBS(N);
    2016                 :            :   int sh;
    2017                 :            : 
    2018         [ +  + ]:    1236390 :   if (ln <= 2)
    2019                 :            :   {
    2020         [ +  + ]:    1228716 :     if (ln == 2) return sqrtispec2_sh(n, r);
    2021         [ +  - ]:    1224993 :     if (ln == 1) return sqrtispec1_sh(n, r);
    2022         [ #  # ]:          0 :     if (r) *r = gen_0;
    2023                 :          0 :     return gen_0;
    2024                 :            :   }
    2025                 :       7674 :   av = avma;
    2026 [ +  + ][ +  + ]:       7674 :   sh = bfffo(n[0]) >> 1;
         [ +  + ][ +  + ]
    2027                 :       7674 :   l2 = (ln + 1) >> 1;
    2028 [ +  + ][ +  + ]:      15285 :   if (sh || (ln & 1)) { /* normalize n, so that n[0] >= 2^BIL / 4 */
         [ +  + ][ +  + ]
    2029                 :       7611 :     GEN s0, t = new_chunk(ln + 1);
    2030                 :       7611 :     t[ln] = 0;
    2031         [ +  + ]:       7611 :     if (sh)
    2032                 :       7446 :       shift_left(t, n, 0,ln-1, 0, sh << 1);
    2033                 :            :     else
    2034                 :        165 :       xmpn_copy(t, n, ln);
    2035                 :       7611 :     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                 :       7611 :     k = sh + (ln & 1) * (BITS_IN_LONG/2);
    2040                 :       7611 :     s0 = remi2n(S, k);
    2041                 :       7611 :     R = addii(shifti(R,-1), mulii(s0, S));
    2042                 :       7611 :     R = shifti(R, 1 - (k<<1));
    2043                 :       7611 :     S = shifti(S, -k);
    2044                 :            :   }
    2045                 :            :   else
    2046                 :         63 :     S = sqrtispec(n, l2, &R);
    2047                 :            : 
    2048         [ +  + ]:       7674 :   if (!r) { avma = (pari_sp)S; return gerepileuptoint(av, S); }
    2049                 :    1236390 :   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                 :   10487391 : sqrtr_abs(GEN x)
    2057                 :            : {
    2058                 :   10487391 :   long l = realprec(x) - 2, e = expo(x), er = e>>1;
    2059                 :   10487391 :   GEN b, c, res = cgetr(2 + l);
    2060                 :   10487391 :   res[1] = evalsigne(1) | evalexpo(er);
    2061         [ +  + ]:   10487391 :   if (e&1) {
    2062                 :    5157471 :     b = new_chunk(l << 1);
    2063                 :    5157471 :     xmpn_copy(b, x+2, l);
    2064                 :    5157471 :     xmpn_zero(b + l,l);
    2065                 :    5157471 :     b = sqrtispec(b, l, &c);
    2066                 :    5157471 :     xmpn_copy(res+2, b+2, l);
    2067         [ +  + ]:    5157471 :     if (cmpii(c, b) > 0) roundr_up_ip(res, l+2);
    2068                 :            :   } else {
    2069                 :            :     ulong u;
    2070                 :    5329920 :     b = new_chunk(2 + (l << 1));
    2071                 :    5329920 :     shift_left(b+1, x+2, 0,l-1, 0, BITS_IN_LONG-1);
    2072                 :    5329920 :     b[0] = uel(x,2)>>1;
    2073                 :    5329920 :     xmpn_zero(b + l+1,l+1);
    2074                 :    5329920 :     b = sqrtispec(b, l+1, &c);
    2075                 :    5329920 :     xmpn_copy(res+2, b+2, l);
    2076                 :    5329920 :     u = uel(b,l+2);
    2077 [ +  + ][ +  + ]:    5329920 :     if ( u&HIGHBIT || (u == ~HIGHBIT && cmpii(c,b) > 0))
                 [ +  - ]
    2078                 :    2587103 :       roundr_up_ip(res, l+2);
    2079                 :            :   }
    2080                 :   10487391 :   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                 :    1033620 : convi_dac(GEN x, ulong l, ulong *res)
    2159                 :            : {
    2160                 :    1033620 :   pari_sp ltop=avma;
    2161                 :            :   ulong m;
    2162                 :            :   GEN x1,x2;
    2163         [ +  + ]:    1506852 :   if (l==1) { *res=itou(x); return; }
    2164                 :     473232 :   m=l>>1;
    2165                 :     473232 :   x1=dvmdii(x,powuu(1000000000UL,m),&x2);
    2166                 :     473232 :   convi_dac(x1,l-m,res+m);
    2167                 :     473232 :   convi_dac(x2,m,res);
    2168                 :     473232 :   avma=ltop;
    2169                 :            : }
    2170                 :            : 
    2171                 :            : /* convert integer --> base 10^9 [not memory clean] */
    2172                 :            : ulong *
    2173                 :     226789 : convi(GEN x, long *l)
    2174                 :            : {
    2175                 :     226789 :   long lz, lx = lgefint(x);
    2176                 :            :   ulong *z;
    2177 [ +  + ][ +  + ]:     226789 :   if (lx == 3 && uel(x,2) < 1000000000UL) {
    2178                 :     139633 :     z = (ulong*)new_chunk(1);
    2179                 :     139633 :     *z = x[2];
    2180                 :     139633 :     *l = 1; return z+1;
    2181                 :            :   }
    2182                 :      87156 :   lz = 1 + (long)bit_accuracy_mul(lx, LOG10_2/9);
    2183                 :      87156 :   z = (ulong*)new_chunk(lz);
    2184                 :      87156 :   convi_dac(x,(ulong)lz,z);
    2185         [ +  + ]:     193761 :   while (z[lz-1]==0) lz--;
    2186                 :     226789 :   *l=lz; return z+lz;
    2187                 :            : }
    2188                 :            : 

Generated by: LCOV version 1.9