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.16.2 lcov report (development 29419-8afb0ed749) Lines: 683 719 95.0 %
Date: 2024-07-02 09:03:41 Functions: 54 56 96.4 %
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; either version 2 of the License, or (at your option) any later
       9             : version. It is distributed in the hope that it will be useful, but WITHOUT
      10             : ANY WARRANTY WHATSOEVER.
      11             : 
      12             : Check the License for details. You should have received a copy of it, along
      13             : with the package; see the file 'COPYING'. If not, write to the Free Software
      14             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      15             : 
      16             : /***********************************************************************/
      17             : /**                                                                   **/
      18             : /**                               GMP KERNEL                          **/
      19             : /** BA2002Sep24                                                       **/
      20             : /***********************************************************************/
      21             : /* GMP t_INT as just like normal t_INT, just the mantissa is the other way
      22             :  * round
      23             :  *
      24             :  *   `How would you like to live in Looking-glass House, Kitty?  I
      25             :  *   wonder if they'd give you milk in there?  Perhaps Looking-glass
      26             :  *   milk isn't good to drink--But oh, Kitty! now we come to the
      27             :  *   passage.  You can just see a little PEEP of the passage in
      28             :  *   Looking-glass House, if you leave the door of our drawing-room
      29             :  *   wide open:  and it's very like our passage as far as you can see,
      30             :  *   only you know it may be quite different on beyond.  Oh, Kitty!
      31             :  *   how nice it would be if we could only get through into Looking-
      32             :  *   glass House!  I'm sure it's got, oh! such beautiful things in it!
      33             :  *
      34             :  *  Through the Looking Glass,  Lewis Carrol
      35             :  *
      36             :  *  (pityful attempt to beat GN code/comments rate)
      37             :  *  */
      38             : 
      39             : #include <gmp.h>
      40             : #include "pari.h"
      41             : #include "paripriv.h"
      42             : #include "../src/kernel/none/tune-gen.h"
      43             : 
      44             : /*We need PARI invmod renamed to invmod_pari*/
      45             : #define INVMOD_PARI
      46             : 
      47           0 : static void *pari_gmp_realloc(void *ptr, size_t old_size, size_t new_size) {
      48           0 :   (void)old_size; return (void *) pari_realloc(ptr,new_size);
      49             : }
      50             : 
      51     1733455 : static void pari_gmp_free(void *ptr, size_t old_size){
      52     1733455 :   (void)old_size; pari_free(ptr);
      53     1733455 : }
      54             : 
      55             : static void *(*old_gmp_malloc)(size_t new_size);
      56             : static void *(*old_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
      57             : static void (*old_gmp_free)(void *ptr, size_t old_size);
      58             : 
      59             : void
      60        1084 : pari_kernel_init(void)
      61             : {
      62        1084 :   mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
      63        1084 :   mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
      64        1084 : }
      65             : 
      66             : const char *
      67           4 : pari_kernel_version(void)
      68             : {
      69             : #ifdef gmp_version
      70           4 :   return gmp_version;
      71             : #else
      72             :   return "";
      73             : #endif
      74             : }
      75             : 
      76             : void
      77        1076 : pari_kernel_close(void)
      78             : {
      79             :   void *(*new_gmp_malloc)(size_t new_size);
      80             :   void *(*new_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
      81             :   void (*new_gmp_free)(void *ptr, size_t old_size);
      82        1076 :   mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
      83        1076 :   if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
      84        1076 :   if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
      85        1076 :   if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
      86        1076 :   mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
      87        1076 : }
      88             : 
      89             : #define LIMBS(x)  ((mp_limb_t *)((x)+2))
      90             : #define NLIMBS(x) (lgefint(x)-2)
      91             : /*This one is for t_REALs to emphasize they are not t_INTs*/
      92             : #define RLIMBS(x)  ((mp_limb_t *)((x)+2))
      93             : #define RNLIMBS(x) (lg(x)-2)
      94             : 
      95             : INLINE void
      96     6679462 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
      97             : {
      98    55891155 :   while (--n >= 0) x[n]=y[n];
      99     6679462 : }
     100             : 
     101             : INLINE void
     102   579713879 : xmpn_mirror(mp_limb_t *x, long n)
     103             : {
     104             :   long i;
     105  3853032311 :   for(i=0;i<(n>>1);i++)
     106             :   {
     107  3273318432 :     ulong m=x[i];
     108  3273318432 :     x[i]=x[n-1-i];
     109  3273318432 :     x[n-1-i]=m;
     110             :   }
     111   579713879 : }
     112             : 
     113             : INLINE void
     114   709155567 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
     115             : {
     116             :   long i;
     117  9764369707 :   for(i=0;i<n;i++)
     118  9055214140 :     z[i]=x[n-1-i];
     119   709155567 : }
     120             : 
     121             : INLINE void
     122   234422454 : xmpn_zero(mp_ptr x, mp_size_t n)
     123             : {
     124  2048764776 :   while (--n >= 0) x[n]=0;
     125   234422454 : }
     126             : 
     127             : INLINE GEN
     128    40502628 : icopy_ef(GEN x, long l)
     129             : {
     130    40502628 :   long lx = lgefint(x);
     131    40502628 :   const GEN y = cgeti(l);
     132             : 
     133   287973933 :   while (--lx > 0) y[lx]=x[lx];
     134    40501283 :   return y;
     135             : }
     136             : 
     137             : /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
     138             :  * GENs but pairs (long *a, long na) representing a list of digits (in basis
     139             :  * BITS_IN_LONG) : a[0], ..., a[na-1]. [ In ordre to facilitate splitting: no
     140             :  * need to reintroduce codewords ]
     141             :  * Use speci(a,na) to visualize the corresponding GEN.
     142             :  */
     143             : 
     144             : /***********************************************************************/
     145             : /**                                                                   **/
     146             : /**                     ADDITION / SUBTRACTION                        **/
     147             : /**                                                                   **/
     148             : /***********************************************************************/
     149             : 
     150             : GEN
     151     2996963 : setloop(GEN a)
     152             : {
     153     2996963 :   pari_sp av = avma - 2 * sizeof(long);
     154     2996963 :   (void)cgetg(lgefint(a) + 3, t_VECSMALL);
     155     2996963 :   return icopy_avma(a, av); /* two cells of extra space after a */
     156             : }
     157             : 
     158             : /* we had a = setloop(?), then some incloops. Reset a to b */
     159             : GEN
     160      174328 : resetloop(GEN a, GEN b) {
     161      174328 :   a[0] = evaltyp(t_INT) | evallg(lgefint(b));
     162      174328 :   affii(b, a); return a;
     163             : }
     164             : 
     165             : /* assume a > 0, initialized by setloop. Do a++ */
     166             : static GEN
     167    99170150 : incpos(GEN a)
     168             : {
     169    99170150 :   long i, l = lgefint(a);
     170    99170155 :   for (i=2; i<l; i++)
     171    99177112 :     if (++uel(a,i)) return a;
     172           3 :   a[l] = 1; l++;
     173           3 :   a[0]=evaltyp(t_INT) | _evallg(l);
     174           3 :   a[1]=evalsigne(1) | evallgefint(l);
     175           3 :   return a;
     176             : }
     177             : 
     178             : /* assume a < 0, initialized by setloop. Do a++ */
     179             : static GEN
     180       66652 : incneg(GEN a)
     181             : {
     182       66652 :   long i, l = lgefint(a)-1;
     183       66652 :   if (uel(a,2)--)
     184             :   {
     185       66648 :     if (!a[l]) /* implies l = 2 */
     186             :     {
     187        1976 :       a[0] = evaltyp(t_INT) | _evallg(2);
     188        1976 :       a[1] = evalsigne(0) | evallgefint(2);
     189             :     }
     190       66648 :     return a;
     191             :   }
     192           5 :   for (i=3; i<=l; i++)
     193           5 :     if (uel(a,i)--) break;
     194           4 :   if (!a[l])
     195             :   {
     196           4 :     a[0] = evaltyp(t_INT) | _evallg(l);
     197           4 :     a[1] = evalsigne(-1) | evallgefint(l);
     198             :   }
     199           4 :   return a;
     200             : }
     201             : 
     202             : /* assume a initialized by setloop. Do a++ */
     203             : GEN
     204    99570324 : incloop(GEN a)
     205             : {
     206    99570324 :   switch(signe(a))
     207             :   {
     208      336504 :     case 0:
     209      336504 :       a[0]=evaltyp(t_INT) | _evallg(3);
     210      336504 :       a[1]=evalsigne(1) | evallgefint(3);
     211      336504 :       a[2]=1; return a;
     212       66652 :     case -1: return incneg(a);
     213    99167168 :     default: return incpos(a);
     214             :   }
     215             : }
     216             : 
     217             : INLINE GEN
     218  2511193964 : adduispec(ulong s, GEN x, long nx)
     219             : {
     220             :   GEN  zd;
     221             :   long lz;
     222             : 
     223  2511193964 :   if (nx == 1) return adduu(uel(x,0), s);
     224   710274275 :   lz = nx+3; zd = cgeti(lz);
     225   713189844 :   if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
     226     5008474 :     zd[lz-1]=1;
     227             :   else
     228   708169189 :     lz--;
     229   713177663 :   zd[1] = evalsigne(1) | evallgefint(lz);
     230   713177663 :   return zd;
     231             : }
     232             : 
     233             : GEN
     234   576036254 : adduispec_offset(ulong s, GEN x, long offset, long nx)
     235             : {
     236   576036254 :   GEN xd=x+2+offset;
     237   799488946 :   while (nx && *(xd+nx-1)==0) nx--;
     238   576036254 :   if (!nx) return utoi(s);
     239   531288886 :   return adduispec(s,xd,nx);
     240             : }
     241             : 
     242             : INLINE GEN
     243  3157498151 : addiispec(GEN x, GEN y, long nx, long ny)
     244             : {
     245             :   GEN zd;
     246             :   long lz;
     247             : 
     248  3157498151 :   if (nx < ny) swapspec(x,y, nx,ny);
     249  3157498151 :   if (ny == 1) return adduispec(*y,x,nx);
     250  1265854608 :   lz = nx+3; zd = cgeti(lz);
     251             : 
     252  1266032402 :   if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
     253    27609842 :     zd[lz-1]=1;
     254             :   else
     255  1238728833 :     lz--;
     256             : 
     257  1266338675 :   zd[1] = evalsigne(1) | evallgefint(lz);
     258  1266338675 :   return zd;
     259             : }
     260             : 
     261             : /* assume x >= y */
     262             : INLINE GEN
     263  1659992489 : subiuspec(GEN x, ulong s, long nx)
     264             : {
     265             :   GEN zd;
     266             :   long lz;
     267             : 
     268  1659992489 :   if (nx == 1) return utoi(x[0] - s);
     269             : 
     270   237407766 :   lz = nx + 2; zd = cgeti(lz);
     271   237868687 :   mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
     272   237868947 :   if (! zd[lz - 1]) { --lz; }
     273             : 
     274   237868947 :   zd[1] = evalsigne(1) | evallgefint(lz);
     275   237868947 :   return zd;
     276             : }
     277             : 
     278             : /* assume x > y */
     279             : INLINE GEN
     280  2861985024 : subiispec(GEN x, GEN y, long nx, long ny)
     281             : {
     282             :   GEN zd;
     283             :   long lz;
     284  2861985024 :   if (ny==1) return subiuspec(x,*y,nx);
     285  1231801979 :   lz = nx+2; zd = cgeti(lz);
     286             : 
     287  1229846807 :   mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
     288  1688082899 :   while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
     289             : 
     290  1229889651 :   zd[1] = evalsigne(1) | evallgefint(lz);
     291  1229889651 :   return zd;
     292             : }
     293             : 
     294             : static void
     295   520712104 : roundr_up_ip(GEN x, long l)
     296             : {
     297   520712104 :   long i = l;
     298             :   for(;;)
     299             :   {
     300   522527064 :     if (++((ulong*)x)[--i]) break;
     301     2160257 :     if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
     302             :   }
     303   520801745 : }
     304             : 
     305             : void
     306   396799942 : affir(GEN x, GEN y)
     307             : {
     308   396799942 :   const long s = signe(x), ly = lg(y);
     309             :   long lx, sh, i;
     310             : 
     311   396799942 :   if (!s)
     312             :   {
     313    40870262 :     y[1] = evalexpo(-bit_accuracy(ly));
     314    40867182 :     return;
     315             :   }
     316   355929680 :   lx = lgefint(x); sh = bfffo(*int_MSW(x));
     317   355929680 :   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
     318   356124577 :   if (sh) {
     319   348660882 :     if (lx <= ly)
     320             :     {
     321   944700884 :       for (i=lx; i<ly; i++) y[i]=0;
     322   245112391 :       mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
     323   245197944 :       xmpn_mirror(LIMBS(y),lx-2);
     324   245531160 :       return;
     325             :     }
     326   103548491 :     mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
     327   103548670 :     uel(y,2) |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
     328   103548670 :     xmpn_mirror(LIMBS(y),ly-2);
     329             :     /* lx > ly: round properly */
     330   103550182 :     if ((uel(x,lx-ly+1)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
     331             :   }
     332             :   else {
     333     7463695 :     GEN xd=int_MSW(x);
     334     7463695 :     if (lx <= ly)
     335             :     {
     336     9328965 :       for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
     337     5112542 :       for (   ; i<ly; i++) y[i]=0;
     338     2334723 :       return;
     339             :     }
     340    14596084 :     for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
     341             :     /* lx > ly: round properly */
     342     5128972 :     if (uel(x,lx-ly+1) & HIGHBIT) roundr_up_ip(y, ly);
     343             :   }
     344             : }
     345             : 
     346             : INLINE GEN
     347   690297741 : shiftispec(GEN x, long nx, long n)
     348             : {
     349             :   long ny,m;
     350             :   GEN yd, y;
     351             : 
     352   690297741 :   if (!n) return icopyspec(x, nx);
     353             : 
     354   662436206 :   if (n > 0)
     355             :   {
     356   389486895 :     long d = dvmdsBIL(n, &m);
     357             :     long i;
     358             : 
     359   389421639 :     ny = nx + d + (m!=0);
     360   389421639 :     y = cgeti(ny + 2); yd = y + 2;
     361   541683957 :     for (i=0; i<d; i++) yd[i] = 0;
     362             : 
     363   388517829 :     if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
     364             :     else
     365             :     {
     366   386991184 :       ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
     367   387072951 :       if (carryd) yd[ny - 1] = carryd;
     368   364041197 :       else ny--;
     369             :     }
     370             :   }
     371             :   else
     372             :   {
     373   272949311 :     long d = dvmdsBIL(-n, &m);
     374             : 
     375   274885644 :     ny = nx - d;
     376   274885644 :     if (ny < 1) return gen_0;
     377   272103143 :     y = cgeti(ny + 2); yd = y + 2;
     378             : 
     379   271776405 :     if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
     380             :     else
     381             :     {
     382   267305639 :       mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
     383   267317428 :       if (yd[ny - 1] == 0)
     384             :       {
     385    59233848 :         if (ny == 1) return gc_const((pari_sp)(yd + 1), gen_0);
     386    49915097 :         ny--;
     387             :       }
     388             :     }
     389             :   }
     390   649311797 :   y[1] = evalsigne(1)|evallgefint(ny + 2);
     391   649311797 :   return y;
     392             : }
     393             : 
     394             : GEN
     395   137152218 : mantissa2nr(GEN x, long n)
     396             : {
     397   137152218 :   long ly, i, m, s = signe(x), lx = lg(x);
     398             :   GEN y;
     399   137152218 :   if (!s) return gen_0;
     400   137150869 :   if (!n)
     401             :   {
     402    30497932 :     y = cgeti(lx);
     403    30492017 :     y[1] = evalsigne(s) | evallgefint(lx);
     404    30492017 :     xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
     405    30491471 :     return y;
     406             :   }
     407   106652937 :   if (n > 0)
     408             :   {
     409      216312 :     GEN z = (GEN)avma;
     410      216312 :     long d = dvmdsBIL(n, &m);
     411             : 
     412      216312 :     ly = lx+d; y = new_chunk(ly);
     413      549392 :     for ( ; d; d--) *--z = 0;
     414      220729 :     if (!m) for (i=2; i<lx; i++) y[i]=x[i];
     415             :     else
     416             :     {
     417      214877 :       const ulong sh = BITS_IN_LONG - m;
     418      214877 :       shift_left(y,x, 2,lx-1, 0,m);
     419      214877 :       i = uel(x,2) >> sh;
     420             :       /* Extend y on the left? */
     421      214877 :       if (i) { ly++; y = new_chunk(1); y[2] = i; }
     422             :     }
     423             :   }
     424             :   else
     425             :   {
     426   106436625 :     ly = lx - dvmdsBIL(-n, &m);
     427   106441716 :     if (ly<3) return gen_0;
     428   106441716 :     y = new_chunk(ly);
     429   106406303 :     if (m) {
     430   106149218 :       shift_right(y,x, 2,ly, 0,m);
     431   106204988 :       if (y[2] == 0)
     432             :       {
     433           0 :         if (ly==3) return gc_const((pari_sp)(y+3), gen_0);
     434           0 :         ly--; set_avma((pari_sp)(++y));
     435             :       }
     436             :     } else {
     437      697924 :       for (i=2; i<ly; i++) y[i]=x[i];
     438             :     }
     439             :   }
     440   106678385 :   xmpn_mirror(LIMBS(y),ly-2);
     441   106719145 :   y[1] = evalsigne(s)|evallgefint(ly);
     442   106719145 :   y[0] = evaltyp(t_INT)|evallg(ly); return y;
     443             : }
     444             : 
     445             : GEN
     446     3473110 : truncr(GEN x)
     447             : {
     448             :   long s, e, d, m, i;
     449             :   GEN y;
     450     3473110 :   if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
     451     1445930 :   d = nbits2lg(e+1); m = remsBIL(e);
     452     1445918 :   if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
     453             : 
     454     1445914 :   y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
     455     1445852 :   if (++m == BITS_IN_LONG)
     456       98937 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     457             :   else
     458             :   {
     459     1396637 :     GEN z=cgeti(d);
     460     2862065 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     461     1396570 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     462     1396574 :     set_avma((pari_sp)y);
     463             :   }
     464     1445778 :   return y;
     465             : }
     466             : 
     467             : /* integral part */
     468             : GEN
     469     6945333 : floorr(GEN x)
     470             : {
     471             :   long e, d, m, i, lx;
     472             :   GEN y;
     473     6945333 :   if (signe(x) >= 0) return truncr(x);
     474     4179475 :   if ((e=expo(x)) < 0) return gen_m1;
     475     3501087 :   d = nbits2lg(e+1); m = remsBIL(e);
     476     3498380 :   lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
     477     3498380 :   y = cgeti(d+1);
     478     3490133 :   if (++m == BITS_IN_LONG)
     479             :   {
     480        2942 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     481        1360 :     i=d; while (i<lx && !x[i]) i++;
     482         876 :     if (i==lx) goto END;
     483             :   }
     484             :   else
     485             :   {
     486     3489257 :     GEN z=cgeti(d);
     487     7807301 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     488     3475693 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     489     3476650 :     if (uel(x,d-1)<<m == 0)
     490             :     {
     491      494236 :       i=d; while (i<lx && !x[i]) i++;
     492      111328 :       if (i==lx) goto END;
     493             :     }
     494             :   }
     495     3405650 :   if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
     496           0 :     y[d++]=1;
     497     3405769 : END:
     498     3477645 :   y[1] = evalsigne(-1) | evallgefint(d);
     499     3477645 :   return y;
     500             : }
     501             : 
     502             : INLINE int
     503  3718405836 : cmpiispec(GEN x, GEN y, long lx, long ly)
     504             : {
     505  3718405836 :   if (lx < ly) return -1;
     506  3433449979 :   if (lx > ly) return  1;
     507  3288049376 :   return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     508             : }
     509             : 
     510             : INLINE int
     511   269941158 : equaliispec(GEN x, GEN y, long lx, long ly)
     512             : {
     513   269941158 :   if (lx != ly) return 0;
     514   269807210 :   return !mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     515             : }
     516             : 
     517             : /***********************************************************************/
     518             : /**                                                                   **/
     519             : /**                          MULTIPLICATION                           **/
     520             : /**                                                                   **/
     521             : /***********************************************************************/
     522             : /* assume ny > 0 */
     523             : INLINE GEN
     524  5322354708 : muluispec(ulong x, GEN y, long ny)
     525             : {
     526  5322354708 :   if (ny == 1)
     527  4243912196 :     return muluu(x, *y);
     528             :   else
     529             :   {
     530  1078442512 :     long lz = ny+3;
     531  1078442512 :     GEN z = cgeti(lz);
     532  1088172482 :     ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
     533  1089396948 :     if (hi) { z[lz - 1] = hi; } else lz--;
     534  1089396948 :     z[1] = evalsigne(1) | evallgefint(lz);
     535  1089396948 :     return z;
     536             :   }
     537             : }
     538             : 
     539             : /* a + b*|y| */
     540             : GEN
     541           0 : addumului(ulong a, ulong b, GEN y)
     542             : {
     543             :   GEN z;
     544             :   long i, lz;
     545             :   ulong hi;
     546           0 :   if (!b || !signe(y)) return utoi(a);
     547           0 :   lz = lgefint(y)+1;
     548           0 :   z = cgeti(lz);
     549           0 :   z[2]=a;
     550           0 :   for(i=3;i<lz;i++) z[i]=0;
     551           0 :   hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
     552           0 :   if (hi) z[lz-1]=hi; else lz--;
     553           0 :   z[1] = evalsigne(1) | evallgefint(lz);
     554           0 :   return gc_const((pari_sp)z, z);
     555             : }
     556             : 
     557             : /***********************************************************************/
     558             : /**                                                                   **/
     559             : /**                          DIVISION                                 **/
     560             : /**                                                                   **/
     561             : /***********************************************************************/
     562             : 
     563             : ulong
     564  1189213311 : umodiu(GEN y, ulong x)
     565             : {
     566  1189213311 :   long sy=signe(y);
     567             :   ulong hi;
     568  1189213311 :   if (!x) pari_err_INV("umodiu",gen_0);
     569  1190045336 :   if (!sy) return 0;
     570   862838500 :   hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
     571   862838500 :   if (!hi) return 0;
     572   848201791 :   return (sy > 0)? hi: x - hi;
     573             : }
     574             : 
     575             : /* return |y| \/ x */
     576             : GEN
     577   109526993 : absdiviu_rem(GEN y, ulong x, ulong *rem)
     578             : {
     579             :   long ly;
     580             :   GEN z;
     581             : 
     582   109526993 :   if (!x) pari_err_INV("absdiviu_rem",gen_0);
     583   109532237 :   if (!signe(y)) { *rem = 0; return gen_0; }
     584             : 
     585    87527899 :   ly = lgefint(y);
     586    87527899 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
     587             : 
     588    73081710 :   z = cgeti(ly);
     589    73081168 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     590    73081409 :   if (z [ly - 1] == 0) ly--;
     591    73081409 :   z[1] = evallgefint(ly) | evalsigne(1);
     592    73081409 :   return z;
     593             : }
     594             : 
     595             : GEN
     596    83917125 : divis_rem(GEN y, long x, long *rem)
     597             : {
     598    83917125 :   long sy=signe(y),ly,s;
     599             :   GEN z;
     600             : 
     601    83917125 :   if (!x) pari_err_INV("divis_rem",gen_0);
     602    83925312 :   if (!sy) { *rem = 0; return gen_0; }
     603    59178922 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     604             : 
     605    59178922 :   ly = lgefint(y);
     606    59178922 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
     607             : 
     608    41129074 :   z = cgeti(ly);
     609    41123367 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     610    41123513 :   if (sy<0) *rem = -  *rem;
     611    41123513 :   if (z[ly - 1] == 0) ly--;
     612    41123513 :   z[1] = evallgefint(ly) | evalsigne(s);
     613    41123513 :   return z;
     614             : }
     615             : 
     616             : GEN
     617      965811 : divis(GEN y, long x)
     618             : {
     619      965811 :   long sy=signe(y),ly,s;
     620             :   GEN z;
     621             : 
     622      965811 :   if (!x) pari_err_INV("divis",gen_0);
     623      965811 :   if (!sy) return gen_0;
     624      965763 :   if (x<0) { s = -sy; x = -x; } else s=sy;
     625             : 
     626      965763 :   ly = lgefint(y);
     627      965763 :   if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
     628             : 
     629      965451 :   z = cgeti(ly);
     630      965448 :   (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     631      965448 :   if (z[ly - 1] == 0) ly--;
     632      965448 :   z[1] = evallgefint(ly) | evalsigne(s);
     633      965448 :   return z;
     634             : }
     635             : 
     636             : /* We keep llx bits of x and lly bits of y*/
     637             : static GEN
     638    74634753 : divrr_with_gmp(GEN x, GEN y)
     639             : {
     640    74634753 :   long lx=RNLIMBS(x),ly=RNLIMBS(y);
     641    74634753 :   long lw=minss(lx,ly);
     642    74635246 :   long lly=minss(lw+1,ly);
     643    74635831 :   GEN  w = cgetg(lw+2, t_REAL);
     644    74615655 :   long lu=lw+lly;
     645    74615655 :   long llx=minss(lu,lx);
     646    74613066 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     647    74595777 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     648             :   mp_limb_t *q,*r;
     649    74568833 :   long e=expo(x)-expo(y);
     650    74568833 :   long sx=signe(x),sy=signe(y);
     651    74568833 :   xmpn_mirrorcopy(z,RLIMBS(y),lly);
     652    74589295 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     653    74628099 :   xmpn_zero(u,lu-llx);
     654    74652430 :   q = (mp_limb_t *)new_chunk(lw+1);
     655    74638435 :   r = (mp_limb_t *)new_chunk(lly);
     656             : 
     657    74617282 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     658             : 
     659             :   /*Round up: This is not exactly correct we should test 2*r>z*/
     660    74658865 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     661    36943639 :     mpn_add_1(q,q,lw+1,1);
     662             : 
     663    74658884 :   xmpn_mirrorcopy(RLIMBS(w),q,lw);
     664             : 
     665    74657472 :   if (q[lw] == 0) e--;
     666    41449361 :   else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
     667           0 :   else { w[2] = HIGHBIT; e++; }
     668    74656743 :   if (sy < 0) sx = -sx;
     669    74656743 :   w[1] = evalsigne(sx) | evalexpo(e);
     670    74651680 :   return gc_const((pari_sp)w, w);
     671             : }
     672             : 
     673             : /* We keep llx bits of x and lly bits of y*/
     674             : static GEN
     675    35069711 : divri_with_gmp(GEN x, GEN y)
     676             : {
     677    35069711 :   long llx=RNLIMBS(x),ly=NLIMBS(y);
     678    35069711 :   long lly=minss(llx+1,ly);
     679    35070115 :   GEN  w = cgetg(llx+2, t_REAL);
     680    35065468 :   long lu=llx+lly, ld=ly-lly;
     681    35065468 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     682    35062514 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     683             :   mp_limb_t *q,*r;
     684    35057299 :   long sh=bfffo(y[ly+1]);
     685    35057299 :   long e=expo(x)-expi(y);
     686    35058551 :   long sx=signe(x),sy=signe(y);
     687    35058551 :   if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
     688      682074 :   else xmpn_copy(z,LIMBS(y)+ld,lly);
     689    35060702 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     690    35067290 :   xmpn_zero(u,lu-llx);
     691    35072493 :   q = (mp_limb_t *)new_chunk(llx+1);
     692    35070407 :   r = (mp_limb_t *)new_chunk(lly);
     693             : 
     694    35065895 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     695             : 
     696             :   /*Round up: This is not exactly correct we should test 2*r>z*/
     697    35074795 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     698    15987494 :     mpn_add_1(q,q,llx+1,1);
     699             : 
     700    35074807 :   xmpn_mirrorcopy(RLIMBS(w),q,llx);
     701             : 
     702    35074418 :   if (q[llx] == 0) e--;
     703    10522066 :   else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
     704           0 :   else { w[2] = HIGHBIT; e++; }
     705    35074333 :   if (sy < 0) sx = -sx;
     706    35074333 :   w[1] = evalsigne(sx) | evalexpo(e);
     707    35073250 :   return gc_const((pari_sp)w, w);
     708             : }
     709             : 
     710             : GEN
     711   150754229 : divri(GEN x, GEN y)
     712             : {
     713   150754229 :   long  s = signe(y);
     714             : 
     715   150754229 :   if (!s) pari_err_INV("divri",gen_0);
     716   150754444 :   if (!signe(x)) return real_0_bit(expo(x) - expi(y));
     717   150525279 :   if (!is_bigint(y)) {
     718   115460556 :     GEN z = divru(x, y[2]);
     719   115459211 :     if (s < 0) togglesign(z);
     720   115459237 :     return z;
     721             :   }
     722    35064098 :   return divri_with_gmp(x,y);
     723             : }
     724             : 
     725             : GEN
     726   140165762 : divrr(GEN x, GEN y)
     727             : {
     728   140165762 :   long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
     729             :   ulong y0,y1;
     730             :   GEN r, r1;
     731             : 
     732   140165762 :   if (!sy) pari_err_INV("divrr",y);
     733   140174013 :   e = expo(x) - expo(y);
     734   140174013 :   if (!sx) return real_0_bit(e);
     735   139677012 :   if (sy<0) sx = -sx;
     736             : 
     737   139677012 :   lx=lg(x); ly=lg(y);
     738   139677012 :   if (ly==3)
     739             :   {
     740    26204921 :     ulong k = x[2], l = (lx>3)? x[3]: 0;
     741             :     LOCAL_HIREMAINDER;
     742    26204921 :     if (k < uel(y,2)) e--;
     743             :     else
     744             :     {
     745     8232577 :       l >>= 1; if (k&1) l |= HIGHBIT;
     746     8232577 :       k >>= 1;
     747             :     }
     748    26204921 :     hiremainder = k; k = divll(l,y[2]);
     749    26204921 :     if (hiremainder > (uel(y,2) >> 1) && !++k) { k = HIGHBIT; e++; }
     750    26204921 :     r = cgetg(3, t_REAL);
     751    26202554 :     r[1] = evalsigne(sx) | evalexpo(e);
     752    26200887 :     r[2] = k; return r;
     753             :   }
     754             : 
     755   113472091 :   if (ly >= prec2lg(DIVRR_GMP_LIMIT))
     756    74634536 :     return divrr_with_gmp(x,y);
     757             : 
     758    38887542 :   lr = minss(lx,ly); r = new_chunk(lr);
     759    38890116 :   r1 = r-1;
     760   131794964 :   r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
     761    38890116 :   r1[lr] = (lx>ly)? x[lr]: 0;
     762    38890116 :   y0 = y[2]; y1 = y[3];
     763   170648657 :   for (i=0; i<lr-1; i++)
     764             :   { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
     765             :     ulong k, qp;
     766             :     LOCAL_HIREMAINDER;
     767             :     LOCAL_OVERFLOW;
     768             : 
     769   131758541 :     if (uel(r1,1) == y0) { qp = ULONG_MAX; k = addll(y0,r1[2]); }
     770             :     else
     771             :     {
     772   131518842 :       if (uel(r1,1) > y0) /* can't happen if i=0 */
     773             :       {
     774           0 :         GEN y1 = y+1;
     775           0 :         j = lr-i; r1[j] = subll(r1[j],y1[j]);
     776           0 :         for (j--; j>0; j--) r1[j] = subllx(r1[j],y1[j]);
     777           0 :         j=i; do uel(r,--j)++; while (j && !r[j]);
     778             :       }
     779   131518842 :       hiremainder = r1[1]; overflow = 0;
     780   131518842 :       qp = divll(r1[2],y0); k = hiremainder;
     781             :     }
     782   131758541 :     j = lr-i+1;
     783   131758541 :     if (!overflow)
     784             :     {
     785             :       long k3, k4;
     786   131625652 :       k3 = mulll(qp,y1);
     787   131625652 :       if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
     788    38946050 :         k4 = subll(hiremainder,k);
     789             :       else
     790             :       {
     791    92679602 :         k3 = subll(k3, r1[3]);
     792    92679602 :         k4 = subllx(hiremainder,k);
     793             :       }
     794   168989642 :       while (!overflow && k4) { qp--; k3=subll(k3,y1); k4=subllx(k4,y0); }
     795             :     }
     796   131758541 :     if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
     797   390433952 :     for (j--; j>1; j--)
     798             :     {
     799   258675411 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     800   258675411 :       hiremainder += overflow;
     801             :     }
     802   131758541 :     if (uel(r1,1) != hiremainder)
     803             :     {
     804      177514 :       if (uel(r1,1) < hiremainder)
     805             :       {
     806      177514 :         qp--;
     807      177514 :         j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
     808      505608 :         for (j--; j>1; j--) r1[j] = addllx(r1[j], y[j]);
     809             :       }
     810             :       else
     811             :       {
     812           0 :         uel(r1,1) -= hiremainder;
     813           0 :         while (r1[1])
     814             :         {
     815           0 :           qp++; if (!qp) { j=i; do uel(r,--j)++; while (j && !r[j]); }
     816           0 :           j = lr-i-(lr-i>=ly); r1[j] = subll(r1[j],y[j]);
     817           0 :           for (j--; j>1; j--) r1[j] = subllx(r1[j],y[j]);
     818           0 :           uel(r1,1) -= overflow;
     819             :         }
     820             :       }
     821             :     }
     822   131758541 :     *++r1 = qp;
     823             :   }
     824             :   /* i = lr-1 */
     825             :   /* round correctly */
     826    38890116 :   if (uel(r1,1) > (y0>>1))
     827             :   {
     828    19659102 :     j=i; do uel(r,--j)++; while (j && !r[j]);
     829             :   }
     830   131972942 :   r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
     831    38890116 :   if (r[0] == 0) e--;
     832    14010036 :   else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
     833             :   else { /* possible only when rounding up to 0x2 0x0 ... */
     834          18 :     r[2] = (long)HIGHBIT; e++;
     835             :   }
     836    38889346 :   r[0] = evaltyp(t_REAL)|evallg(lr);
     837    38945054 :   r[1] = evalsigne(sx) | evalexpo(e);
     838    38935776 :   return r;
     839             : }
     840             : 
     841             : /* Integer division x / y: such that sign(r) = sign(x)
     842             :  *   if z = ONLY_REM return remainder, otherwise return quotient
     843             :  *   if z != NULL set *z to remainder
     844             :  *   *z is the last object on stack (and thus can be disposed of with cgiv
     845             :  *   instead of gerepile)
     846             :  * If *z is zero, we put gen_0 here and no copy.
     847             :  * space needed: lx + ly
     848             :  */
     849             : GEN
     850  1889723474 : dvmdii(GEN x, GEN y, GEN *z)
     851             : {
     852  1889723474 :   long sx=signe(x),sy=signe(y);
     853             :   long lx, ly, lq;
     854             :   pari_sp av;
     855             :   GEN r,q;
     856             : 
     857  1889723474 :   if (!sy) pari_err_INV("dvmdii",y);
     858  1890778450 :   if (!sx)
     859             :   {
     860    66022733 :     if (!z || z == ONLY_REM) return gen_0;
     861    40529670 :     *z=gen_0; return gen_0;
     862             :   }
     863  1824755717 :   lx=lgefint(x);
     864  1824755717 :   ly=lgefint(y); lq=lx-ly;
     865  1824755717 :   if (lq <= 0)
     866             :   {
     867  1135491842 :     if (lq == 0)
     868             :     {
     869  1027992079 :       long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
     870  1027992079 :       if (s>0) goto DIVIDE;
     871   350454389 :       if (s==0)
     872             :       {
     873    30246559 :         if (z == ONLY_REM) return gen_0;
     874    20075599 :         if (z) *z = gen_0;
     875    20075599 :         if (sx < 0) sy = -sy;
     876    20075599 :         return stoi(sy);
     877             :       }
     878             :     }
     879   427707593 :     if (z == ONLY_REM) return icopy(x);
     880    12273232 :     if (z) *z = icopy(x);
     881    12273232 :     return gen_0;
     882             :   }
     883   689263875 : DIVIDE: /* quotient is nonzero */
     884  1366801565 :   av=avma; if (sx<0) sy = -sy;
     885  1366801565 :   if (ly==3)
     886             :   {
     887   557927394 :     ulong lq = lx;
     888             :     ulong si;
     889   557927394 :     q  = cgeti(lq);
     890   557343610 :     si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
     891   558810704 :     if (q[lq - 1] == 0) lq--;
     892   558810704 :     if (z == ONLY_REM)
     893             :     {
     894   324611244 :       if (!si) return gc_const(av, gen_0);
     895   282562203 :       set_avma(av); r = cgeti(3);
     896   281708494 :       r[1] = evalsigne(sx) | evallgefint(3);
     897   281708494 :       r[2] = si; return r;
     898             :     }
     899   234199460 :     q[1] = evalsigne(sy) | evallgefint(lq);
     900   234199460 :     if (!z) return q;
     901   230146102 :     if (!si) { *z=gen_0; return q; }
     902    59652194 :     r=cgeti(3);
     903    59665781 :     r[1] = evalsigne(sx) | evallgefint(3);
     904    59665781 :     r[2] = si; *z=r; return q;
     905             :   }
     906   808874171 :   if (z == ONLY_REM)
     907             :   {
     908   787311578 :     ulong lr = lgefint(y);
     909   787311578 :     ulong lq = lgefint(x)-lgefint(y)+3;
     910   787311578 :     GEN r = cgeti(lr);
     911   782067962 :     GEN q = cgeti(lq);
     912   775493147 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     913   789274757 :     if (!r[lr - 1])
     914             :     {
     915   806211213 :       while(lr>2 && !r[lr - 1]) lr--;
     916   359815976 :       if (lr == 2) return gc_const(av, gen_0); /* exact division */
     917             :     }
     918   774798703 :     r[1] = evalsigne(sx) | evallgefint(lr);
     919   774798703 :     return gc_const((pari_sp)r, r);
     920             :   }
     921             :   else
     922             :   {
     923    21562593 :     ulong lq = lgefint(x)-lgefint(y)+3;
     924    21562593 :     ulong lr = lgefint(y);
     925    21562593 :     GEN q = cgeti(lq);
     926    26564701 :     GEN r = cgeti(lr);
     927    26550559 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     928    26575185 :     if (q[lq - 1] == 0) lq--;
     929    26575185 :     q[1] = evalsigne(sy) | evallgefint(lq);
     930    26575185 :     if (!z) return gc_const((pari_sp)q, q);
     931    22880137 :     if (!r[lr - 1])
     932             :     {
     933    28940005 :       while(lr>2 && !r[lr - 1]) lr--;
     934     5596044 :       if (lr == 2) { *z = gen_0; return gc_const((pari_sp)q, q); } /* exact */
     935             :     }
     936    18786938 :     r[1] = evalsigne(sx) | evallgefint(lr);
     937    18786938 :     *z = gc_const((pari_sp)r, r); return q;
     938             :   }
     939             : }
     940             : 
     941             : /* Montgomery reduction.
     942             :  * N has k words, assume T >= 0 has less than 2k.
     943             :  * Return res := T / B^k mod N, where B = 2^BIL
     944             :  * such that 0 <= res < T/B^k + N  and  res has less than k words */
     945             : GEN
     946    37306967 : red_montgomery(GEN T, GEN N, ulong inv)
     947             : {
     948             :   pari_sp av;
     949             :   GEN Te, Td, Ne, Nd, scratch;
     950    37306967 :   ulong i, j, m, t, d, k = NLIMBS(N);
     951             :   int carry;
     952             :   LOCAL_HIREMAINDER;
     953             :   LOCAL_OVERFLOW;
     954             : 
     955    37306967 :   if (k == 0) return gen_0;
     956    37306967 :   d = NLIMBS(T); /* <= 2*k */
     957    37306967 :   if (d == 0) return gen_0;
     958             : #ifdef DEBUG
     959             :   if (d > 2*k) pari_err_BUG("red_montgomery");
     960             : #endif
     961    37306958 :   if (k == 1)
     962             :   { /* as below, special cased for efficiency */
     963      163341 :     ulong n = uel(N,2);
     964      163341 :     if (d == 1) {
     965      163194 :       hiremainder = uel(T,2);
     966      163194 :       m = hiremainder * inv;
     967      163194 :       (void)addmul(m, n); /* t + m*n = 0 */
     968      163194 :       return utoi(hiremainder);
     969             :     } else { /* d = 2 */
     970         147 :       hiremainder = uel(T,2);
     971         147 :       m = hiremainder * inv;
     972         147 :       (void)addmul(m, n); /* t + m*n = 0 */
     973         147 :       t = addll(hiremainder, uel(T,3));
     974         147 :       if (overflow) t -= n; /* t > n doesn't fit in 1 word */
     975         147 :       return utoi(t);
     976             :     }
     977             :   }
     978             :   /* assume k >= 2 */
     979    37143617 :   av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
     980             : 
     981             :   /* copy T to scratch space (pad with zeroes to 2k words) */
     982    36895637 :   Td = scratch;
     983    36895637 :   Te = T + 2;
     984   532993188 :   for (i=0; i < d     ; i++) *Td++ = *Te++;
     985    66636996 :   for (   ; i < (k<<1); i++) *Td++ = 0;
     986             : 
     987    36895637 :   Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
     988    36895637 :   Ne = N + 1;       /* 1 beyond end of N mantissa */
     989             : 
     990    36895637 :   carry = 0;
     991   280169022 :   for (i=0; i<k; i++) /* set T := T/B nod N, k times */
     992             :   {
     993   243273385 :     Td = Te; /* one beyond end of (new) T mantissa */
     994   243273385 :     Nd = Ne;
     995   243273385 :     hiremainder = *++Td;
     996   243273385 :     m = hiremainder * inv; /* solve T + m N = O(B) */
     997             : 
     998             :     /* set T := (T + mN) / B */
     999   243273385 :     Te = Td;
    1000   243273385 :     (void)addmul(m, *++Nd); /* = 0 */
    1001  2135761109 :     for (j=1; j<k; j++)
    1002             :     {
    1003  1892487724 :       t = addll(addmul(m, *++Nd), *++Td);
    1004  1892487724 :       *Td = t;
    1005  1892487724 :       hiremainder += overflow;
    1006             :     }
    1007   243273385 :     t = addll(hiremainder, *++Td); *Td = t + carry;
    1008   243273385 :     carry = (overflow || (carry && *Td == 0));
    1009             :   }
    1010    36895637 :   if (carry)
    1011             :   { /* Td > N overflows (k+1 words), set Td := Td - N */
    1012       57163 :     GEN NE = N + k+1;
    1013       57163 :     Td = Te;
    1014       57163 :     Nd = Ne;
    1015       57163 :     t = subll(*++Td, *++Nd); *Td = t;
    1016      514333 :     while (Nd < NE) { t = subllx(*++Td, *++Nd); *Td = t; }
    1017             :   }
    1018             : 
    1019             :   /* copy result */
    1020    36895637 :   Td = (GEN)av - 1; /* *Td = high word of final result */
    1021    41032563 :   while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
    1022    36895637 :   k = Td - Te; if (!k) return gc_const(av, gen_0);
    1023    36895637 :   Td = (GEN)av - k; /* will write mantissa there */
    1024    36895637 :   (void)memmove(Td, Te+1, k*sizeof(long));
    1025    36895637 :   Td -= 2;
    1026    36895637 :   Td[0] = evaltyp(t_INT) | evallg(k+2);
    1027    37094752 :   Td[1] = evalsigne(1) | evallgefint(k+2);
    1028             : #ifdef DEBUG
    1029             : {
    1030             :   long l = NLIMBS(N), s = BITS_IN_LONG*l;
    1031             :   GEN R = int2n(s);
    1032             :   GEN res = remii(mulii(T, Fp_inv(R, N)), N);
    1033             :   if (k > lgefint(N)
    1034             :     || !equalii(remii(Td,N),res)
    1035             :     || cmpii(Td, addii(shifti(T, -s), N)) >= 0) pari_err_BUG("red_montgomery");
    1036             : }
    1037             : #endif
    1038    37094752 :   return gc_const((pari_sp)Td, Td);
    1039             : }
    1040             : 
    1041             : /* EXACT INTEGER DIVISION */
    1042             : 
    1043             : /* use undocumented GMP interface */
    1044             : static void
    1045   109826504 : GEN2mpz(mpz_t X, GEN x)
    1046             : {
    1047   109826504 :   long l = lgefint(x)-2;
    1048   109826504 :   X->_mp_alloc = l;
    1049   109826504 :   X->_mp_size = signe(x) > 0? l: -l;
    1050   109826504 :   X->_mp_d = LIMBS(x);
    1051   109826504 : }
    1052             : static void
    1053    54914464 : mpz2GEN(GEN z, mpz_t Z)
    1054             : {
    1055    54914464 :   long l = Z->_mp_size;
    1056    54914464 :   z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
    1057    54914464 : }
    1058             : 
    1059             : #ifdef mpn_divexact_1
    1060             : static GEN
    1061   370547167 : diviuexact_i(GEN x, ulong y)
    1062             : {
    1063   370547167 :   long l = lgefint(x);
    1064   370547167 :   GEN z = cgeti(l);
    1065   370116914 :   mpn_divexact_1(LIMBS(z), LIMBS(x), NLIMBS(x), y);
    1066   370193278 :   if (z[l-1] == 0) l--;
    1067   370193278 :   z[1] = evallgefint(l) | evalsigne(signe(x));
    1068   370193278 :   return z;
    1069             : }
    1070             : #elif 1 && !defined(_WIN64) /* mpz_divexact_ui is not LLP64 friendly   */
    1071             :                             /* assume y != 0 and the division is exact */
    1072             : static GEN
    1073             : diviuexact_i(GEN x, ulong y)
    1074             : {
    1075             :   long l = lgefint(x);
    1076             :   GEN z = cgeti(l);
    1077             :   mpz_t X, Z;
    1078             :   GEN2mpz(X, x);
    1079             :   Z->_mp_alloc = l-2;
    1080             :   Z->_mp_size  = l-2;
    1081             :   Z->_mp_d = LIMBS(z);
    1082             :   mpz_divexact_ui(Z, X, y);
    1083             :   mpz2GEN(z, Z); return z;
    1084             : }
    1085             : #else
    1086             : /* assume y != 0 and the division is exact */
    1087             : static GEN
    1088             : diviuexact_i(GEN x, ulong y)
    1089             : {
    1090             :   /*TODO: implement true exact division.*/
    1091             :   return divii(x,utoi(y));
    1092             : }
    1093             : #endif
    1094             : 
    1095             : GEN
    1096    30747054 : diviuexact(GEN x, ulong y)
    1097             : {
    1098             :   GEN z;
    1099    30747054 :   if (!signe(x)) return gen_0;
    1100    29618482 :   z = diviuexact_i(x, y);
    1101    29616175 :   if (lgefint(z) == 2) pari_err_OP("exact division", x, utoi(y));
    1102    29616099 :   return z;
    1103             : }
    1104             : 
    1105             : /* Find z such that x=y*z, knowing that y | x (unchecked) */
    1106             : GEN
    1107   472196670 : diviiexact(GEN x, GEN y)
    1108             : {
    1109             :   GEN z;
    1110   472196670 :   if (!signe(y)) pari_err_INV("diviiexact",y);
    1111   472357864 :   if (!signe(x)) return gen_0;
    1112   396033574 :   if (lgefint(y) == 3)
    1113             :   {
    1114   341130438 :     z = diviuexact_i(x, y[2]);
    1115   340783915 :     if (signe(y) < 0) togglesign(z);
    1116             :   }
    1117             :   else
    1118             :   {
    1119    54903136 :     long l = lgefint(x);
    1120             :     mpz_t X, Y, Z;
    1121    54903136 :     z = cgeti(l);
    1122    54913380 :     GEN2mpz(X, x);
    1123    54913316 :     GEN2mpz(Y, y);
    1124    54913212 :     Z->_mp_alloc = l-2;
    1125    54913212 :     Z->_mp_size  = l-2;
    1126    54913212 :     Z->_mp_d = LIMBS(z);
    1127    54913212 :     mpz_divexact(Z, X, Y);
    1128    54914475 :     mpz2GEN(z, Z);
    1129             :   }
    1130   395698376 :   if (lgefint(z) == 2) pari_err_OP("exact division", x, y);
    1131   395666465 :   return z;
    1132             : }
    1133             : 
    1134             : /* assume yz != and yz | x */
    1135             : GEN
    1136      199451 : diviuuexact(GEN x, ulong y, ulong z)
    1137             : {
    1138             :   long tmp[4];
    1139             :   ulong t;
    1140             :   LOCAL_HIREMAINDER;
    1141      199451 :   t = mulll(y, z);
    1142      199451 :   if (!hiremainder) return diviuexact(x, t);
    1143           0 :   tmp[0] = evaltyp(t_INT)|_evallg(4);
    1144           0 :   tmp[1] = evalsigne(1)|evallgefint(4);
    1145           0 :   tmp[2] = t;
    1146           0 :   tmp[3] = hiremainder;
    1147           0 :   return diviiexact(x, tmp);
    1148             : }
    1149             : 
    1150             : /********************************************************************/
    1151             : /**                                                                **/
    1152             : /**               INTEGER MULTIPLICATION                           **/
    1153             : /**                                                                **/
    1154             : /********************************************************************/
    1155             : 
    1156             : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
    1157             : GEN
    1158  5662796760 : muliispec(GEN x, GEN y, long nx, long ny)
    1159             : {
    1160             :   GEN zd;
    1161             :   long lz;
    1162             :   ulong hi;
    1163             : 
    1164  5662796760 :   if (nx < ny) swapspec(x,y, nx,ny);
    1165  5662796760 :   if (!ny) return gen_0;
    1166  5662796760 :   if (ny == 1) return muluispec((ulong)*y, x, nx);
    1167             : 
    1168  1024202649 :   lz = nx+ny+2;
    1169  1024202649 :   zd = cgeti(lz);
    1170  1027308022 :   hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
    1171  1034531061 :   if (!hi) lz--;
    1172             :   /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
    1173             : 
    1174  1034531061 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1175  1034531061 :   return zd;
    1176             : }
    1177             : GEN
    1178      221701 : muluui(ulong x, ulong y, GEN z)
    1179             : {
    1180      221701 :   long t, s = signe(z);
    1181             :   GEN r;
    1182             :   LOCAL_HIREMAINDER;
    1183             : 
    1184      221701 :   if (!x || !y || !signe(z)) return gen_0;
    1185      221325 :   t = mulll(x,y);
    1186      221325 :   if (!hiremainder)
    1187      221343 :     r = muluispec(t, z+2, lgefint(z)-2);
    1188             :   else
    1189             :   {
    1190             :     long tmp[2];
    1191           0 :     tmp[1] = hiremainder;
    1192           0 :     tmp[0] = t;
    1193           0 :     r = muliispec(z+2,tmp, lgefint(z)-2, 2);
    1194             :   }
    1195      221284 :   setsigne(r,s); return r;
    1196             : }
    1197             : 
    1198             : GEN
    1199   986449484 : sqrispec(GEN x, long nx)
    1200             : {
    1201             :   GEN zd;
    1202             :   long lz;
    1203             : 
    1204   986449484 :   if (!nx) return gen_0;
    1205   478913959 :   if (nx==1) return sqru(*x);
    1206             : 
    1207   274389767 :   lz = (nx<<1)+2;
    1208   274389767 :   zd = cgeti(lz);
    1209             : #ifdef mpn_sqr
    1210   271613093 :   mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
    1211             : #else
    1212             :   mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
    1213             : #endif
    1214   275214591 :   if (zd[lz-1]==0) lz--;
    1215             : 
    1216   275214591 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1217   275214591 :   return zd;
    1218             : }
    1219             : 
    1220             : INLINE GEN
    1221    41339451 : sqrispec_mirror(GEN x, long nx)
    1222             : {
    1223    41339451 :   GEN cx=new_chunk(nx);
    1224             :   GEN z;
    1225    41300326 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1226    41410784 :   z=sqrispec(cx, nx);
    1227    41471485 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1228    41471299 :   return z;
    1229             : }
    1230             : 
    1231             : /* leaves garbage on the stack. */
    1232             : INLINE GEN
    1233    84034384 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
    1234             : {
    1235             :   GEN cx, cy, z;
    1236    84034384 :   long s = 0;
    1237   113014877 :   while (nx && x[nx-1]==0) { nx--; s++; }
    1238   118292862 :   while (ny && y[ny-1]==0) { ny--; s++; }
    1239    84034384 :   cx=new_chunk(nx); cy=new_chunk(ny);
    1240    83516821 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1241    84492127 :   xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
    1242    85005546 :   z =  nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
    1243    85046917 :   if (s)
    1244             :   {
    1245     7624927 :     long i, lz = lgefint(z) + s;
    1246     7624927 :     (void)new_chunk(s);
    1247     7624923 :     z -= s;
    1248    70863891 :     for (i=0; i<s; i++) z[2+i]=0;
    1249     7624923 :     z[1] = evalsigne(1) | evallgefint(lz);
    1250     7624923 :     z[0] = evaltyp(t_INT) | evallg(lz);
    1251             :   }
    1252    85046912 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1253    85717349 :   return z;
    1254             : }
    1255             : 
    1256             : /* x % (2^n), assuming n >= 0 */
    1257             : GEN
    1258    36392778 : remi2n(GEN x, long n)
    1259             : {
    1260             :   ulong hi;
    1261    36392778 :   long l, k, lx, ly, sx = signe(x);
    1262             :   GEN z, xd, zd;
    1263             : 
    1264    36392778 :   if (!sx || !n) return gen_0;
    1265             : 
    1266    36068736 :   k = dvmdsBIL(n, &l);
    1267    36080959 :   lx = lgefint(x);
    1268    36080959 :   if (lx < k+3) return icopy(x);
    1269             : 
    1270    35312999 :   xd = x + (2 + k);
    1271             :   /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
    1272             :    *              ^--- initial xd  */
    1273    35312999 :   hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
    1274    35312999 :   if (!hi)
    1275             :   { /* strip leading zeroes from result */
    1276     2982500 :     xd--; while (k && !*xd) { k--; xd--; }
    1277     2796937 :     if (!k) return gen_0;
    1278     1857375 :     ly = k+2;
    1279             :   }
    1280             :   else
    1281    32516062 :     ly = k+3;
    1282             : 
    1283    34373437 :   zd = z = cgeti(ly);
    1284    34344349 :   *++zd = evalsigne(sx) | evallgefint(ly);
    1285   499900933 :   xd = x+1; for ( ;k; k--) *++zd = *++xd;
    1286    34344349 :   if (hi) *++zd = hi;
    1287    34344349 :   return z;
    1288             : }
    1289             : 
    1290             : /********************************************************************/
    1291             : /**                                                                **/
    1292             : /**                      INTEGER SQUARE ROOT                       **/
    1293             : /**                                                                **/
    1294             : /********************************************************************/
    1295             : 
    1296             : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
    1297             :  * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
    1298             :  * remainder is 0. R = NULL is allowed. */
    1299             : GEN
    1300     5099913 : sqrtremi(GEN a, GEN *r)
    1301             : {
    1302     5099913 :   long l, na = NLIMBS(a);
    1303             :   mp_size_t nr;
    1304             :   GEN S;
    1305     5099913 :   if (!na) {
    1306         724 :     if (r) *r = gen_0;
    1307         724 :     return gen_0;
    1308             :   }
    1309     5099189 :   l = (na + 5) >> 1; /* 2 + ceil(na/2) */
    1310     5099189 :   S = cgetipos(l);
    1311     5099147 :   if (r) {
    1312     1309127 :     GEN R = cgeti(2 + na);
    1313     1309127 :     nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
    1314     1309127 :     if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
    1315       25558 :     else    { set_avma((pari_sp)S); R = gen_0; }
    1316     1309127 :     *r = R;
    1317             :   }
    1318             :   else
    1319     3790020 :     (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
    1320     5099156 :   return S;
    1321             : }
    1322             : 
    1323             : /* compute sqrt(|a|), assuming a != 0 */
    1324             : GEN
    1325   124951321 : sqrtr_abs(GEN a)
    1326             : {
    1327             :   GEN res;
    1328             :   mp_limb_t *b, *c;
    1329   124951321 :   long l = RNLIMBS(a), e = expo(a), er = e>>1;
    1330             :   long n;
    1331   124951321 :   res = cgetg(2 + l, t_REAL);
    1332   124860617 :   res[1] = evalsigne(1) | evalexpo(er);
    1333   124945879 :   if (e&1)
    1334             :   {
    1335    52527507 :     b = (mp_limb_t *) new_chunk(l<<1);
    1336    52513584 :     xmpn_zero(b,l);
    1337    52514926 :     xmpn_mirrorcopy(b+l, RLIMBS(a), l);
    1338    52534131 :     c = (mp_limb_t *) new_chunk(l);
    1339    52530443 :     n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
    1340    52561201 :     if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
    1341             :   }
    1342             :   else
    1343             :   {
    1344             :     ulong u;
    1345    72418372 :     b = (mp_limb_t *) mantissa2nr(a,-1);
    1346    72460805 :     b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
    1347    72460805 :     b = (mp_limb_t *) new_chunk(l);
    1348    72440020 :     xmpn_zero(b,l+1); /* overwrites the former b[0] */
    1349    72443178 :     c = (mp_limb_t *) new_chunk(l + 1);
    1350    72408691 :     n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
    1351    72499619 :     u = (ulong)*c++;
    1352    72499619 :     if ( u&HIGHBIT || (u == ~HIGHBIT &&
    1353           0 :              (n>l || (n==l && mpn_cmp(b,c,l)>0))))
    1354    35755858 :       mpn_add_1(c,c,l,1);
    1355             :   }
    1356   125068591 :   xmpn_mirrorcopy(RLIMBS(res),c,l);
    1357   125043929 :   return gc_const((pari_sp)res, res);
    1358             : }
    1359             : 
    1360             : /* Normalize a nonnegative integer */
    1361             : GEN
    1362   291626753 : int_normalize(GEN x, long known_zero_words)
    1363             : {
    1364   291626753 :   long i =  lgefint(x) - 1 - known_zero_words;
    1365  2134944192 :   for ( ; i > 1; i--)
    1366  2085606843 :     if (x[i]) { setlgefint(x, i+1); return x; }
    1367    49337349 :   x[1] = evalsigne(0) | evallgefint(2); return x;
    1368             : }
    1369             : 
    1370             : /********************************************************************
    1371             :  **                                                                **
    1372             :  **                           Base Conversion                      **
    1373             :  **                                                                **
    1374             :  ********************************************************************/
    1375             : 
    1376             : ulong *
    1377      427410 : convi(GEN x, long *l)
    1378             : {
    1379      427410 :   long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
    1380      427410 :   GEN str = cgetg(n+1, t_VECSMALL);
    1381      427410 :   unsigned char *res = (unsigned char*) GSTR(str);
    1382      427410 :   long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
    1383             :   long lz;
    1384             :   ulong *z;
    1385             :   long i, j;
    1386             :   unsigned char *t;
    1387      427410 :   while (!*res) {res++; llz--;} /*Strip leading zeros*/
    1388      427410 :   lz  = (8+llz)/9;
    1389      427410 :   z = (ulong*)new_chunk(1+lz);
    1390      427410 :   t=res+llz+9;
    1391      857620 :   for(i=0;i<llz-8;i+=9)
    1392             :   {
    1393             :     ulong s;
    1394      430210 :     t-=18;
    1395      430210 :     s=*t++;
    1396     3871890 :     for (j=1; j<9;j++)
    1397     3441680 :       s=10*s+*t++;
    1398      430210 :     *z++=s;
    1399             :   }
    1400      427410 :   if (i<llz)
    1401             :   {
    1402      423479 :     unsigned char *t = res;
    1403      423479 :     ulong s=*t++;
    1404     1182510 :     for (j=i+1; j<llz;j++)
    1405      759031 :       s=10*s+*t++;
    1406      423479 :     *z++=s;
    1407             :   }
    1408      427410 :   *l = lz;
    1409      427410 :   return z;
    1410             : }

Generated by: LCOV version 1.16