Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - kernel/none - mp_indep.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 20935-1ada2f1) Lines: 521 566 92.0 %
Date: 2017-08-23 06:24:09 Functions: 51 53 96.2 %
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    20686527 : invmod2BIL(ulong b)
      18             : {
      19             :   static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
      20    20686527 :   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    20686527 :   x = x*(2-b*x); /* one more pass necessary */
      25             : #endif
      26    20686527 :   x = x*(2-b*x);
      27    20686527 :   x = x*(2-b*x); return x*(2-b*x);
      28             : }
      29             : 
      30             : void
      31  1004276194 : affrr(GEN x, GEN y)
      32             : {
      33             :   long lx,ly,i;
      34             : 
      35  1004276194 :   y[1] = x[1]; if (!signe(x)) return;
      36             : 
      37   998800814 :   lx=lg(x); ly=lg(y);
      38   998800814 :   if (lx <= ly)
      39             :   {
      40   939369220 :     for (i=2; i<lx; i++) y[i]=x[i];
      41   939369220 :     for (   ; i<ly; i++) y[i]=0;
      42   939369220 :     return;
      43             :   }
      44    59431594 :   for (i=2; i<ly; i++) y[i]=x[i];
      45             :   /* lx > ly: round properly */
      46    59431594 :   if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
      47             : }
      48             : 
      49             : GEN
      50     7050626 : trunc2nr(GEN x, long n)
      51             : {
      52             :   long ex;
      53     7050626 :   if (!signe(x)) return gen_0;
      54     6999632 :   ex = expo(x) + n; if (ex < 0) return gen_0;
      55     6639430 :   return mantissa2nr(x, ex - bit_prec(x) + 1);
      56             : }
      57             : 
      58             : /* x a t_REAL, x = i/2^e, i a t_INT */
      59             : GEN
      60    15127509 : mantissa_real(GEN x, long *e)
      61             : {
      62    15127509 :   *e = bit_prec(x)-1-expo(x);
      63    15127509 :   return mantissa2nr(x, 0);
      64             : }
      65             : 
      66             : GEN
      67   901546604 : mului(ulong x, GEN y)
      68             : {
      69   901546604 :   long s = signe(y);
      70             :   GEN z;
      71             : 
      72   901546604 :   if (!s || !x) return gen_0;
      73   845857939 :   z = muluispec(x, y+2, lgefint(y)-2);
      74   845268276 :   setsigne(z,s); return z;
      75             : }
      76             : 
      77             : GEN
      78   209622374 : mulsi(long x, GEN y)
      79             : {
      80   209622374 :   long s = signe(y);
      81             :   GEN z;
      82             : 
      83   209622374 :   if (!s || !x) return gen_0;
      84   161005308 :   if (x<0) { s = -s; x = -x; }
      85   161005308 :   z = muluispec((ulong)x, y+2, lgefint(y)-2);
      86   161005287 :   setsigne(z,s); return z;
      87             : }
      88             : 
      89             : GEN
      90   177089294 : mulss(long x, long y)
      91             : {
      92             :   long p1;
      93             :   LOCAL_HIREMAINDER;
      94             : 
      95   177089294 :   if (!x || !y) return gen_0;
      96   176991658 :   if (x<0) {
      97       59724 :     x = -x;
      98       59724 :     if (y<0) { y = -y; p1 = mulll(x,y); return uutoi(hiremainder, p1); }
      99       59724 :     p1 = mulll(x,y); return uutoineg(hiremainder, p1);
     100             :   } else {
     101   176931934 :     if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
     102   176931289 :     p1 = mulll(x,y); return uutoi(hiremainder, p1);
     103             :   }
     104             : }
     105             : GEN
     106      888000 : sqrs(long x)
     107             : {
     108             :   long p1;
     109             :   LOCAL_HIREMAINDER;
     110             : 
     111      888000 :   if (!x) return gen_0;
     112      887678 :   if (x<0) x = -x;
     113      887678 :   p1 = mulll(x,x); return uutoi(hiremainder, p1);
     114             : }
     115             : GEN
     116  1973146143 : muluu(ulong x, ulong y)
     117             : {
     118             :   long p1;
     119             :   LOCAL_HIREMAINDER;
     120             : 
     121  1973146143 :   if (!x || !y) return gen_0;
     122  1973099650 :   p1 = mulll(x,y); return uutoi(hiremainder, p1);
     123             : }
     124             : GEN
     125   635679581 : sqru(ulong x)
     126             : {
     127             :   long p1;
     128             :   LOCAL_HIREMAINDER;
     129             : 
     130   635679581 :   if (!x) return gen_0;
     131   635339702 :   p1 = mulll(x,x); return uutoi(hiremainder, p1);
     132             : }
     133             : 
     134             : /* assume x > 1, y != 0. Return u * y with sign s */
     135             : static GEN
     136   477288804 : mulur_2(ulong x, GEN y, long s)
     137             : {
     138   477288804 :   long m, sh, i, lx = lg(y), e = expo(y);
     139   477288804 :   GEN z = cgetr(lx);
     140             :   ulong garde;
     141             :   LOCAL_HIREMAINDER;
     142             : 
     143   477288804 :   y--; garde = mulll(x,y[lx]);
     144   477288804 :   for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
     145   477288804 :   z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
     146   477288804 :   sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
     147   477288804 :   if (sh) shift_left(z,z, 2,lx-1, garde,sh);
     148   477288804 :   z[1] = evalsigne(s) | evalexpo(m+e);
     149   477288804 :   if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
     150   477288804 :   return z;
     151             : }
     152             : 
     153             : INLINE GEN
     154     1168586 : mul0r(GEN x)
     155             : {
     156     1168586 :   long l = lg(x), e = expo(x);
     157     1168586 :   e = (l > 2)? -prec2nbits(l) + e: (e < 0? 2*e: 0);
     158     1168586 :   return real_0_bit(e);
     159             : }
     160             : /* lg(x) > 2 */
     161             : INLINE GEN
     162       65008 : div0r(GEN x) {
     163       65008 :   long l = lg(x), e = expo(x);
     164       65008 :   return real_0_bit(-prec2nbits(l) - e);
     165             : }
     166             : 
     167             : GEN
     168    39257267 : mulsr(long x, GEN y)
     169             : {
     170             :   long s;
     171             : 
     172    39257267 :   if (!x) return mul0r(y);
     173    38239874 :   s = signe(y);
     174    38239874 :   if (!s)
     175             :   {
     176      219059 :     if (x < 0) x = -x;
     177      219059 :     return real_0_bit( expo(y) + expu(x) );
     178             :   }
     179    38020815 :   if (x==1)  return rcopy(y);
     180    37270471 :   if (x==-1) return negr(y);
     181    36279754 :   if (x < 0)
     182    20964269 :     return mulur_2((ulong)-x, y, -s);
     183             :   else
     184    15315485 :     return mulur_2((ulong)x, y, s);
     185             : }
     186             : 
     187             : GEN
     188   487867284 : mulur(ulong x, GEN y)
     189             : {
     190             :   long s;
     191             : 
     192   487867284 :   if (!x) return mul0r(y);
     193   487867284 :   s = signe(y);
     194   487867284 :   if (!s) return real_0_bit( expo(y) + expu(x) );
     195   486076919 :   if (x==1) return rcopy(y);
     196   441009050 :   return mulur_2(x, y, s);
     197             : }
     198             : 
     199             : INLINE void
     200  1022642190 : mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
     201             : {
     202             :   long i;
     203  1022642190 :   if (hi[2] < 0)
     204             :   {
     205   465017476 :     if (z != hi)
     206    30077134 :       for (i=2; i<lz ; i++) z[i] = hi[i];
     207   465017476 :     ez++;
     208             :   }
     209             :   else
     210             :   {
     211   557624714 :     shift_left(z,hi,2,lz-1, garde, 1);
     212   557624714 :     garde <<= 1;
     213             :   }
     214  1022642190 :   if (garde & HIGHBIT)
     215             :   { /* round to nearest */
     216   499409309 :     i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
     217   499409309 :     if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
     218             :   }
     219  1022642190 :   z[1] = evalsigne(sz)|evalexpo(ez);
     220  1022642190 : }
     221             : /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
     222             : INLINE void
     223   731537025 : mulrrz_3end(GEN z, long sz, long ez, ulong garde)
     224             : {
     225   731537025 :   if (z[2] < 0)
     226             :   { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
     227   394978282 :     if (garde & HIGHBIT) z[2]++; /* round properly */
     228   394978282 :     ez++;
     229             :   }
     230             :   else
     231             :   {
     232   336558743 :     uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
     233   336558743 :     if (garde & (1UL<<(BITS_IN_LONG-2)))
     234             :     {
     235   155024797 :       uel(z,2)++; /* round properly, z2+1 can overflow */
     236   155024797 :       if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
     237             :     }
     238             :   }
     239   731537025 :   z[1] = evalsigne(sz)|evalexpo(ez);
     240   731537025 : }
     241             : 
     242             : /* set z <-- x^2 != 0, floating point multiplication.
     243             :  * lz = lg(z) = lg(x) */
     244             : INLINE void
     245    77326050 : sqrz_i(GEN z, GEN x, long lz)
     246             : {
     247    77326050 :   long ez = 2*expo(x);
     248             :   long i, j, lzz, p1;
     249             :   ulong garde;
     250             :   GEN x1;
     251             :   LOCAL_HIREMAINDER;
     252             :   LOCAL_OVERFLOW;
     253             : 
     254    77326050 :   if (lz > MULRR_MULII_LIMIT)
     255             :   {
     256     4808987 :     pari_sp av = avma;
     257     4808987 :     GEN hi = sqrispec_mirror(x+2, lz-2);
     258     4808987 :     mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
     259     4808987 :     avma = av; return;
     260             :   }
     261    72517063 :   if (lz == 3)
     262             :   {
     263    13005368 :     garde = mulll(x[2],x[2]);
     264    13005368 :     z[2] = hiremainder;
     265    13005368 :     mulrrz_3end(z, 1, ez, garde);
     266    13005368 :     return;
     267             :   }
     268             : 
     269    59511695 :   lzz = lz-1; p1 = x[lzz];
     270    59511695 :   if (p1)
     271             :   {
     272    52633827 :     (void)mulll(p1,x[3]);
     273    52633827 :     garde = addmul(p1,x[2]);
     274    52633827 :     z[lzz] = hiremainder;
     275             :   }
     276             :   else
     277             :   {
     278     6877868 :     garde = 0;
     279     6877868 :     z[lzz] = 0;
     280             :   }
     281   200191516 :   for (j=lz-2, x1=x-j; j>=3; j--)
     282             :   {
     283   140679821 :     p1 = x[j]; x1++;
     284   140679821 :     if (p1)
     285             :     {
     286   135336605 :       (void)mulll(p1,x1[lz+1]);
     287   135336605 :       garde = addll(addmul(p1,x1[lz]), garde);
     288   618532239 :       for (i=lzz; i>j; i--)
     289             :       {
     290   483195634 :         hiremainder += overflow;
     291   483195634 :         z[i] = addll(addmul(p1,x1[i]), z[i]);
     292             :       }
     293   135336605 :       z[j] = hiremainder+overflow;
     294             :     }
     295     5343216 :     else z[j]=0;
     296             :   }
     297    59511695 :   p1 = x[2]; x1++;
     298    59511695 :   garde = addll(mulll(p1,x1[lz]), garde);
     299   259703211 :   for (i=lzz; i>2; i--)
     300             :   {
     301   200191516 :     hiremainder += overflow;
     302   200191516 :     z[i] = addll(addmul(p1,x1[i]), z[i]);
     303             :   }
     304    59511695 :   z[2] = hiremainder+overflow;
     305    59511695 :   mulrrz_end(z, z, lz, 1, ez, garde);
     306             : }
     307             : 
     308             : /* lz "large" = lg(y) = lg(z), lg(x) > lz if flag = 1 and >= if flag = 0 */
     309             : INLINE void
     310    60659031 : mulrrz_int(GEN z, GEN x, GEN y, long lz, long flag, long sz)
     311             : {
     312    60659031 :   pari_sp av = avma;
     313    60659031 :   GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
     314    60659031 :   mulrrz_end(z, hi, lz, sz, expo(x)+expo(y), hi[lz]);
     315    60659031 :   avma = av;
     316    60659031 : }
     317             : 
     318             : /* lz = 3 */
     319             : INLINE void
     320   718531657 : mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
     321             : {
     322             :   ulong garde;
     323             :   LOCAL_HIREMAINDER;
     324   718531657 :   if (flag)
     325             :   {
     326    14480229 :     (void)mulll(x[2],y[3]);
     327    14480229 :     garde = addmul(x[2],y[2]);
     328             :   }
     329             :   else
     330   704051428 :     garde = mulll(x[2],y[2]);
     331   718531657 :   z[2] = hiremainder;
     332   718531657 :   mulrrz_3end(z, sz, expo(x)+expo(y), garde);
     333   718531657 : }
     334             : 
     335             : /* set z <-- x*y, floating point multiplication. Trailing 0s for x are
     336             :  * treated efficiently (important application: mulir).
     337             :  * lz = lg(z) = lg(x) <= ly <= lg(y), sz = signe(z). flag = lg(x) < lg(y) */
     338             : INLINE void
     339  1675432003 : mulrrz_i(GEN z, GEN x, GEN y, long lz, long flag, long sz)
     340             : {
     341             :   long ez, i, j, lzz, p1;
     342             :   ulong garde;
     343             :   GEN y1;
     344             :   LOCAL_HIREMAINDER;
     345             :   LOCAL_OVERFLOW;
     346             : 
     347  1675432003 :   if (x == y) { sqrz_i(z,x,lz); return; }
     348  1675432003 :   if (lz > MULRR_MULII_LIMIT) { mulrrz_int(z,x,y,lz,flag,sz); return; }
     349  1614772972 :   if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
     350   896241315 :   ez = expo(x) + expo(y);
     351   896241315 :   if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
     352   896241315 :   lzz=lz-1; p1=x[lzz];
     353   896241315 :   if (p1)
     354             :   {
     355   862182559 :     (void)mulll(p1,y[3]);
     356   862182559 :     garde = addll(addmul(p1,y[2]), garde);
     357   862182559 :     z[lzz] = overflow+hiremainder;
     358             :   }
     359    34058756 :   else z[lzz]=0;
     360  2856770291 :   for (j=lz-2, y1=y-j; j>=3; j--)
     361             :   {
     362  1960528976 :     p1 = x[j]; y1++;
     363  1960528976 :     if (p1)
     364             :     {
     365  1891355383 :       (void)mulll(p1,y1[lz+1]);
     366  1891355383 :       garde = addll(addmul(p1,y1[lz]), garde);
     367  8035428430 :       for (i=lzz; i>j; i--)
     368             :       {
     369  6144073047 :         hiremainder += overflow;
     370  6144073047 :         z[i] = addll(addmul(p1,y1[i]), z[i]);
     371             :       }
     372  1891355383 :       z[j] = hiremainder+overflow;
     373             :     }
     374    69173593 :     else z[j]=0;
     375             :   }
     376   896241315 :   p1 = x[2]; y1++;
     377   896241315 :   garde = addll(mulll(p1,y1[lz]), garde);
     378  3753011606 :   for (i=lzz; i>2; i--)
     379             :   {
     380  2856770291 :     hiremainder += overflow;
     381  2856770291 :     z[i] = addll(addmul(p1,y1[i]), z[i]);
     382             :   }
     383   896241315 :   z[2] = hiremainder+overflow;
     384   896241315 :   mulrrz_end(z, z, lz, sz, ez, garde);
     385             : }
     386             : 
     387             : GEN
     388  1691422723 : mulrr(GEN x, GEN y)
     389             : {
     390             :   long flag, ly, lz, sx, sy;
     391             :   GEN z;
     392             : 
     393  1691422723 :   if (x == y) return sqrr(x);
     394  1691278608 :   sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
     395  1666753841 :   sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
     396  1657438938 :   if (sy < 0) sx = -sx;
     397  1657438938 :   lz = lg(x);
     398  1657438938 :   ly = lg(y);
     399  1657438938 :   if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
     400  1657438938 :   z = cgetr(lz);
     401  1657438938 :   mulrrz_i(z, x,y, lz,flag, sx);
     402  1657438938 :   return z;
     403             : }
     404             : 
     405             : GEN
     406    78139407 : sqrr(GEN x)
     407             : {
     408    78139407 :   long lz, sx = signe(x);
     409             :   GEN z;
     410             : 
     411    78139407 :   if (!sx) return real_0_bit(2*expo(x));
     412    77326050 :   lz = lg(x); z = cgetr(lz);
     413    77326050 :   sqrz_i(z, x, lz);
     414    77326050 :   return z;
     415             : }
     416             : 
     417             : GEN
     418   159490025 : mulir(GEN x, GEN y)
     419             : {
     420   159490025 :   long sx = signe(x), sy;
     421   159490025 :   if (!sx) return mul0r(y);
     422   159338832 :   if (lgefint(x) == 3) {
     423   138859259 :     GEN z = mulur(uel(x,2), y);
     424   138859259 :     if (sx < 0) togglesign(z);
     425   138859259 :     return z;
     426             :   }
     427    20479573 :   sy = signe(y);
     428    20479573 :   if (!sy) return real_0_bit(expi(x) + expo(y));
     429    19414227 :   if (sy < 0) sx = -sx;
     430             :   {
     431    19414227 :     long lz = lg(y), lx = lgefint(x);
     432    19414227 :     GEN hi, z = cgetr(lz);
     433    19414227 :     pari_sp av = avma;
     434    19414227 :     if (lx < (lz>>1) || (lx < lz && lz > MULRR_MULII_LIMIT))
     435             :     { /* size mantissa of x < half size of mantissa z, or lx < lz so large
     436             :        * that mulrr will call mulii anyway: mulii */
     437     1421162 :       x = itor(x, lx);
     438     1421162 :       hi = muliispec_mirror(y+2, x+2, lz-2, lx-2);
     439     1421162 :       mulrrz_end(z, hi, lz, sx, expo(x)+expo(y), hi[lz]);
     440             :     }
     441             :     else /* dubious: complete x with 0s and call mulrr */
     442    17993065 :       mulrrz_i(z, itor(x,lz), y, lz, 0, sx);
     443    19414227 :     avma = av; return z;
     444             :   }
     445             : }
     446             : 
     447             : /* x + y*z, generic. If lgefint(z) <= 3, caller should use faster variants  */
     448             : static GEN
     449    28618240 : addmulii_gen(GEN x, GEN y, GEN z, long lz)
     450             : {
     451    28618240 :   long lx = lgefint(x), ly;
     452             :   pari_sp av;
     453             :   GEN t;
     454    28618240 :   if (lx == 2) return mulii(z,y);
     455    27061531 :   ly = lgefint(y);
     456    27061531 :   if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
     457    27051858 :   av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
     458    27051858 :   t = mulii(z, y);
     459    27051858 :   avma = av; return addii(t,x);
     460             : }
     461             : /* x + y*z, lgefint(z) == 3 */
     462             : static GEN
     463   375309205 : addmulii_lg3(GEN x, GEN y, GEN z)
     464             : {
     465   375309205 :   long s = signe(z), lx, ly;
     466   375309205 :   ulong w = z[2];
     467             :   pari_sp av;
     468             :   GEN t;
     469   375309205 :   if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
     470   293353428 :   lx = lgefint(x);
     471   293353428 :   ly = lgefint(y);
     472   293353428 :   if (lx == 2)
     473             :   { /* x = 0 */
     474    31733121 :     if (ly == 2) return gen_0;
     475    26809378 :     t = muluispec(w, y+2, ly-2);
     476    26809378 :     if (signe(y) < 0) s = -s;
     477    26809378 :     setsigne(t, s); return t;
     478             :   }
     479   261620307 :   if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
     480   256286405 :   av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
     481   256286405 :   t = muluispec(w, y+2, ly-2);
     482   256286405 :   if (signe(y) < 0) s = -s;
     483   256286405 :   setsigne(t, s);
     484   256286405 :   avma = av; return addii(x,t);
     485             : }
     486             : /* x + y*z */
     487             : GEN
     488    29383225 : addmulii(GEN x, GEN y, GEN z)
     489             : {
     490    29383225 :   long lz = lgefint(z);
     491    29383225 :   switch(lz)
     492             :   {
     493      259487 :     case 2: return icopy(x); /* z = 0, wasteful copy */
     494    27047861 :     case 3: return addmulii_lg3(x, y, z);
     495     2075877 :     default:return addmulii_gen(x, y, z, lz);
     496             :   }
     497             : }
     498             : /* x + y*z, returns x itself and not a copy when y*z = 0 */
     499             : GEN
     500  1763246131 : addmulii_inplace(GEN x, GEN y, GEN z)
     501             : {
     502             :   long lz;
     503  1763246131 :   if (lgefint(y) == 2) return x;
     504   374803707 :   lz = lgefint(z);
     505   374803707 :   switch(lz)
     506             :   {
     507           0 :     case 2: return x;
     508   348261344 :     case 3: return addmulii_lg3(x, y, z);
     509    26542363 :     default:return addmulii_gen(x, y, z, lz);
     510             :   }
     511             : }
     512             : 
     513             : /* written by Bruno Haible following an idea of Robert Harley */
     514             : long
     515  2593150130 : vals(ulong z)
     516             : {
     517             :   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};
     518             : #ifdef LONG_IS_64BIT
     519             :   long s;
     520             : #endif
     521             : 
     522  2593150130 :   if (!z) return -1;
     523             : #ifdef LONG_IS_64BIT
     524  2271774872 :   if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
     525             : #endif
     526  2593150130 :   z |= ~z + 1;
     527  2593150130 :   z += z << 4;
     528  2593150130 :   z += z << 6;
     529  2593150130 :   z ^= z << 16; /* or  z -= z<<16 */
     530             : #ifdef LONG_IS_64BIT
     531  2271774872 :   return s + tab[(z&0xffffffff)>>26];
     532             : #else
     533   321375258 :   return tab[z>>26];
     534             : #endif
     535             : }
     536             : 
     537             : GEN
     538          21 : divsi(long x, GEN y)
     539             : {
     540          21 :   long p1, s = signe(y);
     541             :   LOCAL_HIREMAINDER;
     542             : 
     543          21 :   if (!s) pari_err_INV("divsi",gen_0);
     544          21 :   if (!x || lgefint(y)>3 || ((long)y[2])<0) return gen_0;
     545          21 :   hiremainder=0; p1=divll(labs(x),y[2]);
     546          21 :   if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
     547          21 :   if (s<0) p1 = -p1;
     548          21 :   return stoi(p1);
     549             : }
     550             : 
     551             : GEN
     552     2500461 : divir(GEN x, GEN y)
     553             : {
     554             :   GEN z;
     555     2500461 :   long ly = lg(y), lx = lgefint(x);
     556             :   pari_sp av;
     557             : 
     558     2500461 :   if (ly == 2) pari_err_INV("divir",y);
     559     2500461 :   if (lx == 2) return div0r(y);
     560     2435453 :   if (lx == 3) {
     561     1269029 :     z = divur(x[2], y);
     562     1269029 :     if (signe(x) < 0) togglesign(z);
     563     1269029 :     return z;
     564             :   }
     565     1166424 :   z = cgetr(ly); av = avma;
     566     1166424 :   affrr(divrr(itor(x, ly+1), y), z);
     567     1166424 :   avma = av; return z;
     568             : }
     569             : 
     570             : GEN
     571     1302532 : divur(ulong x, GEN y)
     572             : {
     573             :   pari_sp av;
     574     1302532 :   long ly = lg(y);
     575             :   GEN z;
     576             : 
     577     1302532 :   if (ly == 2) pari_err_INV("divur",y);
     578     1302532 :   if (!x) return div0r(y);
     579     1302532 :   if (ly > INVNEWTON_LIMIT) {
     580         571 :     av = avma; z = invr(y);
     581         571 :     if (x == 1) return z;
     582         571 :     return gerepileuptoleaf(av, mulur(x, z));
     583             :   }
     584     1301961 :   z = cgetr(ly); av = avma;
     585     1301961 :   affrr(divrr(utor(x,ly+1), y), z);
     586     1301961 :   avma = av; return z;
     587             : }
     588             : 
     589             : GEN
     590        7324 : divsr(long x, GEN y)
     591             : {
     592             :   pari_sp av;
     593        7324 :   long ly = lg(y);
     594             :   GEN z;
     595             : 
     596        7324 :   if (ly == 2) pari_err_INV("divsr",y);
     597        7324 :   if (!x) return div0r(y);
     598        7324 :   if (ly > INVNEWTON_LIMIT) {
     599           0 :     av = avma; z = invr(y);
     600           0 :     if (x == 1) return z;
     601           0 :     if (x ==-1) { togglesign(z); return z; }
     602           0 :     return gerepileuptoleaf(av, mulsr(x, z));
     603             :   }
     604        7324 :   z = cgetr(ly); av = avma;
     605        7324 :   affrr(divrr(stor(x,ly+1), y), z);
     606        7324 :   avma = av; return z;
     607             : }
     608             : 
     609             : /* returns 1/y, assume y != 0 */
     610             : static GEN
     611    11414103 : invr_basecase(GEN y)
     612             : {
     613    11414103 :   long ly = lg(y);
     614    11414103 :   GEN z = cgetr(ly);
     615    11414103 :   pari_sp av = avma;
     616    11414103 :   affrr(divrr(real_1(ly+1), y), z);
     617    11414103 :   avma = av; return z;
     618             : }
     619             : /* returns 1/b, Newton iteration */
     620             : GEN
     621    11414103 : invr(GEN b)
     622             : {
     623    11414103 :   const long s = 6;
     624    11414103 :   long i, p, l = lg(b);
     625             :   GEN x, a;
     626             :   ulong mask;
     627             : 
     628    11414103 :   if (l <= maxss(INVNEWTON_LIMIT, (1L<<s) + 2)) {
     629    11410609 :     if (l == 2) pari_err_INV("invr",b);
     630    11410609 :     return invr_basecase(b);
     631             :   }
     632        3494 :   mask = quadratic_prec_mask(l-2);
     633        3494 :   for(i=0, p=1; i<s; i++) { p <<= 1; if (mask & 1) p--; mask >>= 1; }
     634        3494 :   x = cgetr(l);
     635        3494 :   a = rcopy(b); a[1] = _evalexpo(0) | evalsigne(1);
     636        3494 :   affrr(invr_basecase(rtor(a, p+2)), x);
     637       12295 :   while (mask > 1)
     638             :   {
     639        5307 :     p <<= 1; if (mask & 1) p--;
     640        5307 :     mask >>= 1;
     641        5307 :     setlg(a, p + 2);
     642        5307 :     setlg(x, p + 2);
     643             :     /* TODO: mulrr(a,x) should be a half product (the higher half is known).
     644             :      * mulrr(x, ) already is */
     645        5307 :     affrr(addrr(x, mulrr(x, subsr(1, mulrr(a,x)))), x);
     646        5307 :     avma = (pari_sp)a;
     647             :   }
     648        3494 :   x[1] = (b[1] & SIGNBITS) | evalexpo(expo(x)-expo(b));
     649        3494 :   avma = (pari_sp)x; return x;
     650             : }
     651             : 
     652             : GEN
     653  1343829598 : modii(GEN x, GEN y)
     654             : {
     655  1343829598 :   switch(signe(x))
     656             :   {
     657   169732373 :     case 0: return gen_0;
     658   799778686 :     case 1: return remii(x,y);
     659             :     default:
     660             :     {
     661   374318539 :       pari_sp av = avma;
     662   374318539 :       (void)new_chunk(lgefint(y));
     663   374315024 :       x = remii(x,y); avma=av;
     664   374319711 :       if (x==gen_0) return x;
     665   352482173 :       return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
     666             :     }
     667             :   }
     668             : }
     669             : 
     670             : void
     671      229877 : modiiz(GEN x, GEN y, GEN z)
     672             : {
     673      229877 :   const pari_sp av = avma;
     674      229877 :   affii(modii(x,y),z); avma=av;
     675      229877 : }
     676             : 
     677             : GEN
     678     3507817 : divrs(GEN x, long y)
     679             : {
     680             :   GEN z;
     681     3507817 :   if (y < 0)
     682             :   {
     683      479394 :     z = divru(x, (ulong)-y);
     684      479394 :     togglesign(z);
     685             :   }
     686             :   else
     687     3028423 :     z = divru(x, (ulong)y);
     688     3507817 :   return z;
     689             : }
     690             : 
     691             : GEN
     692   538649317 : divru(GEN x, ulong y)
     693             : {
     694   538649317 :   long i, lx, sh, e, s = signe(x);
     695             :   ulong garde;
     696             :   GEN z;
     697             :   LOCAL_HIREMAINDER;
     698             : 
     699   538649317 :   if (!y) pari_err_INV("divru",gen_0);
     700   538649317 :   if (!s) return real_0_bit(expo(x) - expu(y));
     701   538585436 :   if (!(y & (y-1))) /* power of 2 */
     702             :   {
     703    36190018 :     if (y == 1) return rcopy(x);
     704    34187275 :     return shiftr(x, -expu(y));
     705             :   }
     706   502395418 :   e = expo(x);
     707   502395418 :   lx = lg(x);
     708   502395418 :   z = cgetr(lx);
     709   502395418 :   if (lx == 3)
     710             :   {
     711   313665902 :     if (y <= uel(x,2))
     712             :     {
     713   313665842 :       hiremainder = 0;
     714   313665842 :       z[2] = divll(x[2],y);
     715             :       /* we may have hiremainder != 0 ==> garde */
     716   313665842 :       garde = divll(0,y);
     717             :     }
     718             :     else
     719             :     {
     720          60 :       hiremainder = x[2];
     721          60 :       z[2] = divll(0,y);
     722          60 :       garde = hiremainder;
     723          60 :       e -= BITS_IN_LONG;
     724             :     }
     725             :   }
     726             :   else
     727             :   {
     728   188729516 :     ulong yp = get_Fl_red(y);
     729   188729516 :     if (y <= uel(x,2))
     730             :     {
     731   188726707 :       hiremainder = 0;
     732   188726707 :       for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
     733             :       /* we may have hiremainder != 0 ==> garde */
     734   188726707 :       garde = divll_pre(0,y,yp);
     735             :     }
     736             :     else
     737             :     {
     738        2809 :       long l = lx-1;
     739        2809 :       hiremainder = x[2];
     740        2809 :       for (i=2; i<l; i++) z[i] = divll_pre(x[i+1],y,yp);
     741        2809 :       z[i] = divll_pre(0,y,yp);
     742        2809 :       garde = hiremainder;
     743        2809 :       e -= BITS_IN_LONG;
     744             :     }
     745             :   }
     746   502395418 :   sh=bfffo(z[2]); /* z[2] != 0 */
     747   502395418 :   if (sh) shift_left(z,z, 2,lx-1, garde,sh);
     748   502395418 :   z[1] = evalsigne(s) | evalexpo(e-sh);
     749   502395418 :   if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
     750   502395418 :   return z;
     751             : }
     752             : 
     753             : GEN
     754   452315100 : truedvmdii(GEN x, GEN y, GEN *z)
     755             : {
     756             :   pari_sp av;
     757             :   GEN r, q, *gptr[2];
     758   452315100 :   if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
     759     1487241 :   if (z == ONLY_REM) return modii(x,y);
     760             : 
     761     1487241 :   av = avma;
     762     1487241 :   q = dvmdii(x,y,&r); /* assume that r is last on stack */
     763     1487241 :   switch(signe(r))
     764             :   {
     765             :     case 0:
     766      750291 :       if (z) *z = gen_0;
     767      750291 :       return q;
     768             :     case 1:
     769      576497 :       if (z) *z = r; else cgiv(r);
     770      576497 :       return q;
     771      160453 :     case -1: break;
     772             :   }
     773      160453 :   q = addis(q, -signe(y));
     774      160453 :   if (!z) return gerepileuptoint(av, q);
     775             : 
     776       19556 :   *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
     777       19556 :   gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(pari_sp)r,gptr,2);
     778       19556 :   return q;
     779             : }
     780             : GEN
     781   450845369 : truedvmdis(GEN x, long y, GEN *z)
     782             : {
     783   450845369 :   pari_sp av = avma;
     784             :   long r;
     785             :   GEN q;
     786             : 
     787   450845369 :   if (z == ONLY_REM) return modis(x, y);
     788   450845369 :   q = divis_rem(x,y,&r);
     789             : 
     790   450845369 :   if (r >= 0)
     791             :   {
     792   444724750 :     if (z) *z = utoi(r);
     793   444724750 :     return q;
     794             :   }
     795     6120619 :   q = gerepileuptoint(av, addis(q, (y < 0)? 1: -1));
     796     6120619 :   if (z) *z = utoi(r + labs(y));
     797     6120619 :   return q;
     798             : }
     799             : GEN
     800     6084595 : truedvmdsi(long x, GEN y, GEN *z)
     801             : {
     802             :   long q, r;
     803     6084595 :   if (z == ONLY_REM) return modsi(x, y);
     804     6084595 :   q = sdivsi_rem(x,y,&r);
     805     6084595 :   if (r >= 0) {
     806     6084567 :     if (z) *z = utoi(r);
     807     6084567 :     return stoi(q);
     808             :   }
     809          28 :   q = q - signe(y);
     810          28 :   if (!z) return stoi(q);
     811             : 
     812           0 :   *z = subiuspec(y+2,(ulong)-r, lgefint(y)-2);
     813           0 :   return stoi(q);
     814             : }
     815             : 
     816             : /* 2^n = shifti(gen_1, n) */
     817             : GEN
     818    74659436 : int2n(long n) {
     819             :   long i, m, l;
     820             :   GEN z;
     821    74659436 :   if (n < 0) return gen_0;
     822    74658668 :   if (n == 0) return gen_1;
     823             : 
     824    73572077 :   l = dvmdsBIL(n, &m) + 3;
     825    73571903 :   z = cgetipos(l);
     826    73571979 :   for (i = 2; i < l; i++) z[i] = 0;
     827    73571979 :   *int_MSW(z) = 1UL << m; return z;
     828             : }
     829             : /* To avoid problems when 2^(BIL-1) < n. Overflow cleanly, where int2n
     830             :  * returns gen_0 */
     831             : GEN
     832     9998801 : int2u(ulong n) {
     833             :   ulong i, m, l;
     834             :   GEN z;
     835     9998801 :   if (n == 0) return gen_1;
     836             : 
     837     9998799 :   l = dvmduBIL(n, &m) + 3;
     838     9998804 :   z = cgetipos(l);
     839     9998793 :   for (i = 2; i < l; i++) z[i] = 0;
     840     9998793 :   *int_MSW(z) = 1UL << m; return z;
     841             : }
     842             : /* 2^n - 1 */
     843             : GEN
     844         476 : int2um1(ulong n) {
     845             :   ulong i, m, l;
     846             :   GEN z;
     847         476 :   if (n == 0) return gen_0;
     848             : 
     849         476 :   l = dvmduBIL(n, &m);
     850         476 :   l += m? 3: 2;
     851         476 :   z = cgetipos(l);
     852         476 :   for (i = 2; i < l; i++) z[i] = ~0UL;
     853         476 :   if (m) *int_MSW(z) = (1UL << m) - 1;
     854         476 :   return z;
     855             : }
     856             : 
     857             : GEN
     858   845299713 : shifti(GEN x, long n)
     859             : {
     860   845299713 :   long s = signe(x);
     861             :   GEN y;
     862             : 
     863   845299713 :   if(s == 0) return gen_0;
     864   769682289 :   y = shiftispec(x + 2, lgefint(x) - 2, n);
     865   769628362 :   if (signe(y)) setsigne(y, s);
     866   769628362 :   return y;
     867             : }
     868             : 
     869             : /* actual operations will take place on a+2 and b+2: we strip the codewords */
     870             : GEN
     871  8486616273 : mulii(GEN a,GEN b)
     872             : {
     873             :   long sa,sb;
     874             :   GEN z;
     875             : 
     876  8486616273 :   sa=signe(a); if (!sa) return gen_0;
     877  6103990916 :   sb=signe(b); if (!sb) return gen_0;
     878  4112208242 :   if (sb<0) sa = -sa;
     879  4112208242 :   z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
     880  4112797986 :   setsigne(z,sa); return z;
     881             : }
     882             : 
     883             : GEN
     884  1456233610 : sqri(GEN a) { return sqrispec(a+2, lgefint(a)-2); }
     885             : 
     886             : /* sqrt()'s result may be off by 1 when a is not representable exactly as a
     887             :  * double [64bit machine] */
     888             : ulong
     889    73667566 : usqrt(ulong a)
     890             : {
     891    73667566 :   ulong x = (ulong)sqrt((double)a);
     892             : #ifdef LONG_IS_64BIT
     893    63113145 :   if (x > LOWMASK || x*x > a) x--;
     894             : #endif
     895    73667566 :   return x;
     896             : }
     897             : 
     898             : /********************************************************************/
     899             : /**                                                                **/
     900             : /**              EXPONENT / CONVERSION t_REAL --> double           **/
     901             : /**                                                                **/
     902             : /********************************************************************/
     903             : 
     904             : #ifdef LONG_IS_64BIT
     905             : long
     906           0 : dblexpo(double x)
     907             : {
     908             :   union { double f; ulong i; } fi;
     909           0 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     910           0 :   const int exp_mid = 0x3ff;/* exponent bias */
     911             : 
     912           0 :   if (x==0.) return -exp_mid;
     913           0 :   fi.f = x;
     914           0 :   return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
     915             : }
     916             : 
     917             : ulong
     918           0 : dblmantissa(double x)
     919             : {
     920             :   union { double f; ulong i; } fi;
     921           0 :   const int expo_len = 11; /* number of bits of exponent */
     922             : 
     923           0 :   if (x==0.) return 0;
     924           0 :   fi.f = x;
     925           0 :   return (fi.i << expo_len) | HIGHBIT;
     926             : }
     927             : 
     928             : GEN
     929     5674328 : dbltor(double x)
     930             : {
     931             :   GEN z;
     932             :   long e;
     933             :   union { double f; ulong i; } fi;
     934     5674328 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     935     5674328 :   const int exp_mid = 0x3ff;/* exponent bias */
     936     5674328 :   const int expo_len = 11; /* number of bits of exponent */
     937             : 
     938     5674328 :   if (x==0.) return real_0_bit(-exp_mid);
     939     5662694 :   fi.f = x; z = cgetr(DEFAULTPREC);
     940             :   {
     941     5662694 :     const ulong a = fi.i;
     942             :     ulong A;
     943     5662694 :     e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
     944     5662694 :     if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
     945     5662694 :     A = a << expo_len;
     946     5662694 :     if (e == -exp_mid)
     947             :     { /* unnormalized values */
     948           0 :       int sh = bfffo(A);
     949           0 :       e -= sh-1;
     950           0 :       z[2] = A << sh;
     951             :     }
     952             :     else
     953     5662694 :       z[2] = HIGHBIT | A;
     954     5662694 :     z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
     955             :   }
     956     5662694 :   return z;
     957             : }
     958             : 
     959             : double
     960    25604676 : rtodbl(GEN x)
     961             : {
     962    25604676 :   long ex,s=signe(x);
     963             :   ulong a;
     964             :   union { double f; ulong i; } fi;
     965    25604676 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     966    25604676 :   const int exp_mid = 0x3ff;/* exponent bias */
     967    25604676 :   const int expo_len = 11; /* number of bits of exponent */
     968             : 
     969    25604676 :   if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
     970             : 
     971             :   /* start by rounding to closest */
     972    25123333 :   a = (x[2] & (HIGHBIT-1)) + 0x400;
     973    25123333 :   if (a & HIGHBIT) { ex++; a=0; }
     974    25123333 :   if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
     975    25123333 :   fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
     976    25123333 :   if (s<0) fi.i |= HIGHBIT;
     977    25123333 :   return fi.f;
     978             : }
     979             : 
     980             : #else /* LONG_IS_64BIT */
     981             : 
     982             : #if   PARI_DOUBLE_FORMAT == 1
     983             : #  define INDEX0 1
     984             : #  define INDEX1 0
     985             : #elif PARI_DOUBLE_FORMAT == 0
     986             : #  define INDEX0 0
     987             : #  define INDEX1 1
     988             : #endif
     989             : 
     990             : long
     991           0 : dblexpo(double x)
     992             : {
     993             :   union { double f; ulong i[2]; } fi;
     994           0 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     995           0 :   const int exp_mid = 0x3ff;/* exponent bias */
     996           0 :   const int shift = mant_len-32;
     997             : 
     998           0 :   if (x==0.) return -exp_mid;
     999           0 :   fi.f = x;
    1000             :   {
    1001           0 :     const ulong a = fi.i[INDEX0];
    1002           0 :     return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
    1003             :   }
    1004             : }
    1005             : 
    1006             : ulong
    1007           0 : dblmantissa(double x)
    1008             : {
    1009             :   union { double f; ulong i[2]; } fi;
    1010           0 :   const int expo_len = 11; /* number of bits of exponent */
    1011             : 
    1012           0 :   if (x==0.) return 0;
    1013           0 :   fi.f = x;
    1014             :   {
    1015           0 :     const ulong a = fi.i[INDEX0];
    1016           0 :     const ulong b = fi.i[INDEX1];
    1017           0 :     return HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
    1018             :   }
    1019             : }
    1020             : 
    1021             : GEN
    1022      944393 : dbltor(double x)
    1023             : {
    1024             :   GEN z;
    1025             :   long e;
    1026             :   union { double f; ulong i[2]; } fi;
    1027      944393 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
    1028      944393 :   const int exp_mid = 0x3ff;/* exponent bias */
    1029      944393 :   const int expo_len = 11; /* number of bits of exponent */
    1030      944393 :   const int shift = mant_len-32;
    1031             : 
    1032      944393 :   if (x==0.) return real_0_bit(-exp_mid);
    1033      943074 :   fi.f = x; z = cgetr(DEFAULTPREC);
    1034             :   {
    1035      943074 :     const ulong a = fi.i[INDEX0];
    1036      943074 :     const ulong b = fi.i[INDEX1];
    1037             :     ulong A, B;
    1038      943074 :     e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
    1039      943074 :     if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
    1040      943074 :     A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
    1041      943074 :     B = b << expo_len;
    1042      943074 :     if (e == -exp_mid)
    1043             :     { /* unnormalized values */
    1044             :       int sh;
    1045           0 :       if (A)
    1046             :       {
    1047           0 :         sh = bfffo(A);
    1048           0 :         e -= sh-1;
    1049           0 :         z[2] = (A << sh) | (B >> (32-sh));
    1050           0 :         z[3] = B << sh;
    1051             :       }
    1052             :       else
    1053             :       {
    1054           0 :         sh = bfffo(B); /* B != 0 */
    1055           0 :         e -= sh-1 + 32;
    1056           0 :         z[2] = B << sh;
    1057           0 :         z[3] = 0;
    1058             :       }
    1059             :     }
    1060             :     else
    1061             :     {
    1062      943074 :       z[3] = B;
    1063      943074 :       z[2] = HIGHBIT | A;
    1064             :     }
    1065      943074 :     z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
    1066             :   }
    1067      943074 :   return z;
    1068             : }
    1069             : 
    1070             : double
    1071     4158195 : rtodbl(GEN x)
    1072             : {
    1073     4158195 :   long ex,s=signe(x),lx=lg(x);
    1074             :   ulong a,b,k;
    1075             :   union { double f; ulong i[2]; } fi;
    1076     4158195 :   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
    1077     4158195 :   const int exp_mid = 0x3ff;/* exponent bias */
    1078     4158195 :   const int expo_len = 11; /* number of bits of exponent */
    1079     4158195 :   const int shift = mant_len-32;
    1080             : 
    1081     4158195 :   if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
    1082             : 
    1083             :   /* start by rounding to closest */
    1084     4078024 :   a = x[2] & (HIGHBIT-1);
    1085     4078024 :   if (lx > 3)
    1086             :   {
    1087     3852156 :     b = x[3] + 0x400UL; if (b < 0x400UL) a++;
    1088     3852156 :     if (a & HIGHBIT) { ex++; a=0; }
    1089             :   }
    1090      225868 :   else b = 0;
    1091     4078024 :   if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
    1092     4078024 :   ex += exp_mid;
    1093     4078024 :   k = (a >> expo_len) | (ex << shift);
    1094     4078024 :   if (s<0) k |= HIGHBIT;
    1095     4078024 :   fi.i[INDEX0] = k;
    1096     4078024 :   fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
    1097     4078024 :   return fi.f;
    1098             : }
    1099             : #endif /* LONG_IS_64BIT */
    1100             : 

Generated by: LCOV version 1.11