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/none - mp_indep.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 29857-ff9b6c6cb1) Lines: 540 572 94.4 %
Date: 2025-01-05 09:10:50 Functions: 52 53 98.1 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : #line 2 "../src/kernel/none/mp_indep.c"
       2             : /* Copyright (C) 2000  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             : /* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */
      17             : ulong
      18    60489644 : invmod2BIL(ulong b)
      19             : {
      20             :   static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
      21    60489644 :   ulong x = b + tab[b & 7]; /* b^(-1) mod 2^4 */
      22             : 
      23             :   /* Newton applied to 1/x - b = 0 */
      24             : #ifdef LONG_IS_64BIT
      25    60489644 :   x = x*(2-b*x); /* one more pass necessary */
      26             : #endif
      27    60489644 :   x = x*(2-b*x);
      28    60489644 :   x = x*(2-b*x); return x*(2-b*x);
      29             : }
      30             : 
      31             : void
      32  1357224113 : affrr(GEN x, GEN y)
      33             : {
      34  1357224113 :   long i, lx, ly = lg(y);
      35  1357224113 :   if (!signe(x))
      36             :   {
      37      126249 :     y[1] = evalexpo(minss(expo(x), -bit_accuracy(ly)));
      38      126249 :     return;
      39             :   }
      40  1357097864 :   y[1] = x[1]; lx = lg(x);
      41  1357097864 :   if (lx <= ly)
      42             :   {
      43  7301758787 :     for (i=2; i<lx; i++) y[i]=x[i];
      44  1300741329 :     for (   ; i<ly; i++) y[i]=0;
      45  1144145704 :     return;
      46             :   }
      47  1140147071 :   for (i=2; i<ly; i++) y[i]=x[i];
      48             :   /* lx > ly: round properly */
      49   212952160 :   if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
      50             : }
      51             : 
      52             : GEN
      53    62466089 : trunc2nr(GEN x, long n)
      54             : {
      55             :   long ex;
      56    62466089 :   if (!signe(x)) return gen_0;
      57    62028237 :   ex = expo(x) + n; if (ex < 0) return gen_0;
      58    59479585 :   return mantissa2nr(x, ex - bit_prec(x) + 1);
      59             : }
      60             : 
      61             : /* x a t_REAL, x = i/2^e, i a t_INT */
      62             : GEN
      63    52623474 : mantissa_real(GEN x, long *e)
      64             : {
      65    52623474 :   *e = bit_prec(x)-1-expo(x);
      66    52623249 :   return mantissa2nr(x, 0);
      67             : }
      68             : 
      69             : GEN
      70   901241030 : mului(ulong x, GEN y)
      71             : {
      72   901241030 :   long s = signe(y);
      73             :   GEN z;
      74             : 
      75   901241030 :   if (!s || !x) return gen_0;
      76   748569990 :   z = muluispec(x, y+2, lgefint(y)-2);
      77   748257389 :   setsigne(z,s); return z;
      78             : }
      79             : 
      80             : GEN
      81   603635258 : mulsi(long x, GEN y)
      82             : {
      83   603635258 :   long s = signe(y);
      84             :   GEN z;
      85             : 
      86   603635258 :   if (!s || !x) return gen_0;
      87   360611205 :   if (x<0) { s = -s; x = -x; }
      88   360611205 :   z = muluispec((ulong)x, y+2, lgefint(y)-2);
      89   360607447 :   setsigne(z,s); return z;
      90             : }
      91             : 
      92             : GEN
      93   201098317 : mulss(long x, long y)
      94             : {
      95             :   ulong p1;
      96             :   LOCAL_HIREMAINDER;
      97             : 
      98   201098317 :   if (!x || !y) return gen_0;
      99   200763017 :   if (x<0) {
     100      414684 :     x = -x;
     101      414684 :     if (y<0) { y = -y; p1 = mulll(x,y); return uutoi(hiremainder, p1); }
     102      299827 :     p1 = mulll(x,y); return uutoineg(hiremainder, p1);
     103             :   } else {
     104   200348333 :     if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
     105   200341999 :     p1 = mulll(x,y); return uutoi(hiremainder, p1);
     106             :   }
     107             : }
     108             : GEN
     109     1497030 : sqrs(long x)
     110             : {
     111             :   ulong p1;
     112             :   LOCAL_HIREMAINDER;
     113             : 
     114     1497030 :   if (!x) return gen_0;
     115     1495259 :   if (x<0) x = -x;
     116     1495259 :   p1 = mulll(x,x); return uutoi(hiremainder, p1);
     117             : }
     118             : GEN
     119  4419846295 : muluu(ulong x, ulong y)
     120             : {
     121             :   ulong p1;
     122             :   LOCAL_HIREMAINDER;
     123             : 
     124  4419846295 :   if (!x || !y) return gen_0;
     125  4420721139 :   p1 = mulll(x,y); return uutoi(hiremainder, p1);
     126             : }
     127             : GEN
     128   605866286 : sqru(ulong x)
     129             : {
     130             :   ulong p1;
     131             :   LOCAL_HIREMAINDER;
     132             : 
     133   605866286 :   if (!x) return gen_0;
     134   605533224 :   p1 = mulll(x,x); return uutoi(hiremainder, p1);
     135             : }
     136             : 
     137             : /* assume x > 1, y != 0. Return u * y with sign s */
     138             : static GEN
     139   701525367 : mulur_2(ulong x, GEN y, long s)
     140             : {
     141   701525367 :   long m, sh, i, lx = lg(y), e = expo(y);
     142   701525367 :   GEN z = cgetg(lx, t_REAL);
     143             :   ulong garde;
     144             :   LOCAL_HIREMAINDER;
     145             : 
     146   700854453 :   y--; garde = mulll(x,y[lx]);
     147  3674397190 :   for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
     148   700854453 :   z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
     149   700854453 :   sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
     150   700854453 :   if (sh) shift_left(z,z, 2,lx-1, garde,sh);
     151   702115102 :   z[1] = evalsigne(s) | evalexpo(m+e);
     152   703518469 :   if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
     153   703307735 :   return z;
     154             : }
     155             : 
     156             : INLINE GEN
     157      473695 : mul0r(GEN x)
     158             : {
     159      473695 :   long l = realprec(x), e = expo(x);
     160      473695 :   e = (l > 0)? e - l: (e < 0? 2*e: 0);
     161      473695 :   return real_0_bit(e);
     162             : }
     163             : /* lg(x) > 2 */
     164             : INLINE GEN
     165      186465 : div0r(GEN x) {
     166      186465 :   long l = realprec(x), e = expo(x);
     167      186465 :   return real_0_bit(-l - e);
     168             : }
     169             : 
     170             : GEN
     171   116512165 : mulsr(long x, GEN y)
     172             : {
     173             :   long s;
     174             : 
     175   116512165 :   if (!x) return mul0r(y);
     176   116512123 :   s = signe(y);
     177   116512123 :   if (!s)
     178             :   {
     179      227461 :     if (x < 0) x = -x;
     180      227461 :     return real_0_bit( expo(y) + expu(x) );
     181             :   }
     182   116284662 :   if (x==1)  return rcopy(y);
     183   107588679 :   if (x==-1) return negr(y);
     184   104473580 :   if (x < 0)
     185    34717706 :     return mulur_2((ulong)-x, y, -s);
     186             :   else
     187    69755874 :     return mulur_2((ulong)x, y, s);
     188             : }
     189             : 
     190             : GEN
     191   919051158 : mulur(ulong x, GEN y)
     192             : {
     193             :   long s;
     194             : 
     195   919051158 :   if (!x) return mul0r(y);
     196   919051151 :   s = signe(y);
     197   919051151 :   if (!s) return real_0_bit( expo(y) + expu(x) );
     198   911757236 :   if (x==1) return rcopy(y);
     199   596811835 :   return mulur_2(x, y, s);
     200             : }
     201             : 
     202             : INLINE void
     203  3550952620 : mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
     204             : {
     205             :   long i;
     206  3550952620 :   if (hi[2] < 0)
     207             :   {
     208  1644444478 :     if (z != hi)
     209  1711571404 :       for (i=2; i<lz ; i++) z[i] = hi[i];
     210  1644444478 :     ez++;
     211             :   }
     212             :   else
     213             :   {
     214  1906508142 :     shift_left(z,hi,2,lz-1, garde, 1);
     215  1909564674 :     garde <<= 1;
     216             :   }
     217  3554009152 :   if (garde & HIGHBIT)
     218             :   { /* round to nearest */
     219  1740891157 :     i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
     220  1723401352 :     if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
     221             :   }
     222  3554009152 :   z[1] = evalsigne(sz)|evalexpo(ez);
     223  3549686315 : }
     224             : /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
     225             : INLINE void
     226   769165848 : mulrrz_3end(GEN z, long sz, long ez, ulong garde)
     227             : {
     228   769165848 :   if (z[2] < 0)
     229             :   { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
     230   360916359 :     if (garde & HIGHBIT) z[2]++; /* round properly */
     231   360916359 :     ez++;
     232             :   }
     233             :   else
     234             :   {
     235   408249489 :     uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
     236   408249489 :     if (garde & (1UL<<(BITS_IN_LONG-2)))
     237             :     {
     238   167522208 :       uel(z,2)++; /* round properly, z2+1 can overflow */
     239   167522208 :       if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
     240             :     }
     241             :   }
     242   769165848 :   z[1] = evalsigne(sz)|evalexpo(ez);
     243   769145928 : }
     244             : 
     245             : /* set z <-- x^2 != 0, floating point multiplication.
     246             :  * lz = lg(z) = lg(x) */
     247             : INLINE void
     248   559140657 : sqrz_i(GEN z, GEN x, long lz)
     249             : {
     250   559140657 :   long ez = 2*expo(x);
     251             :   long i, j, lzz, p1;
     252             :   ulong garde;
     253             :   GEN x1;
     254             :   LOCAL_HIREMAINDER;
     255             :   LOCAL_OVERFLOW;
     256             : 
     257   559140657 :   if (lz > prec2lg(SQRR_SQRI_LIMIT))
     258             :   {
     259    41236943 :     pari_sp av = avma;
     260    41236943 :     GEN hi = sqrispec_mirror(x+2, lz-2);
     261    41516043 :     mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
     262    41498031 :     set_avma(av); return;
     263             :   }
     264   517823388 :   if (lz == 3)
     265             :   {
     266   120784902 :     garde = mulll(x[2],x[2]);
     267   120784902 :     z[2] = hiremainder;
     268   120784902 :     mulrrz_3end(z, 1, ez, garde);
     269   120765435 :     return;
     270             :   }
     271             : 
     272   397038486 :   lzz = lz-1; p1 = x[lzz];
     273   397038486 :   if (p1)
     274             :   {
     275   361054648 :     (void)mulll(p1,x[3]);
     276   361054648 :     garde = addmul(p1,x[2]);
     277   361054648 :     z[lzz] = hiremainder;
     278             :   }
     279             :   else
     280             :   {
     281    35983838 :     garde = 0;
     282    35983838 :     z[lzz] = 0;
     283             :   }
     284  1838077069 :   for (j=lz-2, x1=x-j; j>=3; j--)
     285             :   {
     286  1441038583 :     p1 = x[j]; x1++;
     287  1441038583 :     if (p1)
     288             :     {
     289  1414010235 :       (void)mulll(p1,x1[lz+1]);
     290  1414010235 :       garde = addll(addmul(p1,x1[lz]), garde);
     291  7223779400 :       for (i=lzz; i>j; i--)
     292             :       {
     293  5809769165 :         hiremainder += overflow;
     294  5809769165 :         z[i] = addll(addmul(p1,x1[i]), z[i]);
     295             :       }
     296  1414010235 :       z[j] = hiremainder+overflow;
     297             :     }
     298    27028348 :     else z[j]=0;
     299             :   }
     300   397038486 :   p1 = x[2]; x1++;
     301   397038486 :   garde = addll(mulll(p1,x1[lz]), garde);
     302  2239401908 :   for (i=lzz; i>2; i--)
     303             :   {
     304  1842363422 :     hiremainder += overflow;
     305  1842363422 :     z[i] = addll(addmul(p1,x1[i]), z[i]);
     306             :   }
     307   397038486 :   z[2] = hiremainder+overflow;
     308   397038486 :   mulrrz_end(z, z, lz, 1, ez, garde);
     309             : }
     310             : 
     311             : /* lz "large" = lg(y) = lg(z), lg(x) > lz if flag = 1 and >= if flag = 0 */
     312             : INLINE void
     313    54231483 : mulrrz_int(GEN z, GEN x, GEN y, long lz, long flag, long sz)
     314             : {
     315    54231483 :   pari_sp av = avma;
     316    54231483 :   GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
     317    54231483 :   mulrrz_end(z, hi, lz, sz, expo(x)+expo(y), hi[lz]);
     318    54231483 :   set_avma(av);
     319    54231483 : }
     320             : 
     321             : /* lz = 3 */
     322             : INLINE void
     323   649337639 : mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
     324             : {
     325             :   ulong garde;
     326             :   LOCAL_HIREMAINDER;
     327   649337639 :   if (flag)
     328             :   {
     329    77685423 :     (void)mulll(x[2],y[3]);
     330    77685423 :     garde = addmul(x[2],y[2]);
     331             :   }
     332             :   else
     333   571652216 :     garde = mulll(x[2],y[2]);
     334   649337639 :   z[2] = hiremainder;
     335   649337639 :   mulrrz_3end(z, sz, expo(x)+expo(y), garde);
     336   649425922 : }
     337             : 
     338             : /* set z <-- x*y, floating point multiplication. Trailing 0s for x are
     339             :  * treated efficiently (important application: mulir).
     340             :  * lz = lg(z) = lg(x) <= ly <= lg(y), sz = signe(z). flag = lg(x) < lg(y) */
     341             : INLINE void
     342  3698755541 : mulrrz_i(GEN z, GEN x, GEN y, long lz, long flag, long sz)
     343             : {
     344             :   long ez, i, j, lzz, p1;
     345             :   ulong garde;
     346             :   GEN y1;
     347             :   LOCAL_HIREMAINDER;
     348             :   LOCAL_OVERFLOW;
     349             : 
     350  3698755541 :   if (x == y) { sqrz_i(z,x,lz); return; }
     351  3698755541 :   if (lz > prec2lg(MULRR_MULII_LIMIT)) { mulrrz_int(z,x,y,lz,flag,sz); return; }
     352  3647519252 :   if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
     353  2998183588 :   ez = expo(x) + expo(y);
     354  2998183588 :   if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
     355  2998183588 :   lzz=lz-1; p1=x[lzz];
     356  2998183588 :   if (p1)
     357             :   {
     358  2797055455 :     (void)mulll(p1,y[3]);
     359  2797055455 :     garde = addll(addmul(p1,y[2]), garde);
     360  2797055455 :     z[lzz] = overflow+hiremainder;
     361             :   }
     362   201128133 :   else z[lzz]=0;
     363 13464987202 :   for (j=lz-2, y1=y-j; j>=3; j--)
     364             :   {
     365 10466803614 :     p1 = x[j]; y1++;
     366 10466803614 :     if (p1)
     367             :     {
     368 10168160363 :       (void)mulll(p1,y1[lz+1]);
     369 10168160363 :       garde = addll(addmul(p1,y1[lz]), garde);
     370 61479222559 :       for (i=lzz; i>j; i--)
     371             :       {
     372 51311062196 :         hiremainder += overflow;
     373 51311062196 :         z[i] = addll(addmul(p1,y1[i]), z[i]);
     374             :       }
     375 10168160363 :       z[j] = hiremainder+overflow;
     376             :     }
     377   298643251 :     else z[j]=0;
     378             :   }
     379  2998183588 :   p1 = x[2]; y1++;
     380  2998183588 :   garde = addll(mulll(p1,y1[lz]), garde);
     381 16521900776 :   for (i=lzz; i>2; i--)
     382             :   {
     383 13523717188 :     hiremainder += overflow;
     384 13523717188 :     z[i] = addll(addmul(p1,y1[i]), z[i]);
     385             :   }
     386  2998183588 :   z[2] = hiremainder+overflow;
     387  2998183588 :   mulrrz_end(z, z, lz, sz, ez, garde);
     388             : }
     389             : 
     390             : GEN
     391  3939849426 : mulrr(GEN x, GEN y)
     392             : {
     393             :   long flag, ly, lz, sx, sy;
     394             :   GEN z;
     395             : 
     396  3939849426 :   if (x == y) return sqrr(x);
     397  3938239636 :   sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
     398  3674734422 :   sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
     399  3641332673 :   if (sy < 0) sx = -sx;
     400  3641332673 :   lz = lg(x);
     401  3641332673 :   ly = lg(y);
     402  3641332673 :   if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
     403  3641332673 :   z = cgetg(lz, t_REAL);
     404  3637187891 :   mulrrz_i(z, x,y, lz,flag, sx);
     405  3644383568 :   return z;
     406             : }
     407             : 
     408             : GEN
     409   591363327 : sqrr(GEN x)
     410             : {
     411   591363327 :   long lz, sx = signe(x);
     412             :   GEN z;
     413             : 
     414   591363327 :   if (!sx) return real_0_bit(2*expo(x));
     415   559641126 :   lz = lg(x); z = cgetg(lz, t_REAL);
     416   559156772 :   sqrz_i(z, x, lz);
     417   559840729 :   return z;
     418             : }
     419             : 
     420             : GEN
     421   977795285 : mulir(GEN x, GEN y)
     422             : {
     423   977795285 :   long sx = signe(x), sy;
     424   977795285 :   if (!sx) return mul0r(y);
     425   977321637 :   if (lgefint(x) == 3) {
     426   866869414 :     GEN z = mulur(uel(x,2), y);
     427   868107090 :     if (sx < 0) togglesign(z);
     428   867828713 :     return z;
     429             :   }
     430   110452223 :   sy = signe(y);
     431   110452223 :   if (!sy) return real_0_bit(expi(x) + expo(y));
     432   109767619 :   if (sy < 0) sx = -sx;
     433             :   {
     434   109767619 :     long lz = lg(y), lx = lgefint(x);
     435   109767619 :     GEN hi, z = cgetg(lz, t_REAL);
     436   109397020 :     pari_sp av = avma;
     437   109397020 :     if (lx < (lz>>1) || (lx < lz && lz > prec2lg(MULRR_MULII_LIMIT)))
     438             :     { /* size mantissa of x < half size of mantissa z, or lx < lz so large
     439             :        * that mulrr will call mulii anyway: mulii */
     440    47912192 :       x = itor(x, lg2prec(lx));
     441    47542885 :       hi = muliispec_mirror(y+2, x+2, lz-2, lx-2);
     442    49048669 :       mulrrz_end(z, hi, lz, sx, expo(x)+expo(y), hi[lz]);
     443             :     }
     444             :     else /* dubious: complete x with 0s and call mulrr */
     445    61453197 :       mulrrz_i(z, itor(x, lg2prec(lz)), y, lz, 0, sx);
     446   110879923 :     set_avma(av); return z;
     447             :   }
     448             : }
     449             : 
     450             : /* x + y*z, generic. If lgefint(z) <= 3, caller should use faster variants  */
     451             : static GEN
     452    90310763 : addmulii_gen(GEN x, GEN y, GEN z, long lz)
     453             : {
     454    90310763 :   long lx = lgefint(x), ly;
     455             :   pari_sp av;
     456             :   GEN t;
     457    90310763 :   if (lx == 2) return mulii(z,y);
     458    87542255 :   ly = lgefint(y);
     459    87542255 :   if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
     460    86938465 :   av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
     461    86939424 :   t = mulii(z, y);
     462    86931048 :   set_avma(av); return addii(t,x);
     463             : }
     464             : /* x + y*z, lgefint(z) == 3 */
     465             : static GEN
     466   466161384 : addmulii_lg3(GEN x, GEN y, GEN z)
     467             : {
     468   466161384 :   long s = signe(z), lx, ly;
     469   466161384 :   ulong w = z[2];
     470             :   pari_sp av;
     471             :   GEN t;
     472   466161384 :   if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
     473   384681608 :   lx = lgefint(x);
     474   384681608 :   ly = lgefint(y);
     475   384681608 :   if (lx == 2)
     476             :   { /* x = 0 */
     477    71281689 :     if (ly == 2) return gen_0;
     478    39528096 :     t = muluispec(w, y+2, ly-2);
     479    39529079 :     if (signe(y) < 0) s = -s;
     480    39529079 :     setsigne(t, s); return t;
     481             :   }
     482   313399919 :   if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
     483   272822775 :   av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
     484   272964337 :   t = muluispec(w, y+2, ly-2);
     485   272916145 :   if (signe(y) < 0) s = -s;
     486   272916145 :   setsigne(t, s);
     487   272916145 :   set_avma(av); return addii(x,t);
     488             : }
     489             : /* x + y*z */
     490             : GEN
     491   166876516 : addmulii(GEN x, GEN y, GEN z)
     492             : {
     493   166876516 :   long lz = lgefint(z);
     494   166876516 :   switch(lz)
     495             :   {
     496      740447 :     case 2: return icopy(x); /* z = 0, wasteful copy */
     497   126015782 :     case 3: return addmulii_lg3(x, y, z);
     498    40120287 :     default:return addmulii_gen(x, y, z, lz);
     499             :   }
     500             : }
     501             : /* x + y*z, returns x itself and not a copy when y*z = 0 */
     502             : GEN
     503  1566145347 : addmulii_inplace(GEN x, GEN y, GEN z)
     504             : {
     505             :   long lz;
     506  1566145347 :   if (lgefint(y) == 2) return x;
     507   390157424 :   lz = lgefint(z);
     508   390157424 :   switch(lz)
     509             :   {
     510       17568 :     case 2: return x;
     511   340180190 :     case 3: return addmulii_lg3(x, y, z);
     512    49959666 :     default:return addmulii_gen(x, y, z, lz);
     513             :   }
     514             : }
     515             : 
     516             : /* written by Bruno Haible following an idea of Robert Harley */
     517             : long
     518  2196660075 : vals(ulong z)
     519             : {
     520             :   static char tab[64]={-1,0,1,12,2,6,-1,13,3,-1,7,-1,-1,-1,-1,14,10,4,-1,-1,8,-1,-1,25,-1,-1,-1,-1,-1,21,27,15,31,11,5,-1,-1,-1,-1,-1,9,-1,-1,24,-1,-1,20,26,30,-1,-1,-1,-1,23,-1,19,29,-1,22,18,28,17,16,-1};
     521             : #ifdef LONG_IS_64BIT
     522             :   long s;
     523             : #endif
     524             : 
     525  2196660075 :   if (!z) return -1;
     526             : #ifdef LONG_IS_64BIT
     527  1971541440 :   if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
     528             : #endif
     529  2196660075 :   z |= ~z + 1;
     530  2196660075 :   z += z << 4;
     531  2196660075 :   z += z << 6;
     532  2196660075 :   z ^= z << 16; /* or  z -= z<<16 */
     533             : #ifdef LONG_IS_64BIT
     534  1971541440 :   return s + tab[(z&0xffffffff)>>26];
     535             : #else
     536   225118635 :   return tab[z>>26];
     537             : #endif
     538             : }
     539             : 
     540             : GEN
     541         168 : divsi(long x, GEN y)
     542             : {
     543         168 :   long p1, s = signe(y);
     544             :   LOCAL_HIREMAINDER;
     545             : 
     546         168 :   if (!s) pari_err_INV("divsi",gen_0);
     547         168 :   if (!x || lgefint(y)>3 || ((long)y[2])<0) return gen_0;
     548         168 :   hiremainder=0; p1=divll(labs(x),y[2]);
     549         168 :   if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
     550         168 :   if (s<0) p1 = -p1;
     551         168 :   return stoi(p1);
     552             : }
     553             : 
     554             : GEN
     555     3494754 : divir(GEN x, GEN y)
     556             : {
     557             :   GEN z;
     558     3494754 :   long ly = lg(y), lx = lgefint(x);
     559             :   pari_sp av;
     560             : 
     561     3494754 :   if (ly == 2) pari_err_INV("divir",y);
     562     3494758 :   if (lx == 2) return div0r(y);
     563     3308293 :   if (lx == 3) {
     564     1961486 :     z = divur(x[2], y);
     565     1961484 :     if (signe(x) < 0) togglesign(z);
     566     1961484 :     return z;
     567             :   }
     568     1346807 :   z = cgetg(ly, t_REAL); av = avma;
     569     1346806 :   affrr(divrr(itor(x, lg2prec(ly+1)), y), z);
     570     1346809 :   set_avma(av); return z;
     571             : }
     572             : 
     573             : GEN
     574     2493888 : divur(ulong x, GEN y)
     575             : {
     576             :   pari_sp av;
     577     2493888 :   long p = realprec(y);
     578             :   GEN z;
     579             : 
     580     2493888 :   if (p == 0) pari_err_INV("divur",y);
     581     2493888 :   if (!x) return div0r(y);
     582     2493888 :   if (p > INVNEWTON_LIMIT) {
     583         183 :     av = avma; z = invr(y);
     584         183 :     if (x == 1) return z;
     585         183 :     return gerepileuptoleaf(av, mulur(x, z));
     586             :   }
     587     2493705 :   z = cgetr(p); av = avma;
     588     2493687 :   affrr(divrr(utor(x, p + BITS_IN_LONG), y), z);
     589     2493707 :   set_avma(av); return z;
     590             : }
     591             : 
     592             : GEN
     593         798 : divsr(long x, GEN y)
     594             : {
     595             :   pari_sp av;
     596         798 :   long p = realprec(y);
     597             :   GEN z;
     598             : 
     599         798 :   if (p == 0) pari_err_INV("divsr",y);
     600         798 :   if (!x) return div0r(y);
     601         798 :   if (p > INVNEWTON_LIMIT) {
     602           0 :     av = avma; z = invr(y);
     603           0 :     if (x == 1) return z;
     604           0 :     if (x ==-1) { togglesign(z); return z; }
     605           0 :     return gerepileuptoleaf(av, mulsr(x, z));
     606             :   }
     607         798 :   z = cgetr(p); av = avma;
     608         798 :   affrr(divrr(stor(x,p + BITS_IN_LONG), y), z);
     609         798 :   set_avma(av); return z;
     610             : }
     611             : 
     612             : /* returns 1/y, assume y != 0 */
     613             : static GEN
     614    62373921 : invr_basecase(GEN y)
     615             : {
     616    62373921 :   long p = realprec(y);
     617    62373921 :   GEN z = cgetr(p);
     618    62364245 :   pari_sp av = avma;
     619    62364245 :   affrr(divrr(real_1(p + BITS_IN_LONG), y), z);
     620    62372950 :   set_avma(av); return z;
     621             : }
     622             : /* returns 1/b, Newton iteration */
     623             : GEN
     624    62374874 : invr(GEN b)
     625             : {
     626    62374874 :   const long s = 6;
     627    62374874 :   long i, p, l = lg(b);
     628             :   GEN x, a;
     629             :   ulong mask;
     630             : 
     631    62374874 :   if (l <= maxss(prec2lg(INVNEWTON_LIMIT), (1L<<s) + 2)) {
     632    62370656 :     if (l == 2) pari_err_INV("invr",b);
     633    62370656 :     return invr_basecase(b);
     634             :   }
     635        3675 :   mask = quadratic_prec_mask(l-2);
     636       25725 :   for(i=0, p=1; i<s; i++) { p <<= 1; if (mask & 1) p--; mask >>= 1; }
     637        3675 :   x = cgetg(l, t_REAL);
     638        3675 :   a = rcopy(b); a[1] = _evalexpo(0) | evalsigne(1);
     639        3675 :   affrr(invr_basecase(rtor(a, lg2prec(p+2))), x);
     640       12265 :   while (mask > 1)
     641             :   {
     642        8590 :     p <<= 1; if (mask & 1) p--;
     643        8590 :     mask >>= 1;
     644        8590 :     setlg(a, p + 2);
     645        8590 :     setlg(x, p + 2);
     646             :     /* TODO: mulrr(a,x) should be a half product (the higher half is known).
     647             :      * mulrr(x, ) already is */
     648        8590 :     affrr(addrr(x, mulrr(x, subsr(1, mulrr(a,x)))), x);
     649        8590 :     set_avma((pari_sp)a);
     650             :   }
     651        3675 :   x[1] = (b[1] & SIGNBITS) | evalexpo(expo(x)-expo(b));
     652        3675 :   set_avma((pari_sp)x); return x;
     653             : }
     654             : 
     655             : GEN
     656  2566534889 : modii(GEN x, GEN y)
     657             : {
     658  2566534889 :   switch(signe(x))
     659             :   {
     660   485983008 :     case 0: return gen_0;
     661  1546902465 :     case 1: return remii(x,y);
     662   533649416 :     default:
     663             :     {
     664   533649416 :       pari_sp av = avma;
     665   533649416 :       (void)new_chunk(lgefint(y));
     666   537962812 :       x = remii(x,y); set_avma(av);
     667   538018309 :       if (x==gen_0) return x;
     668   496572045 :       return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
     669             :     }
     670             :   }
     671             : }
     672             : 
     673             : void
     674      782234 : modiiz(GEN x, GEN y, GEN z)
     675             : {
     676      782234 :   const pari_sp av = avma;
     677      782234 :   affii(modii(x,y),z); set_avma(av);
     678      782213 : }
     679             : 
     680             : GEN
     681     9805042 : divrs(GEN x, long y)
     682             : {
     683             :   GEN z;
     684     9805042 :   if (y < 0)
     685             :   {
     686     5753732 :     z = divru(x, (ulong)-y);
     687     5753706 :     togglesign(z);
     688             :   }
     689             :   else
     690     4051310 :     z = divru(x, (ulong)y);
     691     9805016 :   return z;
     692             : }
     693             : 
     694             : GEN
     695  1095335226 : divru(GEN x, ulong y)
     696             : {
     697  1095335226 :   long i, lx, sh, e, s = signe(x);
     698             :   ulong garde;
     699             :   GEN z;
     700             :   LOCAL_HIREMAINDER;
     701             : 
     702  1095335226 :   if (!y) pari_err_INV("divru",gen_0);
     703  1095451112 :   if (!s) return real_0_bit(expo(x) - expu(y));
     704  1094814365 :   if (!(y & (y-1))) /* power of 2 */
     705             :   {
     706    92128301 :     if (y == 1) return rcopy(x);
     707    89414751 :     return shiftr(x, -expu(y));
     708             :   }
     709  1002686064 :   e = expo(x);
     710  1002686064 :   lx = lg(x);
     711  1002686064 :   z = cgetg(lx, t_REAL);
     712  1001987298 :   if (lx == 3)
     713             :   {
     714   160153803 :     if (y <= uel(x,2))
     715             :     {
     716   160153533 :       hiremainder = 0;
     717   160153533 :       z[2] = divll(x[2],y);
     718             :       /* we may have hiremainder != 0 ==> garde */
     719   160153533 :       garde = divll(0,y);
     720             :     }
     721             :     else
     722             :     {
     723         270 :       hiremainder = x[2];
     724         270 :       z[2] = divll(0,y);
     725         270 :       garde = hiremainder;
     726         270 :       e -= BITS_IN_LONG;
     727             :     }
     728             :   }
     729             :   else
     730             :   {
     731   841833495 :     ulong yp = get_Fl_red(y);
     732   842495186 :     if (y <= uel(x,2))
     733             :     {
     734   842514494 :       hiremainder = 0;
     735  7265061827 :       for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
     736             :       /* we may have hiremainder != 0 ==> garde */
     737   846841544 :       garde = divll_pre(0,y,yp);
     738             :     }
     739             :     else
     740             :     {
     741        3231 :       long l = lx-1;
     742        3231 :       hiremainder = x[2];
     743       56821 :       for (i=2; i<l; i++) z[i] = divll_pre(x[i+1],y,yp);
     744        3231 :       z[i] = divll_pre(0,y,yp);
     745        3765 :       garde = hiremainder;
     746        3765 :       e -= BITS_IN_LONG;
     747             :     }
     748             :   }
     749  1005036192 :   sh=bfffo(z[2]); /* z[2] != 0 */
     750  1005036192 :   if (sh) shift_left(z,z, 2,lx-1, garde,sh);
     751  1005186033 :   z[1] = evalsigne(s) | evalexpo(e-sh);
     752  1004270296 :   if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
     753  1004170030 :   return z;
     754             : }
     755             : 
     756             : GEN
     757   132149360 : truedvmdii(GEN x, GEN y, GEN *z)
     758             : {
     759             :   pari_sp av;
     760             :   GEN r, q, *gptr[2];
     761   132149360 :   if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
     762     5271718 :   if (z == ONLY_REM) return modii(x,y);
     763             : 
     764     5271718 :   av = avma;
     765     5271718 :   q = dvmdii(x,y,&r); /* assume that r is last on stack */
     766     5271718 :   switch(signe(r))
     767             :   {
     768      696759 :     case 0:
     769      696759 :       if (z) *z = gen_0;
     770      696759 :       return q;
     771     3549166 :     case 1:
     772     3549166 :       if (z) *z = r; else cgiv(r);
     773     3549166 :       return q;
     774     1025793 :     case -1: break;
     775             :   }
     776     1025793 :   q = addis(q, -signe(y));
     777     1025793 :   if (!z) return gerepileuptoint(av, q);
     778             : 
     779      849555 :   *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
     780      849554 :   gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(pari_sp)r,gptr,2);
     781      849419 :   return q;
     782             : }
     783             : GEN
     784   129033705 : truedvmdis(GEN x, long y, GEN *z)
     785             : {
     786   129033705 :   pari_sp av = avma;
     787             :   long r;
     788             :   GEN q;
     789             : 
     790   129033705 :   if (z == ONLY_REM) return modis(x, y);
     791   129033705 :   q = divis_rem(x,y,&r);
     792             : 
     793   129012525 :   if (r >= 0)
     794             :   {
     795   117462435 :     if (z) *z = utoi(r);
     796   117461247 :     return q;
     797             :   }
     798    11550090 :   q = gerepileuptoint(av, addis(q, (y < 0)? 1: -1));
     799    11562020 :   if (z) *z = utoi(r + labs(y));
     800    11561971 :   return q;
     801             : }
     802             : GEN
     803     6202031 : truedvmdsi(long x, GEN y, GEN *z)
     804             : {
     805             :   long q, r;
     806     6202031 :   if (z == ONLY_REM) return modsi(x, y);
     807     6202031 :   q = sdivsi_rem(x,y,&r);
     808     6202031 :   if (r >= 0) {
     809     6202031 :     if (z) *z = utoi(r);
     810     6202031 :     return stoi(q);
     811             :   }
     812           0 :   q = q - signe(y);
     813           0 :   if (!z) return stoi(q);
     814             : 
     815           0 :   *z = subiuspec(y+2,(ulong)-r, lgefint(y)-2);
     816           0 :   return stoi(q);
     817             : }
     818             : 
     819             : /* 2^n = shifti(gen_1, n) */
     820             : GEN
     821   102093288 : int2n(long n) {
     822             :   long i, m, l;
     823             :   GEN z;
     824   102093288 :   if (n < 0) return gen_0;
     825   102074366 :   if (n == 0) return gen_1;
     826             : 
     827    98430915 :   l = dvmdsBIL(n, &m) + 3;
     828    98429822 :   z = cgetipos(l);
     829   538973223 :   for (i = 2; i < l; i++) z[i] = 0;
     830    98404557 :   *int_MSW(z) = 1UL << m; return z;
     831             : }
     832             : /* To avoid problems when 2^(BIL-1) < n. Overflow cleanly, where int2n
     833             :  * returns gen_0 */
     834             : GEN
     835    25486119 : int2u(ulong n) {
     836             :   ulong i, m, l;
     837             :   GEN z;
     838    25486119 :   if (n == 0) return gen_1;
     839             : 
     840    25486119 :   l = dvmduBIL(n, &m) + 3;
     841    25486657 :   z = cgetipos(l);
     842    51351803 :   for (i = 2; i < l; i++) z[i] = 0;
     843    25484944 :   *int_MSW(z) = 1UL << m; return z;
     844             : }
     845             : /* 2^n - 1 */
     846             : GEN
     847        7644 : int2um1(ulong n) {
     848             :   ulong i, m, l;
     849             :   GEN z;
     850        7644 :   if (n == 0) return gen_0;
     851             : 
     852        7644 :   l = dvmduBIL(n, &m);
     853        7644 :   l += m? 3: 2;
     854        7644 :   z = cgetipos(l);
     855       15427 :   for (i = 2; i < l; i++) z[i] = ~0UL;
     856        7644 :   if (m) *int_MSW(z) = (1UL << m) - 1;
     857        7644 :   return z;
     858             : }
     859             : 
     860             : GEN
     861  2228764523 : shifti(GEN x, long n)
     862             : {
     863  2228764523 :   long s = signe(x);
     864             :   GEN y;
     865             : 
     866  2228764523 :   if(s == 0) return gen_0;
     867  1927056741 :   y = shiftispec(x + 2, lgefint(x) - 2, n);
     868  1926983067 :   if (signe(y)) setsigne(y, s);
     869  1926983067 :   return y;
     870             : }
     871             : 
     872             : /* actual operations will take place on a+2 and b+2: we strip the codewords */
     873             : GEN
     874 20697591830 : mulii(GEN a,GEN b)
     875             : {
     876             :   long sa,sb;
     877             :   GEN z;
     878             : 
     879 20697591830 :   sa=signe(a); if (!sa) return gen_0;
     880 14045190424 :   sb=signe(b); if (!sb) return gen_0;
     881 10328195667 :   if (sb<0) sa = -sa;
     882 10328195667 :   z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
     883 10358423931 :   setsigne(z,sa); return z;
     884             : }
     885             : 
     886             : GEN
     887  1832782288 : sqri(GEN a) { return sqrispec(a+2, lgefint(a)-2); }
     888             : 
     889             : /* sqrt()'s result may be off by 1 when a is not representable exactly as a
     890             :  * double [64bit machine] */
     891             : ulong
     892   228748112 : usqrt(ulong a)
     893             : {
     894   228748112 :   ulong x = (ulong)sqrt((double)a);
     895             : #ifdef LONG_IS_64BIT
     896   195958931 :   if (x > LOWMASK || x*x > a) x--;
     897             : #endif
     898   228748112 :   return x;
     899             : }
     900             : 
     901             : /********************************************************************/
     902             : /**                                                                **/
     903             : /**              EXPONENT / CONVERSION t_REAL --> double           **/
     904             : /**                                                                **/
     905             : /********************************************************************/
     906             : 
     907             : #ifdef LONG_IS_64BIT
     908             : long
     909    22292409 : dblexpo(double x)
     910             : {
     911             :   union { double f; ulong i; } fi;
     912    22292409 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     913    22292409 :   const int exp_mid = 0x3ff;/* exponent bias */
     914             : 
     915    22292409 :   if (x==0.) return -exp_mid;
     916    22292397 :   fi.f = x;
     917    22292397 :   return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
     918             : }
     919             : 
     920             : ulong
     921           0 : dblmantissa(double x)
     922             : {
     923             :   union { double f; ulong i; } fi;
     924           0 :   const int expo_len = 11; /* number of bits of exponent */
     925             : 
     926           0 :   if (x==0.) return 0;
     927           0 :   fi.f = x;
     928           0 :   return (fi.i << expo_len) | HIGHBIT;
     929             : }
     930             : 
     931             : GEN
     932    11493551 : dbltor(double x)
     933             : {
     934             :   GEN z;
     935             :   long e;
     936             :   union { double f; ulong i; } fi;
     937    11493551 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     938    11493551 :   const int exp_mid = 0x3ff;/* exponent bias */
     939    11493551 :   const int expo_len = 11; /* number of bits of exponent */
     940             : 
     941    11493551 :   if (x==0.) return real_0_bit(-exp_mid);
     942    11488231 :   fi.f = x; z = cgetr(DEFAULTPREC);
     943             :   {
     944    11485683 :     const ulong a = fi.i;
     945             :     ulong A;
     946    11485683 :     e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
     947    11485683 :     if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
     948    11485917 :     A = a << expo_len;
     949    11485917 :     if (e == -exp_mid)
     950             :     { /* unnormalized values */
     951           0 :       int sh = bfffo(A);
     952           0 :       e -= sh-1;
     953           0 :       z[2] = A << sh;
     954             :     }
     955             :     else
     956    11485917 :       z[2] = HIGHBIT | A;
     957    11485917 :     z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
     958             :   }
     959    11485917 :   return z;
     960             : }
     961             : 
     962             : double
     963   246817535 : rtodbl(GEN x)
     964             : {
     965   246817535 :   long ex,s=signe(x);
     966             :   ulong a;
     967             :   union { double f; ulong i; } fi;
     968   246817535 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     969   246817535 :   const int exp_mid = 0x3ff;/* exponent bias */
     970   246817535 :   const int expo_len = 11; /* number of bits of exponent */
     971             : 
     972   246817535 :   if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
     973             : 
     974             :   /* start by rounding to closest */
     975   218319627 :   a = (x[2] & (HIGHBIT-1)) + 0x400;
     976   218319627 :   if (a & HIGHBIT) { ex++; a=0; }
     977   218319627 :   if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
     978   218322005 :   fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
     979   218322005 :   if (s<0) fi.i |= HIGHBIT;
     980   218322005 :   return fi.f;
     981             : }
     982             : 
     983             : #else /* LONG_IS_64BIT */
     984             : 
     985             : #if   PARI_DOUBLE_FORMAT == 1
     986             : #  define INDEX0 1
     987             : #  define INDEX1 0
     988             : #elif PARI_DOUBLE_FORMAT == 0
     989             : #  define INDEX0 0
     990             : #  define INDEX1 1
     991             : #endif
     992             : 
     993             : long
     994     3745883 : dblexpo(double x)
     995             : {
     996             :   union { double f; ulong i[2]; } fi;
     997     3745883 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     998     3745883 :   const int exp_mid = 0x3ff;/* exponent bias */
     999     3745883 :   const int shift = mant_len-32;
    1000             : 
    1001     3745883 :   if (x==0.) return -exp_mid;
    1002     3745881 :   fi.f = x;
    1003             :   {
    1004     3745881 :     const ulong a = fi.i[INDEX0];
    1005     3745881 :     return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
    1006             :   }
    1007             : }
    1008             : 
    1009             : ulong
    1010           0 : dblmantissa(double x)
    1011             : {
    1012             :   union { double f; ulong i[2]; } fi;
    1013           0 :   const int expo_len = 11; /* number of bits of exponent */
    1014             : 
    1015           0 :   if (x==0.) return 0;
    1016           0 :   fi.f = x;
    1017             :   {
    1018           0 :     const ulong a = fi.i[INDEX0];
    1019           0 :     const ulong b = fi.i[INDEX1];
    1020           0 :     return HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
    1021             :   }
    1022             : }
    1023             : 
    1024             : GEN
    1025     1877203 : dbltor(double x)
    1026             : {
    1027             :   GEN z;
    1028             :   long e;
    1029             :   union { double f; ulong i[2]; } fi;
    1030     1877203 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
    1031     1877203 :   const int exp_mid = 0x3ff;/* exponent bias */
    1032     1877203 :   const int expo_len = 11; /* number of bits of exponent */
    1033     1877203 :   const int shift = mant_len-32;
    1034             : 
    1035     1877203 :   if (x==0.) return real_0_bit(-exp_mid);
    1036     1876329 :   fi.f = x; z = cgetr(DEFAULTPREC);
    1037             :   {
    1038     1876329 :     const ulong a = fi.i[INDEX0];
    1039     1876329 :     const ulong b = fi.i[INDEX1];
    1040             :     ulong A, B;
    1041     1876329 :     e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
    1042     1876329 :     if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
    1043     1876329 :     A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
    1044     1876329 :     B = b << expo_len;
    1045     1876329 :     if (e == -exp_mid)
    1046             :     { /* unnormalized values */
    1047             :       int sh;
    1048           0 :       if (A)
    1049             :       {
    1050           0 :         sh = bfffo(A);
    1051           0 :         e -= sh-1;
    1052           0 :         z[2] = (A << sh) | (B >> (32-sh));
    1053           0 :         z[3] = B << sh;
    1054             :       }
    1055             :       else
    1056             :       {
    1057           0 :         sh = bfffo(B); /* B != 0 */
    1058           0 :         e -= sh-1 + 32;
    1059           0 :         z[2] = B << sh;
    1060           0 :         z[3] = 0;
    1061             :       }
    1062             :     }
    1063             :     else
    1064             :     {
    1065     1876329 :       z[3] = B;
    1066     1876329 :       z[2] = HIGHBIT | A;
    1067             :     }
    1068     1876329 :     z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
    1069             :   }
    1070     1876329 :   return z;
    1071             : }
    1072             : 
    1073             : double
    1074    40482686 : rtodbl(GEN x)
    1075             : {
    1076    40482686 :   long ex,s=signe(x),lx=lg(x);
    1077             :   ulong a,b,k;
    1078             :   union { double f; ulong i[2]; } fi;
    1079    40482686 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
    1080    40482686 :   const int exp_mid = 0x3ff;/* exponent bias */
    1081    40482686 :   const int expo_len = 11; /* number of bits of exponent */
    1082    40482686 :   const int shift = mant_len-32;
    1083             : 
    1084    40482686 :   if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
    1085             : 
    1086             :   /* start by rounding to closest */
    1087    35665738 :   a = x[2] & (HIGHBIT-1);
    1088    35665738 :   if (lx > 3)
    1089             :   {
    1090    35538709 :     b = x[3] + 0x400UL; if (b < 0x400UL) a++;
    1091    35538709 :     if (a & HIGHBIT) { ex++; a=0; }
    1092             :   }
    1093      127029 :   else b = 0;
    1094    35665738 :   if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
    1095    35665738 :   ex += exp_mid;
    1096    35665738 :   k = (a >> expo_len) | (ex << shift);
    1097    35665738 :   if (s<0) k |= HIGHBIT;
    1098    35665738 :   fi.i[INDEX0] = k;
    1099    35665738 :   fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
    1100    35665738 :   return fi.f;
    1101             : }
    1102             : #endif /* LONG_IS_64BIT */
    1103             : 

Generated by: LCOV version 1.16