Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 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.13.0 lcov report (development 25822-9a3934d752) Lines: 651 719 90.5 %
Date: 2020-09-19 06:09:01 Functions: 55 56 98.2 %
Legend: Lines: hit not hit

          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      374110 : static void pari_gmp_free(void *ptr, size_t old_size){
      51      374110 :   (void)old_size; pari_free(ptr);
      52      374110 : }
      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         992 : pari_kernel_init(void)
      60             : {
      61         992 :   mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
      62         992 :   mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
      63         992 : }
      64             : 
      65             : const char *
      66           4 : pari_kernel_version(void)
      67             : {
      68             : #ifdef gmp_version
      69           4 :   return gmp_version;
      70             : #else
      71             :   return "";
      72             : #endif
      73             : }
      74             : 
      75             : void
      76         984 : pari_kernel_close(void)
      77             : {
      78             :   void *(*new_gmp_malloc)(size_t new_size);
      79             :   void *(*new_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
      80             :   void (*new_gmp_free)(void *ptr, size_t old_size);
      81         984 :   mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
      82         984 :   if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
      83         984 :   if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
      84         984 :   if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
      85         984 :   mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
      86         984 : }
      87             : 
      88             : #define LIMBS(x)  ((mp_limb_t *)((x)+2))
      89             : #define NLIMBS(x) (lgefint(x)-2)
      90             : /*This one is for t_REALs to emphasize they are not t_INTs*/
      91             : #define RLIMBS(x)  ((mp_limb_t *)((x)+2))
      92             : #define RNLIMBS(x) (lg(x)-2)
      93             : 
      94             : INLINE void
      95     3087971 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
      96             : {
      97    32925529 :   while (--n >= 0) x[n]=y[n];
      98     3087971 : }
      99             : 
     100             : INLINE void
     101   308815728 : xmpn_mirror(mp_limb_t *x, long n)
     102             : {
     103             :   long i;
     104  1570356076 :   for(i=0;i<(n>>1);i++)
     105             :   {
     106  1261540348 :     ulong m=x[i];
     107  1261540348 :     x[i]=x[n-1-i];
     108  1261540348 :     x[n-1-i]=m;
     109             :   }
     110   308815728 : }
     111             : 
     112             : INLINE void
     113   378258263 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
     114             : {
     115             :   long i;
     116  3830284506 :   for(i=0;i<n;i++)
     117  3452026243 :     z[i]=x[n-1-i];
     118   378258263 : }
     119             : 
     120             : INLINE void
     121   113583480 : xmpn_zero(mp_ptr x, mp_size_t n)
     122             : {
     123   714456669 :   while (--n >= 0) x[n]=0;
     124   113583480 : }
     125             : 
     126             : INLINE GEN
     127    23186088 : icopy_ef(GEN x, long l)
     128             : {
     129    23186088 :   register long lx = lgefint(x);
     130    23186088 :   const GEN y = cgeti(l);
     131             : 
     132   159720498 :   while (--lx > 0) y[lx]=x[lx];
     133    23184764 :   return y;
     134             : }
     135             : 
     136             : /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
     137             :  * GENs but pairs (long *a, long na) representing a list of digits (in basis
     138             :  * BITS_IN_LONG) : a[0], ..., a[na-1]. [ In ordre to facilitate splitting: no
     139             :  * need to reintroduce codewords ]
     140             :  * Use speci(a,na) to visualize the corresponding GEN.
     141             :  */
     142             : 
     143             : /***********************************************************************/
     144             : /**                                                                   **/
     145             : /**                     ADDITION / SUBTRACTION                        **/
     146             : /**                                                                   **/
     147             : /***********************************************************************/
     148             : 
     149             : GEN
     150     2981282 : setloop(GEN a)
     151             : {
     152     2981282 :   pari_sp av = avma - 2 * sizeof(long);
     153     2981282 :   (void)cgetg(lgefint(a) + 3, t_VECSMALL);
     154     2981286 :   return icopy_avma(a, av); /* two cells of extra space after a */
     155             : }
     156             : 
     157             : /* we had a = setloop(?), then some incloops. Reset a to b */
     158             : GEN
     159      170848 : resetloop(GEN a, GEN b) {
     160      170848 :   a[0] = evaltyp(t_INT) | evallg(lgefint(b));
     161      170848 :   affii(b, a); return a;
     162             : }
     163             : 
     164             : /* assume a > 0, initialized by setloop. Do a++ */
     165             : static GEN
     166    85239514 : incpos(GEN a)
     167             : {
     168    85239514 :   long i, l = lgefint(a);
     169    85239519 :   for (i=2; i<l; i++)
     170    85244277 :     if (++uel(a,i)) return a;
     171           3 :   a[l] = 1; l++;
     172           3 :   a[0]=evaltyp(t_INT) | _evallg(l);
     173           3 :   a[1]=evalsigne(1) | evallgefint(l);
     174           3 :   return a;
     175             : }
     176             : 
     177             : /* assume a < 0, initialized by setloop. Do a++ */
     178             : static GEN
     179       10952 : incneg(GEN a)
     180             : {
     181       10952 :   long i, l = lgefint(a)-1;
     182       10952 :   if (uel(a,2)--)
     183             :   {
     184       10948 :     if (!a[l]) /* implies l = 2 */
     185             :     {
     186         604 :       a[0] = evaltyp(t_INT) | _evallg(2);
     187         604 :       a[1] = evalsigne(0) | evallgefint(2);
     188             :     }
     189       10948 :     return a;
     190             :   }
     191           5 :   for (i=3; i<=l; i++)
     192           5 :     if (uel(a,i)--) break;
     193           4 :   if (!a[l])
     194             :   {
     195           4 :     a[0] = evaltyp(t_INT) | _evallg(l);
     196           4 :     a[1] = evalsigne(-1) | evallgefint(l);
     197             :   }
     198           4 :   return a;
     199             : }
     200             : 
     201             : /* assume a initialized by setloop. Do a++ */
     202             : GEN
     203    85582452 : incloop(GEN a)
     204             : {
     205    85582452 :   switch(signe(a))
     206             :   {
     207      330884 :     case 0:
     208      330884 :       a[0]=evaltyp(t_INT) | _evallg(3);
     209      330884 :       a[1]=evalsigne(1) | evallgefint(3);
     210      330884 :       a[2]=1; return a;
     211       10952 :     case -1: return incneg(a);
     212    85240616 :     default: return incpos(a);
     213             :   }
     214             : }
     215             : 
     216             : INLINE GEN
     217  1553678711 : adduispec(ulong s, GEN x, long nx)
     218             : {
     219             :   GEN  zd;
     220             :   long lz;
     221             : 
     222  1553678711 :   if (nx == 1) return adduu(uel(x,0), s);
     223   264328373 :   lz = nx+3; zd = cgeti(lz);
     224   263796092 :   if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
     225     4971049 :     zd[lz-1]=1;
     226             :   else
     227   258831916 :     lz--;
     228   263802965 :   zd[1] = evalsigne(1) | evallgefint(lz);
     229   263802965 :   return zd;
     230             : }
     231             : 
     232             : GEN
     233   304995423 : adduispec_offset(ulong s, GEN x, long offset, long nx)
     234             : {
     235   304995423 :   GEN xd=x+2+offset;
     236   433984502 :   while (nx && *(xd+nx-1)==0) nx--;
     237   304995423 :   if (!nx) return utoi(s);
     238   281758273 :   return adduispec(s,xd,nx);
     239             : }
     240             : 
     241             : INLINE GEN
     242  2332287812 : addiispec(GEN x, GEN y, long nx, long ny)
     243             : {
     244             :   GEN zd;
     245             :   long lz;
     246             : 
     247  2332287812 :   if (nx < ny) swapspec(x,y, nx,ny);
     248  2332287812 :   if (ny == 1) return adduispec(*y,x,nx);
     249  1137429124 :   lz = nx+3; zd = cgeti(lz);
     250             : 
     251  1133292435 :   if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
     252    15666686 :     zd[lz-1]=1;
     253             :   else
     254  1118622088 :     lz--;
     255             : 
     256  1134288774 :   zd[1] = evalsigne(1) | evallgefint(lz);
     257  1134288774 :   return zd;
     258             : }
     259             : 
     260             : /* assume x >= y */
     261             : INLINE GEN
     262  1155843150 : subiuspec(GEN x, ulong s, long nx)
     263             : {
     264             :   GEN zd;
     265             :   long lz;
     266             : 
     267  1155843150 :   if (nx == 1) return utoi(x[0] - s);
     268             : 
     269   118275157 :   lz = nx + 2; zd = cgeti(lz);
     270   118304524 :   mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
     271   118305373 :   if (! zd[lz - 1]) { --lz; }
     272             : 
     273   118305373 :   zd[1] = evalsigne(1) | evallgefint(lz);
     274   118305373 :   return zd;
     275             : }
     276             : 
     277             : /* assume x > y */
     278             : INLINE GEN
     279  2456412523 : subiispec(GEN x, GEN y, long nx, long ny)
     280             : {
     281             :   GEN zd;
     282             :   long lz;
     283  2456412523 :   if (ny==1) return subiuspec(x,*y,nx);
     284  1325064089 :   lz = nx+2; zd = cgeti(lz);
     285             : 
     286  1314960336 :   mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
     287  1614153627 :   while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
     288             : 
     289  1316382250 :   zd[1] = evalsigne(1) | evallgefint(lz);
     290  1316382250 :   return zd;
     291             : }
     292             : 
     293             : static void
     294   386274995 : roundr_up_ip(GEN x, long l)
     295             : {
     296   386274995 :   long i = l;
     297             :   for(;;)
     298             :   {
     299   388536099 :     if (++((ulong*)x)[--i]) break;
     300     2487975 :     if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
     301             :   }
     302   386279052 : }
     303             : 
     304             : void
     305   255186027 : affir(GEN x, GEN y)
     306             : {
     307   255186027 :   const long s = signe(x), ly = lg(y);
     308             :   long lx, sh, i;
     309             : 
     310   255186027 :   if (!s)
     311             :   {
     312     5531459 :     y[1] = evalexpo(-bit_accuracy(ly));
     313     5531327 :     return;
     314             :   }
     315   249654568 :   lx = lgefint(x); sh = bfffo(*int_MSW(x));
     316   249654568 :   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
     317   249661361 :   if (sh) {
     318   245587630 :     if (lx <= ly)
     319             :     {
     320   582336126 :       for (i=lx; i<ly; i++) y[i]=0;
     321   131323026 :       mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
     322   131323550 :       xmpn_mirror(LIMBS(y),lx-2);
     323   131324649 :       return;
     324             :     }
     325   114264604 :     mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
     326   114264624 :     uel(y,2) |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
     327   114264624 :     xmpn_mirror(LIMBS(y),ly-2);
     328             :     /* lx > ly: round properly */
     329   114264941 :     if ((uel(x,lx-ly+1)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
     330             :   }
     331             :   else {
     332     4073731 :     GEN xd=int_MSW(x);
     333     4073731 :     if (lx <= ly)
     334             :     {
     335     7608801 :       for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
     336     6288921 :       for (   ; i<ly; i++) y[i]=0;
     337     1986686 :       return;
     338             :     }
     339     6861910 :     for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
     340             :     /* lx > ly: round properly */
     341     2087045 :     if (uel(x,lx-ly+1) & HIGHBIT) roundr_up_ip(y, ly);
     342             :   }
     343             : }
     344             : 
     345             : INLINE GEN
     346   512488730 : shiftispec(GEN x, long nx, long n)
     347             : {
     348             :   long ny,m;
     349             :   GEN yd, y;
     350             : 
     351   512488730 :   if (!n) return icopyspec(x, nx);
     352             : 
     353   499377885 :   if (n > 0)
     354             :   {
     355   330290470 :     long d = dvmdsBIL(n, &m);
     356             :     long i;
     357             : 
     358   330291644 :     ny = nx + d + (m!=0);
     359   330291644 :     y = cgeti(ny + 2); yd = y + 2;
     360   904881281 :     for (i=0; i<d; i++) yd[i] = 0;
     361             : 
     362   329574014 :     if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
     363             :     else
     364             :     {
     365   327866768 :       ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
     366   328187465 :       if (carryd) yd[ny - 1] = carryd;
     367   302022672 :       else ny--;
     368             :     }
     369             :   }
     370             :   else
     371             :   {
     372   169087415 :     long d = dvmdsBIL(-n, &m);
     373             : 
     374   169116982 :     ny = nx - d;
     375   169116982 :     if (ny < 1) return gen_0;
     376   167994516 :     y = cgeti(ny + 2); yd = y + 2;
     377             : 
     378   167971837 :     if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
     379             :     else
     380             :     {
     381   166787820 :       mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
     382   166788222 :       if (yd[ny - 1] == 0)
     383             :       {
     384    20981196 :         if (ny == 1) return gc_const((pari_sp)(yd + 1), gen_0);
     385    19062446 :         ny--;
     386             :       }
     387             :     }
     388             :   }
     389   495895672 :   y[1] = evalsigne(1)|evallgefint(ny + 2);
     390   495895672 :   return y;
     391             : }
     392             : 
     393             : GEN
     394    46562993 : mantissa2nr(GEN x, long n)
     395             : {
     396    46562993 :   long ly, i, m, s = signe(x), lx = lg(x);
     397             :   GEN y;
     398    46562993 :   if (!s) return gen_0;
     399    46561692 :   if (!n)
     400             :   {
     401    22966292 :     y = cgeti(lx);
     402    22966026 :     y[1] = evalsigne(s) | evallgefint(lx);
     403    22966026 :     xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
     404    22966222 :     return y;
     405             :   }
     406    23595400 :   if (n > 0)
     407             :   {
     408       66548 :     GEN z = (GEN)avma;
     409       66548 :     long d = dvmdsBIL(n, &m);
     410             : 
     411       66548 :     ly = lx+d; y = new_chunk(ly);
     412      270069 :     for ( ; d; d--) *--z = 0;
     413       68240 :     if (!m) for (i=2; i<lx; i++) y[i]=x[i];
     414             :     else
     415             :     {
     416       65946 :       register const ulong sh = BITS_IN_LONG - m;
     417       65946 :       shift_left(y,x, 2,lx-1, 0,m);
     418       65946 :       i = uel(x,2) >> sh;
     419             :       /* Extend y on the left? */
     420       65946 :       if (i) { ly++; y = new_chunk(1); y[2] = i; }
     421             :     }
     422             :   }
     423             :   else
     424             :   {
     425    23528852 :     ly = lx - dvmdsBIL(-n, &m);
     426    23528697 :     if (ly<3) return gen_0;
     427    23528697 :     y = new_chunk(ly);
     428    23528462 :     if (m) {
     429    23457593 :       shift_right(y,x, 2,ly, 0,m);
     430    23457941 :       if (y[2] == 0)
     431             :       {
     432           0 :         if (ly==3) return gc_const((pari_sp)(y+3), gen_0);
     433           0 :         ly--; set_avma((pari_sp)(++y));
     434             :       }
     435             :     } else {
     436      221167 :       for (i=2; i<ly; i++) y[i]=x[i];
     437             :     }
     438             :   }
     439    23595356 :   xmpn_mirror(LIMBS(y),ly-2);
     440    23595436 :   y[1] = evalsigne(s)|evallgefint(ly);
     441    23595436 :   y[0] = evaltyp(t_INT)|evallg(ly); return y;
     442             : }
     443             : 
     444             : GEN
     445     2056795 : truncr(GEN x)
     446             : {
     447             :   long s, e, d, m, i;
     448             :   GEN y;
     449     2056795 :   if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
     450     1653605 :   d = nbits2lg(e+1); m = remsBIL(e);
     451     1653605 :   if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
     452             : 
     453     1653605 :   y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
     454     1653605 :   if (++m == BITS_IN_LONG)
     455        1316 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     456             :   else
     457             :   {
     458     1652955 :     GEN z=cgeti(d);
     459     3348259 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     460     1652955 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     461     1652955 :     set_avma((pari_sp)y);
     462             :   }
     463     1653605 :   return y;
     464             : }
     465             : 
     466             : /* integral part */
     467             : GEN
     468      867123 : floorr(GEN x)
     469             : {
     470             :   long e, d, m, i, lx;
     471             :   GEN y;
     472      867123 :   if (signe(x) >= 0) return truncr(x);
     473      311207 :   if ((e=expo(x)) < 0) return gen_m1;
     474      128249 :   d = nbits2lg(e+1); m = remsBIL(e);
     475      128249 :   lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
     476      128249 :   y = cgeti(d+1);
     477      128249 :   if (++m == BITS_IN_LONG)
     478             :   {
     479        2766 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     480        8001 :     i=d; while (i<lx && !x[i]) i++;
     481        1109 :     if (i==lx) goto END;
     482             :   }
     483             :   else
     484             :   {
     485      127140 :     GEN z=cgeti(d);
     486      298953 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     487      127140 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     488      127140 :     if (uel(x,d-1)<<m == 0)
     489             :     {
     490      406150 :       i=d; while (i<lx && !x[i]) i++;
     491       35807 :       if (i==lx) goto END;
     492             :     }
     493             :   }
     494      107384 :   if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
     495           0 :     y[d++]=1;
     496      107384 : END:
     497      128249 :   y[1] = evalsigne(-1) | evallgefint(d);
     498      128249 :   return y;
     499             : }
     500             : 
     501             : INLINE int
     502  2944385363 : cmpiispec(GEN x, GEN y, long lx, long ly)
     503             : {
     504  2944385363 :   if (lx < ly) return -1;
     505  2818768536 :   if (lx > ly) return  1;
     506  2694888854 :   return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     507             : }
     508             : 
     509             : INLINE int
     510   189412974 : equaliispec(GEN x, GEN y, long lx, long ly)
     511             : {
     512   189412974 :   if (lx != ly) return 0;
     513   189298972 :   return !mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     514             : }
     515             : 
     516             : /***********************************************************************/
     517             : /**                                                                   **/
     518             : /**                          MULTIPLICATION                           **/
     519             : /**                                                                   **/
     520             : /***********************************************************************/
     521             : /* assume ny > 0 */
     522             : INLINE GEN
     523  3055978736 : muluispec(ulong x, GEN y, long ny)
     524             : {
     525  3055978736 :   if (ny == 1)
     526  2360922852 :     return muluu(x, *y);
     527             :   else
     528             :   {
     529   695055884 :     long lz = ny+3;
     530   695055884 :     GEN z = cgeti(lz);
     531   693800421 :     ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
     532   694155848 :     if (hi) { z[lz - 1] = hi; } else lz--;
     533   694155848 :     z[1] = evalsigne(1) | evallgefint(lz);
     534   694155848 :     return z;
     535             :   }
     536             : }
     537             : 
     538             : /* a + b*|y| */
     539             : GEN
     540      700685 : addumului(ulong a, ulong b, GEN y)
     541             : {
     542             :   GEN z;
     543             :   long i, lz;
     544             :   ulong hi;
     545      700685 :   if (!b || !signe(y)) return utoi(a);
     546      700567 :   lz = lgefint(y)+1;
     547      700567 :   z = cgeti(lz);
     548      700567 :   z[2]=a;
     549     1401134 :   for(i=3;i<lz;i++) z[i]=0;
     550      700567 :   hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
     551      700567 :   if (hi) z[lz-1]=hi; else lz--;
     552      700567 :   z[1] = evalsigne(1) | evallgefint(lz);
     553      700567 :   return gc_const((pari_sp)z, z);
     554             : }
     555             : 
     556             : /***********************************************************************/
     557             : /**                                                                   **/
     558             : /**                          DIVISION                                 **/
     559             : /**                                                                   **/
     560             : /***********************************************************************/
     561             : 
     562             : ulong
     563  1001468062 : umodiu(GEN y, ulong x)
     564             : {
     565  1001468062 :   long sy=signe(y);
     566             :   ulong hi;
     567  1001468062 :   if (!x) pari_err_INV("umodiu",gen_0);
     568  1002791745 :   if (!sy) return 0;
     569   750965150 :   hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
     570   750965150 :   if (!hi) return 0;
     571   740350540 :   return (sy > 0)? hi: x - hi;
     572             : }
     573             : 
     574             : /* return |y| \/ x */
     575             : GEN
     576   187550794 : absdiviu_rem(GEN y, ulong x, ulong *rem)
     577             : {
     578             :   long ly;
     579             :   GEN z;
     580             : 
     581   187550794 :   if (!x) pari_err_INV("absdiviu_rem",gen_0);
     582   187557882 :   if (!signe(y)) { *rem = 0; return gen_0; }
     583             : 
     584   180538719 :   ly = lgefint(y);
     585   180538719 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
     586             : 
     587   173685865 :   z = cgeti(ly);
     588   173619686 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     589   173640797 :   if (z [ly - 1] == 0) ly--;
     590   173640797 :   z[1] = evallgefint(ly) | evalsigne(1);
     591   173640797 :   return z;
     592             : }
     593             : 
     594             : GEN
     595    60768700 : divis_rem(GEN y, long x, long *rem)
     596             : {
     597    60768700 :   long sy=signe(y),ly,s;
     598             :   GEN z;
     599             : 
     600    60768700 :   if (!x) pari_err_INV("divis_rem",gen_0);
     601    60768719 :   if (!sy) { *rem = 0; return gen_0; }
     602    45926198 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     603             : 
     604    45926198 :   ly = lgefint(y);
     605    45926198 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
     606             : 
     607    33216231 :   z = cgeti(ly);
     608    33216234 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     609    33216234 :   if (sy<0) *rem = -  *rem;
     610    33216234 :   if (z[ly - 1] == 0) ly--;
     611    33216234 :   z[1] = evallgefint(ly) | evalsigne(s);
     612    33216234 :   return z;
     613             : }
     614             : 
     615             : GEN
     616     1080974 : divis(GEN y, long x)
     617             : {
     618     1080974 :   long sy=signe(y),ly,s;
     619             :   GEN z;
     620             : 
     621     1080974 :   if (!x) pari_err_INV("divis",gen_0);
     622     1080974 :   if (!sy) return gen_0;
     623     1080966 :   if (x<0) { s = -sy; x = -x; } else s=sy;
     624             : 
     625     1080966 :   ly = lgefint(y);
     626     1080966 :   if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
     627             : 
     628     1080966 :   z = cgeti(ly);
     629     1080966 :   (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     630     1080966 :   if (z[ly - 1] == 0) ly--;
     631     1080966 :   z[1] = evallgefint(ly) | evalsigne(s);
     632     1080966 :   return z;
     633             : }
     634             : 
     635             : /* We keep llx bits of x and lly bits of y*/
     636             : static GEN
     637    74467770 : divrr_with_gmp(GEN x, GEN y)
     638             : {
     639    74467770 :   long lx=RNLIMBS(x),ly=RNLIMBS(y);
     640    74467770 :   long lw=minss(lx,ly);
     641    74464252 :   long lly=minss(lw+1,ly);
     642    74464534 :   GEN  w=cgetr(lw+2);
     643    74441343 :   long lu=lw+lly;
     644    74441343 :   long llx=minss(lu,lx);
     645    74433664 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     646    74399547 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     647             :   mp_limb_t *q,*r;
     648    74360059 :   long e=expo(x)-expo(y);
     649    74360059 :   long sx=signe(x),sy=signe(y);
     650    74360059 :   xmpn_mirrorcopy(z,RLIMBS(y),lly);
     651    74390015 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     652    74411587 :   xmpn_zero(u,lu-llx);
     653    74462422 :   q = (mp_limb_t *)new_chunk(lw+1);
     654    74464632 :   r = (mp_limb_t *)new_chunk(lly);
     655             : 
     656    74422604 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     657             : 
     658             :   /*Round up: This is not exactly correct we should test 2*r>z*/
     659    74547903 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     660    36923968 :     mpn_add_1(q,q,lw+1,1);
     661             : 
     662    74547968 :   xmpn_mirrorcopy(RLIMBS(w),q,lw);
     663             : 
     664    74540455 :   if (q[lw] == 0) e--;
     665    32290018 :   else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
     666          14 :   else { w[2] = HIGHBIT; e++; }
     667    74538899 :   if (sy < 0) sx = -sx;
     668    74538899 :   w[1] = evalsigne(sx) | evalexpo(e);
     669    74516326 :   return gc_const((pari_sp)w, w);
     670             : }
     671             : 
     672             : /* We keep llx bits of x and lly bits of y*/
     673             : static GEN
     674     7056576 : divri_with_gmp(GEN x, GEN y)
     675             : {
     676     7056576 :   long llx=RNLIMBS(x),ly=NLIMBS(y);
     677     7056576 :   long lly=minss(llx+1,ly);
     678     7056647 :   GEN  w=cgetr(llx+2);
     679     7055572 :   long lu=llx+lly, ld=ly-lly;
     680     7055572 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     681     7054770 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     682             :   mp_limb_t *q,*r;
     683     7053642 :   long sh=bfffo(y[ly+1]);
     684     7053642 :   long e=expo(x)-expi(y);
     685     7054951 :   long sx=signe(x),sy=signe(y);
     686     7054951 :   if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
     687      196702 :   else xmpn_copy(z,LIMBS(y)+ld,lly);
     688     7055174 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     689     7056630 :   xmpn_zero(u,lu-llx);
     690     7057441 :   q = (mp_limb_t *)new_chunk(llx+1);
     691     7056789 :   r = (mp_limb_t *)new_chunk(lly);
     692             : 
     693     7055762 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     694             : 
     695             :   /*Round up: This is not exactly correct we should test 2*r>z*/
     696     7058003 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     697     2896138 :     mpn_add_1(q,q,llx+1,1);
     698             : 
     699     7058009 :   xmpn_mirrorcopy(RLIMBS(w),q,llx);
     700             : 
     701     7057920 :   if (q[llx] == 0) e--;
     702     4424231 :   else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
     703           0 :   else { w[2] = HIGHBIT; e++; }
     704     7057871 :   if (sy < 0) sx = -sx;
     705     7057871 :   w[1] = evalsigne(sx) | evalexpo(e);
     706     7057549 :   return gc_const((pari_sp)w, w);
     707             : }
     708             : 
     709             : GEN
     710    34312340 : divri(GEN x, GEN y)
     711             : {
     712    34312340 :   long  s = signe(y);
     713             : 
     714    34312340 :   if (!s) pari_err_INV("divri",gen_0);
     715    34312390 :   if (!signe(x)) return real_0_bit(expo(x) - expi(y));
     716    34163620 :   if (!is_bigint(y)) {
     717    27107541 :     GEN z = divru(x, y[2]);
     718    27107780 :     if (s < 0) togglesign(z);
     719    27107803 :     return z;
     720             :   }
     721     7055914 :   return divri_with_gmp(x,y);
     722             : }
     723             : 
     724             : GEN
     725   175145829 : divrr(GEN x, GEN y)
     726             : {
     727   175145829 :   long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
     728             :   ulong y0,y1;
     729             :   GEN r, r1;
     730             : 
     731   175145829 :   if (!sy) pari_err_INV("divrr",y);
     732   175150314 :   e = expo(x) - expo(y);
     733   175150314 :   if (!sx) return real_0_bit(e);
     734   173277637 :   if (sy<0) sx = -sx;
     735             : 
     736   173277637 :   lx=lg(x); ly=lg(y);
     737   173277637 :   if (ly==3)
     738             :   {
     739    98816924 :     ulong k = x[2], l = (lx>3)? x[3]: 0;
     740             :     LOCAL_HIREMAINDER;
     741    98816924 :     if (k < uel(y,2)) e--;
     742             :     else
     743             :     {
     744    45888829 :       l >>= 1; if (k&1) l |= HIGHBIT;
     745    45888829 :       k >>= 1;
     746             :     }
     747    98816924 :     hiremainder = k; k = divll(l,y[2]);
     748    98816924 :     if (hiremainder > (uel(y,2) >> 1) && !++k) { k = HIGHBIT; e++; }
     749    98816924 :     r = cgetr(3);
     750    98816760 :     r[1] = evalsigne(sx) | evalexpo(e);
     751    98816671 :     r[2] = k; return r;
     752             :   }
     753             : 
     754    74460713 :   if (ly>=DIVRR_GMP_LIMIT)
     755    74470615 :     return divrr_with_gmp(x,y);
     756             : 
     757           0 :   lr = minss(lx,ly); r = new_chunk(lr);
     758           0 :   r1 = r-1;
     759           0 :   r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
     760           0 :   r1[lr] = (lx>ly)? x[lr]: 0;
     761           0 :   y0 = y[2]; y1 = y[3];
     762           0 :   for (i=0; i<lr-1; i++)
     763             :   { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
     764             :     ulong k, qp;
     765             :     LOCAL_HIREMAINDER;
     766             :     LOCAL_OVERFLOW;
     767             : 
     768           0 :     if (uel(r1,1) == y0) { qp = ULONG_MAX; k = addll(y0,r1[2]); }
     769             :     else
     770             :     {
     771           0 :       if (uel(r1,1) > y0) /* can't happen if i=0 */
     772             :       {
     773           0 :         GEN y1 = y+1;
     774           0 :         j = lr-i; r1[j] = subll(r1[j],y1[j]);
     775           0 :         for (j--; j>0; j--) r1[j] = subllx(r1[j],y1[j]);
     776           0 :         j=i; do uel(r,--j)++; while (j && !r[j]);
     777             :       }
     778           0 :       hiremainder = r1[1]; overflow = 0;
     779           0 :       qp = divll(r1[2],y0); k = hiremainder;
     780             :     }
     781           0 :     j = lr-i+1;
     782           0 :     if (!overflow)
     783             :     {
     784             :       long k3, k4;
     785           0 :       k3 = mulll(qp,y1);
     786           0 :       if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
     787           0 :         k4 = subll(hiremainder,k);
     788             :       else
     789             :       {
     790           0 :         k3 = subll(k3, r1[3]);
     791           0 :         k4 = subllx(hiremainder,k);
     792             :       }
     793           0 :       while (!overflow && k4) { qp--; k3=subll(k3,y1); k4=subllx(k4,y0); }
     794             :     }
     795           0 :     if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
     796           0 :     for (j--; j>1; j--)
     797             :     {
     798           0 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     799           0 :       hiremainder += overflow;
     800             :     }
     801           0 :     if (uel(r1,1) != hiremainder)
     802             :     {
     803           0 :       if (uel(r1,1) < hiremainder)
     804             :       {
     805           0 :         qp--;
     806           0 :         j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
     807           0 :         for (j--; j>1; j--) r1[j] = addllx(r1[j], y[j]);
     808             :       }
     809             :       else
     810             :       {
     811           0 :         uel(r1,1) -= hiremainder;
     812           0 :         while (r1[1])
     813             :         {
     814           0 :           qp++; if (!qp) { j=i; do uel(r,--j)++; while (j && !r[j]); }
     815           0 :           j = lr-i-(lr-i>=ly); r1[j] = subll(r1[j],y[j]);
     816           0 :           for (j--; j>1; j--) r1[j] = subllx(r1[j],y[j]);
     817           0 :           uel(r1,1) -= overflow;
     818             :         }
     819             :       }
     820             :     }
     821           0 :     *++r1 = qp;
     822             :   }
     823             :   /* i = lr-1 */
     824             :   /* round correctly */
     825           0 :   if (uel(r1,1) > (y0>>1))
     826             :   {
     827           0 :     j=i; do uel(r,--j)++; while (j && !r[j]);
     828             :   }
     829           0 :   r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
     830           0 :   if (r[0] == 0) e--;
     831           0 :   else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
     832             :   else { /* possible only when rounding up to 0x2 0x0 ... */
     833           0 :     r[2] = (long)HIGHBIT; e++;
     834             :   }
     835           0 :   r[0] = evaltyp(t_REAL)|evallg(lr);
     836           0 :   r[1] = evalsigne(sx) | evalexpo(e);
     837           0 :   return r;
     838             : }
     839             : 
     840             : /* Integer division x / y: such that sign(r) = sign(x)
     841             :  *   if z = ONLY_REM return remainder, otherwise return quotient
     842             :  *   if z != NULL set *z to remainder
     843             :  *   *z is the last object on stack (and thus can be disposed of with cgiv
     844             :  *   instead of gerepile)
     845             :  * If *z is zero, we put gen_0 here and no copy.
     846             :  * space needed: lx + ly
     847             :  */
     848             : GEN
     849  1189532702 : dvmdii(GEN x, GEN y, GEN *z)
     850             : {
     851  1189532702 :   long sx=signe(x),sy=signe(y);
     852             :   long lx, ly, lq;
     853             :   pari_sp av;
     854             :   GEN r,q;
     855             : 
     856  1189532702 :   if (!sy) pari_err_INV("dvmdii",y);
     857  1190316377 :   if (!sx)
     858             :   {
     859    27289104 :     if (!z || z == ONLY_REM) return gen_0;
     860     8302979 :     *z=gen_0; return gen_0;
     861             :   }
     862  1163027273 :   lx=lgefint(x);
     863  1163027273 :   ly=lgefint(y); lq=lx-ly;
     864  1163027273 :   if (lq <= 0)
     865             :   {
     866   659022234 :     if (lq == 0)
     867             :     {
     868   616988562 :       long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
     869   616988562 :       if (s>0) goto DIVIDE;
     870   242940160 :       if (s==0)
     871             :       {
     872    13377741 :         if (z == ONLY_REM) return gen_0;
     873     6966594 :         if (z) *z = gen_0;
     874     6966594 :         if (sx < 0) sy = -sy;
     875     6966594 :         return stoi(sy);
     876             :       }
     877             :     }
     878   271596091 :     if (z == ONLY_REM) return icopy(x);
     879     4817176 :     if (z) *z = icopy(x);
     880     4817176 :     return gen_0;
     881             :   }
     882   504005039 : DIVIDE: /* quotient is non-zero */
     883   878053441 :   av=avma; if (sx<0) sy = -sy;
     884   878053441 :   if (ly==3)
     885             :   {
     886   325890929 :     ulong lq = lx;
     887             :     ulong si;
     888   325890929 :     q  = cgeti(lq);
     889   325655728 :     si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
     890   325680095 :     if (q[lq - 1] == 0) lq--;
     891   325680095 :     if (z == ONLY_REM)
     892             :     {
     893   256706103 :       if (!si) return gc_const(av, gen_0);
     894   226143134 :       set_avma(av); r = cgeti(3);
     895   226102202 :       r[1] = evalsigne(sx) | evallgefint(3);
     896   226102202 :       r[2] = si; return r;
     897             :     }
     898    68973992 :     q[1] = evalsigne(sy) | evallgefint(lq);
     899    68973992 :     if (!z) return q;
     900    66425125 :     if (!si) { *z=gen_0; return q; }
     901    25292442 :     r=cgeti(3);
     902    25293054 :     r[1] = evalsigne(sx) | evallgefint(3);
     903    25293054 :     r[2] = si; *z=r; return q;
     904             :   }
     905   552162512 :   if (z == ONLY_REM)
     906             :   {
     907   541108367 :     ulong lr = lgefint(y);
     908   541108367 :     ulong lq = lgefint(x)-lgefint(y)+3;
     909   541108367 :     GEN r = cgeti(lr);
     910   527629125 :     GEN q = cgeti(lq);
     911   521687077 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     912   549551872 :     if (!r[lr - 1])
     913             :     {
     914   379626616 :       while(lr>2 && !r[lr - 1]) lr--;
     915   157635215 :       if (lr == 2) return gc_const(av, gen_0); /* exact division */
     916             :     }
     917   540817253 :     r[1] = evalsigne(sx) | evallgefint(lr);
     918   540817253 :     return gc_const((pari_sp)r, r);
     919             :   }
     920             :   else
     921             :   {
     922    13080072 :     ulong lq = lgefint(x)-lgefint(y)+3;
     923    13080072 :     ulong lr = lgefint(y);
     924    13080072 :     GEN q = cgeti(lq);
     925    16484933 :     GEN r = cgeti(lr);
     926    16484184 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     927    16486254 :     if (q[lq - 1] == 0) lq--;
     928    16486254 :     q[1] = evalsigne(sy) | evallgefint(lq);
     929    16486254 :     if (!z) return gc_const((pari_sp)q, q);
     930    13650014 :     if (!r[lr - 1])
     931             :     {
     932    22473520 :       while(lr>2 && !r[lr - 1]) lr--;
     933     3375817 :       if (lr == 2) { *z = gen_0; return gc_const((pari_sp)q, q); } /* exact */
     934             :     }
     935    11106384 :     r[1] = evalsigne(sx) | evallgefint(lr);
     936    11106384 :     *z = gc_const((pari_sp)r, r); return q;
     937             :   }
     938             : }
     939             : 
     940             : /* Montgomery reduction.
     941             :  * N has k words, assume T >= 0 has less than 2k.
     942             :  * Return res := T / B^k mod N, where B = 2^BIL
     943             :  * such that 0 <= res < T/B^k + N  and  res has less than k words */
     944             : GEN
     945    22103781 : red_montgomery(GEN T, GEN N, ulong inv)
     946             : {
     947             :   pari_sp av;
     948             :   GEN Te, Td, Ne, Nd, scratch;
     949    22103781 :   ulong i, j, m, t, d, k = NLIMBS(N);
     950             :   int carry;
     951             :   LOCAL_HIREMAINDER;
     952             :   LOCAL_OVERFLOW;
     953             : 
     954    22103781 :   if (k == 0) return gen_0;
     955    22103781 :   d = NLIMBS(T); /* <= 2*k */
     956    22103781 :   if (d == 0) return gen_0;
     957             : #ifdef DEBUG
     958             :   if (d > 2*k) pari_err_BUG("red_montgomery");
     959             : #endif
     960    22103772 :   if (k == 1)
     961             :   { /* as below, special cased for efficiency */
     962         447 :     ulong n = uel(N,2);
     963         447 :     if (d == 1) {
     964         447 :       hiremainder = uel(T,2);
     965         447 :       m = hiremainder * inv;
     966         447 :       (void)addmul(m, n); /* t + m*n = 0 */
     967         447 :       return utoi(hiremainder);
     968             :     } else { /* d = 2 */
     969           0 :       hiremainder = uel(T,2);
     970           0 :       m = hiremainder * inv;
     971           0 :       (void)addmul(m, n); /* t + m*n = 0 */
     972           0 :       t = addll(hiremainder, uel(T,3));
     973           0 :       if (overflow) t -= n; /* t > n doesn't fit in 1 word */
     974           0 :       return utoi(t);
     975             :     }
     976             :   }
     977             :   /* assume k >= 2 */
     978    22103325 :   av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
     979             : 
     980             :   /* copy T to scratch space (pad with zeroes to 2k words) */
     981    21573072 :   Td = scratch;
     982    21573072 :   Te = T + 2;
     983   308544033 :   for (i=0; i < d     ; i++) *Td++ = *Te++;
     984    43534186 :   for (   ; i < (k<<1); i++) *Td++ = 0;
     985             : 
     986    21573072 :   Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
     987    21573072 :   Ne = N + 1;       /* 1 beyond end of N mantissa */
     988             : 
     989    21573072 :   carry = 0;
     990   166917585 :   for (i=0; i<k; i++) /* set T := T/B nod N, k times */
     991             :   {
     992   145344513 :     Td = Te; /* one beyond end of (new) T mantissa */
     993   145344513 :     Nd = Ne;
     994   145344513 :     hiremainder = *++Td;
     995   145344513 :     m = hiremainder * inv; /* solve T + m N = O(B) */
     996             : 
     997             :     /* set T := (T + mN) / B */
     998   145344513 :     Te = Td;
     999   145344513 :     (void)addmul(m, *++Nd); /* = 0 */
    1000  1292639028 :     for (j=1; j<k; j++)
    1001             :     {
    1002  1147294515 :       t = addll(addmul(m, *++Nd), *++Td);
    1003  1147294515 :       *Td = t;
    1004  1147294515 :       hiremainder += overflow;
    1005             :     }
    1006   145344513 :     t = addll(hiremainder, *++Td); *Td = t + carry;
    1007   145344513 :     carry = (overflow || (carry && *Td == 0));
    1008             :   }
    1009    21573072 :   if (carry)
    1010             :   { /* Td > N overflows (k+1 words), set Td := Td - N */
    1011       41812 :     GEN NE = N + k+1;
    1012       41812 :     Td = Te;
    1013       41812 :     Nd = Ne;
    1014       41812 :     t = subll(*++Td, *++Nd); *Td = t;
    1015      356148 :     while (Nd < NE) { t = subllx(*++Td, *++Nd); *Td = t; }
    1016             :   }
    1017             : 
    1018             :   /* copy result */
    1019    21573072 :   Td = (GEN)av - 1; /* *Td = high word of final result */
    1020    25138919 :   while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
    1021    21573072 :   k = Td - Te; if (!k) return gc_const(av, gen_0);
    1022    21573072 :   Td = (GEN)av - k; /* will write mantissa there */
    1023    21573072 :   (void)memmove(Td, Te+1, k*sizeof(long));
    1024    21573072 :   Td -= 2;
    1025    21573072 :   Td[0] = evaltyp(t_INT) | evallg(k+2);
    1026    21749330 :   Td[1] = evalsigne(1) | evallgefint(k+2);
    1027             : #ifdef DEBUG
    1028             : {
    1029             :   long l = NLIMBS(N), s = BITS_IN_LONG*l;
    1030             :   GEN R = int2n(s);
    1031             :   GEN res = remii(mulii(T, Fp_inv(R, N)), N);
    1032             :   if (k > lgefint(N)
    1033             :     || !equalii(remii(Td,N),res)
    1034             :     || cmpii(Td, addii(shifti(T, -s), N)) >= 0) pari_err_BUG("red_montgomery");
    1035             : }
    1036             : #endif
    1037    21749330 :   return gc_const((pari_sp)Td, Td);
    1038             : }
    1039             : 
    1040             : /* EXACT INTEGER DIVISION */
    1041             : 
    1042             : /* use undocumented GMP interface */
    1043             : static void
    1044    49993886 : GEN2mpz(mpz_t X, GEN x)
    1045             : {
    1046    49993886 :   long l = lgefint(x)-2;
    1047    49993886 :   X->_mp_alloc = l;
    1048    49993886 :   X->_mp_size = signe(x) > 0? l: -l;
    1049    49993886 :   X->_mp_d = LIMBS(x);
    1050    49993886 : }
    1051             : static void
    1052    24996945 : mpz2GEN(GEN z, mpz_t Z)
    1053             : {
    1054    24996945 :   long l = Z->_mp_size;
    1055    24996945 :   z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
    1056    24996945 : }
    1057             : 
    1058             : #ifdef mpn_divexact_1
    1059             : static GEN
    1060   283050029 : diviuexact_i(GEN x, ulong y)
    1061             : {
    1062   283050029 :   long l = lgefint(x);
    1063   283050029 :   GEN z = cgeti(l);
    1064   283018667 :   mpn_divexact_1(LIMBS(z), LIMBS(x), NLIMBS(x), y);
    1065   283019018 :   if (z[l-1] == 0) l--;
    1066   283019018 :   z[1] = evallgefint(l) | evalsigne(signe(x));
    1067   283019018 :   return z;
    1068             : }
    1069             : #elif 1 && !defined(_WIN64) /* mpz_divexact_ui is not LLP64 friendly   */
    1070             :                             /* assume y != 0 and the division is exact */
    1071             : static GEN
    1072             : diviuexact_i(GEN x, ulong y)
    1073             : {
    1074             :   long l = lgefint(x);
    1075             :   GEN z = cgeti(l);
    1076             :   mpz_t X, Z;
    1077             :   GEN2mpz(X, x);
    1078             :   Z->_mp_alloc = l-2;
    1079             :   Z->_mp_size  = l-2;
    1080             :   Z->_mp_d = LIMBS(z);
    1081             :   mpz_divexact_ui(Z, X, y);
    1082             :   mpz2GEN(z, Z); return z;
    1083             : }
    1084             : #else
    1085             : /* assume y != 0 and the division is exact */
    1086             : static GEN
    1087             : diviuexact_i(GEN x, ulong y)
    1088             : {
    1089             :   /*TODO: implement true exact division.*/
    1090             :   return divii(x,utoi(y));
    1091             : }
    1092             : #endif
    1093             : 
    1094             : GEN
    1095    40615755 : diviuexact(GEN x, ulong y)
    1096             : {
    1097             :   GEN z;
    1098    40615755 :   if (!signe(x)) return gen_0;
    1099    39505243 :   z = diviuexact_i(x, y);
    1100    39501291 :   if (lgefint(z) == 2) pari_err_OP("exact division", x, utoi(y));
    1101    39501320 :   return z;
    1102             : }
    1103             : 
    1104             : /* Find z such that x=y*z, knowing that y | x (unchecked) */
    1105             : GEN
    1106   317893453 : diviiexact(GEN x, GEN y)
    1107             : {
    1108             :   GEN z;
    1109   317893453 :   if (!signe(y)) pari_err_INV("diviiexact",y);
    1110   317907284 :   if (!signe(x)) return gen_0;
    1111   268524720 :   if (lgefint(y) == 3)
    1112             :   {
    1113   243530962 :     z = diviuexact_i(x, y[2]);
    1114   243517611 :     if (signe(y) < 0) togglesign(z);
    1115             :   }
    1116             :   else
    1117             :   {
    1118    24993758 :     long l = lgefint(x);
    1119             :     mpz_t X, Y, Z;
    1120    24993758 :     z = cgeti(l);
    1121    24996943 :     GEN2mpz(X, x);
    1122    24996943 :     GEN2mpz(Y, y);
    1123    24996943 :     Z->_mp_alloc = l-2;
    1124    24996943 :     Z->_mp_size  = l-2;
    1125    24996943 :     Z->_mp_d = LIMBS(z);
    1126    24996943 :     mpz_divexact(Z, X, Y);
    1127    24996945 :     mpz2GEN(z, Z);
    1128             :   }
    1129   268514746 :   if (lgefint(z) == 2) pari_err_OP("exact division", x, y);
    1130   268514478 :   return z;
    1131             : }
    1132             : 
    1133             : /* assume yz != and yz | x */
    1134             : GEN
    1135      174160 : diviuuexact(GEN x, ulong y, ulong z)
    1136             : {
    1137             :   long tmp[4];
    1138             :   ulong t;
    1139             :   LOCAL_HIREMAINDER;
    1140      174160 :   t = mulll(y, z);
    1141      174160 :   if (!hiremainder) return diviuexact(x, t);
    1142           0 :   tmp[0] = evaltyp(t_INT)|_evallg(4);
    1143           0 :   tmp[1] = evalsigne(1)|evallgefint(4);
    1144           0 :   tmp[2] = t;
    1145           0 :   tmp[3] = hiremainder;
    1146           0 :   return diviiexact(x, tmp);
    1147             : }
    1148             : 
    1149             : /********************************************************************/
    1150             : /**                                                                **/
    1151             : /**               INTEGER MULTIPLICATION                           **/
    1152             : /**                                                                **/
    1153             : /********************************************************************/
    1154             : 
    1155             : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
    1156             : GEN
    1157  3245083074 : muliispec(GEN x, GEN y, long nx, long ny)
    1158             : {
    1159             :   GEN zd;
    1160             :   long lz;
    1161             :   ulong hi;
    1162             : 
    1163  3245083074 :   if (nx < ny) swapspec(x,y, nx,ny);
    1164  3245083074 :   if (!ny) return gen_0;
    1165  3245083074 :   if (ny == 1) return muluispec((ulong)*y, x, nx);
    1166             : 
    1167   814304676 :   lz = nx+ny+2;
    1168   814304676 :   zd = cgeti(lz);
    1169   811871644 :   hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
    1170   820083861 :   if (!hi) lz--;
    1171             :   /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
    1172             : 
    1173   820083861 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1174   820083861 :   return zd;
    1175             : }
    1176             : GEN
    1177      196384 : muluui(ulong x, ulong y, GEN z)
    1178             : {
    1179      196384 :   long t, s = signe(z);
    1180             :   GEN r;
    1181             :   LOCAL_HIREMAINDER;
    1182             : 
    1183      196384 :   if (!x || !y || !signe(z)) return gen_0;
    1184      195988 :   t = mulll(x,y);
    1185      195988 :   if (!hiremainder)
    1186      195988 :     r = muluispec(t, z+2, lgefint(z)-2);
    1187             :   else
    1188             :   {
    1189             :     long tmp[2];
    1190           0 :     tmp[1] = hiremainder;
    1191           0 :     tmp[0] = t;
    1192           0 :     r = muliispec(z+2,tmp, lgefint(z)-2, 2);
    1193             :   }
    1194      195988 :   setsigne(r,s); return r;
    1195             : }
    1196             : 
    1197             : GEN
    1198   768910511 : sqrispec(GEN x, long nx)
    1199             : {
    1200             :   GEN zd;
    1201             :   long lz;
    1202             : 
    1203   768910511 :   if (!nx) return gen_0;
    1204   261907167 :   if (nx==1) return sqru(*x);
    1205             : 
    1206   176614714 :   lz = (nx<<1)+2;
    1207   176614714 :   zd = cgeti(lz);
    1208             : #ifdef mpn_sqr
    1209   171230240 :   mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
    1210             : #else
    1211             :   mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
    1212             : #endif
    1213   175360310 :   if (zd[lz-1]==0) lz--;
    1214             : 
    1215   175360310 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1216   175360310 :   return zd;
    1217             : }
    1218             : 
    1219             : INLINE GEN
    1220     8532020 : sqrispec_mirror(GEN x, long nx)
    1221             : {
    1222     8532020 :   GEN cx=new_chunk(nx);
    1223             :   GEN z;
    1224     8532025 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1225     8532027 :   z=sqrispec(cx, nx);
    1226     8532027 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1227     8532027 :   return z;
    1228             : }
    1229             : 
    1230             : /* leaves garbage on the stack. */
    1231             : INLINE GEN
    1232    31105664 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
    1233             : {
    1234             :   GEN cx, cy, z;
    1235    31105664 :   long s = 0;
    1236    52976084 :   while (nx && x[nx-1]==0) { nx--; s++; }
    1237    49360813 :   while (ny && y[ny-1]==0) { ny--; s++; }
    1238    31105664 :   cx=new_chunk(nx); cy=new_chunk(ny);
    1239    31105664 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1240    31105664 :   xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
    1241    31105664 :   z =  nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
    1242    31105664 :   if (s)
    1243             :   {
    1244     2412441 :     long i, lz = lgefint(z) + s;
    1245     2412441 :     (void)new_chunk(s);
    1246     2412441 :     z -= s;
    1247    42538010 :     for (i=0; i<s; i++) z[2+i]=0;
    1248     2412441 :     z[1] = evalsigne(1) | evallgefint(lz);
    1249     2412441 :     z[0] = evaltyp(t_INT) | evallg(lz);
    1250             :   }
    1251    31105664 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1252    31105664 :   return z;
    1253             : }
    1254             : 
    1255             : /* x % (2^n), assuming n >= 0 */
    1256             : GEN
    1257    31277136 : remi2n(GEN x, long n)
    1258             : {
    1259             :   ulong hi;
    1260    31277136 :   long l, k, lx, ly, sx = signe(x);
    1261             :   GEN z, xd, zd;
    1262             : 
    1263    31277136 :   if (!sx || !n) return gen_0;
    1264             : 
    1265    30995498 :   k = dvmdsBIL(n, &l);
    1266    31417869 :   lx = lgefint(x);
    1267    31417869 :   if (lx < k+3) return icopy(x);
    1268             : 
    1269    31013226 :   xd = x + (2 + k);
    1270             :   /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
    1271             :    *              ^--- initial xd  */
    1272    31013226 :   hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
    1273    31013226 :   if (!hi)
    1274             :   { /* strip leading zeroes from result */
    1275     2360013 :     xd--; while (k && !*xd) { k--; xd--; }
    1276     2316396 :     if (!k) return gen_0;
    1277     1234838 :     ly = k+2;
    1278             :   }
    1279             :   else
    1280    28696830 :     ly = k+3;
    1281             : 
    1282    29931668 :   zd = z = cgeti(ly);
    1283    29525142 :   *++zd = evalsigne(sx) | evallgefint(ly);
    1284   277220411 :   xd = x+1; for ( ;k; k--) *++zd = *++xd;
    1285    29525142 :   if (hi) *++zd = hi;
    1286    29525142 :   return z;
    1287             : }
    1288             : 
    1289             : /********************************************************************/
    1290             : /**                                                                **/
    1291             : /**                      INTEGER SQUARE ROOT                       **/
    1292             : /**                                                                **/
    1293             : /********************************************************************/
    1294             : 
    1295             : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
    1296             :  * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
    1297             :  * remainder is 0. R = NULL is allowed. */
    1298             : GEN
    1299     2861177 : sqrtremi(GEN a, GEN *r)
    1300             : {
    1301     2861177 :   long l, na = NLIMBS(a);
    1302             :   mp_size_t nr;
    1303             :   GEN S;
    1304     2861177 :   if (!na) {
    1305          60 :     if (r) *r = gen_0;
    1306          60 :     return gen_0;
    1307             :   }
    1308     2861117 :   l = (na + 5) >> 1; /* 2 + ceil(na/2) */
    1309     2861117 :   S = cgetipos(l);
    1310     2861117 :   if (r) {
    1311      350187 :     GEN R = cgeti(2 + na);
    1312      350187 :     nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
    1313      350187 :     if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
    1314       11252 :     else    { set_avma((pari_sp)S); R = gen_0; }
    1315      350187 :     *r = R;
    1316             :   }
    1317             :   else
    1318     2510930 :     (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
    1319     2861117 :   return S;
    1320             : }
    1321             : 
    1322             : /* compute sqrt(|a|), assuming a != 0 */
    1323             : GEN
    1324    32117466 : sqrtr_abs(GEN a)
    1325             : {
    1326             :   GEN res;
    1327             :   mp_limb_t *b, *c;
    1328    32117466 :   long l = RNLIMBS(a), e = expo(a), er = e>>1;
    1329             :   long n;
    1330    32117466 :   res = cgetr(2 + l);
    1331    32115634 :   res[1] = evalsigne(1) | evalexpo(er);
    1332    32117115 :   if (e&1)
    1333             :   {
    1334    15250196 :     b = (mp_limb_t *) new_chunk(l<<1);
    1335    15249982 :     xmpn_zero(b,l);
    1336    15250521 :     xmpn_mirrorcopy(b+l, RLIMBS(a), l);
    1337    15251132 :     c = (mp_limb_t *) new_chunk(l);
    1338    15250682 :     n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
    1339    15251802 :     if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
    1340             :   }
    1341             :   else
    1342             :   {
    1343             :     ulong u;
    1344    16866919 :     b = (mp_limb_t *) mantissa2nr(a,-1);
    1345    16866844 :     b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
    1346    16866844 :     b = (mp_limb_t *) new_chunk(l);
    1347    16866800 :     xmpn_zero(b,l+1); /* overwrites the former b[0] */
    1348    16866765 :     c = (mp_limb_t *) new_chunk(l + 1);
    1349    16866371 :     n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
    1350    16868533 :     u = (ulong)*c++;
    1351    16868533 :     if ( u&HIGHBIT || (u == ~HIGHBIT &&
    1352           0 :              (n>l || (n==l && mpn_cmp(b,c,l)>0))))
    1353     8277731 :       mpn_add_1(c,c,l,1);
    1354             :   }
    1355    32120395 :   xmpn_mirrorcopy(RLIMBS(res),c,l);
    1356    32120005 :   return gc_const((pari_sp)res, res);
    1357             : }
    1358             : 
    1359             : /* Normalize a non-negative integer */
    1360             : GEN
    1361   180955352 : int_normalize(GEN x, long known_zero_words)
    1362             : {
    1363   180955352 :   long i =  lgefint(x) - 1 - known_zero_words;
    1364  1235018024 :   for ( ; i > 1; i--)
    1365  1204819704 :     if (x[i]) { setlgefint(x, i+1); return x; }
    1366    30198320 :   x[1] = evalsigne(0) | evallgefint(2); return x;
    1367             : }
    1368             : 
    1369             : /********************************************************************
    1370             :  **                                                                **
    1371             :  **                           Base Conversion                      **
    1372             :  **                                                                **
    1373             :  ********************************************************************/
    1374             : 
    1375             : ulong *
    1376      396383 : convi(GEN x, long *l)
    1377             : {
    1378      396383 :   long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
    1379      396383 :   GEN str = cgetg(n+1, t_VECSMALL);
    1380      396383 :   unsigned char *res = (unsigned char*) GSTR(str);
    1381      396383 :   long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
    1382             :   long lz;
    1383             :   ulong *z;
    1384             :   long i, j;
    1385             :   unsigned char *t;
    1386      396383 :   while (!*res) {res++; llz--;} /*Strip leading zeros*/
    1387      396383 :   lz  = (8+llz)/9;
    1388      396383 :   z = (ulong*)new_chunk(1+lz);
    1389      396383 :   t=res+llz+9;
    1390      800107 :   for(i=0;i<llz-8;i+=9)
    1391             :   {
    1392             :     ulong s;
    1393      403724 :     t-=18;
    1394      403724 :     s=*t++;
    1395     3633516 :     for (j=1; j<9;j++)
    1396     3229792 :       s=10*s+*t++;
    1397      403724 :     *z++=s;
    1398             :   }
    1399      396383 :   if (i<llz)
    1400             :   {
    1401      392748 :     unsigned char *t = res;
    1402      392748 :     ulong s=*t++;
    1403     1113932 :     for (j=i+1; j<llz;j++)
    1404      721184 :       s=10*s+*t++;
    1405      392748 :     *z++=s;
    1406             :   }
    1407      396383 :   *l = lz;
    1408      396383 :   return z;
    1409             : }

Generated by: LCOV version 1.13