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 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.10.0 lcov report (development 19823-d80e022) Lines: 639 716 89.2 %
Date: 2016-12-03 05:49:13 Functions: 53 54 98.1 %
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      202793 : static void pari_gmp_free(void *ptr, size_t old_size){
      51      202793 :   (void)old_size; pari_free(ptr);
      52      202793 : }
      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         776 : pari_kernel_init(void)
      60             : {
      61         776 :   mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
      62         776 :   mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
      63         776 : }
      64             : 
      65             : void
      66        1284 : 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        1284 :   mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
      72        1284 :   if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
      73        1284 :   if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
      74        1284 :   if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
      75        1284 :   mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
      76        1284 : }
      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     1277613 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
      86             : {
      87     1277613 :   while (--n >= 0) x[n]=y[n];
      88     1277613 : }
      89             : 
      90             : INLINE void
      91   179546301 : xmpn_mirror(mp_limb_t *x, long n)
      92             : {
      93             :   long i;
      94  1037630380 :   for(i=0;i<(n>>1);i++)
      95             :   {
      96   858084079 :     ulong m=x[i];
      97   858084079 :     x[i]=x[n-1-i];
      98   858084079 :     x[n-1-i]=m;
      99             :   }
     100   179546301 : }
     101             : 
     102             : INLINE void
     103   273911998 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
     104             : {
     105             :   long i;
     106  2577444515 :   for(i=0;i<n;i++)
     107  2303532517 :     z[i]=x[n-1-i];
     108   273911998 : }
     109             : 
     110             : INLINE void
     111    57037877 : xmpn_zero(mp_ptr x, mp_size_t n)
     112             : {
     113    57037877 :   while (--n >= 0) x[n]=0;
     114    57037877 : }
     115             : 
     116             : INLINE GEN
     117    19943831 : icopy_ef(GEN x, long l)
     118             : {
     119    19943831 :   register long lx = lgefint(x);
     120    19943831 :   const GEN y = cgeti(l);
     121             : 
     122    19943812 :   while (--lx > 0) y[lx]=x[lx];
     123    19943812 :   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     2412749 : setloop(GEN a)
     141             : {
     142     2412749 :   pari_sp av = avma - 2 * sizeof(long);
     143     2412749 :   (void)cgetg(lgefint(a) + 3, t_VECSMALL);
     144     2412757 :   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      156748 : resetloop(GEN a, GEN b) {
     150      156748 :   a[0] = evaltyp(t_INT) | evallg(lgefint(b));
     151      156748 :   affii(b, a); return a;
     152             : }
     153             : 
     154             : /* assume a > 0, initialized by setloop. Do a++ */
     155             : static GEN
     156    32665445 : incpos(GEN a)
     157             : {
     158    32665445 :   long i, l = lgefint(a);
     159    32665450 :   for (i=2; i<l; i++)
     160    32665446 :     if (++uel(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           4 :   return a;
     165             : }
     166             : 
     167             : /* assume a < 0, initialized by setloop. Do a++ */
     168             : static GEN
     169        6184 : incneg(GEN a)
     170             : {
     171        6184 :   long i, l = lgefint(a)-1;
     172        6184 :   if (uel(a,2)--)
     173             :   {
     174        6180 :     if (!a[l]) /* implies l = 2 */
     175             :     {
     176         532 :       a[0] = evaltyp(t_INT) | _evallg(2);
     177         532 :       a[1] = evalsigne(0) | evallgefint(2);
     178             :     }
     179        6180 :     return a;
     180             :   }
     181           5 :   for (i=3; i<=l; i++)
     182           5 :     if (uel(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           4 :   return a;
     189             : }
     190             : 
     191             : /* assume a initialized by setloop. Do a++ */
     192             : GEN
     193    32671699 : incloop(GEN a)
     194             : {
     195    32671699 :   switch(signe(a))
     196             :   {
     197             :     case 0:
     198       20708 :       a[0]=evaltyp(t_INT) | _evallg(3);
     199       20708 :       a[1]=evalsigne(1) | evallgefint(3);
     200       20708 :       a[2]=1; return a;
     201        6184 :     case -1: return incneg(a);
     202    32644807 :     default: return incpos(a);
     203             :   }
     204             : }
     205             : 
     206             : INLINE GEN
     207  1080556221 : adduispec(ulong s, GEN x, long nx)
     208             : {
     209             :   GEN  zd;
     210             :   long lz;
     211             : 
     212  1080556221 :   if (nx == 1) return adduu(uel(x,0), s);
     213   285984949 :   lz = nx+3; zd = cgeti(lz);
     214   285368644 :   if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
     215     4916480 :     zd[lz-1]=1;
     216             :   else
     217   280572027 :     lz--;
     218   285488507 :   zd[1] = evalsigne(1) | evallgefint(lz);
     219   285488507 :   return zd;
     220             : }
     221             : 
     222             : GEN
     223   318608687 : adduispec_offset(ulong s, GEN x, long offset, long nx)
     224             : {
     225   318608687 :   GEN xd=x+2+offset;
     226   318608687 :   while (nx && *(xd+nx-1)==0) nx--;
     227   318608687 :   if (!nx) return utoi(s);
     228   300071892 :   return adduispec(s,xd,nx);
     229             : }
     230             : 
     231             : INLINE GEN
     232  1411799021 : addiispec(GEN x, GEN y, long nx, long ny)
     233             : {
     234             :   GEN zd;
     235             :   long lz;
     236             : 
     237  1411799021 :   if (nx < ny) swapspec(x,y, nx,ny);
     238  1411799021 :   if (ny == 1) return adduispec(*y,x,nx);
     239   670014160 :   lz = nx+3; zd = cgeti(lz);
     240             : 
     241   669386641 :   if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
     242    10902653 :     zd[lz-1]=1;
     243             :   else
     244   658966364 :     lz--;
     245             : 
     246   669869017 :   zd[1] = evalsigne(1) | evallgefint(lz);
     247   669869017 :   return zd;
     248             : }
     249             : 
     250             : /* assume x >= y */
     251             : INLINE GEN
     252   814334731 : subiuspec(GEN x, ulong s, long nx)
     253             : {
     254             :   GEN zd;
     255             :   long lz;
     256             : 
     257   814334731 :   if (nx == 1) return utoi(x[0] - s);
     258             : 
     259    74294345 :   lz = nx + 2; zd = cgeti(lz);
     260    74281694 :   mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
     261    74282156 :   if (! zd[lz - 1]) { --lz; }
     262             : 
     263    74282156 :   zd[1] = evalsigne(1) | evallgefint(lz);
     264    74282156 :   return zd;
     265             : }
     266             : 
     267             : /* assume x > y */
     268             : INLINE GEN
     269  1418819808 : subiispec(GEN x, GEN y, long nx, long ny)
     270             : {
     271             :   GEN zd;
     272             :   long lz;
     273  1418819808 :   if (ny==1) return subiuspec(x,*y,nx);
     274   613878784 :   lz = nx+2; zd = cgeti(lz);
     275             : 
     276   610726366 :   mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
     277   612510494 :   while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
     278             : 
     279   612510494 :   zd[1] = evalsigne(1) | evallgefint(lz);
     280   612510494 :   return zd;
     281             : }
     282             : 
     283             : static void
     284   266169333 : roundr_up_ip(GEN x, long l)
     285             : {
     286   266169333 :   long i = l;
     287             :   for(;;)
     288             :   {
     289   266430931 :     if (++((ulong*)x)[--i]) break;
     290      337403 :     if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
     291      261598 :   }
     292   266169333 : }
     293             : 
     294             : void
     295   100452939 : affir(GEN x, GEN y)
     296             : {
     297   100452939 :   const long s = signe(x), ly = lg(y);
     298             :   long lx, sh, i;
     299             : 
     300   100452939 :   if (!s)
     301             :   {
     302     1145709 :     y[1] = evalexpo(-bit_accuracy(ly));
     303     1145709 :     return;
     304             :   }
     305    99307230 :   lx = lgefint(x); sh = bfffo(*int_MSW(x));
     306    99307230 :   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
     307    99307230 :   if (sh) {
     308    98210780 :     if (lx <= ly)
     309             :     {
     310    65145821 :       for (i=lx; i<ly; i++) y[i]=0;
     311    65145821 :       mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
     312    65145821 :       xmpn_mirror(LIMBS(y),lx-2);
     313    65145821 :       return;
     314             :     }
     315    33064959 :     mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
     316    33064959 :     uel(y,2) |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
     317    33064959 :     xmpn_mirror(LIMBS(y),ly-2);
     318             :     /* lx > ly: round properly */
     319    33064959 :     if ((uel(x,lx-ly+1)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
     320             :   }
     321             :   else {
     322     1096450 :     GEN xd=int_MSW(x);
     323     1096450 :     if (lx <= ly)
     324             :     {
     325      640519 :       for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
     326      640519 :       for (   ; i<ly; i++) y[i]=0;
     327      640519 :       return;
     328             :     }
     329      455931 :     for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
     330             :     /* lx > ly: round properly */
     331      455931 :     if (uel(x,lx-ly+1) & HIGHBIT) roundr_up_ip(y, ly);
     332             :   }
     333             : }
     334             : 
     335             : INLINE GEN
     336   340352746 : shiftispec(GEN x, long nx, long n)
     337             : {
     338             :   long ny,m;
     339             :   GEN yd, y;
     340             : 
     341   340352746 :   if (!n) return icopyspec(x, nx);
     342             : 
     343   323413000 :   if (n > 0)
     344             :   {
     345   215317800 :     long d = dvmdsBIL(n, &m);
     346             :     long i;
     347             : 
     348   215303488 :     ny = nx + d + (m!=0);
     349   215303488 :     y = cgeti(ny + 2); yd = y + 2;
     350   215161531 :     for (i=0; i<d; i++) yd[i] = 0;
     351             : 
     352   215161531 :     if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
     353             :     else
     354             :     {
     355   214432814 :       ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
     356   214515557 :       if (carryd) yd[ny - 1] = carryd;
     357   202721752 :       else ny--;
     358             :     }
     359             :   }
     360             :   else
     361             :   {
     362   108095200 :     long d = dvmdsBIL(-n, &m);
     363             : 
     364   108095204 :     ny = nx - d;
     365   108095204 :     if (ny < 1) return gen_0;
     366   107754230 :     y = cgeti(ny + 2); yd = y + 2;
     367             : 
     368   107753414 :     if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
     369             :     else
     370             :     {
     371   107291848 :       mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
     372   107291907 :       if (yd[ny - 1] == 0)
     373             :       {
     374     8536937 :         if (ny == 1) { avma = (pari_sp)(yd + 1); return gen_0; }
     375     7616913 :         ny--;
     376             :       }
     377             :     }
     378             :   }
     379   322048861 :   y[1] = evalsigne(1)|evallgefint(ny + 2);
     380   322048861 :   return y;
     381             : }
     382             : 
     383             : GEN
     384    21959132 : mantissa2nr(GEN x, long n)
     385             : {
     386    21959132 :   long ly, i, m, s = signe(x), lx = lg(x);
     387             :   GEN y;
     388    21959132 :   if (!s) return gen_0;
     389    21957824 :   if (!n)
     390             :   {
     391     8254331 :     y = cgeti(lx);
     392     8254331 :     y[1] = evalsigne(s) | evallgefint(lx);
     393     8254331 :     xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
     394     8254331 :     return y;
     395             :   }
     396    13703493 :   if (n > 0)
     397             :   {
     398      296739 :     GEN z = (GEN)avma;
     399      296739 :     long d = dvmdsBIL(n, &m);
     400             : 
     401      296739 :     ly = lx+d; y = new_chunk(ly);
     402      296739 :     for ( ; d; d--) *--z = 0;
     403      296739 :     if (!m) for (i=2; i<lx; i++) y[i]=x[i];
     404             :     else
     405             :     {
     406      296196 :       register const ulong sh = BITS_IN_LONG - m;
     407      296196 :       shift_left(y,x, 2,lx-1, 0,m);
     408      296196 :       i = uel(x,2) >> sh;
     409             :       /* Extend y on the left? */
     410      296196 :       if (i) { ly++; y = new_chunk(1); y[2] = i; }
     411             :     }
     412             :   }
     413             :   else
     414             :   {
     415    13406754 :     ly = lx - dvmdsBIL(-n, &m);
     416    13406754 :     if (ly<3) return gen_0;
     417    13406754 :     y = new_chunk(ly);
     418    13406754 :     if (m) {
     419    13375505 :       shift_right(y,x, 2,ly, 0,m);
     420    13375505 :       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       31249 :       for (i=2; i<ly; i++) y[i]=x[i];
     427             :     }
     428             :   }
     429    13703493 :   xmpn_mirror(LIMBS(y),ly-2);
     430    13703493 :   y[1] = evalsigne(s)|evallgefint(ly);
     431    13703493 :   y[0] = evaltyp(t_INT)|evallg(ly); return y;
     432             : }
     433             : 
     434             : GEN
     435     1485907 : truncr(GEN x)
     436             : {
     437             :   long s, e, d, m, i;
     438             :   GEN y;
     439     1485907 :   if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
     440     1303626 :   d = nbits2prec(e+1); m = remsBIL(e);
     441     1303626 :   if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
     442             : 
     443     1303626 :   y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
     444     1303626 :   if (++m == BITS_IN_LONG)
     445          43 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     446             :   else
     447             :   {
     448     1303583 :     GEN z=cgeti(d);
     449     1303583 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     450     1303583 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     451     1303583 :     avma=(pari_sp)y;
     452             :   }
     453     1303626 :   return y;
     454             : }
     455             : 
     456             : /* integral part */
     457             : GEN
     458      441499 : floorr(GEN x)
     459             : {
     460             :   long e, d, m, i, lx;
     461             :   GEN y;
     462      441499 :   if (signe(x) >= 0) return truncr(x);
     463      181657 :   if ((e=expo(x)) < 0) return gen_m1;
     464       75332 :   d = nbits2prec(e+1); m = remsBIL(e);
     465       75332 :   lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
     466       75332 :   y = cgeti(d+1);
     467       75332 :   if (++m == BITS_IN_LONG)
     468             :   {
     469          19 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     470          19 :     i=d; while (i<lx && !x[i]) i++;
     471          19 :     if (i==lx) goto END;
     472             :   }
     473             :   else
     474             :   {
     475       75313 :     GEN z=cgeti(d);
     476       75313 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     477       75313 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     478       75313 :     if (uel(x,d-1)<<m == 0)
     479             :     {
     480       19034 :       i=d; while (i<lx && !x[i]) i++;
     481       19034 :       if (i==lx) goto END;
     482             :     }
     483             :   }
     484       68559 :   if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
     485           0 :     y[d++]=1;
     486             : END:
     487       75332 :   y[1] = evalsigne(-1) | evallgefint(d);
     488       75332 :   return y;
     489             : }
     490             : 
     491             : INLINE int
     492  1753333240 : cmpiispec(GEN x, GEN y, long lx, long ly)
     493             : {
     494  1753333240 :   if (lx < ly) return -1;
     495  1673937014 :   if (lx > ly) return  1;
     496  1596862790 :   return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     497             : }
     498             : 
     499             : INLINE int
     500    55330780 : equaliispec(GEN x, GEN y, long lx, long ly)
     501             : {
     502    55330780 :   if (lx != ly) return 0;
     503    55328271 :   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  1950660381 : muluispec(ulong x, GEN y, long ny)
     514             : {
     515  1950660381 :   if (ny == 1)
     516  1505851550 :     return muluu(x, *y);
     517             :   else
     518             :   {
     519   444808831 :     long lz = ny+3;
     520   444808831 :     GEN z = cgeti(lz);
     521   442431648 :     ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
     522   443080754 :     if (hi) { z[lz - 1] = hi; } else lz--;
     523   443080754 :     z[1] = evalsigne(1) | evallgefint(lz);
     524   443080754 :     return z;
     525             :   }
     526             : }
     527             : 
     528             : /* a + b*|y| */
     529             : GEN
     530      759665 : addumului(ulong a, ulong b, GEN y)
     531             : {
     532             :   GEN z;
     533             :   long i, lz;
     534             :   ulong hi;
     535      759665 :   if (!b || !signe(y)) return utoi(a);
     536      759599 :   lz = lgefint(y)+1;
     537      759599 :   z = cgeti(lz);
     538      759599 :   z[2]=a;
     539      759599 :   for(i=3;i<lz;i++) z[i]=0;
     540      759599 :   hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
     541      759599 :   if (hi) z[lz-1]=hi; else lz--;
     542      759599 :   z[1] = evalsigne(1) | evallgefint(lz);
     543      759599 :   avma=(pari_sp)z; return z;
     544             : }
     545             : 
     546             : /***********************************************************************/
     547             : /**                                                                   **/
     548             : /**                          DIVISION                                 **/
     549             : /**                                                                   **/
     550             : /***********************************************************************/
     551             : 
     552             : ulong
     553   774334220 : umodiu(GEN y, ulong x)
     554             : {
     555   774334220 :   long sy=signe(y);
     556             :   ulong hi;
     557   774334220 :   if (!x) pari_err_INV("umodiu",gen_0);
     558   774361042 :   if (!sy) return 0;
     559   574501580 :   hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
     560   574511022 :   if (!hi) return 0;
     561   566377452 :   return (sy > 0)? hi: x - hi;
     562             : }
     563             : 
     564             : /* return |y| \/ x */
     565             : GEN
     566   125634714 : diviu_rem(GEN y, ulong x, ulong *rem)
     567             : {
     568             :   long ly;
     569             :   GEN z;
     570             : 
     571   125634714 :   if (!x) pari_err_INV("diviu_rem",gen_0);
     572   125635019 :   if (!signe(y)) { *rem = 0; return gen_0; }
     573             : 
     574   125096853 :   ly = lgefint(y);
     575   125096853 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
     576             : 
     577   125037363 :   z = cgeti(ly);
     578   125016544 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     579   125039472 :   if (z [ly - 1] == 0) ly--;
     580   125039472 :   z[1] = evallgefint(ly) | evalsigne(1);
     581   125039472 :   return z;
     582             : }
     583             : 
     584             : GEN
     585    29942181 : divis_rem(GEN y, long x, long *rem)
     586             : {
     587    29942181 :   long sy=signe(y),ly,s;
     588             :   GEN z;
     589             : 
     590    29942181 :   if (!x) pari_err_INV("divis_rem",gen_0);
     591    29942181 :   if (!sy) { *rem = 0; return gen_0; }
     592    24218306 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     593             : 
     594    24218306 :   ly = lgefint(y);
     595    24218306 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
     596             : 
     597    19321869 :   z = cgeti(ly);
     598    19321869 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     599    19321869 :   if (sy<0) *rem = -  *rem;
     600    19321869 :   if (z[ly - 1] == 0) ly--;
     601    19321869 :   z[1] = evallgefint(ly) | evalsigne(s);
     602    19321869 :   return z;
     603             : }
     604             : 
     605             : GEN
     606      471143 : divis(GEN y, long x)
     607             : {
     608      471143 :   long sy=signe(y),ly,s;
     609             :   GEN z;
     610             : 
     611      471143 :   if (!x) pari_err_INV("divis",gen_0);
     612      471143 :   if (!sy) return gen_0;
     613      471143 :   if (x<0) { s = -sy; x = -x; } else s=sy;
     614             : 
     615      471143 :   ly = lgefint(y);
     616      471143 :   if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
     617             : 
     618      471143 :   z = cgeti(ly);
     619      471143 :   (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     620      471143 :   if (z[ly - 1] == 0) ly--;
     621      471143 :   z[1] = evallgefint(ly) | evalsigne(s);
     622      471143 :   return z;
     623             : }
     624             : 
     625             : /* We keep llx bits of x and lly bits of y*/
     626             : static GEN
     627    31731208 : divrr_with_gmp(GEN x, GEN y)
     628             : {
     629    31731208 :   long lx=RNLIMBS(x),ly=RNLIMBS(y);
     630    31731208 :   long lw=minss(lx,ly);
     631    31731208 :   long lly=minss(lw+1,ly);
     632    31731208 :   GEN  w=cgetr(lw+2);
     633    31731208 :   long lu=lw+lly;
     634    31731208 :   long llx=minss(lu,lx);
     635    31731208 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     636    31731208 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     637             :   mp_limb_t *q,*r;
     638    31731208 :   long e=expo(x)-expo(y);
     639    31731208 :   long sx=signe(x),sy=signe(y);
     640    31731208 :   xmpn_mirrorcopy(z,RLIMBS(y),lly);
     641    31731208 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     642    31731208 :   xmpn_zero(u,lu-llx);
     643    31731208 :   q = (mp_limb_t *)new_chunk(lw+1);
     644    31731208 :   r = (mp_limb_t *)new_chunk(lly);
     645             : 
     646    31731208 :   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    31731208 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     650    15694023 :     mpn_add_1(q,q,lw+1,1);
     651             : 
     652    31731208 :   xmpn_mirrorcopy(RLIMBS(w),q,lw);
     653             : 
     654    31731208 :   if (q[lw] == 0) e--;
     655    12002402 :   else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
     656          42 :   else { w[2] = HIGHBIT; e++; }
     657    31731208 :   if (sy < 0) sx = -sx;
     658    31731208 :   w[1] = evalsigne(sx) | evalexpo(e);
     659    31731208 :   avma=(pari_sp) w;
     660    31731208 :   return w;
     661             : }
     662             : 
     663             : /* We keep llx bits of x and lly bits of y*/
     664             : static GEN
     665     3593978 : divri_with_gmp(GEN x, GEN y)
     666             : {
     667     3593978 :   long llx=RNLIMBS(x),ly=NLIMBS(y);
     668     3593978 :   long lly=minss(llx+1,ly);
     669     3593978 :   GEN  w=cgetr(llx+2);
     670     3593978 :   long lu=llx+lly, ld=ly-lly;
     671     3593978 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     672     3593978 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     673             :   mp_limb_t *q,*r;
     674     3593978 :   long sh=bfffo(y[ly+1]);
     675     3593978 :   long e=expo(x)-expi(y);
     676     3593978 :   long sx=signe(x),sy=signe(y);
     677     3593978 :   if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
     678       87329 :   else xmpn_copy(z,LIMBS(y)+ld,lly);
     679     3593978 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     680     3593978 :   xmpn_zero(u,lu-llx);
     681     3593978 :   q = (mp_limb_t *)new_chunk(llx+1);
     682     3593978 :   r = (mp_limb_t *)new_chunk(lly);
     683             : 
     684     3593978 :   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     3593978 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     688     1606272 :     mpn_add_1(q,q,llx+1,1);
     689             : 
     690     3593978 :   xmpn_mirrorcopy(RLIMBS(w),q,llx);
     691             : 
     692     3593978 :   if (q[llx] == 0) e--;
     693     2872964 :   else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
     694           0 :   else { w[2] = HIGHBIT; e++; }
     695     3593978 :   if (sy < 0) sx = -sx;
     696     3593978 :   w[1] = evalsigne(sx) | evalexpo(e);
     697     3593978 :   avma=(pari_sp) w;
     698     3593978 :   return w;
     699             : }
     700             : 
     701             : GEN
     702    17750481 : divri(GEN x, GEN y)
     703             : {
     704    17750481 :   long  s = signe(y);
     705             : 
     706    17750481 :   if (!s) pari_err_INV("divri",gen_0);
     707    17750481 :   if (!signe(x)) return real_0_bit(expo(x) - expi(y));
     708    17717042 :   if (!is_bigint(y)) {
     709    14123064 :     GEN z = divru(x, y[2]);
     710    14123064 :     if (s < 0) togglesign(z);
     711    14123064 :     return z;
     712             :   }
     713     3593978 :   return divri_with_gmp(x,y);
     714             : }
     715             : 
     716             : GEN
     717    82456562 : divrr(GEN x, GEN y)
     718             : {
     719    82456562 :   long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
     720             :   ulong y0,y1;
     721             :   GEN r, r1;
     722             : 
     723    82456562 :   if (!sy) pari_err_INV("divrr",y);
     724    82456562 :   e = expo(x) - expo(y);
     725    82456562 :   if (!sx) return real_0_bit(e);
     726    81023208 :   if (sy<0) sx = -sx;
     727             : 
     728    81023208 :   lx=lg(x); ly=lg(y);
     729    81023208 :   if (ly==3)
     730             :   {
     731    49292000 :     ulong k = x[2], l = (lx>3)? x[3]: 0;
     732             :     LOCAL_HIREMAINDER;
     733    49292000 :     if (k < uel(y,2)) e--;
     734             :     else
     735             :     {
     736    22663880 :       l >>= 1; if (k&1) l |= HIGHBIT;
     737    22663880 :       k >>= 1;
     738             :     }
     739    49292000 :     hiremainder = k; k = divll(l,y[2]);
     740    49292000 :     if (hiremainder & HIGHBIT)
     741             :     {
     742    13668211 :       k++;
     743    13668211 :       if (!k) { k = HIGHBIT; e++; }
     744             :     }
     745    49292000 :     r = cgetr(3);
     746    49292000 :     r[1] = evalsigne(sx) | evalexpo(e);
     747    49292000 :     r[2] = k; return r;
     748             :   }
     749             : 
     750    31731208 :   if (ly>=DIVRR_GMP_LIMIT)
     751    31731208 :     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 uel(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 :         uel(r1,1) -= hiremainder;
     811           0 :         while (r1[1])
     812             :         {
     813           0 :           qp++; if (!qp) { j=i; do uel(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 :           uel(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 uel(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           0 :   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   805345783 : dvmdii(GEN x, GEN y, GEN *z)
     849             : {
     850   805345783 :   long sx=signe(x),sy=signe(y);
     851             :   long lx, ly, lq;
     852             :   pari_sp av;
     853             :   GEN r,q;
     854             : 
     855   805345783 :   if (!sy) pari_err_INV("dvmdii",y);
     856   806911821 :   if (!sx)
     857             :   {
     858    15930014 :     if (!z || z == ONLY_REM) return gen_0;
     859     7685322 :     *z=gen_0; return gen_0;
     860             :   }
     861   790981807 :   lx=lgefint(x);
     862   790981807 :   ly=lgefint(y); lq=lx-ly;
     863   790981807 :   if (lq <= 0)
     864             :   {
     865   440598925 :     if (lq == 0)
     866             :     {
     867   423297363 :       long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
     868   423703506 :       if (s>0) goto DIVIDE;
     869   137370438 :       if (s==0)
     870             :       {
     871     6466320 :         if (z == ONLY_REM) return gen_0;
     872     3252282 :         if (z) *z = gen_0;
     873     3252282 :         if (sx < 0) sy = -sy;
     874     3252282 :         return stoi(sy);
     875             :       }
     876             :     }
     877   148205680 :     if (z == ONLY_REM) return icopy(x);
     878     8012381 :     if (z) *z = icopy(x);
     879     8012381 :     return gen_0;
     880             :   }
     881             : DIVIDE: /* quotient is non-zero */
     882   636715950 :   av=avma; if (sx<0) sy = -sy;
     883   636715950 :   if (ly==3)
     884             :   {
     885   335673732 :     ulong lq = lx;
     886             :     ulong si;
     887   335673732 :     q  = cgeti(lq);
     888   335481715 :     si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
     889   335709794 :     if (q[lq - 1] == 0) lq--;
     890   335709794 :     if (z == ONLY_REM)
     891             :     {
     892   294145564 :       avma=av; if (!si) return gen_0;
     893   269707419 :       r=cgeti(3);
     894   269543583 :       r[1] = evalsigne(sx) | evallgefint(3);
     895   269543583 :       r[2] = si; return r;
     896             :     }
     897    41564230 :     q[1] = evalsigne(sy) | evallgefint(lq);
     898    41564230 :     if (!z) return q;
     899    37616385 :     if (!si) { *z=gen_0; return q; }
     900    17082278 :     r=cgeti(3);
     901    17082266 :     r[1] = evalsigne(sx) | evallgefint(3);
     902    17082266 :     r[2] = si; *z=r; return q;
     903             :   }
     904   301042218 :   if (z == ONLY_REM)
     905             :   {
     906   291885822 :     ulong lr = lgefint(y);
     907   291885822 :     ulong lq = lgefint(x)-lgefint(y)+3;
     908   291885822 :     GEN r = cgeti(lr);
     909   289572752 :     GEN q = cgeti(lq);
     910   288595507 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     911   292516708 :     if (!r[lr - 1])
     912             :     {
     913    83978325 :       while(lr>2 && !r[lr - 1]) lr--;
     914    83978325 :       if (lr == 2) {avma=av; return gen_0;} /* exact division */
     915             :     }
     916   288841142 :     r[1] = evalsigne(sx) | evallgefint(lr);
     917   288841142 :     avma = (pari_sp) r; return r;
     918             :   }
     919             :   else
     920             :   {
     921     9156396 :     ulong lq = lgefint(x)-lgefint(y)+3;
     922     9156396 :     ulong lr = lgefint(y);
     923     9156396 :     GEN q = cgeti(lq);
     924     9156396 :     GEN r = cgeti(lr);
     925     9156396 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     926     9156396 :     if (q[lq - 1] == 0) lq--;
     927     9156396 :     q[1] = evalsigne(sy) | evallgefint(lq);
     928     9156396 :     if (!z) { avma = (pari_sp)q; return q; }
     929     8554084 :     if (!r[lr - 1])
     930             :     {
     931     3631583 :       while(lr>2 && !r[lr - 1]) lr--;
     932     3631583 :       if (lr == 2) {avma=(pari_sp) q; *z=gen_0; return q;} /* exact division */
     933             :     }
     934     5895520 :     r[1] = evalsigne(sx) | evallgefint(lr);
     935     5895520 :     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     6711901 : red_montgomery(GEN T, GEN N, ulong inv)
     945             : {
     946             :   pari_sp av;
     947             :   GEN Te, Td, Ne, Nd, scratch;
     948     6711901 :   ulong i, j, m, t, d, k = NLIMBS(N);
     949             :   int carry;
     950             :   LOCAL_HIREMAINDER;
     951             :   LOCAL_OVERFLOW;
     952             : 
     953     6711901 :   if (k == 0) return gen_0;
     954     6711901 :   d = NLIMBS(T); /* <= 2*k */
     955     6711901 :   if (d == 0) return gen_0;
     956             : #ifdef DEBUG
     957             :   if (d > 2*k) pari_err_BUG("red_montgomery");
     958             : #endif
     959     6711892 :   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     6711892 :   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     6720976 :   Td = scratch;
     981     6720976 :   Te = T + 2;
     982     6720976 :   for (i=0; i < d     ; i++) *Td++ = *Te++;
     983     6720976 :   for (   ; i < (k<<1); i++) *Td++ = 0;
     984             : 
     985     6720976 :   Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
     986     6720976 :   Ne = N + 1;       /* 1 beyond end of N mantissa */
     987             : 
     988     6720976 :   carry = 0;
     989    27682600 :   for (i=0; i<k; i++) /* set T := T/B nod N, k times */
     990             :   {
     991    20961624 :     Td = Te; /* one beyond end of (new) T mantissa */
     992    20961624 :     Nd = Ne;
     993    20961624 :     hiremainder = *++Td;
     994    20961624 :     m = hiremainder * inv; /* solve T + m N = O(B) */
     995             : 
     996             :     /* set T := (T + mN) / B */
     997    20961624 :     Te = Td;
     998    20961624 :     (void)addmul(m, *++Nd); /* = 0 */
     999    81611952 :     for (j=1; j<k; j++)
    1000             :     {
    1001    60650328 :       t = addll(addmul(m, *++Nd), *++Td);
    1002    60650328 :       *Td = t;
    1003    60650328 :       hiremainder += overflow;
    1004             :     }
    1005    20961624 :     t = addll(hiremainder, *++Td); *Td = t + carry;
    1006    20961624 :     carry = (overflow || (carry && *Td == 0));
    1007             :   }
    1008     6720976 :   if (carry)
    1009             :   { /* Td > N overflows (k+1 words), set Td := Td - N */
    1010       13411 :     Td = Te;
    1011       13411 :     Nd = Ne;
    1012       13411 :     t = subll(*++Td, *++Nd); *Td = t;
    1013       13411 :     while (Td < (GEN)av) { t = subllx(*++Td, *++Nd); *Td = t; }
    1014             :   }
    1015             : 
    1016             :   /* copy result */
    1017     6720976 :   Td = (GEN)av - 1; /* *Td = high word of final result */
    1018     6720976 :   while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
    1019     6720976 :   k = Td - Te; if (!k) { avma = av; return gen_0; }
    1020     6720976 :   Td = (GEN)av - k; /* will write mantissa there */
    1021     6720976 :   (void)memmove(Td, Te+1, k*sizeof(long));
    1022     6720976 :   Td -= 2;
    1023     6720976 :   Td[0] = evaltyp(t_INT) | evallg(k+2);
    1024     6718746 :   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     6718746 :   avma = (pari_sp)Td; return Td;
    1036             : }
    1037             : 
    1038             : /* EXACT INTEGER DIVISION */
    1039             : 
    1040             : #if 1 && !defined(_WIN64) /* use undocumented GMP interface */
    1041             :                           /* mpz_divexact_ui is not LLP64 friendly */
    1042             : static void
    1043   223529041 : GEN2mpz(mpz_t X, GEN x)
    1044             : {
    1045   223529041 :   long l = lgefint(x)-2;
    1046   223529041 :   X->_mp_alloc = l;
    1047   223529041 :   X->_mp_size = signe(x) > 0? l: -l;
    1048   223529041 :   X->_mp_d = LIMBS(x);
    1049   223529041 : }
    1050             : static void
    1051   212761387 : mpz2GEN(GEN z, mpz_t Z)
    1052             : {
    1053   212761387 :   long l = Z->_mp_size;
    1054   212761387 :   z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
    1055   212761387 : }
    1056             : 
    1057             : /* Find z such that x=y*z, knowing that y | x (unchecked) */
    1058             : GEN
    1059   231420318 : diviiexact(GEN x, GEN y)
    1060             : {
    1061   231420318 :   if (!signe(y)) pari_err_INV("diviiexact",y);
    1062   231420615 :   if (lgefint(y) == 3)
    1063             :   {
    1064   217787706 :     GEN z = diviuexact(x, y[2]);
    1065   217786937 :     if (signe(y) < 0) togglesign(z);
    1066   217787585 :     return z;
    1067             :   }
    1068    13632909 :   if (!signe(x)) return gen_0;
    1069             :   {
    1070    10776397 :     long l = lgefint(x);
    1071             :     mpz_t X, Y, Z;
    1072    10776397 :     GEN z = cgeti(l);
    1073    10776397 :     GEN2mpz(X, x);
    1074    10776397 :     GEN2mpz(Y, y);
    1075    10776397 :     Z->_mp_alloc = l-2;
    1076    10776397 :     Z->_mp_size  = l-2;
    1077    10776397 :     Z->_mp_d = LIMBS(z);
    1078    10776397 :     mpz_divexact(Z, X, Y);
    1079    10776397 :     mpz2GEN(z, Z); return z;
    1080             :   }
    1081             : }
    1082             : #else
    1083             : /* Find z such that x=y*z, knowing that y | x (unchecked)
    1084             :  * Method: y0 z0 = x0 mod B = 2^BITS_IN_LONG ==> z0 = 1/y0 mod B.
    1085             :  *    Set x := (x - z0 y) / B, updating only relevant words, and repeat */
    1086             : GEN
    1087             : diviiexact(GEN x, GEN y)
    1088             : {
    1089             :   /*TODO: use mpn_bdivmod instead*/
    1090             :   if (!signe(y)) pari_err_INV("diviiexact",y);
    1091             :   return divii(x,y);
    1092             : }
    1093             : #endif
    1094             : 
    1095             : #ifdef mpn_divexact_1
    1096             : GEN
    1097             : diviuexact(GEN y, ulong x)
    1098             : {
    1099             :   long ly;
    1100             :   GEN z;
    1101             :   if (!signe(y)) return gen_0;
    1102             :   ly = lgefint(y);
    1103             :   z = cgeti(ly);
    1104             :   mpn_divexact_1(LIMBS(z), LIMBS(y), NLIMBS(y), x);
    1105             :   if (z [ly - 1] == 0) ly--;
    1106             :   z[1] = evallgefint(ly) | evalsigne(signe(y));
    1107             :   return z;
    1108             : }
    1109             : #elif 1 && !defined(_WIN64) /* mpz_divexact_ui is not LLP64 friendly   */
    1110             :                             /* assume y != 0 and the division is exact */
    1111             : GEN
    1112   244193737 : diviuexact(GEN x, ulong y)
    1113             : {
    1114   244193737 :   if (!signe(x)) return gen_0;
    1115             :   {
    1116   201985173 :     long l = lgefint(x);
    1117             :     mpz_t X, Z;
    1118   201985173 :     GEN z = cgeti(l);
    1119   201971875 :     GEN2mpz(X, x);
    1120   201972703 :     Z->_mp_alloc = l-2;
    1121   201972703 :     Z->_mp_size  = l-2;
    1122   201972703 :     Z->_mp_d = LIMBS(z);
    1123   201972703 :     mpz_divexact_ui(Z, X, y);
    1124   201984475 :     mpz2GEN(z, Z); return z;
    1125             :   }
    1126             : }
    1127             : #else
    1128             : /* assume y != 0 and the division is exact */
    1129             : GEN
    1130             : diviuexact(GEN x, ulong y)
    1131             : {
    1132             :   /*TODO: implement true exact division.*/
    1133             :   return divii(x,utoi(y));
    1134             : }
    1135             : #endif
    1136             : 
    1137             : /* assume yz != and yz | x */
    1138             : GEN
    1139      513237 : diviuuexact(GEN x, ulong y, ulong z)
    1140             : {
    1141             :   long tmp[4];
    1142             :   ulong t;
    1143             :   LOCAL_HIREMAINDER;
    1144      513237 :   t = mulll(y, z);
    1145      513237 :   if (!hiremainder) return diviuexact(x, t);
    1146           0 :   tmp[0] = evaltyp(t_INT)|_evallg(4);
    1147           0 :   tmp[1] = evalsigne(1)|evallgefint(4);
    1148           0 :   tmp[2] = t;
    1149           0 :   tmp[3] = hiremainder;
    1150           0 :   return diviiexact(x, tmp);
    1151             : }
    1152             : 
    1153             : 
    1154             : /********************************************************************/
    1155             : /**                                                                **/
    1156             : /**               INTEGER MULTIPLICATION                           **/
    1157             : /**                                                                **/
    1158             : /********************************************************************/
    1159             : 
    1160             : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
    1161             : GEN
    1162  1832496249 : muliispec(GEN x, GEN y, long nx, long ny)
    1163             : {
    1164             :   GEN zd;
    1165             :   long lz;
    1166             :   ulong hi;
    1167             : 
    1168  1832496249 :   if (nx < ny) swapspec(x,y, nx,ny);
    1169  1832496249 :   if (!ny) return gen_0;
    1170  1832496249 :   if (ny == 1) return muluispec((ulong)*y, x, nx);
    1171             : 
    1172   474340381 :   lz = nx+ny+2;
    1173   474340381 :   zd = cgeti(lz);
    1174   473183932 :   hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
    1175   475466260 :   if (!hi) lz--;
    1176             :   /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
    1177             : 
    1178   475466260 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1179   475466260 :   return zd;
    1180             : }
    1181             : GEN
    1182      535445 : muluui(ulong x, ulong y, GEN z)
    1183             : {
    1184      535445 :   long t, s = signe(z);
    1185             :   GEN r;
    1186             :   LOCAL_HIREMAINDER;
    1187             : 
    1188      535445 :   if (!x || !y || !signe(z)) return gen_0;
    1189      535049 :   t = mulll(x,y);
    1190      535049 :   if (!hiremainder)
    1191      535049 :     r = muluispec(t, z+2, lgefint(z)-2);
    1192             :   else
    1193             :   {
    1194             :     long tmp[2];
    1195           0 :     tmp[1] = hiremainder;
    1196           0 :     tmp[0] = t;
    1197           0 :     r = muliispec(z+2,tmp, lgefint(z)-2, 2);
    1198             :   }
    1199      535049 :   setsigne(r,s); return r;
    1200             : }
    1201             : 
    1202             : GEN
    1203   640219712 : sqrispec(GEN x, long nx)
    1204             : {
    1205             :   GEN zd;
    1206             :   long lz;
    1207             : 
    1208   640219712 :   if (!nx) return gen_0;
    1209   158142123 :   if (nx==1) return sqru(*x);
    1210             : 
    1211    96553190 :   lz = (nx<<1)+2;
    1212    96553190 :   zd = cgeti(lz);
    1213             : #ifdef mpn_sqr
    1214    96200792 :   mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
    1215             : #else
    1216             :   mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
    1217             : #endif
    1218    97056270 :   if (zd[lz-1]==0) lz--;
    1219             : 
    1220    97056270 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1221    97056270 :   return zd;
    1222             : }
    1223             : 
    1224             : INLINE GEN
    1225     4837347 : sqrispec_mirror(GEN x, long nx)
    1226             : {
    1227     4837347 :   GEN cx=new_chunk(nx);
    1228             :   GEN z;
    1229     4837347 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1230     4837347 :   z=sqrispec(cx, nx);
    1231     4837347 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1232     4837347 :   return z;
    1233             : }
    1234             : 
    1235             : /* leaves garbage on the stack. */
    1236             : INLINE GEN
    1237    62794681 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
    1238             : {
    1239             :   GEN cx, cy, z;
    1240    62794681 :   long s = 0;
    1241    62794681 :   while (nx && x[nx-1]==0) { nx--; s++; }
    1242    62794681 :   while (ny && y[ny-1]==0) { ny--; s++; }
    1243    62794681 :   cx=new_chunk(nx); cy=new_chunk(ny);
    1244    62794681 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1245    62794681 :   xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
    1246    62794681 :   z =  nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
    1247    62794681 :   if (s)
    1248             :   {
    1249     7810561 :     long i, lz = lgefint(z) + s;
    1250     7810561 :     (void)new_chunk(s);
    1251     7810561 :     z -= s;
    1252     7810561 :     for (i=0; i<s; i++) z[2+i]=0;
    1253     7810561 :     z[1] = evalsigne(1) | evallgefint(lz);
    1254     7810561 :     z[0] = evaltyp(t_INT) | evallg(lz);
    1255             :   }
    1256    62794681 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1257    62794681 :   return z;
    1258             : }
    1259             : 
    1260             : /* x % (2^n), assuming n >= 0 */
    1261             : GEN
    1262    10056411 : remi2n(GEN x, long n)
    1263             : {
    1264             :   ulong hi;
    1265    10056411 :   long l, k, lx, ly, sx = signe(x);
    1266             :   GEN z, xd, zd;
    1267             : 
    1268    10056411 :   if (!sx || !n) return gen_0;
    1269             : 
    1270     9908970 :   k = dvmdsBIL(n, &l);
    1271     9931891 :   lx = lgefint(x);
    1272     9931891 :   if (lx < k+3) return icopy(x);
    1273             : 
    1274     9841741 :   xd = x + (2 + k);
    1275             :   /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
    1276             :    *              ^--- initial xd  */
    1277     9841741 :   hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
    1278     9841741 :   if (!hi)
    1279             :   { /* strip leading zeroes from result */
    1280     1089234 :     xd--; while (k && !*xd) { k--; xd--; }
    1281     1089234 :     if (!k) return gen_0;
    1282      483959 :     ly = k+2;
    1283             :   }
    1284             :   else
    1285     8752507 :     ly = k+3;
    1286             : 
    1287     9236466 :   zd = z = cgeti(ly);
    1288     9116234 :   *++zd = evalsigne(sx) | evallgefint(ly);
    1289     9116234 :   xd = x+1; for ( ;k; k--) *++zd = *++xd;
    1290     9116234 :   if (hi) *++zd = hi;
    1291     9116234 :   return z;
    1292             : }
    1293             : 
    1294             : /********************************************************************/
    1295             : /**                                                                **/
    1296             : /**                      INTEGER SQUARE ROOT                       **/
    1297             : /**                                                                **/
    1298             : /********************************************************************/
    1299             : 
    1300             : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
    1301             :  * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
    1302             :  * remainder is 0. R = NULL is allowed. */
    1303             : GEN
    1304     2658920 : sqrtremi(GEN a, GEN *r)
    1305             : {
    1306     2658920 :   long l, na = NLIMBS(a);
    1307             :   mp_size_t nr;
    1308             :   GEN S;
    1309     2658920 :   if (!na) {
    1310           0 :     if (r) *r = gen_0;
    1311           0 :     return gen_0;
    1312             :   }
    1313     2658920 :   l = (na + 5) >> 1; /* 2 + ceil(na/2) */
    1314     2658920 :   S = cgetipos(l);
    1315     2658920 :   if (r) {
    1316        3775 :     GEN R = cgeti(2 + na);
    1317        3775 :     nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
    1318        3775 :     if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
    1319        3394 :     else    { avma = (pari_sp)S; R = gen_0; }
    1320        3775 :     *r = R;
    1321             :   }
    1322             :   else
    1323     2655145 :     (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
    1324     2658920 :   return S;
    1325             : }
    1326             : 
    1327             : /* compute sqrt(|a|), assuming a != 0 */
    1328             : GEN
    1329    21712691 : sqrtr_abs(GEN a)
    1330             : {
    1331             :   GEN res;
    1332             :   mp_limb_t *b, *c;
    1333    21712691 :   long l = RNLIMBS(a), e = expo(a), er = e>>1;
    1334             :   long n;
    1335    21712691 :   res = cgetr(2 + l);
    1336    21712691 :   res[1] = evalsigne(1) | evalexpo(er);
    1337    21712691 :   if (e&1)
    1338             :   {
    1339    11136687 :     b = (mp_limb_t *) new_chunk(l<<1);
    1340    11136687 :     xmpn_zero(b,l);
    1341    11136687 :     xmpn_mirrorcopy(b+l, RLIMBS(a), l);
    1342    11136687 :     c = (mp_limb_t *) new_chunk(l);
    1343    11136687 :     n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
    1344    11136687 :     if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
    1345             :   }
    1346             :   else
    1347             :   {
    1348             :     ulong u;
    1349    10576004 :     b = (mp_limb_t *) mantissa2nr(a,-1);
    1350    10576004 :     b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
    1351    10576004 :     b = (mp_limb_t *) new_chunk(l);
    1352    10576004 :     xmpn_zero(b,l+1); /* overwrites the former b[0] */
    1353    10576004 :     c = (mp_limb_t *) new_chunk(l + 1);
    1354    10576004 :     n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
    1355    10576004 :     u = (ulong)*c++;
    1356    10576004 :     if ( u&HIGHBIT || (u == ~HIGHBIT &&
    1357           0 :              (n>l || (n==l && mpn_cmp(b,c,l)>0))))
    1358     5144668 :       mpn_add_1(c,c,l,1);
    1359             :   }
    1360    21712691 :   xmpn_mirrorcopy(RLIMBS(res),c,l);
    1361    21712691 :   avma = (pari_sp)res; return res;
    1362             : }
    1363             : 
    1364             : /* Normalize a non-negative integer */
    1365             : GEN
    1366   145157348 : int_normalize(GEN x, long known_zero_words)
    1367             : {
    1368   145157348 :   long i =  lgefint(x) - 1 - known_zero_words;
    1369  1199093324 :   for ( ; i > 1; i--)
    1370  1138533473 :     if (x[i]) { setlgefint(x, i+1); return x; }
    1371    60559851 :   x[1] = evalsigne(0) | evallgefint(2); return x;
    1372             : }
    1373             : 
    1374             : /********************************************************************
    1375             :  **                                                                **
    1376             :  **                           Base Conversion                      **
    1377             :  **                                                                **
    1378             :  ********************************************************************/
    1379             : 
    1380             : ulong *
    1381      486112 : convi(GEN x, long *l)
    1382             : {
    1383      486112 :   long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
    1384      486112 :   GEN str = cgetg(n+1, t_VECSMALL);
    1385      486112 :   unsigned char *res = (unsigned char*) GSTR(str);
    1386      486112 :   long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
    1387             :   long lz;
    1388             :   ulong *z;
    1389             :   long i, j;
    1390             :   unsigned char *t;
    1391      486112 :   while (!*res) {res++; llz--;} /*Strip leading zeros*/
    1392      486112 :   lz  = (8+llz)/9;
    1393      486112 :   z = (ulong*)new_chunk(1+lz);
    1394      486112 :   t=res+llz+9;
    1395     1214538 :   for(i=0;i<llz-8;i+=9)
    1396             :   {
    1397             :     ulong s;
    1398      728426 :     t-=18;
    1399      728426 :     s=*t++;
    1400     6555834 :     for (j=1; j<9;j++)
    1401     5827408 :       s=10*s+*t++;
    1402      728426 :     *z++=s;
    1403             :   }
    1404      486112 :   if (i<llz)
    1405             :   {
    1406      470983 :     unsigned char *t = res;
    1407      470983 :     ulong s=*t++;
    1408     1983835 :     for (j=i+1; j<llz;j++)
    1409     1512852 :       s=10*s+*t++;
    1410      470983 :     *z++=s;
    1411             :   }
    1412      486112 :   *l = lz;
    1413      486112 :   return z;
    1414             : }

Generated by: LCOV version 1.11