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.18.1 lcov report (development 30387-12623f1693) Lines: 683 719 95.0 %
Date: 2025-07-11 09:22:46 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     1755987 : static void pari_gmp_free(void *ptr, size_t old_size){
      52     1755987 :   (void)old_size; pari_free(ptr);
      53     1755987 : }
      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        1112 : pari_kernel_init(void)
      61             : {
      62        1112 :   mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
      63        1112 :   mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
      64        1112 : }
      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        1104 : 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        1104 :   mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
      83        1104 :   if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
      84        1104 :   if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
      85        1104 :   if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
      86        1104 :   mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
      87        1104 : }
      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     6890517 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
      97             : {
      98    57047568 :   while (--n >= 0) x[n]=y[n];
      99     6890517 : }
     100             : 
     101             : INLINE void
     102   587305463 : xmpn_mirror(mp_limb_t *x, long n)
     103             : {
     104             :   long i;
     105  3922097916 :   for(i=0;i<(n>>1);i++)
     106             :   {
     107  3334792453 :     ulong m=x[i];
     108  3334792453 :     x[i]=x[n-1-i];
     109  3334792453 :     x[n-1-i]=m;
     110             :   }
     111   587305463 : }
     112             : 
     113             : INLINE void
     114   718729771 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
     115             : {
     116             :   long i;
     117  9936132870 :   for(i=0;i<n;i++)
     118  9217403099 :     z[i]=x[n-1-i];
     119   718729771 : }
     120             : 
     121             : INLINE void
     122   237525521 : xmpn_zero(mp_ptr x, mp_size_t n)
     123             : {
     124  2075766673 :   while (--n >= 0) x[n]=0;
     125   237525521 : }
     126             : 
     127             : INLINE GEN
     128    41271209 : icopy_ef(GEN x, long l)
     129             : {
     130    41271209 :   long lx = lgefint(x);
     131    41271209 :   const GEN y = cgeti(l);
     132             : 
     133   294283785 :   while (--lx > 0) y[lx]=x[lx];
     134    41269865 :   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 order to facilitate splitting: no
     140             :  * need to reintroduce codewords. */
     141             : 
     142             : /***********************************************************************/
     143             : /**                                                                   **/
     144             : /**                     ADDITION / SUBTRACTION                        **/
     145             : /**                                                                   **/
     146             : /***********************************************************************/
     147             : 
     148             : GEN
     149     2997894 : setloop(GEN a)
     150             : {
     151     2997894 :   pari_sp av = avma - 2 * sizeof(long);
     152     2997894 :   (void)cgetg(lgefint(a) + 3, t_VECSMALL);
     153     2997894 :   return icopy_avma(a, av); /* two cells of extra space after a */
     154             : }
     155             : 
     156             : /* we had a = setloop(?), then some incloops. Reset a to b */
     157             : GEN
     158      174328 : resetloop(GEN a, GEN b) {
     159      174328 :   a[0] = evaltyp(t_INT) | evallg(lgefint(b));
     160      174328 :   affii(b, a); return a;
     161             : }
     162             : 
     163             : /* assume a > 0, initialized by setloop. Do a++ */
     164             : static GEN
     165   103265964 : incpos(GEN a)
     166             : {
     167   103265964 :   long i, l = lgefint(a);
     168   103265969 :   for (i=2; i<l; i++)
     169   103270283 :     if (++uel(a,i)) return a;
     170           3 :   a[l] = 1; l++;
     171           3 :   a[0]=evaltyp(t_INT) | _evallg(l);
     172           3 :   a[1]=evalsigne(1) | evallgefint(l);
     173           3 :   return a;
     174             : }
     175             : 
     176             : /* assume a < 0, initialized by setloop. Do a++ */
     177             : static GEN
     178       66684 : incneg(GEN a)
     179             : {
     180       66684 :   long i, l = lgefint(a)-1;
     181       66684 :   if (uel(a,2)--)
     182             :   {
     183       66680 :     if (!a[l]) /* implies l = 2 */
     184             :     {
     185        1980 :       a[0] = evaltyp(t_INT) | _evallg(2);
     186        1980 :       a[1] = evalsigne(0) | evallgefint(2);
     187             :     }
     188       66680 :     return a;
     189             :   }
     190           5 :   for (i=3; i<=l; i++)
     191           5 :     if (uel(a,i)--) break;
     192           4 :   if (!a[l])
     193             :   {
     194           4 :     a[0] = evaltyp(t_INT) | _evallg(l);
     195           4 :     a[1] = evalsigne(-1) | evallgefint(l);
     196             :   }
     197           4 :   return a;
     198             : }
     199             : 
     200             : /* assume a initialized by setloop. Do a++ */
     201             : GEN
     202   103669911 : incloop(GEN a)
     203             : {
     204   103669911 :   switch(signe(a))
     205             :   {
     206      336652 :     case 0:
     207      336652 :       a[0]=evaltyp(t_INT) | _evallg(3);
     208      336652 :       a[1]=evalsigne(1) | evallgefint(3);
     209      336652 :       a[2]=1; return a;
     210       66684 :     case -1: return incneg(a);
     211   103266575 :     default: return incpos(a);
     212             :   }
     213             : }
     214             : 
     215             : INLINE GEN
     216  2783279784 : adduispec(ulong s, GEN x, long nx)
     217             : {
     218             :   GEN  zd;
     219             :   long lz;
     220             : 
     221  2783279784 :   if (nx == 1) return adduu(uel(x,0), s);
     222   900789471 :   lz = nx+3; zd = cgeti(lz);
     223   903757658 :   if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
     224     5022064 :     zd[lz-1]=1;
     225             :   else
     226   898740256 :     lz--;
     227   903762320 :   zd[1] = evalsigne(1) | evallgefint(lz);
     228   903762320 :   return zd;
     229             : }
     230             : 
     231             : GEN
     232   764425805 : adduispec_offset(ulong s, GEN x, long offset, long nx)
     233             : {
     234   764425805 :   GEN xd=x+2+offset;
     235   996168912 :   while (nx && *(xd+nx-1)==0) nx--;
     236   764425805 :   if (!nx) return utoi(s);
     237   719714488 :   return adduispec(s,xd,nx);
     238             : }
     239             : 
     240             : INLINE GEN
     241  3332638220 : addiispec(GEN x, GEN y, long nx, long ny)
     242             : {
     243             :   GEN zd;
     244             :   long lz;
     245             : 
     246  3332638220 :   if (nx < ny) swapspec(x,y, nx,ny);
     247  3332638220 :   if (ny == 1) return adduispec(*y,x,nx);
     248  1355814126 :   lz = nx+3; zd = cgeti(lz);
     249             : 
     250  1355897165 :   if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
     251    30430363 :     zd[lz-1]=1;
     252             :   else
     253  1325748106 :     lz--;
     254             : 
     255  1356178469 :   zd[1] = evalsigne(1) | evallgefint(lz);
     256  1356178469 :   return zd;
     257             : }
     258             : 
     259             : /* assume x >= y */
     260             : INLINE GEN
     261  1774416858 : subiuspec(GEN x, ulong s, long nx)
     262             : {
     263             :   GEN zd;
     264             :   long lz;
     265             : 
     266  1774416858 :   if (nx == 1) return utoi(x[0] - s);
     267             : 
     268   251813601 :   lz = nx + 2; zd = cgeti(lz);
     269   252293708 :   mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
     270   252294328 :   if (! zd[lz - 1]) { --lz; }
     271             : 
     272   252294328 :   zd[1] = evalsigne(1) | evallgefint(lz);
     273   252294328 :   return zd;
     274             : }
     275             : 
     276             : /* assume x > y */
     277             : INLINE GEN
     278  3019012421 : subiispec(GEN x, GEN y, long nx, long ny)
     279             : {
     280             :   GEN zd;
     281             :   long lz;
     282  3019012421 :   if (ny==1) return subiuspec(x,*y,nx);
     283  1274789543 :   lz = nx+2; zd = cgeti(lz);
     284             : 
     285  1272683510 :   mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
     286  1761115911 :   while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
     287             : 
     288  1272944419 :   zd[1] = evalsigne(1) | evallgefint(lz);
     289  1272944419 :   return zd;
     290             : }
     291             : 
     292             : static void
     293   527272399 : roundr_up_ip(GEN x, long l)
     294             : {
     295   527272399 :   long i = l;
     296             :   for(;;)
     297             :   {
     298   529067706 :     if (++((ulong*)x)[--i]) break;
     299     2190298 :     if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
     300             :   }
     301   527320324 : }
     302             : 
     303             : void
     304   403965023 : affir(GEN x, GEN y)
     305             : {
     306   403965023 :   const long s = signe(x), ly = lg(y);
     307             :   long lx, sh, i;
     308             : 
     309   403965023 :   if (!s)
     310             :   {
     311    42152147 :     y[1] = evalexpo(-bit_accuracy(ly));
     312    42148106 :     return;
     313             :   }
     314   361812876 :   lx = lgefint(x); sh = bfffo(*int_MSW(x));
     315   361812876 :   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
     316   362023884 :   if (sh) {
     317   354369849 :     if (lx <= ly)
     318             :     {
     319   961785974 :       for (i=lx; i<ly; i++) y[i]=0;
     320   249391126 :       mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
     321   249460471 :       xmpn_mirror(LIMBS(y),lx-2);
     322   249729635 :       return;
     323             :     }
     324   104978723 :     mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
     325   104979516 :     uel(y,2) |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
     326   104979516 :     xmpn_mirror(LIMBS(y),ly-2);
     327             :     /* lx > ly: round properly */
     328   104977023 :     if ((uel(x,lx-ly+1)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
     329             :   }
     330             :   else {
     331     7654035 :     GEN xd=int_MSW(x);
     332     7654035 :     if (lx <= ly)
     333             :     {
     334     9516859 :       for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
     335     5236357 :       for (   ; i<ly; i++) y[i]=0;
     336     2363422 :       return;
     337             :     }
     338    15121651 :     for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
     339             :     /* lx > ly: round properly */
     340     5290613 :     if (uel(x,lx-ly+1) & HIGHBIT) roundr_up_ip(y, ly);
     341             :   }
     342             : }
     343             : 
     344             : INLINE GEN
     345   715156434 : shiftispec(GEN x, long nx, long n)
     346             : {
     347             :   long ny,m;
     348             :   GEN yd, y;
     349             : 
     350   715156434 :   if (!n) return icopyspec(x, nx);
     351             : 
     352   685863458 :   if (n > 0)
     353             :   {
     354   399617224 :     long d = dvmdsBIL(n, &m);
     355             :     long i;
     356             : 
     357   399523267 :     ny = nx + d + (m!=0);
     358   399523267 :     y = cgeti(ny + 2); yd = y + 2;
     359   553248262 :     for (i=0; i<d; i++) yd[i] = 0;
     360             : 
     361   398595744 :     if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
     362             :     else
     363             :     {
     364   397054073 :       ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
     365   397113208 :       if (carryd) yd[ny - 1] = carryd;
     366   373831259 :       else ny--;
     367             :     }
     368             :   }
     369             :   else
     370             :   {
     371   286246234 :     long d = dvmdsBIL(-n, &m);
     372             : 
     373   288164797 :     ny = nx - d;
     374   288164797 :     if (ny < 1) return gen_0;
     375   285156557 :     y = cgeti(ny + 2); yd = y + 2;
     376             : 
     377   284838431 :     if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
     378             :     else
     379             :     {
     380   280174736 :       mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
     381   280188443 :       if (yd[ny - 1] == 0)
     382             :       {
     383    61013406 :         if (ny == 1) return gc_const((pari_sp)(yd + 1), gen_0);
     384    51398242 :         ny--;
     385             :       }
     386             :     }
     387             :   }
     388   672050691 :   y[1] = evalsigne(1)|evallgefint(ny + 2);
     389   672050691 :   return y;
     390             : }
     391             : 
     392             : GEN
     393   138203462 : mantissa2nr(GEN x, long n)
     394             : {
     395   138203462 :   long ly, i, m, s = signe(x), lx = lg(x);
     396             :   GEN y;
     397   138203462 :   if (!s) return gen_0;
     398   138202113 :   if (!n)
     399             :   {
     400    30690684 :     y = cgeti(lx);
     401    30687432 :     y[1] = evalsigne(s) | evallgefint(lx);
     402    30687432 :     xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
     403    30687188 :     return y;
     404             :   }
     405   107511429 :   if (n > 0)
     406             :   {
     407      218441 :     GEN z = (GEN)avma;
     408      218441 :     long d = dvmdsBIL(n, &m);
     409             : 
     410      218442 :     ly = lx+d; y = new_chunk(ly);
     411      557570 :     for ( ; d; d--) *--z = 0;
     412      222873 :     if (!m) for (i=2; i<lx; i++) y[i]=x[i];
     413             :     else
     414             :     {
     415      217002 :       const ulong sh = BITS_IN_LONG - m;
     416      217002 :       shift_left(y,x, 2,lx-1, 0,m);
     417      217002 :       i = uel(x,2) >> sh;
     418             :       /* Extend y on the left? */
     419      217002 :       if (i) { ly++; y = new_chunk(1); y[2] = i; }
     420             :     }
     421             :   }
     422             :   else
     423             :   {
     424   107292988 :     ly = lx - dvmdsBIL(-n, &m);
     425   107289960 :     if (ly<3) return gen_0;
     426   107289960 :     y = new_chunk(ly);
     427   107251295 :     if (m) {
     428   106992879 :       shift_right(y,x, 2,ly, 0,m);
     429   107030400 :       if (y[2] == 0)
     430             :       {
     431           0 :         if (ly==3) return gc_const((pari_sp)(y+3), gen_0);
     432           0 :         ly--; set_avma((pari_sp)(++y));
     433             :       }
     434             :     } else {
     435      701560 :       for (i=2; i<ly; i++) y[i]=x[i];
     436             :     }
     437             :   }
     438   107507258 :   xmpn_mirror(LIMBS(y),ly-2);
     439   107555993 :   y[1] = evalsigne(s)|evallgefint(ly);
     440   107555993 :   y[0] = evaltyp(t_INT)|evallg(ly); return y;
     441             : }
     442             : 
     443             : GEN
     444     3649588 : truncr(GEN x)
     445             : {
     446             :   long s, e, d, m, i;
     447             :   GEN y;
     448     3649588 :   if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
     449     1533410 :   d = nbits2lg(e+1); m = remsBIL(e);
     450     1533395 :   if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
     451             : 
     452     1533391 :   y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
     453     1533325 :   if (++m == BITS_IN_LONG)
     454       95989 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     455             :   else
     456             :   {
     457     1485589 :     GEN z=cgeti(d);
     458     3037803 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     459     1485503 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     460     1485503 :     set_avma((pari_sp)y);
     461             :   }
     462     1533221 :   return y;
     463             : }
     464             : 
     465             : /* integral part */
     466             : GEN
     467     7137395 : floorr(GEN x)
     468             : {
     469             :   long e, d, m, i, lx;
     470             :   GEN y;
     471     7137395 :   if (signe(x) >= 0) return truncr(x);
     472     4245868 :   if ((e=expo(x)) < 0) return gen_m1;
     473     3525683 :   d = nbits2lg(e+1); m = remsBIL(e);
     474     3524118 :   lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
     475     3524118 :   y = cgeti(d+1);
     476     3513612 :   if (++m == BITS_IN_LONG)
     477             :   {
     478        3053 :     for (i=2; i<d; i++) y[d-i+1]=x[i];
     479        1450 :     i=d; while (i<lx && !x[i]) i++;
     480         906 :     if (i==lx) goto END;
     481             :   }
     482             :   else
     483             :   {
     484     3512706 :     GEN z=cgeti(d);
     485     7855216 :     for (i=2; i<d; i++) z[d-i+1]=x[i];
     486     3499482 :     mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
     487     3500078 :     if (uel(x,d-1)<<m == 0)
     488             :     {
     489      521540 :       i=d; while (i<lx && !x[i]) i++;
     490      121321 :       if (i==lx) goto END;
     491             :     }
     492             :   }
     493     3423065 :   if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
     494           0 :     y[d++]=1;
     495     3422862 : END:
     496     3500781 :   y[1] = evalsigne(-1) | evallgefint(d);
     497     3500781 :   return y;
     498             : }
     499             : 
     500             : INLINE int
     501  3976437370 : cmpiispec(GEN x, GEN y, long lx, long ly)
     502             : {
     503  3976437370 :   if (lx < ly) return -1;
     504  3653704482 :   if (lx > ly) return  1;
     505  3504584417 :   return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     506             : }
     507             : 
     508             : INLINE int
     509   269590177 : equaliispec(GEN x, GEN y, long lx, long ly)
     510             : {
     511   269590177 :   if (lx != ly) return 0;
     512   269448929 :   return !mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
     513             : }
     514             : 
     515             : /***********************************************************************/
     516             : /**                                                                   **/
     517             : /**                          MULTIPLICATION                           **/
     518             : /**                                                                   **/
     519             : /***********************************************************************/
     520             : /* assume ny > 0 */
     521             : INLINE GEN
     522  5604865397 : muluispec(ulong x, GEN y, long ny)
     523             : {
     524  5604865397 :   if (ny == 1)
     525  4507074680 :     return muluu(x, *y);
     526             :   else
     527             :   {
     528  1097790717 :     long lz = ny+3;
     529  1097790717 :     GEN z = cgeti(lz);
     530  1107107086 :     ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
     531  1108479340 :     if (hi) { z[lz - 1] = hi; } else lz--;
     532  1108479340 :     z[1] = evalsigne(1) | evallgefint(lz);
     533  1108479340 :     return z;
     534             :   }
     535             : }
     536             : 
     537             : /* a + b*|y| */
     538             : GEN
     539           0 : addumului(ulong a, ulong b, GEN y)
     540             : {
     541             :   GEN z;
     542             :   long i, lz;
     543             :   ulong hi;
     544           0 :   if (!b || !signe(y)) return utoi(a);
     545           0 :   lz = lgefint(y)+1;
     546           0 :   z = cgeti(lz);
     547           0 :   z[2]=a;
     548           0 :   for(i=3;i<lz;i++) z[i]=0;
     549           0 :   hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
     550           0 :   if (hi) z[lz-1]=hi; else lz--;
     551           0 :   z[1] = evalsigne(1) | evallgefint(lz);
     552           0 :   return gc_const((pari_sp)z, z);
     553             : }
     554             : 
     555             : /***********************************************************************/
     556             : /**                                                                   **/
     557             : /**                          DIVISION                                 **/
     558             : /**                                                                   **/
     559             : /***********************************************************************/
     560             : 
     561             : ulong
     562  1324895234 : umodiu(GEN y, ulong x)
     563             : {
     564  1324895234 :   long sy=signe(y);
     565             :   ulong hi;
     566  1324895234 :   if (!x) pari_err_INV("umodiu",gen_0);
     567  1325848334 :   if (!sy) return 0;
     568   914545050 :   hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
     569   914545050 :   if (!hi) return 0;
     570   897946938 :   return (sy > 0)? hi: x - hi;
     571             : }
     572             : 
     573             : /* return |y| \/ x */
     574             : GEN
     575   110609613 : absdiviu_rem(GEN y, ulong x, ulong *rem)
     576             : {
     577             :   long ly;
     578             :   GEN z;
     579             : 
     580   110609613 :   if (!x) pari_err_INV("absdiviu_rem",gen_0);
     581   110615702 :   if (!signe(y)) { *rem = 0; return gen_0; }
     582             : 
     583    83682579 :   ly = lgefint(y);
     584    83682579 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
     585             : 
     586    68980877 :   z = cgeti(ly);
     587    68980526 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     588    68980759 :   if (z [ly - 1] == 0) ly--;
     589    68980759 :   z[1] = evallgefint(ly) | evalsigne(1);
     590    68980759 :   return z;
     591             : }
     592             : 
     593             : GEN
     594    86196561 : divis_rem(GEN y, long x, long *rem)
     595             : {
     596    86196561 :   long sy=signe(y),ly,s;
     597             :   GEN z;
     598             : 
     599    86196561 :   if (!x) pari_err_INV("divis_rem",gen_0);
     600    86207033 :   if (!sy) { *rem = 0; return gen_0; }
     601    60367167 :   if (x<0) { s = -sy; x = -x; } else s = sy;
     602             : 
     603    60367167 :   ly = lgefint(y);
     604    60367167 :   if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
     605             : 
     606    41575171 :   z = cgeti(ly);
     607    41570634 :   *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     608    41570721 :   if (sy<0) *rem = -  *rem;
     609    41570721 :   if (z[ly - 1] == 0) ly--;
     610    41570721 :   z[1] = evallgefint(ly) | evalsigne(s);
     611    41570721 :   return z;
     612             : }
     613             : 
     614             : GEN
     615     1006043 : divis(GEN y, long x)
     616             : {
     617     1006043 :   long sy=signe(y),ly,s;
     618             :   GEN z;
     619             : 
     620     1006043 :   if (!x) pari_err_INV("divis",gen_0);
     621     1006043 :   if (!sy) return gen_0;
     622     1005995 :   if (x<0) { s = -sy; x = -x; } else s=sy;
     623             : 
     624     1005995 :   ly = lgefint(y);
     625     1005995 :   if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
     626             : 
     627     1005683 :   z = cgeti(ly);
     628     1005682 :   (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
     629     1005682 :   if (z[ly - 1] == 0) ly--;
     630     1005682 :   z[1] = evallgefint(ly) | evalsigne(s);
     631     1005682 :   return z;
     632             : }
     633             : 
     634             : /* We keep llx bits of x and lly bits of y*/
     635             : static GEN
     636    76154676 : divrr_with_gmp(GEN x, GEN y)
     637             : {
     638    76154676 :   long lx=RNLIMBS(x),ly=RNLIMBS(y);
     639    76154676 :   long lw=minss(lx,ly);
     640    76155547 :   long lly=minss(lw+1,ly);
     641    76156711 :   GEN  w = cgetg(lw+2, t_REAL);
     642    76139333 :   long lu=lw+lly;
     643    76139333 :   long llx=minss(lu,lx);
     644    76136784 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     645    76116242 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     646             :   mp_limb_t *q,*r;
     647    76083631 :   long e=expo(x)-expo(y);
     648    76083631 :   long sx=signe(x),sy=signe(y);
     649    76083631 :   xmpn_mirrorcopy(z,RLIMBS(y),lly);
     650    76097204 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     651    76144156 :   xmpn_zero(u,lu-llx);
     652    76172522 :   q = (mp_limb_t *)new_chunk(lw+1);
     653    76155517 :   r = (mp_limb_t *)new_chunk(lly);
     654             : 
     655    76131720 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     656             : 
     657             :   /*Round up: This is not exactly correct we should test 2*r>z*/
     658    76180447 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     659    37780020 :     mpn_add_1(q,q,lw+1,1);
     660             : 
     661    76180484 :   xmpn_mirrorcopy(RLIMBS(w),q,lw);
     662             : 
     663    76179203 :   if (q[lw] == 0) e--;
     664    42026566 :   else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
     665           0 :   else { w[2] = HIGHBIT; e++; }
     666    76178164 :   if (sy < 0) sx = -sx;
     667    76178164 :   w[1] = evalsigne(sx) | evalexpo(e);
     668    76170783 :   return gc_const((pari_sp)w, w);
     669             : }
     670             : 
     671             : /* We keep llx bits of x and lly bits of y*/
     672             : static GEN
     673    35250418 : divri_with_gmp(GEN x, GEN y)
     674             : {
     675    35250418 :   long llx=RNLIMBS(x),ly=NLIMBS(y);
     676    35250418 :   long lly=minss(llx+1,ly);
     677    35250693 :   GEN  w = cgetg(llx+2, t_REAL);
     678    35246435 :   long lu=llx+lly, ld=ly-lly;
     679    35246435 :   mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
     680    35242812 :   mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
     681             :   mp_limb_t *q,*r;
     682    35236885 :   long sh=bfffo(y[ly+1]);
     683    35236885 :   long e=expo(x)-expi(y);
     684    35237641 :   long sx=signe(x),sy=signe(y);
     685    35237641 :   if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
     686      685168 :   else xmpn_copy(z,LIMBS(y)+ld,lly);
     687    35239401 :   xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
     688    35247013 :   xmpn_zero(u,lu-llx);
     689    35253280 :   q = (mp_limb_t *)new_chunk(llx+1);
     690    35250647 :   r = (mp_limb_t *)new_chunk(lly);
     691             : 
     692    35245921 :   mpn_tdiv_qr(q,r,0,u,lu,z,lly);
     693             : 
     694             :   /*Round up: This is not exactly correct we should test 2*r>z*/
     695    35255811 :   if (uel(r,lly-1) > (uel(z,lly-1)>>1))
     696    16045013 :     mpn_add_1(q,q,llx+1,1);
     697             : 
     698    35255819 :   xmpn_mirrorcopy(RLIMBS(w),q,llx);
     699             : 
     700    35255485 :   if (q[llx] == 0) e--;
     701    10627954 :   else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
     702           0 :   else { w[2] = HIGHBIT; e++; }
     703    35255337 :   if (sy < 0) sx = -sx;
     704    35255337 :   w[1] = evalsigne(sx) | evalexpo(e);
     705    35254098 :   return gc_const((pari_sp)w, w);
     706             : }
     707             : 
     708             : GEN
     709   151402426 : divri(GEN x, GEN y)
     710             : {
     711   151402426 :   long  s = signe(y);
     712             : 
     713   151402426 :   if (!s) pari_err_INV("divri",gen_0);
     714   151402830 :   if (!signe(x)) return real_0_bit(expo(x) - expi(y));
     715   151166258 :   if (!is_bigint(y)) {
     716   115921313 :     GEN z = divru(x, y[2]);
     717   115919291 :     if (s < 0) togglesign(z);
     718   115919337 :     return z;
     719             :   }
     720    35243857 :   return divri_with_gmp(x,y);
     721             : }
     722             : 
     723             : GEN
     724   142802202 : divrr(GEN x, GEN y)
     725             : {
     726   142802202 :   long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
     727             :   ulong y0,y1;
     728             :   GEN r, r1;
     729             : 
     730   142802202 :   if (!sy) pari_err_INV("divrr",y);
     731   142810051 :   e = expo(x) - expo(y);
     732   142810051 :   if (!sx) return real_0_bit(e);
     733   142303905 :   if (sy<0) sx = -sx;
     734             : 
     735   142303905 :   lx=lg(x); ly=lg(y);
     736   142303905 :   if (ly==3)
     737             :   {
     738    26342191 :     ulong k = x[2], l = (lx>3)? x[3]: 0;
     739             :     LOCAL_HIREMAINDER;
     740    26342191 :     if (k < uel(y,2)) e--;
     741             :     else
     742             :     {
     743     8290171 :       l >>= 1; if (k&1) l |= HIGHBIT;
     744     8290171 :       k >>= 1;
     745             :     }
     746    26342191 :     hiremainder = k; k = divll(l,y[2]);
     747    26342191 :     if (hiremainder > (uel(y,2) >> 1) && !++k) { k = HIGHBIT; e++; }
     748    26342191 :     r = cgetg(3, t_REAL);
     749    26340910 :     r[1] = evalsigne(sx) | evalexpo(e);
     750    26338695 :     r[2] = k; return r;
     751             :   }
     752             : 
     753   115961714 :   if (ly >= prec2lg(DIVRR_GMP_LIMIT))
     754    76154355 :     return divrr_with_gmp(x,y);
     755             : 
     756    39854711 :   lr = minss(lx,ly); r = new_chunk(lr);
     757    39856538 :   r1 = r-1;
     758   135495937 :   r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
     759    39856538 :   r1[lr] = (lx>ly)? x[lr]: 0;
     760    39856538 :   y0 = y[2]; y1 = y[3];
     761   175314755 :   for (i=0; i<lr-1; i++)
     762             :   { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
     763             :     ulong k, qp;
     764             :     LOCAL_HIREMAINDER;
     765             :     LOCAL_OVERFLOW;
     766             : 
     767   135458217 :     if (uel(r1,1) == y0) { qp = ULONG_MAX; k = addll(y0,r1[2]); }
     768             :     else
     769             :     {
     770   135205215 :       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   135205215 :       hiremainder = r1[1]; overflow = 0;
     778   135205215 :       qp = divll(r1[2],y0); k = hiremainder;
     779             :     }
     780   135458217 :     j = lr-i+1;
     781   135458217 :     if (!overflow)
     782             :     {
     783             :       long k3, k4;
     784   135300271 :       k3 = mulll(qp,y1);
     785   135300271 :       if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
     786    39909888 :         k4 = subll(hiremainder,k);
     787             :       else
     788             :       {
     789    95390383 :         k3 = subll(k3, r1[3]);
     790    95390383 :         k4 = subllx(hiremainder,k);
     791             :       }
     792   173712642 :       while (!overflow && k4) { qp--; k3=subll(k3,y1); k4=subllx(k4,y0); }
     793             :     }
     794   135458217 :     if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
     795   402146268 :     for (j--; j>1; j--)
     796             :     {
     797   266688051 :       r1[j] = subll(r1[j], addmul(qp,y[j]));
     798   266688051 :       hiremainder += overflow;
     799             :     }
     800   135458217 :     if (uel(r1,1) != hiremainder)
     801             :     {
     802      187478 :       if (uel(r1,1) < hiremainder)
     803             :       {
     804      187478 :         qp--;
     805      187478 :         j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
     806      534281 :         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   135458217 :     *++r1 = qp;
     821             :   }
     822             :   /* i = lr-1 */
     823             :   /* round correctly */
     824    39856538 :   if (uel(r1,1) > (y0>>1))
     825             :   {
     826    20134089 :     j=i; do uel(r,--j)++; while (j && !r[j]);
     827             :   }
     828   135679040 :   r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
     829    39856538 :   if (r[0] == 0) e--;
     830    14281576 :   else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
     831             :   else { /* possible only when rounding up to 0x2 0x0 ... */
     832          18 :     r[2] = (long)HIGHBIT; e++;
     833             :   }
     834    39855696 :   r[0] = evaltyp(t_REAL)|evallg(lr);
     835    39914128 :   r[1] = evalsigne(sx) | evalexpo(e);
     836    39907206 :   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 can be disposed of with cgiv
     843             :  * If *z is zero, we put gen_0 here and no copy.
     844             :  * space needed: lx + ly
     845             :  */
     846             : GEN
     847  2196026897 : dvmdii(GEN x, GEN y, GEN *z)
     848             : {
     849  2196026897 :   long sx=signe(x),sy=signe(y);
     850             :   long lx, ly, lq;
     851             :   pari_sp av;
     852             :   GEN r,q;
     853             : 
     854  2196026897 :   if (!sy) pari_err_INV("dvmdii",y);
     855  2197120880 :   if (!sx)
     856             :   {
     857    69919530 :     if (!z || z == ONLY_REM) return gen_0;
     858    41102861 :     *z=gen_0; return gen_0;
     859             :   }
     860  2127201350 :   lx=lgefint(x);
     861  2127201350 :   ly=lgefint(y); lq=lx-ly;
     862  2127201350 :   if (lq <= 0)
     863             :   {
     864  1260682015 :     if (lq == 0)
     865             :     {
     866  1071849645 :       long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
     867  1071849645 :       if (s>0) goto DIVIDE;
     868   369581737 :       if (s==0)
     869             :       {
     870    33439563 :         if (z == ONLY_REM) return gen_0;
     871    22327706 :         if (z) *z = gen_0;
     872    22327706 :         if (sx < 0) sy = -sy;
     873    22327706 :         return stoi(sy);
     874             :       }
     875             :     }
     876   524974544 :     if (z == ONLY_REM) return icopy(x);
     877    14332338 :     if (z) *z = icopy(x);
     878    14332338 :     return gen_0;
     879             :   }
     880   866519335 : DIVIDE: /* quotient is nonzero */
     881  1568787243 :   av=avma; if (sx<0) sy = -sy;
     882  1568787243 :   if (ly==3)
     883             :   {
     884   574061258 :     ulong lq = lx;
     885             :     ulong si;
     886   574061258 :     q  = cgeti(lq);
     887   573615232 :     si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
     888   574140128 :     if (q[lq - 1] == 0) lq--;
     889   574140128 :     if (z == ONLY_REM)
     890             :     {
     891   333962548 :       if (!si) return gc_const(av, gen_0);
     892   287638082 :       set_avma(av); r = cgeti(3);
     893   287099839 :       r[1] = evalsigne(sx) | evallgefint(3);
     894   287099839 :       r[2] = si; return r;
     895             :     }
     896   240177580 :     q[1] = evalsigne(sy) | evallgefint(lq);
     897   240177580 :     if (!z) return q;
     898   236064192 :     if (!si) { *z=gen_0; return q; }
     899    64202540 :     r=cgeti(3);
     900    64229916 :     r[1] = evalsigne(sx) | evallgefint(3);
     901    64229916 :     r[2] = si; *z=r; return q;
     902             :   }
     903   994725985 :   if (z == ONLY_REM)
     904             :   {
     905   971573299 :     ulong lr = lgefint(y);
     906   971573299 :     ulong lq = lgefint(x)-lgefint(y)+3;
     907   971573299 :     GEN r = cgeti(lr);
     908   965105019 :     GEN q = cgeti(lq);
     909   958636061 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     910   972277172 :     if (!r[lr - 1])
     911             :     {
     912  1182253658 :       while(lr>2 && !r[lr - 1]) lr--;
     913   545231634 :       if (lr == 2) return gc_const(av, gen_0); /* exact division */
     914             :     }
     915   957337498 :     r[1] = evalsigne(sx) | evallgefint(lr);
     916   957337498 :     return gc_const((pari_sp)r, r);
     917             :   }
     918             :   else
     919             :   {
     920    23152686 :     ulong lq = lgefint(x)-lgefint(y)+3;
     921    23152686 :     ulong lr = lgefint(y);
     922    23152686 :     GEN q = cgeti(lq);
     923    28641665 :     GEN r = cgeti(lr);
     924    28627064 :     mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
     925    28653897 :     if (q[lq - 1] == 0) lq--;
     926    28653897 :     q[1] = evalsigne(sy) | evallgefint(lq);
     927    28653897 :     if (!z) return gc_const((pari_sp)q, q);
     928    24949961 :     if (!r[lr - 1])
     929             :     {
     930    38510623 :       while(lr>2 && !r[lr - 1]) lr--;
     931     6343219 :       if (lr == 2) { *z = gen_0; return gc_const((pari_sp)q, q); } /* exact */
     932             :     }
     933    20169050 :     r[1] = evalsigne(sx) | evallgefint(lr);
     934    20169050 :     *z = gc_const((pari_sp)r, r); return q;
     935             :   }
     936             : }
     937             : 
     938             : /* Montgomery reduction.
     939             :  * N has k words, assume T >= 0 has less than 2k.
     940             :  * Return res := T / B^k mod N, where B = 2^BIL
     941             :  * such that 0 <= res < T/B^k + N  and  res has less than k words */
     942             : GEN
     943    35025180 : red_montgomery(GEN T, GEN N, ulong inv)
     944             : {
     945             :   pari_sp av;
     946             :   GEN Te, Td, Ne, Nd, scratch;
     947    35025180 :   ulong i, j, m, t, d, k = NLIMBS(N);
     948             :   int carry;
     949             :   LOCAL_HIREMAINDER;
     950             :   LOCAL_OVERFLOW;
     951             : 
     952    35025180 :   if (k == 0) return gen_0;
     953    35025180 :   d = NLIMBS(T); /* <= 2*k */
     954    35025180 :   if (d == 0) return gen_0;
     955             : #ifdef DEBUG
     956             :   if (d > 2*k) pari_err_BUG("red_montgomery");
     957             : #endif
     958    35025171 :   if (k == 1)
     959             :   { /* as below, special cased for efficiency */
     960      163341 :     ulong n = uel(N,2);
     961      163341 :     if (d == 1) {
     962      163194 :       hiremainder = uel(T,2);
     963      163194 :       m = hiremainder * inv;
     964      163194 :       (void)addmul(m, n); /* t + m*n = 0 */
     965      163194 :       return utoi(hiremainder);
     966             :     } else { /* d = 2 */
     967         147 :       hiremainder = uel(T,2);
     968         147 :       m = hiremainder * inv;
     969         147 :       (void)addmul(m, n); /* t + m*n = 0 */
     970         147 :       t = addll(hiremainder, uel(T,3));
     971         147 :       if (overflow) t -= n; /* t > n doesn't fit in 1 word */
     972         147 :       return utoi(t);
     973             :     }
     974             :   }
     975             :   /* assume k >= 2 */
     976    34861830 :   av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
     977             : 
     978             :   /* copy T to scratch space (pad with zeroes to 2k words) */
     979    34457974 :   Td = scratch;
     980    34457974 :   Te = T + 2;
     981   482588050 :   for (i=0; i < d     ; i++) *Td++ = *Te++;
     982    61629692 :   for (   ; i < (k<<1); i++) *Td++ = 0;
     983             : 
     984    34457974 :   Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
     985    34457974 :   Ne = N + 1;       /* 1 beyond end of N mantissa */
     986             : 
     987    34457974 :   carry = 0;
     988   257372207 :   for (i=0; i<k; i++) /* set T := T/B nod N, k times */
     989             :   {
     990   222914233 :     Td = Te; /* one beyond end of (new) T mantissa */
     991   222914233 :     Nd = Ne;
     992   222914233 :     hiremainder = *++Td;
     993   222914233 :     m = hiremainder * inv; /* solve T + m N = O(B) */
     994             : 
     995             :     /* set T := (T + mN) / B */
     996   222914233 :     Te = Td;
     997   222914233 :     (void)addmul(m, *++Nd); /* = 0 */
     998  1934927276 :     for (j=1; j<k; j++)
     999             :     {
    1000  1712013043 :       t = addll(addmul(m, *++Nd), *++Td);
    1001  1712013043 :       *Td = t;
    1002  1712013043 :       hiremainder += overflow;
    1003             :     }
    1004   222914233 :     t = addll(hiremainder, *++Td); *Td = t + carry;
    1005   222914233 :     carry = (overflow || (carry && *Td == 0));
    1006             :   }
    1007    34457974 :   if (carry)
    1008             :   { /* Td > N overflows (k+1 words), set Td := Td - N */
    1009       75961 :     GEN NE = N + k+1;
    1010       75961 :     Td = Te;
    1011       75961 :     Nd = Ne;
    1012       75961 :     t = subll(*++Td, *++Nd); *Td = t;
    1013      764666 :     while (Nd < NE) { t = subllx(*++Td, *++Nd); *Td = t; }
    1014             :   }
    1015             : 
    1016             :   /* copy result */
    1017    34457974 :   Td = (GEN)av - 1; /* *Td = high word of final result */
    1018    38066542 :   while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
    1019    34457974 :   k = Td - Te; if (!k) return gc_const(av, gen_0);
    1020    34457974 :   Td = (GEN)av - k; /* will write mantissa there */
    1021    34457974 :   (void)memmove(Td, Te+1, k*sizeof(long));
    1022    34457974 :   Td -= 2;
    1023    34457974 :   Td[0] = evaltyp(t_INT) | evallg(k+2);
    1024    34779064 :   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    34779064 :   return gc_const((pari_sp)Td, Td);
    1036             : }
    1037             : 
    1038             : /* EXACT INTEGER DIVISION */
    1039             : 
    1040             : /* use undocumented GMP interface */
    1041             : static void
    1042   115382511 : GEN2mpz(mpz_t X, GEN x)
    1043             : {
    1044   115382511 :   long l = lgefint(x)-2;
    1045   115382511 :   X->_mp_alloc = l;
    1046   115382511 :   X->_mp_size = signe(x) > 0? l: -l;
    1047   115382511 :   X->_mp_d = LIMBS(x);
    1048   115382511 : }
    1049             : static void
    1050    57692615 : mpz2GEN(GEN z, mpz_t Z)
    1051             : {
    1052    57692615 :   long l = Z->_mp_size;
    1053    57692615 :   z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
    1054    57692615 : }
    1055             : 
    1056             : #ifdef mpn_divexact_1
    1057             : static GEN
    1058   419123013 : diviuexact_i(GEN x, ulong y)
    1059             : {
    1060   419123013 :   long l = lgefint(x);
    1061   419123013 :   GEN z = cgeti(l);
    1062   418648841 :   mpn_divexact_1(LIMBS(z), LIMBS(x), NLIMBS(x), y);
    1063   418693404 :   if (z[l-1] == 0) l--;
    1064   418693404 :   z[1] = evallgefint(l) | evalsigne(signe(x));
    1065   418693404 :   return z;
    1066             : }
    1067             : #elif 1 && !defined(_WIN64) /* mpz_divexact_ui is not LLP64 friendly   */
    1068             :                             /* assume y != 0 and the division is exact */
    1069             : static GEN
    1070             : diviuexact_i(GEN x, ulong y)
    1071             : {
    1072             :   long l = lgefint(x);
    1073             :   GEN z = cgeti(l);
    1074             :   mpz_t X, Z;
    1075             :   GEN2mpz(X, x);
    1076             :   Z->_mp_alloc = l-2;
    1077             :   Z->_mp_size  = l-2;
    1078             :   Z->_mp_d = LIMBS(z);
    1079             :   mpz_divexact_ui(Z, X, y);
    1080             :   mpz2GEN(z, Z); return z;
    1081             : }
    1082             : #else
    1083             : /* assume y != 0 and the division is exact */
    1084             : static GEN
    1085             : diviuexact_i(GEN x, ulong y)
    1086             : {
    1087             :   /*TODO: implement true exact division.*/
    1088             :   return divii(x,utoi(y));
    1089             : }
    1090             : #endif
    1091             : 
    1092             : GEN
    1093    39544437 : diviuexact(GEN x, ulong y)
    1094             : {
    1095             :   GEN z;
    1096    39544437 :   if (!signe(x)) return gen_0;
    1097    38415861 :   z = diviuexact_i(x, y);
    1098    38412397 :   if (lgefint(z) == 2) pari_err_OP("exact division", x, utoi(y));
    1099    38412447 :   return z;
    1100             : }
    1101             : 
    1102             : /* Find z such that x=y*z, knowing that y | x (unchecked) */
    1103             : GEN
    1104   531624323 : diviiexact(GEN x, GEN y)
    1105             : {
    1106             :   GEN z;
    1107   531624323 :   if (!signe(y)) pari_err_INV("diviiexact",y);
    1108   531796115 :   if (!signe(x)) return gen_0;
    1109   438589457 :   if (lgefint(y) == 3)
    1110             :   {
    1111   380913855 :     z = diviuexact_i(x, y[2]);
    1112   380479799 :     if (signe(y) < 0) togglesign(z);
    1113             :   }
    1114             :   else
    1115             :   {
    1116    57675602 :     long l = lgefint(x);
    1117             :     mpz_t X, Y, Z;
    1118    57675602 :     z = cgeti(l);
    1119    57691404 :     GEN2mpz(X, x);
    1120    57691335 :     GEN2mpz(Y, y);
    1121    57691217 :     Z->_mp_alloc = l-2;
    1122    57691217 :     Z->_mp_size  = l-2;
    1123    57691217 :     Z->_mp_d = LIMBS(z);
    1124    57691217 :     mpz_divexact(Z, X, Y);
    1125    57692637 :     mpz2GEN(z, Z);
    1126             :   }
    1127   438172404 :   if (lgefint(z) == 2) pari_err_OP("exact division", x, y);
    1128   438147808 :   return z;
    1129             : }
    1130             : 
    1131             : /* assume yz != and yz | x */
    1132             : GEN
    1133      200495 : diviuuexact(GEN x, ulong y, ulong z)
    1134             : {
    1135             :   long tmp[4];
    1136             :   ulong t;
    1137             :   LOCAL_HIREMAINDER;
    1138      200495 :   t = mulll(y, z);
    1139      200495 :   if (!hiremainder) return diviuexact(x, t);
    1140           0 :   tmp[0] = evaltyp(t_INT)|_evallg(4);
    1141           0 :   tmp[1] = evalsigne(1)|evallgefint(4);
    1142           0 :   tmp[2] = t;
    1143           0 :   tmp[3] = hiremainder;
    1144           0 :   return diviiexact(x, tmp);
    1145             : }
    1146             : 
    1147             : /********************************************************************/
    1148             : /**                                                                **/
    1149             : /**               INTEGER MULTIPLICATION                           **/
    1150             : /**                                                                **/
    1151             : /********************************************************************/
    1152             : 
    1153             : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
    1154             : GEN
    1155  5940234325 : muliispec(GEN x, GEN y, long nx, long ny)
    1156             : {
    1157             :   GEN zd;
    1158             :   long lz;
    1159             :   ulong hi;
    1160             : 
    1161  5940234325 :   if (nx < ny) swapspec(x,y, nx,ny);
    1162  5940234325 :   if (!ny) return gen_0;
    1163  5940234325 :   if (ny == 1) return muluispec((ulong)*y, x, nx);
    1164             : 
    1165  1057578953 :   lz = nx+ny+2;
    1166  1057578953 :   zd = cgeti(lz);
    1167  1060488039 :   hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
    1168  1067676867 :   if (!hi) lz--;
    1169             :   /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
    1170             : 
    1171  1067676867 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1172  1067676867 :   return zd;
    1173             : }
    1174             : GEN
    1175      222731 : muluui(ulong x, ulong y, GEN z)
    1176             : {
    1177      222731 :   long t, s = signe(z);
    1178             :   GEN r;
    1179             :   LOCAL_HIREMAINDER;
    1180             : 
    1181      222731 :   if (!x || !y || !signe(z)) return gen_0;
    1182      222350 :   t = mulll(x,y);
    1183      222350 :   if (!hiremainder)
    1184      222357 :     r = muluispec(t, z+2, lgefint(z)-2);
    1185             :   else
    1186             :   {
    1187             :     long tmp[2];
    1188           0 :     tmp[1] = hiremainder;
    1189           0 :     tmp[0] = t;
    1190           0 :     r = muliispec(z+2,tmp, lgefint(z)-2, 2);
    1191             :   }
    1192      222323 :   setsigne(r,s); return r;
    1193             : }
    1194             : 
    1195             : GEN
    1196  1029912551 : sqrispec(GEN x, long nx)
    1197             : {
    1198             :   GEN zd;
    1199             :   long lz;
    1200             : 
    1201  1029912551 :   if (!nx) return gen_0;
    1202   492857097 :   if (nx==1) return sqru(*x);
    1203             : 
    1204   279012447 :   lz = (nx<<1)+2;
    1205   279012447 :   zd = cgeti(lz);
    1206             : #ifdef mpn_sqr
    1207   275783575 :   mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
    1208             : #else
    1209             :   mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
    1210             : #endif
    1211   279514922 :   if (zd[lz-1]==0) lz--;
    1212             : 
    1213   279514922 :   zd[1] = evalsigne(1) | evallgefint(lz);
    1214   279514922 :   return zd;
    1215             : }
    1216             : 
    1217             : INLINE GEN
    1218    41404768 : sqrispec_mirror(GEN x, long nx)
    1219             : {
    1220    41404768 :   GEN cx=new_chunk(nx);
    1221             :   GEN z;
    1222    41358629 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1223    41463516 :   z=sqrispec(cx, nx);
    1224    41536283 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1225    41536390 :   return z;
    1226             : }
    1227             : 
    1228             : /* leaves garbage on the stack. */
    1229             : INLINE GEN
    1230    85101315 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
    1231             : {
    1232             :   GEN cx, cy, z;
    1233    85101315 :   long s = 0;
    1234   116240350 :   while (nx && x[nx-1]==0) { nx--; s++; }
    1235   122857072 :   while (ny && y[ny-1]==0) { ny--; s++; }
    1236    85101315 :   cx=new_chunk(nx); cy=new_chunk(ny);
    1237    84644232 :   xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
    1238    85570670 :   xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
    1239    86125552 :   z =  nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
    1240    86114288 :   if (s)
    1241             :   {
    1242     7856748 :     long i, lz = lgefint(z) + s;
    1243     7856748 :     (void)new_chunk(s);
    1244     7856747 :     z -= s;
    1245    76751538 :     for (i=0; i<s; i++) z[2+i]=0;
    1246     7856747 :     z[1] = evalsigne(1) | evallgefint(lz);
    1247     7856747 :     z[0] = evaltyp(t_INT) | evallg(lz);
    1248             :   }
    1249    86114286 :   xmpn_mirror(LIMBS(z), NLIMBS(z));
    1250    86792119 :   return z;
    1251             : }
    1252             : 
    1253             : /* x % (2^n), assuming n >= 0 */
    1254             : GEN
    1255    39007509 : remi2n(GEN x, long n)
    1256             : {
    1257             :   ulong hi;
    1258    39007509 :   long l, k, lx, ly, sx = signe(x);
    1259             :   GEN z, xd, zd;
    1260             : 
    1261    39007509 :   if (!sx || !n) return gen_0;
    1262             : 
    1263    38685852 :   k = dvmdsBIL(n, &l);
    1264    38697044 :   lx = lgefint(x);
    1265    38697044 :   if (lx < k+3) return icopy(x);
    1266             : 
    1267    37750517 :   xd = x + (2 + k);
    1268             :   /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
    1269             :    *              ^--- initial xd  */
    1270    37750517 :   hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
    1271    37750517 :   if (!hi)
    1272             :   { /* strip leading zeroes from result */
    1273     4272013 :     xd--; while (k && !*xd) { k--; xd--; }
    1274     4072570 :     if (!k) return gen_0;
    1275     2190646 :     ly = k+2;
    1276             :   }
    1277             :   else
    1278    33677947 :     ly = k+3;
    1279             : 
    1280    35868593 :   zd = z = cgeti(ly);
    1281    35826712 :   *++zd = evalsigne(sx) | evallgefint(ly);
    1282   535875736 :   xd = x+1; for ( ;k; k--) *++zd = *++xd;
    1283    35826712 :   if (hi) *++zd = hi;
    1284    35826712 :   return z;
    1285             : }
    1286             : 
    1287             : /********************************************************************/
    1288             : /**                                                                **/
    1289             : /**                      INTEGER SQUARE ROOT                       **/
    1290             : /**                                                                **/
    1291             : /********************************************************************/
    1292             : 
    1293             : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
    1294             :  * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
    1295             :  * remainder is 0. R = NULL is allowed. */
    1296             : GEN
    1297     5182386 : sqrtremi(GEN a, GEN *r)
    1298             : {
    1299     5182386 :   long l, na = NLIMBS(a);
    1300             :   mp_size_t nr;
    1301             :   GEN S;
    1302     5182386 :   if (!na) {
    1303         724 :     if (r) *r = gen_0;
    1304         724 :     return gen_0;
    1305             :   }
    1306     5181662 :   l = (na + 5) >> 1; /* 2 + ceil(na/2) */
    1307     5181662 :   S = cgetipos(l);
    1308     5181622 :   if (r) {
    1309     1310193 :     GEN R = cgeti(2 + na);
    1310     1310193 :     nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
    1311     1310193 :     if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
    1312       26515 :     else    { set_avma((pari_sp)S); R = gen_0; }
    1313     1310193 :     *r = R;
    1314             :   }
    1315             :   else
    1316     3871429 :     (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
    1317     5181632 :   return S;
    1318             : }
    1319             : 
    1320             : /* compute sqrt(|a|), assuming a != 0 */
    1321             : GEN
    1322   126391400 : sqrtr_abs(GEN a)
    1323             : {
    1324             :   GEN res;
    1325             :   mp_limb_t *b, *c;
    1326   126391400 :   long l = RNLIMBS(a), e = expo(a), er = e>>1;
    1327             :   long n;
    1328   126391400 :   res = cgetg(2 + l, t_REAL);
    1329   126306369 :   res[1] = evalsigne(1) | evalexpo(er);
    1330   126383123 :   if (e&1)
    1331             :   {
    1332    53331449 :     b = (mp_limb_t *) new_chunk(l<<1);
    1333    53316607 :     xmpn_zero(b,l);
    1334    53317351 :     xmpn_mirrorcopy(b+l, RLIMBS(a), l);
    1335    53330996 :     c = (mp_limb_t *) new_chunk(l);
    1336    53324209 :     n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
    1337    53360359 :     if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
    1338             :   }
    1339             :   else
    1340             :   {
    1341             :     ulong u;
    1342    73051674 :     b = (mp_limb_t *) mantissa2nr(a,-1);
    1343    73086000 :     b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
    1344    73086000 :     b = (mp_limb_t *) new_chunk(l);
    1345    73062659 :     xmpn_zero(b,l+1); /* overwrites the former b[0] */
    1346    73065429 :     c = (mp_limb_t *) new_chunk(l + 1);
    1347    73026878 :     n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
    1348    73126242 :     u = (ulong)*c++;
    1349    73126242 :     if ( u&HIGHBIT || (u == ~HIGHBIT &&
    1350           0 :              (n>l || (n==l && mpn_cmp(b,c,l)>0))))
    1351    36026045 :       mpn_add_1(c,c,l,1);
    1352             :   }
    1353   126495658 :   xmpn_mirrorcopy(RLIMBS(res),c,l);
    1354   126470557 :   return gc_const((pari_sp)res, res);
    1355             : }
    1356             : 
    1357             : /* Normalize a nonnegative integer */
    1358             : GEN
    1359   307844128 : int_normalize(GEN x, long known_zero_words)
    1360             : {
    1361   307844128 :   long i =  lgefint(x) - 1 - known_zero_words;
    1362  2775400321 :   for ( ; i > 1; i--)
    1363  2723507060 :     if (x[i]) { setlgefint(x, i+1); return x; }
    1364    51893261 :   x[1] = evalsigne(0) | evallgefint(2); return x;
    1365             : }
    1366             : 
    1367             : /********************************************************************
    1368             :  **                                                                **
    1369             :  **                           Base Conversion                      **
    1370             :  **                                                                **
    1371             :  ********************************************************************/
    1372             : 
    1373             : ulong *
    1374      439099 : convi(GEN x, long *l)
    1375             : {
    1376      439099 :   long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
    1377      439099 :   GEN str = cgetg(n+1, t_VECSMALL);
    1378      439099 :   unsigned char *res = (unsigned char*) GSTR(str);
    1379      439099 :   long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
    1380             :   long lz;
    1381             :   ulong *z;
    1382             :   long i, j;
    1383             :   unsigned char *t;
    1384      439099 :   while (!*res) {res++; llz--;} /*Strip leading zeros*/
    1385      439099 :   lz  = (8+llz)/9;
    1386      439099 :   z = (ulong*)new_chunk(1+lz);
    1387      439099 :   t=res+llz+9;
    1388      870517 :   for(i=0;i<llz-8;i+=9)
    1389             :   {
    1390             :     ulong s;
    1391      431418 :     t-=18;
    1392      431418 :     s=*t++;
    1393     3882762 :     for (j=1; j<9;j++)
    1394     3451344 :       s=10*s+*t++;
    1395      431418 :     *z++=s;
    1396             :   }
    1397      439099 :   if (i<llz)
    1398             :   {
    1399      435116 :     unsigned char *t = res;
    1400      435116 :     ulong s=*t++;
    1401     1229599 :     for (j=i+1; j<llz;j++)
    1402      794483 :       s=10*s+*t++;
    1403      435116 :     *z++=s;
    1404             :   }
    1405      439099 :   *l = lz;
    1406      439099 :   return z;
    1407             : }

Generated by: LCOV version 1.16