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.12.1 lcov report (development 24988-2584e74448) Lines: 539 569 94.7 %
Date: 2020-01-26 05:57:03 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. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */
      16             : ulong
      17    33699194 : invmod2BIL(ulong b)
      18             : {
      19             :   static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
      20    33699194 :   ulong x = b + tab[b & 7]; /* b^(-1) mod 2^4 */
      21             : 
      22             :   /* Newton applied to 1/x - b = 0 */
      23             : #ifdef LONG_IS_64BIT
      24    33699194 :   x = x*(2-b*x); /* one more pass necessary */
      25             : #endif
      26    33699194 :   x = x*(2-b*x);
      27    33699194 :   x = x*(2-b*x); return x*(2-b*x);
      28             : }
      29             : 
      30             : void
      31  1322548058 : affrr(GEN x, GEN y)
      32             : {
      33  1322548058 :   long i, lx, ly = lg(y);
      34  1322548058 :   if (!signe(x))
      35             :   {
      36     9061141 :     y[1] = evalexpo(minss(expo(x), -bit_accuracy(ly)));
      37     9061141 :     return;
      38             :   }
      39  1313486917 :   y[1] = x[1]; lx = lg(x);
      40  1313486917 :   if (lx <= ly)
      41             :   {
      42  1236722003 :     for (i=2; i<lx; i++) y[i]=x[i];
      43  1236722003 :     for (   ; i<ly; i++) y[i]=0;
      44  1236722003 :     return;
      45             :   }
      46    76764914 :   for (i=2; i<ly; i++) y[i]=x[i];
      47             :   /* lx > ly: round properly */
      48    76764914 :   if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
      49             : }
      50             : 
      51             : GEN
      52    11073722 : trunc2nr(GEN x, long n)
      53             : {
      54             :   long ex;
      55    11073722 :   if (!signe(x)) return gen_0;
      56    10962451 :   ex = expo(x) + n; if (ex < 0) return gen_0;
      57    10420465 :   return mantissa2nr(x, ex - bit_prec(x) + 1);
      58             : }
      59             : 
      60             : /* x a t_REAL, x = i/2^e, i a t_INT */
      61             : GEN
      62    36232744 : mantissa_real(GEN x, long *e)
      63             : {
      64    36232744 :   *e = bit_prec(x)-1-expo(x);
      65    36232744 :   return mantissa2nr(x, 0);
      66             : }
      67             : 
      68             : GEN
      69   715062704 : mului(ulong x, GEN y)
      70             : {
      71   715062704 :   long s = signe(y);
      72             :   GEN z;
      73             : 
      74   715062704 :   if (!s || !x) return gen_0;
      75   583551109 :   z = muluispec(x, y+2, lgefint(y)-2);
      76   581590496 :   setsigne(z,s); return z;
      77             : }
      78             : 
      79             : GEN
      80   274993138 : mulsi(long x, GEN y)
      81             : {
      82   274993138 :   long s = signe(y);
      83             :   GEN z;
      84             : 
      85   274993138 :   if (!s || !x) return gen_0;
      86   201893840 :   if (x<0) { s = -s; x = -x; }
      87   201893840 :   z = muluispec((ulong)x, y+2, lgefint(y)-2);
      88   201892635 :   setsigne(z,s); return z;
      89             : }
      90             : 
      91             : GEN
      92   194550747 : mulss(long x, long y)
      93             : {
      94             :   long p1;
      95             :   LOCAL_HIREMAINDER;
      96             : 
      97   194550747 :   if (!x || !y) return gen_0;
      98   194366507 :   if (x<0) {
      99      116585 :     x = -x;
     100      116585 :     if (y<0) { y = -y; p1 = mulll(x,y); return uutoi(hiremainder, p1); }
     101      116585 :     p1 = mulll(x,y); return uutoineg(hiremainder, p1);
     102             :   } else {
     103   194249922 :     if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
     104   194245963 :     p1 = mulll(x,y); return uutoi(hiremainder, p1);
     105             :   }
     106             : }
     107             : GEN
     108      361264 : sqrs(long x)
     109             : {
     110             :   long p1;
     111             :   LOCAL_HIREMAINDER;
     112             : 
     113      361264 :   if (!x) return gen_0;
     114      360942 :   if (x<0) x = -x;
     115      360942 :   p1 = mulll(x,x); return uutoi(hiremainder, p1);
     116             : }
     117             : GEN
     118  2044495544 : muluu(ulong x, ulong y)
     119             : {
     120             :   long p1;
     121             :   LOCAL_HIREMAINDER;
     122             : 
     123  2044495544 :   if (!x || !y) return gen_0;
     124  2044448105 :   p1 = mulll(x,y); return uutoi(hiremainder, p1);
     125             : }
     126             : GEN
     127   246353231 : sqru(ulong x)
     128             : {
     129             :   long p1;
     130             :   LOCAL_HIREMAINDER;
     131             : 
     132   246353231 :   if (!x) return gen_0;
     133   246072664 :   p1 = mulll(x,x); return uutoi(hiremainder, p1);
     134             : }
     135             : 
     136             : /* assume x > 1, y != 0. Return u * y with sign s */
     137             : static GEN
     138   585746261 : mulur_2(ulong x, GEN y, long s)
     139             : {
     140   585746261 :   long m, sh, i, lx = lg(y), e = expo(y);
     141   585746261 :   GEN z = cgetr(lx);
     142             :   ulong garde;
     143             :   LOCAL_HIREMAINDER;
     144             : 
     145   585741572 :   y--; garde = mulll(x,y[lx]);
     146   585741572 :   for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
     147   585741572 :   z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
     148   585741572 :   sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
     149   585741572 :   if (sh) shift_left(z,z, 2,lx-1, garde,sh);
     150   585765395 :   z[1] = evalsigne(s) | evalexpo(m+e);
     151   585750454 :   if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
     152   585749445 :   return z;
     153             : }
     154             : 
     155             : INLINE GEN
     156     1532477 : mul0r(GEN x)
     157             : {
     158     1532477 :   long l = lg(x), e = expo(x);
     159     1532477 :   e = (l > 2)? -prec2nbits(l) + e: (e < 0? 2*e: 0);
     160     1532477 :   return real_0_bit(e);
     161             : }
     162             : /* lg(x) > 2 */
     163             : INLINE GEN
     164      105294 : div0r(GEN x) {
     165      105294 :   long l = lg(x), e = expo(x);
     166      105294 :   return real_0_bit(-prec2nbits(l) - e);
     167             : }
     168             : 
     169             : GEN
     170    70661304 : mulsr(long x, GEN y)
     171             : {
     172             :   long s;
     173             : 
     174    70661304 :   if (!x) return mul0r(y);
     175    69643911 :   s = signe(y);
     176    69643911 :   if (!s)
     177             :   {
     178      225890 :     if (x < 0) x = -x;
     179      225890 :     return real_0_bit( expo(y) + expu(x) );
     180             :   }
     181    69418021 :   if (x==1)  return rcopy(y);
     182    63109844 :   if (x==-1) return negr(y);
     183    61780905 :   if (x < 0)
     184    29606301 :     return mulur_2((ulong)-x, y, -s);
     185             :   else
     186    32174604 :     return mulur_2((ulong)x, y, s);
     187             : }
     188             : 
     189             : GEN
     190   606542912 : mulur(ulong x, GEN y)
     191             : {
     192             :   long s;
     193             : 
     194   606542912 :   if (!x) return mul0r(y);
     195   606542912 :   s = signe(y);
     196   606542912 :   if (!s) return real_0_bit( expo(y) + expu(x) );
     197   602010582 :   if (x==1) return rcopy(y);
     198   523979962 :   return mulur_2(x, y, s);
     199             : }
     200             : 
     201             : INLINE void
     202  1440036211 : mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
     203             : {
     204             :   long i;
     205  1440036211 :   if (hi[2] < 0)
     206             :   {
     207   661738602 :     if (z != hi)
     208    19235541 :       for (i=2; i<lz ; i++) z[i] = hi[i];
     209   661738602 :     ez++;
     210             :   }
     211             :   else
     212             :   {
     213   778297609 :     shift_left(z,hi,2,lz-1, garde, 1);
     214   778201218 :     garde <<= 1;
     215             :   }
     216  1439939820 :   if (garde & HIGHBIT)
     217             :   { /* round to nearest */
     218   702730086 :     i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
     219   702730086 :     if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
     220             :   }
     221  1439939820 :   z[1] = evalsigne(sz)|evalexpo(ez);
     222  1439073593 : }
     223             : /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
     224             : INLINE void
     225  1110944516 : mulrrz_3end(GEN z, long sz, long ez, ulong garde)
     226             : {
     227  1110944516 :   if (z[2] < 0)
     228             :   { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
     229   602331680 :     if (garde & HIGHBIT) z[2]++; /* round properly */
     230   602331680 :     ez++;
     231             :   }
     232             :   else
     233             :   {
     234   508612836 :     uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
     235   508612836 :     if (garde & (1UL<<(BITS_IN_LONG-2)))
     236             :     {
     237   233660261 :       uel(z,2)++; /* round properly, z2+1 can overflow */
     238   233660261 :       if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
     239             :     }
     240             :   }
     241  1110944516 :   z[1] = evalsigne(sz)|evalexpo(ez);
     242  1110933866 : }
     243             : 
     244             : /* set z <-- x^2 != 0, floating point multiplication.
     245             :  * lz = lg(z) = lg(x) */
     246             : INLINE void
     247   107058988 : sqrz_i(GEN z, GEN x, long lz)
     248             : {
     249   107058988 :   long ez = 2*expo(x);
     250             :   long i, j, lzz, p1;
     251             :   ulong garde;
     252             :   GEN x1;
     253             :   LOCAL_HIREMAINDER;
     254             :   LOCAL_OVERFLOW;
     255             : 
     256   107058988 :   if (lz > SQRR_SQRI_LIMIT)
     257             :   {
     258     7621274 :     pari_sp av = avma;
     259     7621274 :     GEN hi = sqrispec_mirror(x+2, lz-2);
     260     7621274 :     mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
     261     7621274 :     set_avma(av); return;
     262             :   }
     263    99437714 :   if (lz == 3)
     264             :   {
     265    22838646 :     garde = mulll(x[2],x[2]);
     266    22838646 :     z[2] = hiremainder;
     267    22838646 :     mulrrz_3end(z, 1, ez, garde);
     268    22838633 :     return;
     269             :   }
     270             : 
     271    76599068 :   lzz = lz-1; p1 = x[lzz];
     272    76599068 :   if (p1)
     273             :   {
     274    68171651 :     (void)mulll(p1,x[3]);
     275    68171651 :     garde = addmul(p1,x[2]);
     276    68171651 :     z[lzz] = hiremainder;
     277             :   }
     278             :   else
     279             :   {
     280     8427417 :     garde = 0;
     281     8427417 :     z[lzz] = 0;
     282             :   }
     283   258653766 :   for (j=lz-2, x1=x-j; j>=3; j--)
     284             :   {
     285   182054698 :     p1 = x[j]; x1++;
     286   182054698 :     if (p1)
     287             :     {
     288   174390984 :       (void)mulll(p1,x1[lz+1]);
     289   174390984 :       garde = addll(addmul(p1,x1[lz]), garde);
     290   757892569 :       for (i=lzz; i>j; i--)
     291             :       {
     292   583501585 :         hiremainder += overflow;
     293   583501585 :         z[i] = addll(addmul(p1,x1[i]), z[i]);
     294             :       }
     295   174390984 :       z[j] = hiremainder+overflow;
     296             :     }
     297     7663714 :     else z[j]=0;
     298             :   }
     299    76599068 :   p1 = x[2]; x1++;
     300    76599068 :   garde = addll(mulll(p1,x1[lz]), garde);
     301   335253175 :   for (i=lzz; i>2; i--)
     302             :   {
     303   258654107 :     hiremainder += overflow;
     304   258654107 :     z[i] = addll(addmul(p1,x1[i]), z[i]);
     305             :   }
     306    76599068 :   z[2] = hiremainder+overflow;
     307    76599068 :   mulrrz_end(z, z, lz, 1, ez, garde);
     308             : }
     309             : 
     310             : /* lz "large" = lg(y) = lg(z), lg(x) > lz if flag = 1 and >= if flag = 0 */
     311             : INLINE void
     312    18575112 : mulrrz_int(GEN z, GEN x, GEN y, long lz, long flag, long sz)
     313             : {
     314    18575112 :   pari_sp av = avma;
     315    18575112 :   GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
     316    18575112 :   mulrrz_end(z, hi, lz, sz, expo(x)+expo(y), hi[lz]);
     317    18575112 :   set_avma(av);
     318    18575112 : }
     319             : 
     320             : /* lz = 3 */
     321             : INLINE void
     322  1088100852 : mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
     323             : {
     324             :   ulong garde;
     325             :   LOCAL_HIREMAINDER;
     326  1088100852 :   if (flag)
     327             :   {
     328    38154992 :     (void)mulll(x[2],y[3]);
     329    38154992 :     garde = addmul(x[2],y[2]);
     330             :   }
     331             :   else
     332  1049945860 :     garde = mulll(x[2],y[2]);
     333  1088100852 :   z[2] = hiremainder;
     334  1088100852 :   mulrrz_3end(z, sz, expo(x)+expo(y), garde);
     335  1088097481 : }
     336             : 
     337             : /* set z <-- x*y, floating point multiplication. Trailing 0s for x are
     338             :  * treated efficiently (important application: mulir).
     339             :  * lz = lg(z) = lg(x) <= ly <= lg(y), sz = signe(z). flag = lg(x) < lg(y) */
     340             : INLINE void
     341  2430558197 : mulrrz_i(GEN z, GEN x, GEN y, long lz, long flag, long sz)
     342             : {
     343             :   long ez, i, j, lzz, p1;
     344             :   ulong garde;
     345             :   GEN y1;
     346             :   LOCAL_HIREMAINDER;
     347             :   LOCAL_OVERFLOW;
     348             : 
     349  2430558197 :   if (x == y) { sqrz_i(z,x,lz); return; }
     350  2430558197 :   if (lz > MULRR_MULII_LIMIT) { mulrrz_int(z,x,y,lz,flag,sz); return; }
     351  2411983085 :   if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
     352  1323886897 :   ez = expo(x) + expo(y);
     353  1323886897 :   if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
     354  1323886897 :   lzz=lz-1; p1=x[lzz];
     355  1323886897 :   if (p1)
     356             :   {
     357  1267664304 :     (void)mulll(p1,y[3]);
     358  1267664304 :     garde = addll(addmul(p1,y[2]), garde);
     359  1267664304 :     z[lzz] = overflow+hiremainder;
     360             :   }
     361    56222593 :   else z[lzz]=0;
     362  5218116458 :   for (j=lz-2, y1=y-j; j>=3; j--)
     363             :   {
     364  3894229561 :     p1 = x[j]; y1++;
     365  3894229561 :     if (p1)
     366             :     {
     367  3764134969 :       (void)mulll(p1,y1[lz+1]);
     368  3764134969 :       garde = addll(addmul(p1,y1[lz]), garde);
     369 26144642061 :       for (i=lzz; i>j; i--)
     370             :       {
     371 22380507092 :         hiremainder += overflow;
     372 22380507092 :         z[i] = addll(addmul(p1,y1[i]), z[i]);
     373             :       }
     374  3764134969 :       z[j] = hiremainder+overflow;
     375             :     }
     376   130094592 :     else z[j]=0;
     377             :   }
     378  1323886897 :   p1 = x[2]; y1++;
     379  1323886897 :   garde = addll(mulll(p1,y1[lz]), garde);
     380  6543639050 :   for (i=lzz; i>2; i--)
     381             :   {
     382  5219752153 :     hiremainder += overflow;
     383  5219752153 :     z[i] = addll(addmul(p1,y1[i]), z[i]);
     384             :   }
     385  1323886897 :   z[2] = hiremainder+overflow;
     386  1323886897 :   mulrrz_end(z, z, lz, sz, ez, garde);
     387             : }
     388             : 
     389             : GEN
     390  2467689390 : mulrr(GEN x, GEN y)
     391             : {
     392             :   long flag, ly, lz, sx, sy;
     393             :   GEN z;
     394             : 
     395  2467689390 :   if (x == y) return sqrr(x);
     396  2467475113 :   sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
     397  2418953613 :   sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
     398  2401100014 :   if (sy < 0) sx = -sx;
     399  2401100014 :   lz = lg(x);
     400  2401100014 :   ly = lg(y);
     401  2401100014 :   if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
     402  2401100014 :   z = cgetr(lz);
     403  2399170194 :   mulrrz_i(z, x,y, lz,flag, sx);
     404  2401545058 :   return z;
     405             : }
     406             : 
     407             : GEN
     408   115575320 : sqrr(GEN x)
     409             : {
     410   115575320 :   long lz, sx = signe(x);
     411             :   GEN z;
     412             : 
     413   115575320 :   if (!sx) return real_0_bit(2*expo(x));
     414   107058381 :   lz = lg(x); z = cgetr(lz);
     415   107057410 :   sqrz_i(z, x, lz);
     416   107058640 :   return z;
     417             : }
     418             : 
     419             : GEN
     420   280245395 : mulir(GEN x, GEN y)
     421             : {
     422   280245395 :   long sx = signe(x), sy;
     423   280245395 :   if (!sx) return mul0r(y);
     424   279730311 :   if (lgefint(x) == 3) {
     425   234257539 :     GEN z = mulur(uel(x,2), y);
     426   234254643 :     if (sx < 0) togglesign(z);
     427   234254092 :     return z;
     428             :   }
     429    45472772 :   sy = signe(y);
     430    45472772 :   if (!sy) return real_0_bit(expi(x) + expo(y));
     431    40809229 :   if (sy < 0) sx = -sx;
     432             :   {
     433    40809229 :     long lz = lg(y), lx = lgefint(x);
     434    40809229 :     GEN hi, z = cgetr(lz);
     435    40809229 :     pari_sp av = avma;
     436    40809229 :     if (lx < (lz>>1) || (lx < lz && lz > MULRR_MULII_LIMIT))
     437             :     { /* size mantissa of x < half size of mantissa z, or lx < lz so large
     438             :        * that mulrr will call mulii anyway: mulii */
     439    12603941 :       x = itor(x, lx);
     440    12603941 :       hi = muliispec_mirror(y+2, x+2, lz-2, lx-2);
     441    12603941 :       mulrrz_end(z, hi, lz, sx, expo(x)+expo(y), hi[lz]);
     442             :     }
     443             :     else /* dubious: complete x with 0s and call mulrr */
     444    28205288 :       mulrrz_i(z, itor(x,lz), y, lz, 0, sx);
     445    40809229 :     set_avma(av); return z;
     446             :   }
     447             : }
     448             : 
     449             : /* x + y*z, generic. If lgefint(z) <= 3, caller should use faster variants  */
     450             : static GEN
     451    37384749 : addmulii_gen(GEN x, GEN y, GEN z, long lz)
     452             : {
     453    37384749 :   long lx = lgefint(x), ly;
     454             :   pari_sp av;
     455             :   GEN t;
     456    37384749 :   if (lx == 2) return mulii(z,y);
     457    35264122 :   ly = lgefint(y);
     458    35264122 :   if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
     459    35250272 :   av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
     460    35250272 :   t = mulii(z, y);
     461    35250272 :   set_avma(av); return addii(t,x);
     462             : }
     463             : /* x + y*z, lgefint(z) == 3 */
     464             : static GEN
     465   347424655 : addmulii_lg3(GEN x, GEN y, GEN z)
     466             : {
     467   347424655 :   long s = signe(z), lx, ly;
     468   347424655 :   ulong w = z[2];
     469             :   pari_sp av;
     470             :   GEN t;
     471   347424655 :   if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
     472   286902770 :   lx = lgefint(x);
     473   286902770 :   ly = lgefint(y);
     474   286902770 :   if (lx == 2)
     475             :   { /* x = 0 */
     476    58107900 :     if (ly == 2) return gen_0;
     477    25402231 :     t = muluispec(w, y+2, ly-2);
     478    25401789 :     if (signe(y) < 0) s = -s;
     479    25401789 :     setsigne(t, s); return t;
     480             :   }
     481   228794870 :   if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
     482   201475180 :   av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
     483   201474426 :   t = muluispec(w, y+2, ly-2);
     484   201472170 :   if (signe(y) < 0) s = -s;
     485   201472170 :   setsigne(t, s);
     486   201472170 :   set_avma(av); return addii(x,t);
     487             : }
     488             : /* x + y*z */
     489             : GEN
     490   114841577 : addmulii(GEN x, GEN y, GEN z)
     491             : {
     492   114841577 :   long lz = lgefint(z);
     493   114841577 :   switch(lz)
     494             :   {
     495      469404 :     case 2: return icopy(x); /* z = 0, wasteful copy */
     496   107856985 :     case 3: return addmulii_lg3(x, y, z);
     497     6515188 :     default:return addmulii_gen(x, y, z, lz);
     498             :   }
     499             : }
     500             : /* x + y*z, returns x itself and not a copy when y*z = 0 */
     501             : GEN
     502  1159910060 : addmulii_inplace(GEN x, GEN y, GEN z)
     503             : {
     504             :   long lz;
     505  1159910060 :   if (lgefint(y) == 2) return x;
     506   270433183 :   lz = lgefint(z);
     507   270433183 :   switch(lz)
     508             :   {
     509        1225 :     case 2: return x;
     510   239562397 :     case 3: return addmulii_lg3(x, y, z);
     511    30869561 :     default:return addmulii_gen(x, y, z, lz);
     512             :   }
     513             : }
     514             : 
     515             : /* written by Bruno Haible following an idea of Robert Harley */
     516             : long
     517  2559880107 : vals(ulong z)
     518             : {
     519             :   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};
     520             : #ifdef LONG_IS_64BIT
     521             :   long s;
     522             : #endif
     523             : 
     524  2559880107 :   if (!z) return -1;
     525             : #ifdef LONG_IS_64BIT
     526  2247703371 :   if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
     527             : #endif
     528  2559880107 :   z |= ~z + 1;
     529  2559880107 :   z += z << 4;
     530  2559880107 :   z += z << 6;
     531  2559880107 :   z ^= z << 16; /* or  z -= z<<16 */
     532             : #ifdef LONG_IS_64BIT
     533  2247703371 :   return s + tab[(z&0xffffffff)>>26];
     534             : #else
     535   312176736 :   return tab[z>>26];
     536             : #endif
     537             : }
     538             : 
     539             : GEN
     540          21 : divsi(long x, GEN y)
     541             : {
     542          21 :   long p1, s = signe(y);
     543             :   LOCAL_HIREMAINDER;
     544             : 
     545          21 :   if (!s) pari_err_INV("divsi",gen_0);
     546          21 :   if (!x || lgefint(y)>3 || ((long)y[2])<0) return gen_0;
     547          21 :   hiremainder=0; p1=divll(labs(x),y[2]);
     548          21 :   if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
     549          21 :   if (s<0) p1 = -p1;
     550          21 :   return stoi(p1);
     551             : }
     552             : 
     553             : GEN
     554     2629481 : divir(GEN x, GEN y)
     555             : {
     556             :   GEN z;
     557     2629481 :   long ly = lg(y), lx = lgefint(x);
     558             :   pari_sp av;
     559             : 
     560     2629481 :   if (ly == 2) pari_err_INV("divir",y);
     561     2629481 :   if (lx == 2) return div0r(y);
     562     2524187 :   if (lx == 3) {
     563     1346507 :     z = divur(x[2], y);
     564     1346507 :     if (signe(x) < 0) togglesign(z);
     565     1346507 :     return z;
     566             :   }
     567     1177680 :   z = cgetr(ly); av = avma;
     568     1177680 :   affrr(divrr(itor(x, ly+1), y), z);
     569     1177680 :   set_avma(av); return z;
     570             : }
     571             : 
     572             : GEN
     573     1397906 : divur(ulong x, GEN y)
     574             : {
     575             :   pari_sp av;
     576     1397906 :   long ly = lg(y);
     577             :   GEN z;
     578             : 
     579     1397906 :   if (ly == 2) pari_err_INV("divur",y);
     580     1397906 :   if (!x) return div0r(y);
     581     1397906 :   if (ly > INVNEWTON_LIMIT) {
     582         254 :     av = avma; z = invr(y);
     583         254 :     if (x == 1) return z;
     584         254 :     return gerepileuptoleaf(av, mulur(x, z));
     585             :   }
     586     1397652 :   z = cgetr(ly); av = avma;
     587     1397652 :   affrr(divrr(utor(x,ly+1), y), z);
     588     1397652 :   set_avma(av); return z;
     589             : }
     590             : 
     591             : GEN
     592         182 : divsr(long x, GEN y)
     593             : {
     594             :   pari_sp av;
     595         182 :   long ly = lg(y);
     596             :   GEN z;
     597             : 
     598         182 :   if (ly == 2) pari_err_INV("divsr",y);
     599         182 :   if (!x) return div0r(y);
     600         182 :   if (ly > INVNEWTON_LIMIT) {
     601           0 :     av = avma; z = invr(y);
     602           0 :     if (x == 1) return z;
     603           0 :     if (x ==-1) { togglesign(z); return z; }
     604           0 :     return gerepileuptoleaf(av, mulsr(x, z));
     605             :   }
     606         182 :   z = cgetr(ly); av = avma;
     607         182 :   affrr(divrr(stor(x,ly+1), y), z);
     608         182 :   set_avma(av); return z;
     609             : }
     610             : 
     611             : /* returns 1/y, assume y != 0 */
     612             : static GEN
     613    18016666 : invr_basecase(GEN y)
     614             : {
     615    18016666 :   long ly = lg(y);
     616    18016666 :   GEN z = cgetr(ly);
     617    18014533 :   pari_sp av = avma;
     618    18014533 :   affrr(divrr(real_1(ly+1), y), z);
     619    18017507 :   set_avma(av); return z;
     620             : }
     621             : /* returns 1/b, Newton iteration */
     622             : GEN
     623    18016603 : invr(GEN b)
     624             : {
     625    18016603 :   const long s = 6;
     626    18016603 :   long i, p, l = lg(b);
     627             :   GEN x, a;
     628             :   ulong mask;
     629             : 
     630    18016603 :   if (l <= maxss(INVNEWTON_LIMIT, (1L<<s) + 2)) {
     631    17689445 :     if (l == 2) pari_err_INV("invr",b);
     632    17689445 :     return invr_basecase(b);
     633             :   }
     634      327246 :   mask = quadratic_prec_mask(l-2);
     635      327246 :   for(i=0, p=1; i<s; i++) { p <<= 1; if (mask & 1) p--; mask >>= 1; }
     636      327246 :   x = cgetr(l);
     637      327246 :   a = rcopy(b); a[1] = _evalexpo(0) | evalsigne(1);
     638      327246 :   affrr(invr_basecase(rtor(a, p+2)), x);
     639     1130373 :   while (mask > 1)
     640             :   {
     641      475881 :     p <<= 1; if (mask & 1) p--;
     642      475881 :     mask >>= 1;
     643      475881 :     setlg(a, p + 2);
     644      475881 :     setlg(x, p + 2);
     645             :     /* TODO: mulrr(a,x) should be a half product (the higher half is known).
     646             :      * mulrr(x, ) already is */
     647      475881 :     affrr(addrr(x, mulrr(x, subsr(1, mulrr(a,x)))), x);
     648      475881 :     set_avma((pari_sp)a);
     649             :   }
     650      327246 :   x[1] = (b[1] & SIGNBITS) | evalexpo(expo(x)-expo(b));
     651      327246 :   set_avma((pari_sp)x); return x;
     652             : }
     653             : 
     654             : GEN
     655  1469862604 : modii(GEN x, GEN y)
     656             : {
     657  1469862604 :   switch(signe(x))
     658             :   {
     659   283979497 :     case 0: return gen_0;
     660   845447222 :     case 1: return remii(x,y);
     661             :     default:
     662             :     {
     663   340435885 :       pari_sp av = avma;
     664   340435885 :       (void)new_chunk(lgefint(y));
     665   340442919 :       x = remii(x,y); set_avma(av);
     666   340429774 :       if (x==gen_0) return x;
     667   313268271 :       return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
     668             :     }
     669             :   }
     670             : }
     671             : 
     672             : void
     673      265830 : modiiz(GEN x, GEN y, GEN z)
     674             : {
     675      265830 :   const pari_sp av = avma;
     676      265830 :   affii(modii(x,y),z); set_avma(av);
     677      265830 : }
     678             : 
     679             : GEN
     680     8815883 : divrs(GEN x, long y)
     681             : {
     682             :   GEN z;
     683     8815883 :   if (y < 0)
     684             :   {
     685      838424 :     z = divru(x, (ulong)-y);
     686      838424 :     togglesign(z);
     687             :   }
     688             :   else
     689     7977459 :     z = divru(x, (ulong)y);
     690     8815883 :   return z;
     691             : }
     692             : 
     693             : GEN
     694   594881518 : divru(GEN x, ulong y)
     695             : {
     696   594881518 :   long i, lx, sh, e, s = signe(x);
     697             :   ulong garde;
     698             :   GEN z;
     699             :   LOCAL_HIREMAINDER;
     700             : 
     701   594881518 :   if (!y) pari_err_INV("divru",gen_0);
     702   595044162 :   if (!s) return real_0_bit(expo(x) - expu(y));
     703   594912307 :   if (!(y & (y-1))) /* power of 2 */
     704             :   {
     705    45221702 :     if (y == 1) return rcopy(x);
     706    42951158 :     return shiftr(x, -expu(y));
     707             :   }
     708   549690605 :   e = expo(x);
     709   549690605 :   lx = lg(x);
     710   549690605 :   z = cgetr(lx);
     711   549750660 :   if (lx == 3)
     712             :   {
     713   334738664 :     if (y <= uel(x,2))
     714             :     {
     715   334738604 :       hiremainder = 0;
     716   334738604 :       z[2] = divll(x[2],y);
     717             :       /* we may have hiremainder != 0 ==> garde */
     718   334738604 :       garde = divll(0,y);
     719             :     }
     720             :     else
     721             :     {
     722          60 :       hiremainder = x[2];
     723          60 :       z[2] = divll(0,y);
     724          60 :       garde = hiremainder;
     725          60 :       e -= BITS_IN_LONG;
     726             :     }
     727             :   }
     728             :   else
     729             :   {
     730   215011996 :     ulong yp = get_Fl_red(y);
     731   214960898 :     if (y <= uel(x,2))
     732             :     {
     733   214958075 :       hiremainder = 0;
     734   214958075 :       for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
     735             :       /* we may have hiremainder != 0 ==> garde */
     736   215115029 :       garde = divll_pre(0,y,yp);
     737             :     }
     738             :     else
     739             :     {
     740        2823 :       long l = lx-1;
     741        2823 :       hiremainder = x[2];
     742        2823 :       for (i=2; i<l; i++) z[i] = divll_pre(x[i+1],y,yp);
     743        2823 :       z[i] = divll_pre(0,y,yp);
     744        2823 :       garde = hiremainder;
     745        2823 :       e -= BITS_IN_LONG;
     746             :     }
     747             :   }
     748   549858073 :   sh=bfffo(z[2]); /* z[2] != 0 */
     749   549858073 :   if (sh) shift_left(z,z, 2,lx-1, garde,sh);
     750   549895343 :   z[1] = evalsigne(s) | evalexpo(e-sh);
     751   549770103 :   if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
     752   549766524 :   return z;
     753             : }
     754             : 
     755             : GEN
     756    79300174 : truedvmdii(GEN x, GEN y, GEN *z)
     757             : {
     758             :   pari_sp av;
     759             :   GEN r, q, *gptr[2];
     760    79300174 :   if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
     761     2714673 :   if (z == ONLY_REM) return modii(x,y);
     762             : 
     763     2714673 :   av = avma;
     764     2714673 :   q = dvmdii(x,y,&r); /* assume that r is last on stack */
     765     2714673 :   switch(signe(r))
     766             :   {
     767             :     case 0:
     768      427821 :       if (z) *z = gen_0;
     769      427821 :       return q;
     770             :     case 1:
     771     1585478 :       if (z) *z = r; else cgiv(r);
     772     1585478 :       return q;
     773      701374 :     case -1: break;
     774             :   }
     775      701374 :   q = addis(q, -signe(y));
     776      701374 :   if (!z) return gerepileuptoint(av, q);
     777             : 
     778      293224 :   *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
     779      293224 :   gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(pari_sp)r,gptr,2);
     780      293224 :   return q;
     781             : }
     782             : GEN
     783    76715735 : truedvmdis(GEN x, long y, GEN *z)
     784             : {
     785    76715735 :   pari_sp av = avma;
     786             :   long r;
     787             :   GEN q;
     788             : 
     789    76715735 :   if (z == ONLY_REM) return modis(x, y);
     790    76715735 :   q = divis_rem(x,y,&r);
     791             : 
     792    76715755 :   if (r >= 0)
     793             :   {
     794    72063239 :     if (z) *z = utoi(r);
     795    72063156 :     return q;
     796             :   }
     797     4652516 :   q = gerepileuptoint(av, addis(q, (y < 0)? 1: -1));
     798     4652516 :   if (z) *z = utoi(r + labs(y));
     799     4652516 :   return q;
     800             : }
     801             : GEN
     802     6168301 : truedvmdsi(long x, GEN y, GEN *z)
     803             : {
     804             :   long q, r;
     805     6168301 :   if (z == ONLY_REM) return modsi(x, y);
     806     6168301 :   q = sdivsi_rem(x,y,&r);
     807     6168301 :   if (r >= 0) {
     808     6167993 :     if (z) *z = utoi(r);
     809     6167993 :     return stoi(q);
     810             :   }
     811         308 :   q = q - signe(y);
     812         308 :   if (!z) return stoi(q);
     813             : 
     814           0 :   *z = subiuspec(y+2,(ulong)-r, lgefint(y)-2);
     815           0 :   return stoi(q);
     816             : }
     817             : 
     818             : /* 2^n = shifti(gen_1, n) */
     819             : GEN
     820    64630434 : int2n(long n) {
     821             :   long i, m, l;
     822             :   GEN z;
     823    64630434 :   if (n < 0) return gen_0;
     824    64628179 :   if (n == 0) return gen_1;
     825             : 
     826    64007807 :   l = dvmdsBIL(n, &m) + 3;
     827    64007155 :   z = cgetipos(l);
     828    64003879 :   for (i = 2; i < l; i++) z[i] = 0;
     829    64003879 :   *int_MSW(z) = 1UL << m; return z;
     830             : }
     831             : /* To avoid problems when 2^(BIL-1) < n. Overflow cleanly, where int2n
     832             :  * returns gen_0 */
     833             : GEN
     834    12551441 : int2u(ulong n) {
     835             :   ulong i, m, l;
     836             :   GEN z;
     837    12551441 :   if (n == 0) return gen_1;
     838             : 
     839    12551441 :   l = dvmduBIL(n, &m) + 3;
     840    12551434 :   z = cgetipos(l);
     841    12551341 :   for (i = 2; i < l; i++) z[i] = 0;
     842    12551341 :   *int_MSW(z) = 1UL << m; return z;
     843             : }
     844             : /* 2^n - 1 */
     845             : GEN
     846        9772 : int2um1(ulong n) {
     847             :   ulong i, m, l;
     848             :   GEN z;
     849        9772 :   if (n == 0) return gen_0;
     850             : 
     851        9772 :   l = dvmduBIL(n, &m);
     852        9772 :   l += m? 3: 2;
     853        9772 :   z = cgetipos(l);
     854        9772 :   for (i = 2; i < l; i++) z[i] = ~0UL;
     855        9772 :   if (m) *int_MSW(z) = (1UL << m) - 1;
     856        9772 :   return z;
     857             : }
     858             : 
     859             : GEN
     860  1087243063 : shifti(GEN x, long n)
     861             : {
     862  1087243063 :   long s = signe(x);
     863             :   GEN y;
     864             : 
     865  1087243063 :   if(s == 0) return gen_0;
     866   948980147 :   y = shiftispec(x + 2, lgefint(x) - 2, n);
     867   948738902 :   if (signe(y)) setsigne(y, s);
     868   948738902 :   return y;
     869             : }
     870             : 
     871             : /* actual operations will take place on a+2 and b+2: we strip the codewords */
     872             : GEN
     873  8947021187 : mulii(GEN a,GEN b)
     874             : {
     875             :   long sa,sb;
     876             :   GEN z;
     877             : 
     878  8947021187 :   sa=signe(a); if (!sa) return gen_0;
     879  6915495967 :   sb=signe(b); if (!sb) return gen_0;
     880  4882493685 :   if (sb<0) sa = -sa;
     881  4882493685 :   z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
     882  4884520953 :   setsigne(z,sa); return z;
     883             : }
     884             : 
     885             : GEN
     886  1181485082 : sqri(GEN a) { return sqrispec(a+2, lgefint(a)-2); }
     887             : 
     888             : /* sqrt()'s result may be off by 1 when a is not representable exactly as a
     889             :  * double [64bit machine] */
     890             : ulong
     891   103641972 : usqrt(ulong a)
     892             : {
     893   103641972 :   ulong x = (ulong)sqrt((double)a);
     894             : #ifdef LONG_IS_64BIT
     895    88801054 :   if (x > LOWMASK || x*x > a) x--;
     896             : #endif
     897   103641972 :   return x;
     898             : }
     899             : 
     900             : /********************************************************************/
     901             : /**                                                                **/
     902             : /**              EXPONENT / CONVERSION t_REAL --> double           **/
     903             : /**                                                                **/
     904             : /********************************************************************/
     905             : 
     906             : #ifdef LONG_IS_64BIT
     907             : long
     908         894 : dblexpo(double x)
     909             : {
     910             :   union { double f; ulong i; } fi;
     911         894 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     912         894 :   const int exp_mid = 0x3ff;/* exponent bias */
     913             : 
     914         894 :   if (x==0.) return -exp_mid;
     915         894 :   fi.f = x;
     916         894 :   return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
     917             : }
     918             : 
     919             : ulong
     920           0 : dblmantissa(double x)
     921             : {
     922             :   union { double f; ulong i; } fi;
     923           0 :   const int expo_len = 11; /* number of bits of exponent */
     924             : 
     925           0 :   if (x==0.) return 0;
     926           0 :   fi.f = x;
     927           0 :   return (fi.i << expo_len) | HIGHBIT;
     928             : }
     929             : 
     930             : GEN
     931     9442988 : dbltor(double x)
     932             : {
     933             :   GEN z;
     934             :   long e;
     935             :   union { double f; ulong i; } fi;
     936     9442988 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     937     9442988 :   const int exp_mid = 0x3ff;/* exponent bias */
     938     9442988 :   const int expo_len = 11; /* number of bits of exponent */
     939             : 
     940     9442988 :   if (x==0.) return real_0_bit(-exp_mid);
     941     9438908 :   fi.f = x; z = cgetr(DEFAULTPREC);
     942             :   {
     943     9438908 :     const ulong a = fi.i;
     944             :     ulong A;
     945     9438908 :     e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
     946     9438908 :     if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
     947     9438908 :     A = a << expo_len;
     948     9438908 :     if (e == -exp_mid)
     949             :     { /* unnormalized values */
     950           0 :       int sh = bfffo(A);
     951           0 :       e -= sh-1;
     952           0 :       z[2] = A << sh;
     953             :     }
     954             :     else
     955     9438908 :       z[2] = HIGHBIT | A;
     956     9438908 :     z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
     957             :   }
     958     9438908 :   return z;
     959             : }
     960             : 
     961             : double
     962    39716111 : rtodbl(GEN x)
     963             : {
     964    39716111 :   long ex,s=signe(x);
     965             :   ulong a;
     966             :   union { double f; ulong i; } fi;
     967    39716111 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     968    39716111 :   const int exp_mid = 0x3ff;/* exponent bias */
     969    39716111 :   const int expo_len = 11; /* number of bits of exponent */
     970             : 
     971    39716111 :   if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
     972             : 
     973             :   /* start by rounding to closest */
     974    39092905 :   a = (x[2] & (HIGHBIT-1)) + 0x400;
     975    39092905 :   if (a & HIGHBIT) { ex++; a=0; }
     976    39092905 :   if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
     977    39093244 :   fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
     978    39093244 :   if (s<0) fi.i |= HIGHBIT;
     979    39093244 :   return fi.f;
     980             : }
     981             : 
     982             : #else /* LONG_IS_64BIT */
     983             : 
     984             : #if   PARI_DOUBLE_FORMAT == 1
     985             : #  define INDEX0 1
     986             : #  define INDEX1 0
     987             : #elif PARI_DOUBLE_FORMAT == 0
     988             : #  define INDEX0 0
     989             : #  define INDEX1 1
     990             : #endif
     991             : 
     992             : long
     993         149 : dblexpo(double x)
     994             : {
     995             :   union { double f; ulong i[2]; } fi;
     996         149 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     997         149 :   const int exp_mid = 0x3ff;/* exponent bias */
     998         149 :   const int shift = mant_len-32;
     999             : 
    1000         149 :   if (x==0.) return -exp_mid;
    1001         149 :   fi.f = x;
    1002             :   {
    1003         149 :     const ulong a = fi.i[INDEX0];
    1004         149 :     return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
    1005             :   }
    1006             : }
    1007             : 
    1008             : ulong
    1009           0 : dblmantissa(double x)
    1010             : {
    1011             :   union { double f; ulong i[2]; } fi;
    1012           0 :   const int expo_len = 11; /* number of bits of exponent */
    1013             : 
    1014           0 :   if (x==0.) return 0;
    1015           0 :   fi.f = x;
    1016             :   {
    1017           0 :     const ulong a = fi.i[INDEX0];
    1018           0 :     const ulong b = fi.i[INDEX1];
    1019           0 :     return HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
    1020             :   }
    1021             : }
    1022             : 
    1023             : GEN
    1024     1572660 : dbltor(double x)
    1025             : {
    1026             :   GEN z;
    1027             :   long e;
    1028             :   union { double f; ulong i[2]; } fi;
    1029     1572660 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
    1030     1572660 :   const int exp_mid = 0x3ff;/* exponent bias */
    1031     1572660 :   const int expo_len = 11; /* number of bits of exponent */
    1032     1572660 :   const int shift = mant_len-32;
    1033             : 
    1034     1572660 :   if (x==0.) return real_0_bit(-exp_mid);
    1035     1571992 :   fi.f = x; z = cgetr(DEFAULTPREC);
    1036             :   {
    1037     1571992 :     const ulong a = fi.i[INDEX0];
    1038     1571992 :     const ulong b = fi.i[INDEX1];
    1039             :     ulong A, B;
    1040     1571992 :     e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
    1041     1571992 :     if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
    1042     1571992 :     A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
    1043     1571992 :     B = b << expo_len;
    1044     1571992 :     if (e == -exp_mid)
    1045             :     { /* unnormalized values */
    1046             :       int sh;
    1047           0 :       if (A)
    1048             :       {
    1049           0 :         sh = bfffo(A);
    1050           0 :         e -= sh-1;
    1051           0 :         z[2] = (A << sh) | (B >> (32-sh));
    1052           0 :         z[3] = B << sh;
    1053             :       }
    1054             :       else
    1055             :       {
    1056           0 :         sh = bfffo(B); /* B != 0 */
    1057           0 :         e -= sh-1 + 32;
    1058           0 :         z[2] = B << sh;
    1059           0 :         z[3] = 0;
    1060             :       }
    1061             :     }
    1062             :     else
    1063             :     {
    1064     1571992 :       z[3] = B;
    1065     1571992 :       z[2] = HIGHBIT | A;
    1066             :     }
    1067     1571992 :     z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
    1068             :   }
    1069     1571992 :   return z;
    1070             : }
    1071             : 
    1072             : double
    1073     6401749 : rtodbl(GEN x)
    1074             : {
    1075     6401749 :   long ex,s=signe(x),lx=lg(x);
    1076             :   ulong a,b,k;
    1077             :   union { double f; ulong i[2]; } fi;
    1078     6401749 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
    1079     6401749 :   const int exp_mid = 0x3ff;/* exponent bias */
    1080     6401749 :   const int expo_len = 11; /* number of bits of exponent */
    1081     6401749 :   const int shift = mant_len-32;
    1082             : 
    1083     6401749 :   if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
    1084             : 
    1085             :   /* start by rounding to closest */
    1086     6298168 :   a = x[2] & (HIGHBIT-1);
    1087     6298168 :   if (lx > 3)
    1088             :   {
    1089     6168231 :     b = x[3] + 0x400UL; if (b < 0x400UL) a++;
    1090     6168231 :     if (a & HIGHBIT) { ex++; a=0; }
    1091             :   }
    1092      129937 :   else b = 0;
    1093     6298168 :   if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
    1094     6298168 :   ex += exp_mid;
    1095     6298168 :   k = (a >> expo_len) | (ex << shift);
    1096     6298168 :   if (s<0) k |= HIGHBIT;
    1097     6298168 :   fi.i[INDEX0] = k;
    1098     6298168 :   fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
    1099     6298168 :   return fi.f;
    1100             : }
    1101             : #endif /* LONG_IS_64BIT */
    1102             : 

Generated by: LCOV version 1.13