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/gmp - mp.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16861-9aeb453) Lines: 641 716 89.5 %
Date: 2014-10-08 Functions: 53 54 98.1 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 421 533 79.0 %

           Branch data     Line data    Source code
       1                 :            : #line 2 "../src/kernel/gmp/mp.c"
       2                 :            : /* Copyright (C) 2002-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                 :            : /**                               GMP KERNEL                          **/
      18                 :            : /** BA2002Sep24                                                       **/
      19                 :            : /***********************************************************************/
      20                 :            : /* GMP t_INT as just like normal t_INT, just the mantissa is the other way
      21                 :            :  * round
      22                 :            :  *
      23                 :            :  *   `How would you like to live in Looking-glass House, Kitty?  I
      24                 :            :  *   wonder if they'd give you milk in there?  Perhaps Looking-glass
      25                 :            :  *   milk isn't good to drink--But oh, Kitty! now we come to the
      26                 :            :  *   passage.  You can just see a little PEEP of the passage in
      27                 :            :  *   Looking-glass House, if you leave the door of our drawing-room
      28                 :            :  *   wide open:  and it's very like our passage as far as you can see,
      29                 :            :  *   only you know it may be quite different on beyond.  Oh, Kitty!
      30                 :            :  *   how nice it would be if we could only get through into Looking-
      31                 :            :  *   glass House!  I'm sure it's got, oh! such beautiful things in it!
      32                 :            :  *
      33                 :            :  *  Through the Looking Glass,  Lewis Carrol
      34                 :            :  *
      35                 :            :  *  (pityful attempt to beat GN code/comments rate)
      36                 :            :  *  */
      37                 :            : 
      38                 :            : #include <gmp.h>
      39                 :            : #include "pari.h"
      40                 :            : #include "paripriv.h"
      41                 :            : #include "../src/kernel/none/tune-gen.h"
      42                 :            : 
      43                 :            : /*We need PARI invmod renamed to invmod_pari*/
      44                 :            : #define INVMOD_PARI
      45                 :            : 
      46                 :          0 : static void *pari_gmp_realloc(void *ptr, size_t old_size, size_t new_size) {
      47                 :          0 :   (void)old_size; return (void *) pari_realloc(ptr,new_size);
      48                 :            : }
      49                 :            : 
      50                 :       3002 : static void pari_gmp_free(void *ptr, size_t old_size){
      51                 :       3002 :   (void)old_size; pari_free(ptr);
      52                 :       3002 : }
      53                 :            : 
      54                 :            : static void *(*old_gmp_malloc)(size_t new_size);
      55                 :            : static void *(*old_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
      56                 :            : static void (*old_gmp_free)(void *ptr, size_t old_size);
      57                 :            : 
      58                 :            : void
      59                 :        696 : pari_kernel_init(void)
      60                 :            : {
      61                 :        696 :   mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
      62                 :        696 :   mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
      63                 :        696 : }
      64                 :            : 
      65                 :            : void
      66                 :        684 : pari_kernel_close(void)
      67                 :            : {
      68                 :            :   void *(*new_gmp_malloc)(size_t new_size);
      69                 :            :   void *(*new_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
      70                 :            :   void (*new_gmp_free)(void *ptr, size_t old_size);
      71                 :        684 :   mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
      72         [ +  - ]:        684 :   if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
      73         [ +  - ]:        684 :   if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
      74         [ +  - ]:        684 :   if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
      75                 :        684 :   mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
      76                 :        684 : }
      77                 :            : 
      78                 :            : #define LIMBS(x)  ((mp_limb_t *)((x)+2))
      79                 :            : #define NLIMBS(x) (lgefint(x)-2)
      80                 :            : /*This one is for t_REALs to emphasize they are not t_INTs*/
      81                 :            : #define RLIMBS(x)  ((mp_limb_t *)((x)+2))
      82                 :            : #define RNLIMBS(x) (lg(x)-2)
      83                 :            : 
      84                 :            : INLINE void
      85                 :    1659731 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
      86                 :            : {
      87         [ +  + ]:   25766868 :   while (--n >= 0) x[n]=y[n];
      88                 :    1659731 : }
      89                 :            : 
      90                 :            : INLINE void
      91                 :  167354613 : xmpn_mirror(mp_limb_t *x, long n)
      92                 :            : {
      93                 :            :   long i;
      94         [ +  + ]:  945934755 :   for(i=0;i<(n>>1);i++)
      95                 :            :   {
      96                 :  778580142 :     ulong m=x[i];
      97                 :  778580142 :     x[i]=x[n-1-i];
      98                 :  778580142 :     x[n-1-i]=m;
      99                 :            :   }
     100                 :  167354613 : }
     101                 :            : 
     102                 :            : INLINE void
     103                 :  191287791 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
     104                 :            : {
     105                 :            :   long i;
     106         [ +  + ]: 2007114226 :   for(i=0;i<n;i++)
     107                 : 1815826435 :     z[i]=x[n-1-i];
     108                 :  191287791 : }
     109                 :            : 
     110                 :            : INLINE void
     111                 :   29540663 : xmpn_zero(mp_ptr x, mp_size_t n)
     112                 :            : {
     113         [ +  + ]:  213397903 :   while (--n >= 0) x[n]=0;
     114                 :   29540663 : }
     115                 :            : 
     116                 :            : INLINE GEN
     117                 :    7141756 : icopy_ef(GEN x, long l)
     118                 :            : {
     119                 :    7141756 :   register long lx = lgefint(x);
     120                 :    7141756 :   const GEN y = cgeti(l);
     121                 :            : 
     122         [ +  + ]:   41649104 :   while (--lx > 0) y[lx]=x[lx];
     123                 :    7141756 :   return y;
     124                 :            : }
     125                 :            : 
     126                 :            : /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
     127                 :            :  * GENs but pairs (long *a, long na) representing a list of digits (in basis
     128                 :            :  * BITS_IN_LONG) : a[0], ..., a[na-1]. [ In ordre to facilitate splitting: no
     129                 :            :  * need to reintroduce codewords ]
     130                 :            :  * Use speci(a,na) to visualize the corresponding GEN.
     131                 :            :  */
     132                 :            : 
     133                 :            : /***********************************************************************/
     134                 :            : /**                                                                   **/
     135                 :            : /**                     ADDITION / SUBTRACTION                        **/
     136                 :            : /**                                                                   **/
     137                 :            : /***********************************************************************/
     138                 :            : 
     139                 :            : GEN
     140                 :    1464759 : setloop(GEN a)
     141                 :            : {
     142                 :    1464759 :   pari_sp av = avma - 2 * sizeof(long);
     143                 :    1464759 :   (void)cgetg(lgefint(a) + 3, t_VECSMALL);
     144                 :    1464759 :   return icopy_avma(a, av); /* two cells of extra space after a */
     145                 :            : }
     146                 :            : 
     147                 :            : /* we had a = setloop(?), then some incloops. Reset a to b */
     148                 :            : GEN
     149                 :     145484 : resetloop(GEN a, GEN b) {
     150                 :     145484 :   a[0] = evaltyp(t_INT) | evallg(lgefint(b));
     151                 :     145484 :   affii(b, a); return a;
     152                 :            : }
     153                 :            : 
     154                 :            : /* assume a > 0, initialized by setloop. Do a++ */
     155                 :            : static GEN
     156                 :   17149718 : incpos(GEN a)
     157                 :            : {
     158                 :   17149718 :   long i, l = lgefint(a);
     159         [ +  + ]:   17149723 :   for (i=2; i<l; i++)
     160         [ +  + ]:   17149719 :     if (++a[i]) return a;
     161                 :          4 :   a[l] = 1; l++;
     162                 :          4 :   a[0]=evaltyp(t_INT) | _evallg(l);
     163                 :          4 :   a[1]=evalsigne(1) | evallgefint(l);
     164                 :   17149718 :   return a;
     165                 :            : }
     166                 :            : 
     167                 :            : /* assume a < 0, initialized by setloop. Do a++ */
     168                 :            : static GEN
     169                 :       5408 : incneg(GEN a)
     170                 :            : {
     171                 :       5408 :   long i, l = lgefint(a)-1;
     172         [ +  + ]:       5408 :   if (a[2]--)
     173                 :            :   {
     174         [ +  + ]:       5404 :     if (!a[l]) /* implies l = 2 */
     175                 :            :     {
     176                 :        420 :       a[0] = evaltyp(t_INT) | _evallg(2);
     177                 :        420 :       a[1] = evalsigne(0) | evallgefint(2);
     178                 :            :     }
     179                 :       5404 :     return a;
     180                 :            :   }
     181         [ +  - ]:          5 :   for (i=3; i<=l; i++)
     182         [ +  + ]:          5 :     if (a[i]--) break;
     183         [ +  - ]:          4 :   if (!a[l])
     184                 :            :   {
     185                 :          4 :     a[0] = evaltyp(t_INT) | _evallg(l);
     186                 :          4 :     a[1] = evalsigne(-1) | evallgefint(l);
     187                 :            :   }
     188                 :       5408 :   return a;
     189                 :            : }
     190                 :            : 
     191                 :            : /* assume a initialized by setloop. Do a++ */
     192                 :            : GEN
     193                 :   17163890 : incloop(GEN a)
     194                 :            : {
     195      [ +  +  + ]:   17163890 :   switch(signe(a))
     196                 :            :   {
     197                 :            :     case 0:
     198                 :       8764 :       a[0]=evaltyp(t_INT) | _evallg(3);
     199                 :       8764 :       a[1]=evalsigne(1) | evallgefint(3);
     200                 :       8764 :       a[2]=1; return a;
     201                 :       5408 :     case -1: return incneg(a);
     202                 :   17163890 :     default: return incpos(a);
     203                 :            :   }
     204                 :            : }
     205                 :            : 
     206                 :            : INLINE GEN
     207                 :  925532616 : adduispec(ulong s, GEN x, long nx)
     208                 :            : {
     209                 :            :   GEN  zd;
     210                 :            :   long lz;
     211                 :            : 
     212         [ +  + ]:  925532616 :   if (nx == 1) return adduu(uel(x,0), s);
     213                 :  248108079 :   lz = nx+3; zd = cgeti(lz);
     214         [ +  + ]:  248108079 :   if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
     215                 :    4860602 :     zd[lz-1]=1;
     216                 :            :   else
     217                 :  243247477 :     lz--;
     218                 :  248108079 :   zd[1] = evalsigne(1) | evallgefint(lz);
     219                 :  925532616 :   return zd;
     220                 :            : }
     221                 :            : 
     222                 :            : GEN
     223                 :  285598861 : adduispec_offset(ulong s, GEN x, long offset, long nx)
     224                 :            : {
     225                 :  285598861 :   GEN xd=x+2+offset;
     226 [ +  + ][ +  + ]:  374819079 :   while (nx && *(xd+nx-1)==0) nx--;
     227         [ +  + ]:  285598861 :   if (!nx) return utoi(s);
     228                 :  285598861 :   return adduispec(s,xd,nx);
     229                 :            : }
     230                 :            : 
     231                 :            : INLINE GEN
     232                 : 1079088644 : addiispec(GEN x, GEN y, long nx, long ny)
     233                 :            : {
     234                 :            :   GEN zd;
     235                 :            :   long lz;
     236                 :            : 
     237         [ +  + ]: 1079088644 :   if (nx < ny) swapspec(x,y, nx,ny);
     238         [ +  + ]: 1079088644 :   if (ny == 1) return adduispec(*y,x,nx);
     239                 :  447433249 :   lz = nx+3; zd = cgeti(lz);
     240                 :            : 
     241         [ +  + ]:  447433249 :   if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
     242                 :    6269690 :     zd[lz-1]=1;
     243                 :            :   else
     244                 :  441163559 :     lz--;
     245                 :            : 
     246                 :  447433249 :   zd[1] = evalsigne(1) | evallgefint(lz);
     247                 : 1079088644 :   return zd;
     248                 :            : }
     249                 :            : 
     250                 :            : /* assume x >= y */
     251                 :            : INLINE GEN
     252                 :  832692702 : subiuspec(GEN x, ulong s, long nx)
     253                 :            : {
     254                 :            :   GEN zd;
     255                 :            :   long lz;
     256                 :            : 
     257         [ +  + ]:  832692702 :   if (nx == 1) return utoi(x[0] - s);
     258                 :            : 
     259                 :   54140679 :   lz = nx + 2; zd = cgeti(lz);
     260                 :   54140679 :   mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
     261         [ +  + ]:   54140679 :   if (! zd[lz - 1]) { --lz; }
     262                 :            : 
     263                 :   54140679 :   zd[1] = evalsigne(1) | evallgefint(lz);
     264                 :  832692702 :   return zd;
     265                 :            : }
     266                 :            : 
     267                 :            : /* assume x > y */
     268                 :            : INLINE GEN
     269                 : 1301830723 : subiispec(GEN x, GEN y, long nx, long ny)
     270                 :            : {
     271                 :            :   GEN zd;
     272                 :            :   long lz;
     273         [ +  + ]: 1301830723 :   if (ny==1) return subiuspec(x,*y,nx);
     274                 :  471611294 :   lz = nx+2; zd = cgeti(lz);
     275                 :            : 
     276                 :  471611294 :   mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
     277 [ +  - ][ +  + ]:  537628073 :   while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
     278                 :            : 
     279                 :  471611294 :   zd[1] = evalsigne(1) | evallgefint(lz);
     280                 : 1301830723 :   return zd;
     281                 :            : }
     282                 :            : 
     283                 :            : static void
     284                 :   84021413 : roundr_up_ip(GEN x, long l)
     285                 :            : {
     286                 :   84021413 :   long i = l;
     287                 :            :   for(;;)
     288                 :            :   {
     289         [ +  + ]:   84195805 :     if (++((ulong*)x)[--i]) break;
     290         [ +  + ]:     205535 :     if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
     291                 :     174392 :   }
     292                 :   84021413 : }
     293                 :            : 
     294                 :            : void
     295                 :  100266345 : affir(GEN x, GEN y)
     296                 :            : {
     297                 :  100266345 :   const long s = signe(x), ly = lg(y);
     298                 :            :   long lx, sh, i;
     299                 :            : 
     300         [ +  + ]:  100266345 :   if (!s)
     301                 :            :   {
     302                 :     794588 :     y[1] = evalexpo(-bit_accuracy(ly));
     303                 :     794588 :     return;
     304                 :            :   }
     305                 :   99471757 :   lx = lgefint(x); sh = bfffo(*int_MSW(x));
     306                 :   99471757 :   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
     307         [ +  + ]:   99471757 :   if (sh) {
     308         [ +  + ]:   98114107 :     if (lx <= ly)
     309                 :            :     {
     310         [ +  + ]:  214681653 :       for (i=lx; i<ly; i++) y[i]=0;
     311                 :   65173881 :       mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
     312                 :   65173881 :       xmpn_mirror(LIMBS(y),lx-2);
     313                 :   65173881 :       return;
     314                 :            :     }
     315                 :   32940226 :     mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
     316                 :   32940226 :     y[2] |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
     317                 :   32940226 :     xmpn_mirror(LIMBS(y),ly-2);
     318                 :            :     /* lx > ly: round properly */
     319         [ +  + ]:   32940226 :     if ((x[lx-ly+1]<<sh) & HIGHBIT) roundr_up_ip(y, ly);
     320                 :            :   }
     321                 :            :   else {
     322                 :    1357650 :     GEN xd=int_MSW(x);
     323         [ +  + ]:    1357650 :     if (lx <= ly)
     324                 :            :     {
     325         [ +  + ]:    2381086 :       for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
     326         [ +  + ]:    2387449 :       for (   ; i<ly; i++) y[i]=0;
     327                 :     841163 :       return;
     328                 :            :     }
     329         [ +  + ]:    1569765 :     for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
     330                 :            :     /* lx > ly: round properly */
     331         [ +  + ]:  100266345 :     if (x[lx-ly+1] & HIGHBIT) roundr_up_ip(y, ly);
     332                 :            :   }
     333                 :            : }
     334                 :            : 
     335                 :            : INLINE GEN
     336                 :  202240135 : shiftispec(GEN x, long nx, long n)
     337                 :            : {
     338                 :            :   long ny,m;
     339                 :            :   GEN yd, y;
     340                 :            : 
     341         [ +  + ]:  202240135 :   if (!n) return icopyspec(x, nx);
     342                 :            : 
     343         [ +  + ]:  197642124 :   if (n > 0)
     344                 :            :   {
     345                 :  145636325 :     long d = dvmdsBIL(n, &m);
     346                 :            :     long i;
     347                 :            : 
     348                 :  145636325 :     ny = nx + d + (m!=0);
     349                 :  145636325 :     y = cgeti(ny + 2); yd = y + 2;
     350         [ +  + ]:  154971578 :     for (i=0; i<d; i++) yd[i] = 0;
     351                 :            : 
     352         [ +  + ]:  145636325 :     if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
     353                 :            :     else
     354                 :            :     {
     355                 :  144942854 :       ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
     356         [ +  + ]:  144942854 :       if (carryd) yd[ny - 1] = carryd;
     357                 :  142405805 :       else ny--;
     358                 :            :     }
     359                 :            :   }
     360                 :            :   else
     361                 :            :   {
     362                 :   52005799 :     long d = dvmdsBIL(-n, &m);
     363                 :            : 
     364                 :   52005799 :     ny = nx - d;
     365         [ +  + ]:   52005799 :     if (ny < 1) return gen_0;
     366                 :   51751121 :     y = cgeti(ny + 2); yd = y + 2;
     367                 :            : 
     368         [ +  + ]:   51751121 :     if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
     369                 :            :     else
     370                 :            :     {
     371                 :   50861198 :       mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
     372         [ +  + ]:   50861198 :       if (yd[ny - 1] == 0)
     373                 :            :       {
     374         [ +  + ]:    4966257 :         if (ny == 1) { avma = (pari_sp)(yd + 1); return gen_0; }
     375                 :    4341705 :         ny--;
     376                 :            :       }
     377                 :            :     }
     378                 :            :   }
     379                 :  196762894 :   y[1] = evalsigne(1)|evallgefint(ny + 2);
     380                 :  202240135 :   return y;
     381                 :            : }
     382                 :            : 
     383                 :            : GEN
     384                 :   15428113 : mantissa2nr(GEN x, long n)
     385                 :            : {
     386                 :   15428113 :   long ly, i, m, s = signe(x), lx = lg(x);
     387                 :            :   GEN y;
     388         [ +  + ]:   15428113 :   if (!s) return gen_0;
     389         [ +  + ]:   15426804 :   if (!n)
     390                 :            :   {
     391                 :    7411361 :     y = cgeti(lx);
     392                 :    7411361 :     y[1] = evalsigne(s) | evallgefint(lx);
     393                 :    7411361 :     xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
     394                 :    7411361 :     return y;
     395                 :            :   }
     396         [ +  + ]:    8015443 :   if (n > 0)
     397                 :            :   {
     398                 :     405672 :     GEN z = (GEN)avma;
     399                 :     405672 :     long d = dvmdsBIL(n, &m);
     400                 :            : 
     401                 :     405672 :     ly = lx+d; y = new_chunk(ly);
     402         [ +  + ]:     457103 :     for ( ; d; d--) *--z = 0;
     403 [ +  + ][ +  + ]:     406441 :     if (!m) for (i=2; i<lx; i++) y[i]=x[i];
     404                 :            :     else
     405                 :            :     {
     406                 :     405443 :       register const ulong sh = BITS_IN_LONG - m;
     407                 :     405443 :       shift_left(y,x, 2,lx-1, 0,m);
     408                 :     405443 :       i = uel(x,2) >> sh;
     409                 :            :       /* Extend y on the left? */
     410         [ +  - ]:     405443 :       if (i) { ly++; y = new_chunk(1); y[2] = i; }
     411                 :            :     }
     412                 :            :   }
     413                 :            :   else
     414                 :            :   {
     415                 :    7609771 :     ly = lx - dvmdsBIL(-n, &m);
     416         [ -  + ]:    7609771 :     if (ly<3) return gen_0;
     417                 :    7609771 :     y = new_chunk(ly);
     418         [ +  + ]:    7609771 :     if (m) {
     419                 :    7586507 :       shift_right(y,x, 2,ly, 0,m);
     420         [ -  + ]:    7586507 :       if (y[2] == 0)
     421                 :            :       {
     422         [ #  # ]:          0 :         if (ly==3) { avma = (pari_sp)(y+3); return gen_0; }
     423                 :          0 :         ly--; avma = (pari_sp)(++y);
     424                 :            :       }
     425                 :            :     } else {
     426         [ +  + ]:      73502 :       for (i=2; i<ly; i++) y[i]=x[i];
     427                 :            :     }
     428                 :            :   }
     429                 :    8015443 :   xmpn_mirror(LIMBS(y),ly-2);
     430                 :    8015443 :   y[1] = evalsigne(s)|evallgefint(ly);
     431                 :   15428113 :   y[0] = evaltyp(t_INT)|evallg(ly); return y;
     432                 :            : }
     433                 :            : 
     434                 :            : GEN
     435                 :     387006 : truncr(GEN x)
     436                 :            : {
     437                 :            :   long s, e, d, m, i;
     438                 :            :   GEN y;
     439 [ +  + ][ +  + ]:     387006 :   if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
     440                 :     208131 :   d = nbits2prec(e+1); m = remsBIL(e);
     441         [ -  + ]:     208131 :   if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
     442                 :            : 
     443                 :     208131 :   y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
     444         [ +  + ]:     208131 :   if (++m == BITS_IN_LONG)
     445         [ +  + ]:         41 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     446                 :            :   else
     447                 :            :   {
     448                 :     208111 :     GEN z=cgeti(d);
     449         [ +  + ]:     416687 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     450                 :     208111 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     451                 :     208111 :     avma=(pari_sp)y;
     452                 :            :   }
     453                 :     387006 :   return y;
     454                 :            : }
     455                 :            : 
     456                 :            : /* integral part */
     457                 :            : GEN
     458                 :     376516 : floorr(GEN x)
     459                 :            : {
     460                 :            :   long e, d, m, i, lx;
     461                 :            :   GEN y;
     462         [ +  + ]:     376516 :   if (signe(x) >= 0) return truncr(x);
     463         [ +  + ]:     137705 :   if ((e=expo(x)) < 0) return gen_m1;
     464                 :      60683 :   d = nbits2prec(e+1); m = remsBIL(e);
     465         [ -  + ]:      60683 :   lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
     466                 :      60683 :   y = cgeti(d+1);
     467         [ +  + ]:      60683 :   if (++m == BITS_IN_LONG)
     468                 :            :   {
     469         [ +  + ]:         40 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     470 [ +  - ][ -  + ]:         20 :     i=d; while (i<lx && !x[i]) i++;
     471         [ -  + ]:         20 :     if (i==lx) goto END;
     472                 :            :   }
     473                 :            :   else
     474                 :            :   {
     475                 :      60663 :     GEN z=cgeti(d);
     476         [ +  + ]:     122695 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     477                 :      60663 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     478         [ +  + ]:      60663 :     if (x[d-1]<<m == 0)
     479                 :            :     {
     480 [ +  + ][ +  + ]:      96451 :       i=d; while (i<lx && !x[i]) i++;
     481         [ +  + ]:      18417 :       if (i==lx) goto END;
     482                 :            :     }
     483                 :            :   }
     484         [ -  + ]:      49197 :   if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
     485                 :          0 :     y[d++]=1;
     486                 :            : END:
     487                 :      60683 :   y[1] = evalsigne(-1) | evallgefint(d);
     488                 :     376516 :   return y;
     489                 :            : }
     490                 :            : 
     491                 :            : INLINE int
     492                 : 1470336639 : cmpiispec(GEN x, GEN y, long lx, long ly)
     493                 :            : {
     494         [ +  + ]: 1470336639 :   if (lx < ly) return -1;
     495         [ +  + ]: 1378798503 :   if (lx > ly) return  1;
     496                 : 1470336639 :   return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     497                 :            : }
     498                 :            : 
     499                 :            : INLINE int
     500                 :   21513573 : equaliispec(GEN x, GEN y, long lx, long ly)
     501                 :            : {
     502         [ +  + ]:   21513573 :   if (lx != ly) return 0;
     503                 :   21513573 :   return !mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     504                 :            : }
     505                 :            : 
     506                 :            : /***********************************************************************/
     507                 :            : /**                                                                   **/
     508                 :            : /**                          MULTIPLICATION                           **/
     509                 :            : /**                                                                   **/
     510                 :            : /***********************************************************************/
     511                 :            : /* assume ny > 0 */
     512                 :            : INLINE GEN
     513                 : 1171443555 : muluispec(ulong x, GEN y, long ny)
     514                 :            : {
     515         [ +  + ]: 1171443555 :   if (ny == 1)
     516                 :  905670169 :     return muluu(x, *y);
     517                 :            :   else
     518                 :            :   {
     519                 :  265773386 :     long lz = ny+3;
     520                 :  265773386 :     GEN z = cgeti(lz);
     521                 :  265773386 :     ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
     522         [ +  + ]:  265773386 :     if (hi) { z[lz - 1] = hi; } else lz--;
     523                 :  265773386 :     z[1] = evalsigne(1) | evallgefint(lz);
     524                 : 1171443555 :     return z;
     525                 :            :   }
     526                 :            : }
     527                 :            : 
     528                 :            : /* a + b*|y| */
     529                 :            : GEN
     530                 :     985952 : addumului(ulong a, ulong b, GEN y)
     531                 :            : {
     532                 :            :   GEN z;
     533                 :            :   long i, lz;
     534                 :            :   ulong hi;
     535         [ +  + ]:     985952 :   if (!signe(y)) return utoi(a);
     536                 :     985296 :   lz = lgefint(y)+1;
     537                 :     985296 :   z = cgeti(lz);
     538                 :     985296 :   z[2]=a;
     539         [ +  + ]:    3887029 :   for(i=3;i<lz;i++) z[i]=0;
     540                 :     985296 :   hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
     541         [ +  + ]:     985296 :   if (hi) z[lz-1]=hi; else lz--;
     542                 :     985296 :   z[1] = evalsigne(1) | evallgefint(lz);
     543                 :     985952 :   avma=(pari_sp)z; return z;
     544                 :            : }
     545                 :            : 
     546                 :            : /***********************************************************************/
     547                 :            : /**                                                                   **/
     548                 :            : /**                          DIVISION                                 **/
     549                 :            : /**                                                                   **/
     550                 :            : /***********************************************************************/
     551                 :            : 
     552                 :            : ulong
     553                 :  411714533 : umodiu(GEN y, ulong x)
     554                 :            : {
     555                 :  411714533 :   long sy=signe(y);
     556                 :            :   ulong hi;
     557         [ -  + ]:  411714533 :   if (!x) pari_err_INV("umodiu",gen_0);
     558         [ +  + ]:  411714533 :   if (!sy) return 0;
     559                 :  353004684 :   hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
     560         [ +  + ]:  353004684 :   if (!hi) return 0;
     561         [ +  + ]:  411714533 :   return (sy > 0)? hi: x - hi;
     562                 :            : }
     563                 :            : 
     564                 :            : /* return |y| \/ x */
     565                 :            : GEN
     566                 :  108927917 : diviu_rem(GEN y, ulong x, ulong *rem)
     567                 :            : {
     568                 :            :   long ly;
     569                 :            :   GEN z;
     570                 :            : 
     571         [ -  + ]:  108927917 :   if (!x) pari_err_INV("diviu_rem",gen_0);
     572         [ +  + ]:  108927917 :   if (!signe(y)) { *rem = 0; return gen_0; }
     573                 :            : 
     574                 :  108399473 :   ly = lgefint(y);
     575 [ +  + ][ +  + ]:  108399473 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
     576                 :            : 
     577                 :  108358321 :   z = cgeti(ly);
     578                 :  108358321 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     579         [ +  + ]:  108358321 :   if (z [ly - 1] == 0) ly--;
     580                 :  108358321 :   z[1] = evallgefint(ly) | evalsigne(1);
     581                 :  108927917 :   return z;
     582                 :            : }
     583                 :            : 
     584                 :            : GEN
     585                 :   35345893 : divis_rem(GEN y, long x, long *rem)
     586                 :            : {
     587                 :   35345893 :   long sy=signe(y),ly,s;
     588                 :            :   GEN z;
     589                 :            : 
     590         [ -  + ]:   35345893 :   if (!x) pari_err_INV("divis_rem",gen_0);
     591         [ +  + ]:   35345893 :   if (!sy) { *rem = 0; return gen_0; }
     592         [ +  + ]:   29236173 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     593                 :            : 
     594                 :   29236173 :   ly = lgefint(y);
     595 [ +  + ][ +  + ]:   29236173 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
     596                 :            : 
     597                 :   24185773 :   z = cgeti(ly);
     598                 :   24185773 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     599         [ +  + ]:   24185773 :   if (sy<0) *rem = -  *rem;
     600         [ +  + ]:   24185773 :   if (z[ly - 1] == 0) ly--;
     601                 :   24185773 :   z[1] = evallgefint(ly) | evalsigne(s);
     602                 :   35345893 :   return z;
     603                 :            : }
     604                 :            : 
     605                 :            : GEN
     606                 :      24293 : divis(GEN y, long x)
     607                 :            : {
     608                 :      24293 :   long sy=signe(y),ly,s;
     609                 :            :   GEN z;
     610                 :            : 
     611         [ -  + ]:      24293 :   if (!x) pari_err_INV("divis",gen_0);
     612         [ -  + ]:      24293 :   if (!sy) return gen_0;
     613         [ -  + ]:      24293 :   if (x<0) { s = -sy; x = -x; } else s=sy;
     614                 :            : 
     615                 :      24293 :   ly = lgefint(y);
     616 [ +  + ][ -  + ]:      24293 :   if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
     617                 :            : 
     618                 :      24293 :   z = cgeti(ly);
     619                 :      24293 :   (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     620         [ +  + ]:      24293 :   if (z[ly - 1] == 0) ly--;
     621                 :      24293 :   z[1] = evallgefint(ly) | evalsigne(s);
     622                 :      24293 :   return z;
     623                 :            : }
     624                 :            : 
     625                 :            : /* We keep llx bits of x and lly bits of y*/
     626                 :            : static GEN
     627                 :   14501205 : divrr_with_gmp(GEN x, GEN y)
     628                 :            : {
     629                 :   14501205 :   long lx=RNLIMBS(x),ly=RNLIMBS(y);
     630                 :   14501205 :   long lw=minss(lx,ly);
     631                 :   14501205 :   long lly=minss(lw+1,ly);
     632                 :   14501205 :   GEN  w=cgetr(lw+2);
     633                 :   14501205 :   long lu=lw+lly;
     634                 :   14501205 :   long llx=minss(lu,lx);
     635                 :   14501205 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     636                 :   14501205 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     637                 :            :   mp_limb_t *q,*r;
     638                 :   14501205 :   long e=expo(x)-expo(y);
     639                 :   14501205 :   long sx=signe(x),sy=signe(y);
     640                 :   14501205 :   xmpn_mirrorcopy(z,RLIMBS(y),lly);
     641                 :   14501205 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     642                 :   14501205 :   xmpn_zero(u,lu-llx);
     643                 :   14501205 :   q = (mp_limb_t *)new_chunk(lw+1);
     644                 :   14501205 :   r = (mp_limb_t *)new_chunk(lly);
     645                 :            : 
     646                 :   14501205 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     647                 :            : 
     648                 :            :   /*Round up: This is not exactly correct we should test 2*r>z*/
     649         [ +  + ]:   14501205 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     650                 :    7086677 :     mpn_add_1(q,q,lw+1,1);
     651                 :            : 
     652                 :   14501205 :   xmpn_mirrorcopy(RLIMBS(w),q,lw);
     653                 :            : 
     654         [ +  + ]:   14501205 :   if (q[lw] == 0) e--;
     655         [ +  + ]:    4449346 :   else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
     656                 :          4 :   else { w[2] = HIGHBIT; e++; }
     657         [ +  + ]:   14501205 :   if (sy < 0) sx = -sx;
     658                 :   14501205 :   w[1] = evalsigne(sx) | evalexpo(e);
     659                 :   14501205 :   avma=(pari_sp) w;
     660                 :   14501205 :   return w;
     661                 :            : }
     662                 :            : 
     663                 :            : /* We keep llx bits of x and lly bits of y*/
     664                 :            : static GEN
     665                 :    3072049 : divri_with_gmp(GEN x, GEN y)
     666                 :            : {
     667                 :    3072049 :   long llx=RNLIMBS(x),ly=NLIMBS(y);
     668                 :    3072049 :   long lly=minss(llx+1,ly);
     669                 :    3072049 :   GEN  w=cgetr(llx+2);
     670                 :    3072049 :   long lu=llx+lly, ld=ly-lly;
     671                 :    3072049 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     672                 :    3072049 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     673                 :            :   mp_limb_t *q,*r;
     674                 :    3072049 :   long sh=bfffo(y[ly+1]);
     675                 :    3072049 :   long e=expo(x)-expi(y);
     676                 :    3072049 :   long sx=signe(x),sy=signe(y);
     677         [ +  + ]:    3072049 :   if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
     678                 :      76337 :   else xmpn_copy(z,LIMBS(y)+ld,lly);
     679                 :    3072049 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     680                 :    3072049 :   xmpn_zero(u,lu-llx);
     681                 :    3072049 :   q = (mp_limb_t *)new_chunk(llx+1);
     682                 :    3072049 :   r = (mp_limb_t *)new_chunk(lly);
     683                 :            : 
     684                 :    3072049 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     685                 :            : 
     686                 :            :   /*Round up: This is not exactly correct we should test 2*r>z*/
     687         [ +  + ]:    3072049 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     688                 :    1382489 :     mpn_add_1(q,q,llx+1,1);
     689                 :            : 
     690                 :    3072049 :   xmpn_mirrorcopy(RLIMBS(w),q,llx);
     691                 :            : 
     692         [ +  + ]:    3072049 :   if (q[llx] == 0) e--;
     693         [ +  - ]:    2574252 :   else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
     694                 :          0 :   else { w[2] = HIGHBIT; e++; }
     695         [ +  + ]:    3072049 :   if (sy < 0) sx = -sx;
     696                 :    3072049 :   w[1] = evalsigne(sx) | evalexpo(e);
     697                 :    3072049 :   avma=(pari_sp) w;
     698                 :    3072049 :   return w;
     699                 :            : }
     700                 :            : 
     701                 :            : GEN
     702                 :   13603131 : divri(GEN x, GEN y)
     703                 :            : {
     704                 :   13603131 :   long  s = signe(y);
     705                 :            : 
     706         [ -  + ]:   13603131 :   if (!s) pari_err_INV("divri",gen_0);
     707         [ +  + ]:   13603131 :   if (!signe(x)) return real_0_bit(expo(x) - expi(y));
     708         [ +  + ]:   13593443 :   if (!is_bigint(y)) {
     709                 :   10521394 :     GEN z = divru(x, y[2]);
     710         [ +  + ]:   10521394 :     if (s < 0) togglesign(z);
     711                 :   10521394 :     return z;
     712                 :            :   }
     713                 :   13603131 :   return divri_with_gmp(x,y);
     714                 :            : }
     715                 :            : 
     716                 :            : GEN
     717                 :   74350881 : divrr(GEN x, GEN y)
     718                 :            : {
     719                 :   74350881 :   long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
     720                 :            :   ulong y0,y1;
     721                 :            :   GEN r, r1;
     722                 :            : 
     723         [ -  + ]:   74350881 :   if (!sy) pari_err_INV("divrr",y);
     724                 :   74350881 :   e = expo(x) - expo(y);
     725         [ +  + ]:   74350881 :   if (!sx) return real_0_bit(e);
     726         [ +  + ]:   72786099 :   if (sy<0) sx = -sx;
     727                 :            : 
     728                 :   72786099 :   lx=lg(x); ly=lg(y);
     729         [ +  + ]:   72786099 :   if (ly==3)
     730                 :            :   {
     731         [ +  + ]:   58284894 :     ulong k = x[2], l = (lx>3)? x[3]: 0;
     732                 :            :     LOCAL_HIREMAINDER;
     733         [ +  + ]:   58284894 :     if (k < uel(y,2)) e--;
     734                 :            :     else
     735                 :            :     {
     736         [ +  + ]:   27061390 :       l >>= 1; if (k&1) l |= HIGHBIT;
     737                 :   27061390 :       k >>= 1;
     738                 :            :     }
     739                 :   58284894 :     hiremainder = k; k = divll(l,y[2]);
     740         [ +  + ]:   58284894 :     if (hiremainder & HIGHBIT)
     741                 :            :     {
     742                 :   16126829 :       k++;
     743         [ -  + ]:   16126829 :       if (!k) { k = HIGHBIT; e++; }
     744                 :            :     }
     745                 :   58284894 :     r = cgetr(3);
     746                 :   58284894 :     r[1] = evalsigne(sx) | evalexpo(e);
     747                 :   58284894 :     r[2] = k; return r;
     748                 :            :   }
     749                 :            : 
     750         [ +  - ]:   14501205 :   if (ly>=DIVRR_GMP_LIMIT)
     751                 :   14501205 :     return divrr_with_gmp(x,y);
     752                 :            : 
     753                 :          0 :   lr = minss(lx,ly); r = new_chunk(lr);
     754                 :          0 :   r1 = r-1;
     755         [ #  # ]:          0 :   r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
     756         [ #  # ]:          0 :   r1[lr] = (lx>ly)? x[lr]: 0;
     757                 :          0 :   y0 = y[2]; y1 = y[3];
     758         [ #  # ]:          0 :   for (i=0; i<lr-1; i++)
     759                 :            :   { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
     760                 :            :     ulong k, qp;
     761                 :            :     LOCAL_HIREMAINDER;
     762                 :            :     LOCAL_OVERFLOW;
     763                 :            : 
     764         [ #  # ]:          0 :     if (uel(r1,1) == y0)
     765                 :            :     {
     766                 :          0 :       qp = ULONG_MAX; k = addll(y0,r1[2]);
     767                 :            :     }
     768                 :            :     else
     769                 :            :     {
     770         [ #  # ]:          0 :       if (uel(r1,1) > y0) /* can't happen if i=0 */
     771                 :            :       {
     772                 :          0 :         GEN y1 = y+1;
     773                 :          0 :         j = lr-i; r1[j] = subll(r1[j],y1[j]);
     774         [ #  # ]:          0 :         for (j--; j>0; j--) r1[j] = subllx(r1[j],y1[j]);
     775 [ #  # ][ #  # ]:          0 :         j=i; do ((ulong*)r)[--j]++; while (j && !r[j]);
     776                 :            :       }
     777                 :          0 :       hiremainder = r1[1]; overflow = 0;
     778                 :          0 :       qp = divll(r1[2],y0); k = hiremainder;
     779                 :            :     }
     780                 :          0 :     j = lr-i+1;
     781         [ #  # ]:          0 :     if (!overflow)
     782                 :            :     {
     783                 :            :       long k3, k4;
     784                 :          0 :       k3 = mulll(qp,y1);
     785         [ #  # ]:          0 :       if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
     786                 :          0 :         k4 = subll(hiremainder,k);
     787                 :            :       else
     788                 :            :       {
     789                 :          0 :         k3 = subll(k3, r1[3]);
     790                 :          0 :         k4 = subllx(hiremainder,k);
     791                 :            :       }
     792 [ #  # ][ #  # ]:          0 :       while (!overflow && k4) { qp--; k3=subll(k3,y1); k4=subllx(k4,y0); }
     793                 :            :     }
     794         [ #  # ]:          0 :     if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
     795         [ #  # ]:          0 :     for (j--; j>1; j--)
     796                 :            :     {
     797                 :          0 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     798                 :          0 :       hiremainder += overflow;
     799                 :            :     }
     800         [ #  # ]:          0 :     if (uel(r1,1) != hiremainder)
     801                 :            :     {
     802         [ #  # ]:          0 :       if (uel(r1,1) < hiremainder)
     803                 :            :       {
     804                 :          0 :         qp--;
     805                 :          0 :         j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
     806         [ #  # ]:          0 :         for (j--; j>1; j--) r1[j] = addllx(r1[j], y[j]);
     807                 :            :       }
     808                 :            :       else
     809                 :            :       {
     810                 :          0 :         r1[1] -= hiremainder;
     811         [ #  # ]:          0 :         while (r1[1])
     812                 :            :         {
     813 [ #  # ][ #  # ]:          0 :           qp++; if (!qp) { j=i; do ((ulong*)r)[--j]++; while (j && !r[j]); }
                 [ #  # ]
     814                 :          0 :           j = lr-i-(lr-i>=ly); r1[j] = subll(r1[j],y[j]);
     815         [ #  # ]:          0 :           for (j--; j>1; j--) r1[j] = subllx(r1[j],y[j]);
     816                 :          0 :           r1[1] -= overflow;
     817                 :            :         }
     818                 :            :       }
     819                 :            :     }
     820                 :          0 :     *++r1 = qp;
     821                 :            :   }
     822                 :            :   /* i = lr-1 */
     823                 :            :   /* round correctly */
     824         [ #  # ]:          0 :   if (uel(r1,1) > (y0>>1))
     825                 :            :   {
     826 [ #  # ][ #  # ]:          0 :     j=i; do r[--j]++; while (j && !r[j]);
     827                 :            :   }
     828         [ #  # ]:          0 :   r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
     829         [ #  # ]:          0 :   if (r[0] == 0) e--;
     830         [ #  # ]:          0 :   else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
     831                 :            :   else { /* possible only when rounding up to 0x2 0x0 ... */
     832                 :          0 :     r[2] = (long)HIGHBIT; e++;
     833                 :            :   }
     834                 :          0 :   r[0] = evaltyp(t_REAL)|evallg(lr);
     835                 :          0 :   r[1] = evalsigne(sx) | evalexpo(e);
     836                 :   74350881 :   return r;
     837                 :            : }
     838                 :            : 
     839                 :            : /* Integer division x / y: such that sign(r) = sign(x)
     840                 :            :  *   if z = ONLY_REM return remainder, otherwise return quotient
     841                 :            :  *   if z != NULL set *z to remainder
     842                 :            :  *   *z is the last object on stack (and thus can be disposed of with cgiv
     843                 :            :  *   instead of gerepile)
     844                 :            :  * If *z is zero, we put gen_0 here and no copy.
     845                 :            :  * space needed: lx + ly
     846                 :            :  */
     847                 :            : GEN
     848                 :  528696484 : dvmdii(GEN x, GEN y, GEN *z)
     849                 :            : {
     850                 :  528696484 :   long sx=signe(x),sy=signe(y);
     851                 :            :   long lx, ly, lq;
     852                 :            :   pari_sp av;
     853                 :            :   GEN r,q;
     854                 :            : 
     855         [ +  + ]:  528696484 :   if (!sy) pari_err_INV("dvmdii",y);
     856         [ +  + ]:  528696464 :   if (!sx)
     857                 :            :   {
     858 [ +  + ][ +  + ]:   12562324 :     if (!z || z == ONLY_REM) return gen_0;
     859                 :    5214279 :     *z=gen_0; return gen_0;
     860                 :            :   }
     861                 :  516134140 :   lx=lgefint(x);
     862                 :  516134140 :   ly=lgefint(y); lq=lx-ly;
     863         [ +  + ]:  516134140 :   if (lq <= 0)
     864                 :            :   {
     865         [ +  + ]:  241072911 :     if (lq == 0)
     866                 :            :     {
     867                 :  226887676 :       long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
     868         [ +  + ]:  226887676 :       if (s>0) goto DIVIDE;
     869         [ +  + ]:   79614700 :       if (s==0)
     870                 :            :       {
     871         [ +  + ]:    5314873 :         if (z == ONLY_REM) return gen_0;
     872         [ +  + ]:    1456338 :         if (z) *z = gen_0;
     873         [ +  + ]:    1456338 :         if (sx < 0) sy = -sy;
     874                 :    1456338 :         return stoi(sy);
     875                 :            :       }
     876                 :            :     }
     877         [ +  + ]:   88485062 :     if (z == ONLY_REM) return icopy(x);
     878         [ +  + ]:    4189812 :     if (z) *z = icopy(x);
     879                 :    4189812 :     return gen_0;
     880                 :            :   }
     881                 :            : DIVIDE: /* quotient is non-zero */
     882         [ +  + ]:  422334205 :   av=avma; if (sx<0) sy = -sy;
     883         [ +  + ]:  422334205 :   if (ly==3)
     884                 :            :   {
     885                 :  215648757 :     ulong lq = lx;
     886                 :            :     ulong si;
     887                 :  215648757 :     q  = cgeti(lq);
     888                 :  215648757 :     si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
     889         [ +  + ]:  215648757 :     if (q[lq - 1] == 0) lq--;
     890         [ +  + ]:  215648757 :     if (z == ONLY_REM)
     891                 :            :     {
     892         [ +  + ]:  194573318 :       avma=av; if (!si) return gen_0;
     893                 :  179465609 :       r=cgeti(3);
     894                 :  179465609 :       r[1] = evalsigne(sx) | evallgefint(3);
     895                 :  179465609 :       r[2] = si; return r;
     896                 :            :     }
     897                 :   21075439 :     q[1] = evalsigne(sy) | evallgefint(lq);
     898         [ +  + ]:   21075439 :     if (!z) return q;
     899         [ +  + ]:   20571468 :     if (!si) { *z=gen_0; return q; }
     900                 :    8054783 :     r=cgeti(3);
     901                 :    8054783 :     r[1] = evalsigne(sx) | evallgefint(3);
     902                 :    8054783 :     r[2] = si; *z=r; return q;
     903                 :            :   }
     904         [ +  + ]:  206685448 :   if (z == ONLY_REM)
     905                 :            :   {
     906                 :  201739819 :     ulong lr = lgefint(y);
     907                 :  201739819 :     ulong lq = lgefint(x)-lgefint(y)+3;
     908                 :  201739819 :     GEN r = cgeti(lr);
     909                 :  201739819 :     GEN q = cgeti(lq);
     910                 :  201739819 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     911         [ +  + ]:  201739819 :     if (!r[lr - 1])
     912                 :            :     {
     913 [ +  + ][ +  + ]:   99777857 :       while(lr>2 && !r[lr - 1]) lr--;
     914         [ +  + ]:   44175061 :       if (lr == 2) {avma=av; return gen_0;} /* exact division */
     915                 :            :     }
     916                 :  198983680 :     r[1] = evalsigne(sx) | evallgefint(lr);
     917                 :  198983680 :     avma = (pari_sp) r; return r;
     918                 :            :   }
     919                 :            :   else
     920                 :            :   {
     921                 :    4945629 :     ulong lq = lgefint(x)-lgefint(y)+3;
     922                 :    4945629 :     ulong lr = lgefint(y);
     923                 :    4945629 :     GEN q = cgeti(lq);
     924                 :    4945629 :     GEN r = cgeti(lr);
     925                 :    4945629 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     926         [ +  + ]:    4945629 :     if (q[lq - 1] == 0) lq--;
     927                 :    4945629 :     q[1] = evalsigne(sy) | evallgefint(lq);
     928         [ +  + ]:    4945629 :     if (!z) { avma = (pari_sp)q; return q; }
     929         [ +  + ]:    4397421 :     if (!r[lr - 1])
     930                 :            :     {
     931 [ +  + ][ +  + ]:    7088537 :       while(lr>2 && !r[lr - 1]) lr--;
     932         [ +  + ]:    1513234 :       if (lr == 2) {avma=(pari_sp) q; *z=gen_0; return q;} /* exact division */
     933                 :            :     }
     934                 :    3206661 :     r[1] = evalsigne(sx) | evallgefint(lr);
     935                 :  528696464 :     avma = (pari_sp) r; *z = r; return q;
     936                 :            :   }
     937                 :            : }
     938                 :            : 
     939                 :            : /* Montgomery reduction.
     940                 :            :  * N has k words, assume T >= 0 has less than 2k.
     941                 :            :  * Return res := T / B^k mod N, where B = 2^BIL
     942                 :            :  * such that 0 <= res < T/B^k + N  and  res has less than k words */
     943                 :            : GEN
     944                 :    5326029 : red_montgomery(GEN T, GEN N, ulong inv)
     945                 :            : {
     946                 :            :   pari_sp av;
     947                 :            :   GEN Te, Td, Ne, Nd, scratch;
     948                 :    5326029 :   ulong i, j, m, t, d, k = NLIMBS(N);
     949                 :            :   int carry;
     950                 :            :   LOCAL_HIREMAINDER;
     951                 :            :   LOCAL_OVERFLOW;
     952                 :            : 
     953         [ -  + ]:    5326029 :   if (k == 0) return gen_0;
     954                 :    5326029 :   d = NLIMBS(T); /* <= 2*k */
     955         [ +  + ]:    5326029 :   if (d == 0) return gen_0;
     956                 :            : #ifdef DEBUG
     957                 :            :   if (d > 2*k) pari_err_BUG("red_montgomery");
     958                 :            : #endif
     959         [ -  + ]:    5325750 :   if (k == 1)
     960                 :            :   { /* as below, special cased for efficiency */
     961                 :          0 :     ulong n = uel(N,2);
     962         [ #  # ]:          0 :     if (d == 1) {
     963                 :          0 :       hiremainder = uel(T,2);
     964                 :          0 :       m = hiremainder * inv;
     965                 :          0 :       (void)addmul(m, n); /* t + m*n = 0 */
     966                 :          0 :       return utoi(hiremainder);
     967                 :            :     } else { /* d = 2 */
     968                 :          0 :       hiremainder = uel(T,2);
     969                 :          0 :       m = hiremainder * inv;
     970                 :          0 :       (void)addmul(m, n); /* t + m*n = 0 */
     971                 :          0 :       t = addll(hiremainder, uel(T,3));
     972         [ #  # ]:          0 :       if (overflow) t -= n; /* t > n doesn't fit in 1 word */
     973                 :          0 :       return utoi(t);
     974                 :            :     }
     975                 :            :   }
     976                 :            :   /* assume k >= 2 */
     977                 :    5325750 :   av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
     978                 :            : 
     979                 :            :   /* copy T to scratch space (pad with zeroes to 2k words) */
     980                 :    5325750 :   Td = scratch;
     981                 :    5325750 :   Te = T + 2;
     982         [ +  + ]:   22536519 :   for (i=0; i < d     ; i++) *Td++ = *Te++;
     983         [ +  + ]:   14200809 :   for (   ; i < (k<<1); i++) *Td++ = 0;
     984                 :            : 
     985                 :    5325750 :   Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
     986                 :    5325750 :   Ne = N + 1;       /* 1 beyond end of N mantissa */
     987                 :            : 
     988                 :    5325750 :   carry = 0;
     989         [ +  + ]:   18368664 :   for (i=0; i<k; i++) /* set T := T/B nod N, k times */
     990                 :            :   {
     991                 :   13042914 :     Td = Te; /* one beyond end of (new) T mantissa */
     992                 :   13042914 :     Nd = Ne;
     993                 :   13042914 :     hiremainder = *++Td;
     994                 :   13042914 :     m = hiremainder * inv; /* solve T + m N = O(B) */
     995                 :            : 
     996                 :            :     /* set T := (T + mN) / B */
     997                 :   13042914 :     Te = Td;
     998                 :   13042914 :     (void)addmul(m, *++Nd); /* = 0 */
     999         [ +  + ]:   35241954 :     for (j=1; j<k; j++)
    1000                 :            :     {
    1001                 :   22199040 :       t = addll(addmul(m, *++Nd), *++Td);
    1002                 :   22199040 :       *Td = t;
    1003                 :   22199040 :       hiremainder += overflow;
    1004                 :            :     }
    1005                 :   13042914 :     t = addll(hiremainder, *++Td); *Td = t + carry;
    1006 [ +  + ][ +  + ]:   13042914 :     carry = (overflow || (carry && *Td == 0));
                 [ -  + ]
    1007                 :            :   }
    1008         [ +  + ]:    5325750 :   if (carry)
    1009                 :            :   { /* Td > N overflows (k+1 words), set Td := Td - N */
    1010                 :      10830 :     Td = Te;
    1011                 :      10830 :     Nd = Ne;
    1012                 :      10830 :     t = subll(*++Td, *++Nd); *Td = t;
    1013         [ +  + ]:      41565 :     while (Td < (GEN)av) { t = subllx(*++Td, *++Nd); *Td = t; }
    1014                 :            :   }
    1015                 :            : 
    1016                 :            :   /* copy result */
    1017                 :    5325750 :   Td = (GEN)av - 1; /* *Td = high word of final result */
    1018 [ +  + ][ +  - ]:    8906037 :   while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
    1019         [ -  + ]:    5325750 :   k = Td - Te; if (!k) { avma = av; return gen_0; }
    1020                 :    5325750 :   Td = (GEN)av - k; /* will write mantissa there */
    1021                 :    5325750 :   (void)memmove(Td, Te+1, k*sizeof(long));
    1022                 :    5325750 :   Td -= 2;
    1023                 :    5325750 :   Td[0] = evaltyp(t_INT) | evallg(k+2);
    1024                 :    5325750 :   Td[1] = evalsigne(1) | evallgefint(k+2);
    1025                 :            : #ifdef DEBUG
    1026                 :            : {
    1027                 :            :   long l = NLIMBS(N), s = BITS_IN_LONG*l;
    1028                 :            :   GEN R = int2n(s);
    1029                 :            :   GEN res = remii(mulii(T, Fp_inv(R, N)), N);
    1030                 :            :   if (k > lgefint(N)
    1031                 :            :     || !equalii(remii(Td,N),res)
    1032                 :            :     || cmpii(Td, addii(shifti(T, -s), N)) >= 0) pari_err_BUG("red_montgomery");
    1033                 :            : }
    1034                 :            : #endif
    1035                 :    5326029 :   avma = (pari_sp)Td; return Td;
    1036                 :            : }
    1037                 :            : 
    1038                 :            : /* EXACT INTEGER DIVISION */
    1039                 :            : 
    1040                 :            : #if 1 /* use undocumented GMP interface */
    1041                 :            : static void
    1042                 :   98343465 : GEN2mpz(mpz_t X, GEN x)
    1043                 :            : {
    1044                 :   98343465 :   long l = lgefint(x)-2;
    1045                 :   98343465 :   X->_mp_alloc = l;
    1046         [ +  + ]:   98343465 :   X->_mp_size = signe(x) > 0? l: -l;
    1047                 :   98343465 :   X->_mp_d = LIMBS(x);
    1048                 :   98343465 : }
    1049                 :            : static void
    1050                 :   89706516 : mpz2GEN(GEN z, mpz_t Z)
    1051                 :            : {
    1052                 :   89706516 :   long l = Z->_mp_size;
    1053         [ +  + ]:   89706516 :   z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
    1054                 :   89706516 : }
    1055                 :            : 
    1056                 :            : /* assume y != 0 and the division is exact */
    1057                 :            : GEN
    1058                 :  209090562 : diviuexact(GEN x, ulong y)
    1059                 :            : {
    1060         [ +  + ]:  209090562 :   if (!signe(x)) return gen_0;
    1061                 :            :   {
    1062                 :   81069567 :     long l = lgefint(x);
    1063                 :            :     mpz_t X, Z;
    1064                 :   81069567 :     GEN z = cgeti(l);
    1065                 :   81069567 :     GEN2mpz(X, x);
    1066                 :   81069567 :     Z->_mp_alloc = l-2;
    1067                 :   81069567 :     Z->_mp_size  = l-2;
    1068                 :   81069567 :     Z->_mp_d = LIMBS(z);
    1069                 :   81069567 :     mpz_divexact_ui(Z, X, y);
    1070                 :  209090562 :     mpz2GEN(z, Z); return z;
    1071                 :            :   }
    1072                 :            : }
    1073                 :            : 
    1074                 :            : /* Find z such that x=y*z, knowing that y | x (unchecked) */
    1075                 :            : GEN
    1076                 :  217567585 : diviiexact(GEN x, GEN y)
    1077                 :            : {
    1078         [ -  + ]:  217567585 :   if (!signe(y)) pari_err_INV("diviiexact",y);
    1079         [ +  + ]:  217567585 :   if (lgefint(y) == 3)
    1080                 :            :   {
    1081                 :  204402411 :     GEN z = diviuexact(x, y[2]);
    1082         [ +  + ]:  204402411 :     if (signe(y) < 0) togglesign(z);
    1083                 :  204402411 :     return z;
    1084                 :            :   }
    1085         [ +  + ]:   13165174 :   if (!signe(x)) return gen_0;
    1086                 :            :   {
    1087                 :    8636949 :     long l = lgefint(x);
    1088                 :            :     mpz_t X, Y, Z;
    1089                 :    8636949 :     GEN z = cgeti(l);
    1090                 :    8636949 :     GEN2mpz(X, x);
    1091                 :    8636949 :     GEN2mpz(Y, y);
    1092                 :    8636949 :     Z->_mp_alloc = l-2;
    1093                 :    8636949 :     Z->_mp_size  = l-2;
    1094                 :    8636949 :     Z->_mp_d = LIMBS(z);
    1095                 :    8636949 :     mpz_divexact(Z, X, Y);
    1096                 :  217567585 :     mpz2GEN(z, Z); return z;
    1097                 :            :   }
    1098                 :            : }
    1099                 :            : #else
    1100                 :            : /* assume y != 0 and the division is exact */
    1101                 :            : GEN
    1102                 :            : diviuexact(GEN x, ulong y)
    1103                 :            : {
    1104                 :            :   /*TODO: implement true exact division.*/
    1105                 :            :   return divii(x,utoi(y));
    1106                 :            : }
    1107                 :            : 
    1108                 :            : /* Find z such that x=y*z, knowing that y | x (unchecked)
    1109                 :            :  * Method: y0 z0 = x0 mod B = 2^BITS_IN_LONG ==> z0 = 1/y0 mod B.
    1110                 :            :  *    Set x := (x - z0 y) / B, updating only relevant words, and repeat */
    1111                 :            : GEN
    1112                 :            : diviiexact(GEN x, GEN y)
    1113                 :            : {
    1114                 :            :   /*TODO: use mpn_bdivmod instead*/
    1115                 :            :   return divii(x,y);
    1116                 :            : }
    1117                 :            : #endif
    1118                 :            : 
    1119                 :            : /* assume yz != and yz | x */
    1120                 :            : GEN
    1121                 :      53196 : diviuuexact(GEN x, ulong y, ulong z)
    1122                 :            : {
    1123                 :            :   long tmp[4];
    1124                 :            :   ulong t;
    1125                 :            :   LOCAL_HIREMAINDER;
    1126                 :      53196 :   t = mulll(y, z);
    1127         [ +  - ]:      53196 :   if (!hiremainder) return diviuexact(x, t);
    1128                 :          0 :   tmp[0] = evaltyp(t_INT)|_evallg(4);
    1129                 :          0 :   tmp[1] = evalsigne(1)|evallgefint(4);
    1130                 :          0 :   tmp[2] = t;
    1131                 :          0 :   tmp[3] = hiremainder;
    1132                 :      53196 :   return diviiexact(x, tmp);
    1133                 :            : }
    1134                 :            : 
    1135                 :            : 
    1136                 :            : /********************************************************************/
    1137                 :            : /**                                                                **/
    1138                 :            : /**               INTEGER MULTIPLICATION                           **/
    1139                 :            : /**                                                                **/
    1140                 :            : /********************************************************************/
    1141                 :            : 
    1142                 :            : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
    1143                 :            : GEN
    1144                 : 1058422087 : muliispec(GEN x, GEN y, long nx, long ny)
    1145                 :            : {
    1146                 :            :   GEN zd;
    1147                 :            :   long lz;
    1148                 :            :   ulong hi;
    1149                 :            : 
    1150         [ +  + ]: 1058422087 :   if (nx < ny) swapspec(x,y, nx,ny);
    1151         [ -  + ]: 1058422087 :   if (!ny) return gen_0;
    1152         [ +  + ]: 1058422087 :   if (ny == 1) return muluispec((ulong)*y, x, nx);
    1153                 :            : 
    1154                 :  306192498 :   lz = nx+ny+2;
    1155                 :  306192498 :   zd = cgeti(lz);
    1156                 :  306192498 :   hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
    1157         [ +  + ]:  306192498 :   if (!hi) lz--;
    1158                 :            :   /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
    1159                 :            : 
    1160                 :  306192498 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1161                 : 1058422087 :   return zd;
    1162                 :            : }
    1163                 :            : GEN
    1164                 :      75404 : muluui(ulong x, ulong y, GEN z)
    1165                 :            : {
    1166                 :      75404 :   long t, s = signe(z);
    1167                 :            :   GEN r;
    1168                 :            :   LOCAL_HIREMAINDER;
    1169                 :            : 
    1170 [ +  - ][ +  - ]:      75404 :   if (!x || !y || !signe(z)) return gen_0;
                 [ +  + ]
    1171                 :      75008 :   t = mulll(x,y);
    1172         [ +  - ]:      75008 :   if (!hiremainder)
    1173                 :      75008 :     r = muluispec(t, z+2, lgefint(z)-2);
    1174                 :            :   else
    1175                 :            :   {
    1176                 :            :     long tmp[2];
    1177                 :          0 :     tmp[1] = hiremainder;
    1178                 :          0 :     tmp[0] = t;
    1179                 :          0 :     r = muliispec(z+2,tmp, lgefint(z)-2, 2);
    1180                 :            :   }
    1181                 :      75404 :   setsigne(r,s); return r;
    1182                 :            : }
    1183                 :            : 
    1184                 :            : GEN
    1185                 :  567010964 : sqrispec(GEN x, long nx)
    1186                 :            : {
    1187                 :            :   GEN zd;
    1188                 :            :   long lz;
    1189                 :            : 
    1190         [ +  + ]:  567010964 :   if (!nx) return gen_0;
    1191         [ +  + ]:   86825643 :   if (nx==1) return sqru(*x);
    1192                 :            : 
    1193                 :   40782545 :   lz = (nx<<1)+2;
    1194                 :   40782545 :   zd = cgeti(lz);
    1195                 :            : #ifdef mpn_sqr
    1196                 :   40782545 :   mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
    1197                 :            : #else
    1198                 :            :   mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
    1199                 :            : #endif
    1200         [ +  + ]:   40782545 :   if (zd[lz-1]==0) lz--;
    1201                 :            : 
    1202                 :   40782545 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1203                 :  567010964 :   return zd;
    1204                 :            : }
    1205                 :            : 
    1206                 :            : INLINE GEN
    1207                 :    6700958 : sqrispec_mirror(GEN x, long nx)
    1208                 :            : {
    1209                 :    6700958 :   GEN cx=new_chunk(nx);
    1210                 :            :   GEN z;
    1211                 :    6700958 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1212                 :    6700958 :   z=sqrispec(cx, nx);
    1213                 :    6700958 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1214                 :    6700958 :   return z;
    1215                 :            : }
    1216                 :            : 
    1217                 :            : /* leaves garbage on the stack. */
    1218                 :            : INLINE GEN
    1219                 :   54524105 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
    1220                 :            : {
    1221                 :            :   GEN cx, cy, z;
    1222                 :   54524105 :   long s = 0;
    1223 [ +  - ][ +  + ]:   63920032 :   while (nx && x[nx-1]==0) { nx--; s++; }
    1224 [ +  - ][ +  + ]:   65483788 :   while (ny && y[ny-1]==0) { ny--; s++; }
    1225                 :   54524105 :   cx=new_chunk(nx); cy=new_chunk(ny);
    1226                 :   54524105 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1227                 :   54524105 :   xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
    1228         [ +  + ]:   54524105 :   z =  nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
    1229         [ +  + ]:   54524105 :   if (s)
    1230                 :            :   {
    1231                 :    3750764 :     long i, lz = lgefint(z) + s;
    1232                 :    3750764 :     (void)new_chunk(s);
    1233                 :    3750764 :     z -= s;
    1234         [ +  + ]:   24106374 :     for (i=0; i<s; i++) z[2+i]=0;
    1235                 :    3750764 :     z[1] = evalsigne(1) | evallgefint(lz);
    1236                 :    3750764 :     z[0] = evaltyp(t_INT) | evallg(lz);
    1237                 :            :   }
    1238                 :   54524105 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1239                 :   54524105 :   return z;
    1240                 :            : }
    1241                 :            : 
    1242                 :            : /* x % (2^n), assuming n >= 0 */
    1243                 :            : GEN
    1244                 :    7194382 : remi2n(GEN x, long n)
    1245                 :            : {
    1246                 :            :   ulong hi;
    1247                 :    7194382 :   long l, k, lx, ly, sx = signe(x);
    1248                 :            :   GEN z, xd, zd;
    1249                 :            : 
    1250 [ +  + ][ -  + ]:    7194382 :   if (!sx || !n) return gen_0;
    1251                 :            : 
    1252                 :    7181062 :   k = dvmdsBIL(n, &l);
    1253                 :    7181062 :   lx = lgefint(x);
    1254         [ +  + ]:    7181062 :   if (lx < k+3) return icopy(x);
    1255                 :            : 
    1256                 :    7129666 :   xd = x + (2 + k);
    1257                 :            :   /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
    1258                 :            :    *              ^--- initial xd  */
    1259                 :    7129666 :   hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
    1260         [ +  + ]:    7129666 :   if (!hi)
    1261                 :            :   { /* strip leading zeroes from result */
    1262 [ +  + ][ +  + ]:     488885 :     xd--; while (k && !*xd) { k--; xd--; }
    1263         [ +  + ]:     483245 :     if (!k) return gen_0;
    1264                 :     336257 :     ly = k+2;
    1265                 :            :   }
    1266                 :            :   else
    1267                 :    6646421 :     ly = k+3;
    1268                 :            : 
    1269                 :    6982678 :   zd = z = cgeti(ly);
    1270                 :    6982678 :   *++zd = evalsigne(sx) | evallgefint(ly);
    1271         [ +  + ]:   28119644 :   xd = x+1; for ( ;k; k--) *++zd = *++xd;
    1272         [ +  + ]:    6982678 :   if (hi) *++zd = hi;
    1273                 :    7194382 :   return z;
    1274                 :            : }
    1275                 :            : 
    1276                 :            : /********************************************************************/
    1277                 :            : /**                                                                **/
    1278                 :            : /**                      INTEGER SQUARE ROOT                       **/
    1279                 :            : /**                                                                **/
    1280                 :            : /********************************************************************/
    1281                 :            : 
    1282                 :            : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
    1283                 :            :  * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
    1284                 :            :  * remainder is 0. R = NULL is allowed. */
    1285                 :            : GEN
    1286                 :    1397073 : sqrtremi(GEN a, GEN *r)
    1287                 :            : {
    1288                 :    1397073 :   long l, na = NLIMBS(a);
    1289                 :            :   mp_size_t nr;
    1290                 :            :   GEN S;
    1291         [ -  + ]:    1397073 :   if (!na) {
    1292         [ #  # ]:          0 :     if (r) *r = gen_0;
    1293                 :          0 :     return gen_0;
    1294                 :            :   }
    1295                 :    1397073 :   l = (na + 5) >> 1; /* 2 + ceil(na/2) */
    1296                 :    1397073 :   S = cgetipos(l);
    1297         [ +  + ]:    1397073 :   if (r) {
    1298                 :       3398 :     GEN R = cgeti(2 + na);
    1299                 :       3398 :     nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
    1300         [ +  + ]:       3398 :     if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
    1301                 :       3034 :     else    { avma = (pari_sp)S; R = gen_0; }
    1302                 :       3398 :     *r = R;
    1303                 :            :   }
    1304                 :            :   else
    1305                 :    1393675 :     (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
    1306                 :    1397073 :   return S;
    1307                 :            : }
    1308                 :            : 
    1309                 :            : /* compute sqrt(|a|), assuming a != 0 */
    1310                 :            : GEN
    1311                 :   11967409 : sqrtr_abs(GEN a)
    1312                 :            : {
    1313                 :            :   GEN res;
    1314                 :            :   mp_limb_t *b, *c;
    1315                 :   11967409 :   long l = RNLIMBS(a), e = expo(a), er = e>>1;
    1316                 :            :   long n;
    1317                 :   11967409 :   res = cgetr(2 + l);
    1318                 :   11967409 :   res[1] = evalsigne(1) | evalexpo(er);
    1319         [ +  + ]:   11967409 :   if (e&1)
    1320                 :            :   {
    1321                 :    6512140 :     b = (mp_limb_t *) new_chunk(l<<1);
    1322                 :    6512140 :     xmpn_zero(b,l);
    1323                 :    6512140 :     xmpn_mirrorcopy(b+l, RLIMBS(a), l);
    1324                 :    6512140 :     c = (mp_limb_t *) new_chunk(l);
    1325                 :    6512140 :     n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
    1326 [ +  + ][ +  + ]:    6512140 :     if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
                 [ +  + ]
    1327                 :            :   }
    1328                 :            :   else
    1329                 :            :   {
    1330                 :            :     ulong u;
    1331                 :    5455269 :     b = (mp_limb_t *) mantissa2nr(a,-1);
    1332                 :    5455269 :     b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
    1333                 :    5455269 :     b = (mp_limb_t *) new_chunk(l);
    1334                 :    5455269 :     xmpn_zero(b,l+1); /* overwrites the former b[0] */
    1335                 :    5455269 :     c = (mp_limb_t *) new_chunk(l + 1);
    1336                 :    5455269 :     n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
    1337                 :    5455269 :     u = (ulong)*c++;
    1338 [ +  + ][ +  + ]:    5455269 :     if ( u&HIGHBIT || (u == ~HIGHBIT &&
                 [ -  + ]
    1339 [ #  # ][ #  # ]:          0 :              (n>l || (n==l && mpn_cmp(b,c,l)>0))))
    1340                 :    2549138 :       mpn_add_1(c,c,l,1);
    1341                 :            :   }
    1342                 :   11967409 :   xmpn_mirrorcopy(RLIMBS(res),c,l);
    1343                 :   11967409 :   avma = (pari_sp)res; return res;
    1344                 :            : }
    1345                 :            : 
    1346                 :            : /* Normalize a non-negative integer */
    1347                 :            : GEN
    1348                 :  129857078 : int_normalize(GEN x, long known_zero_words)
    1349                 :            : {
    1350                 :  129857078 :   long i =  lgefint(x) - 1 - known_zero_words;
    1351         [ +  + ]: 1096918201 :   for ( ; i > 1; i--)
    1352         [ +  + ]: 1034156669 :     if (x[i]) { setlgefint(x, i+1); return x; }
    1353                 :  129857078 :   x[1] = evalsigne(0) | evallgefint(2); return x;
    1354                 :            : }
    1355                 :            : 
    1356                 :            : /********************************************************************
    1357                 :            :  **                                                                **
    1358                 :            :  **                           Base Conversion                      **
    1359                 :            :  **                                                                **
    1360                 :            :  ********************************************************************/
    1361                 :            : 
    1362                 :            : ulong *
    1363                 :     349378 : convi(GEN x, long *l)
    1364                 :            : {
    1365                 :     349378 :   long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
    1366                 :     349378 :   GEN str = cgetg(n+1, t_VECSMALL);
    1367                 :     349378 :   unsigned char *res = (unsigned char*) GSTR(str);
    1368                 :     349378 :   long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
    1369                 :            :   long lz;
    1370                 :            :   ulong *z;
    1371                 :            :   long i, j;
    1372                 :            :   unsigned char *t;
    1373         [ -  + ]:     349378 :   while (!*res) {res++; llz--;} /*Strip leading zeros*/
    1374                 :     349378 :   lz  = (8+llz)/9;
    1375                 :     349378 :   z = (ulong*)new_chunk(1+lz);
    1376                 :     349378 :   t=res+llz+9;
    1377         [ +  + ]:     905051 :   for(i=0;i<llz-8;i+=9)
    1378                 :            :   {
    1379                 :            :     ulong s;
    1380                 :     555673 :     t-=18;
    1381                 :     555673 :     s=*t++;
    1382         [ +  + ]:    5001057 :     for (j=1; j<9;j++)
    1383                 :    4445384 :       s=10*s+*t++;
    1384                 :     555673 :     *z++=s;
    1385                 :            :   }
    1386         [ +  + ]:     349378 :   if (i<llz)
    1387                 :            :   {
    1388                 :     335481 :     unsigned char *t = res;
    1389                 :     335481 :     ulong s=*t++;
    1390         [ +  + ]:    1089192 :     for (j=i+1; j<llz;j++)
    1391                 :     753711 :       s=10*s+*t++;
    1392                 :     335481 :     *z++=s;
    1393                 :            :   }
    1394                 :     349378 :   *l = lz;
    1395                 :     349378 :   return z;
    1396                 :            : }

Generated by: LCOV version 1.9