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-bordeaux1.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 - ratlift.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16624-25b9976) Lines: 88 94 93.6 %
Date: 2014-06-24 Functions: 2 2 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 78 90 86.7 %

           Branch data     Line data    Source code
       1                 :            : #line 2 "../src/kernel/none/ratlift.c"
       2                 :            : /* Copyright (C) 2003  The PARI group.
       3                 :            : 
       4                 :            : This file is part of the PARI/GP package.
       5                 :            : 
       6                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       7                 :            : terms of the GNU General Public License as published by the Free Software
       8                 :            : Foundation. 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                 :            : /*==========================================================
      16                 :            :  * Fp_ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax)
      17                 :            :  *==========================================================
      18                 :            :  * Reconstruct rational number from its residue x mod m
      19                 :            :  *    Given t_INT x, m, amax>=0, bmax>0 such that
      20                 :            :  *         0 <= x < m;  2*amax*bmax < m
      21                 :            :  *    attempts to find t_INT a, b such that
      22                 :            :  *         (1) a = b*x (mod m)
      23                 :            :  *         (2) |a| <= amax, 0 < b <= bmax
      24                 :            :  *         (3) gcd(m, b) = gcd(a, b)
      25                 :            :  *    If unsuccessful, it will return 0 and leave a,b unchanged  (and
      26                 :            :  *    caller may deduce no such a,b exist).  If successful, sets a,b
      27                 :            :  *    and returns 1.  If there exist a,b satisfying (1), (2), and
      28                 :            :  *         (3') gcd(m, b) = 1
      29                 :            :  *    then they are uniquely determined subject to (1),(2) and
      30                 :            :  *         (3'') gcd(a, b) = 1,
      31                 :            :  *    and will be returned by the routine.  (The caller may wish to
      32                 :            :  *    check gcd(a,b)==1, either directly or based on known prime
      33                 :            :  *    divisors of m, depending on the application.)
      34                 :            :  * Reference:
      35                 :            :  @article {MR97c:11116,
      36                 :            :      AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.},
      37                 :            :       TITLE = {Efficient rational number reconstruction},
      38                 :            :     JOURNAL = {J. Symbolic Comput.},
      39                 :            :      VOLUME = {20},
      40                 :            :        YEAR = {1995},
      41                 :            :      NUMBER = {3},
      42                 :            :       PAGES = {287--297},
      43                 :            :  }
      44                 :            :  * Preprint available from:
      45                 :            :  * ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz */
      46                 :            : 
      47                 :            : /* #define DEBUG_RATLIFT */
      48                 :            : static ulong
      49                 :      99169 : get_vmax(GEN r, long lb, long lbb)
      50                 :            : {
      51                 :      99169 :   long lr = lb - lgefint(r);
      52                 :            :   ulong vmax;
      53         [ +  + ]:      99169 :   if (lr > 1)        /* still more than a word's worth to go */
      54                 :       5833 :     vmax = ULONG_MAX;        /* (cannot in fact happen) */
      55                 :            :   else
      56                 :            :   { /* take difference of bit lengths */
      57 [ +  + ][ +  + ]:      93336 :     long lbr = bfffo(*int_MSW(r));
         [ +  + ][ +  + ]
      58                 :      93336 :     lr = lr*BITS_IN_LONG - lbb + lbr;
      59 [ +  + ][ +  + ]:      93336 :     if ((ulong)lr > BITS_IN_LONG)
      60                 :      16798 :       vmax = ULONG_MAX;
      61         [ +  + ]:      76538 :     else if (lr == 0)
      62                 :        628 :       vmax = 1UL;
      63                 :            :     else
      64                 :      75910 :       vmax = 1UL << (lr-1); /* pessimistic but faster than a division */
      65                 :            :   }
      66                 :            : #ifdef DEBUG_RATLIFT
      67                 :            :   err_printf("rl-fs: vmax=%lu\n", vmax);
      68                 :            : #endif
      69                 :      99169 :   return vmax;
      70                 :            : }
      71                 :            : 
      72                 :            : /* Assume x,m,amax >= 0,bmax > 0 are t_INTs, 0 <= x < m, 2 amax * bmax < m */
      73                 :            : int
      74                 :     105802 : Fp_ratlift(GEN x, GEN m, GEN amax, GEN bmax, GEN *a, GEN *b)
      75                 :            : {
      76                 :            :   GEN d, d1, v, v1, q, r;
      77                 :     105802 :   pari_sp av = avma, av1, lim;
      78                 :            :   long lb, lbb, s, s0;
      79                 :            :   ulong vmax;
      80                 :            :   ulong xu, xu1, xv, xv1; /* Lehmer stage recurrence matrix */
      81                 :            :   int lhmres;             /* Lehmer stage return value */
      82                 :            : 
      83                 :            :   /* special cases x=0 and/or amax=0 */
      84         [ +  + ]:     105802 :   if (!signe(x)) { *a = gen_0; *b = gen_1; return 1; }
      85         [ -  + ]:      72090 :   if (!signe(amax)) return 0;
      86                 :            :   /* assert: m > x > 0, amax > 0 */
      87                 :            : 
      88                 :            :   /* check whether a=x, b=1 is a solution */
      89         [ +  + ]:      72090 :   if (cmpii(x,amax) <= 0) { *a = icopy(x); *b = gen_1; return 1; }
      90                 :            : 
      91                 :            :   /* There is no special case for single-word numbers since this is
      92                 :            :    * mainly meant to be used with large moduli. */
      93                 :      63093 :   (void)new_chunk(lgefint(bmax) + lgefint(amax)); /* room for a,b */
      94                 :      63093 :   d = m; d1 = x;
      95                 :      63093 :   v = gen_0; v1 = gen_1;
      96                 :            :   /* assert d1 > amax, v1 <= bmax here */
      97                 :      63093 :   lb = lgefint(bmax);
      98 [ +  + ][ +  + ]:      63093 :   lbb = bfffo(*int_MSW(bmax));
         [ +  + ][ +  + ]
      99                 :      63093 :   s = 1;
     100                 :      63093 :   av1 = avma; lim = stack_lim(av, 1);
     101                 :            : 
     102                 :            :   /* General case: Euclidean division chain starting with m div x, and
     103                 :            :    * with bounds on the sequence of convergents' denoms v_j.
     104                 :            :    * Just to be different from what invmod and bezout are doing, we work
     105                 :            :    * here with the all-nonnegative matrices [u,u1;v,v1]=prod_j([0,1;1,q_j]).
     106                 :            :    * Loop invariants:
     107                 :            :    * (a) (sign)*[-v,v1]*x = [d,d1] (mod m)  (componentwise)
     108                 :            :    * (sign initially +1, changes with each Euclidean step)
     109                 :            :    * so [a,b] will be obtained in the form [-+d,v] or [+-d1,v1];
     110                 :            :    * this congruence is a consequence of
     111                 :            :    *
     112                 :            :    * (b) [x,m]~ = [u,u1;v,v1]*[d1,d]~,
     113                 :            :    * where u,u1 is the usual numerator sequence starting with 1,0
     114                 :            :    * instead of 0,1  (just multiply the eqn on the left by the inverse
     115                 :            :    * matrix, which is det*[v1,-u1;-v,u], where "det" is the same as the
     116                 :            :    * "(sign)" in (a)).  From m = v*d1 + v1*d and
     117                 :            :    *
     118                 :            :    * (c) d > d1 >= 0, 0 <= v < v1,
     119                 :            :    * we have d >= m/(2*v1), so while v1 remains smaller than m/(2*amax),
     120                 :            :    * the pair [-(sign)*d,v] satisfies (1) but violates (2) (d > amax).
     121                 :            :    * Conversely, v1 > bmax indicates that no further solutions will be
     122                 :            :    * forthcoming;  [-(sign)*d,v] will be the last, and first, candidate.
     123                 :            :    * Thus there's at most one point in the chain division where a solution
     124                 :            :    * can live:  v < bmax, v1 >= m/(2*amax) > bmax,  and this is acceptable
     125                 :            :    * iff in fact d <= amax  (e.g. m=221, x=34 or 35, amax=bmax=10 fail on
     126                 :            :    * this count while x=32,33,36,37 succeed).  However, a division may leave
     127                 :            :    * a zero residue before we ever reach this point  (consider m=210, x=35,
     128                 :            :    * amax=bmax=10),  and our caller may find that gcd(d,v) > 1  (Examples:
     129                 :            :    * keep m=210 and consider any of x=29,31,32,33,34,36,37,38,39,40,41).
     130                 :            :    * Furthermore, at the start of the loop body we have in fact
     131                 :            :    *
     132                 :            :    * (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0,
     133                 :            :    * (and are never done already).
     134                 :            :    *
     135                 :            :    * Main loop is similar to those of invmod() and bezout(), except for
     136                 :            :    * having to determine appropriate vmax bounds, and checking termination
     137                 :            :    * conditions.  The signe(d1) condition is only for paranoia */
     138         [ +  + ]:     101130 :   while (lgefint(d) > 3 && signe(d1))
           [ +  -  +  + ]
                 [ +  - ]
     139                 :            :   {
     140                 :            :     /* determine vmax for lgcdii so as to ensure v won't overshoot.
     141                 :            :      * If v+v1 > bmax, the next step would take v1 beyond the limit, so
     142                 :            :      * since [+-d1,v1] is not a solution, we give up.  Otherwise if v+v1
     143                 :            :      * is way shorter than bmax, use vmax=MAXULUNG.  Otherwise, set vmax
     144                 :            :      * to a crude lower approximation of bmax/(v+v1), or to 1, which will
     145                 :            :      * allow the inner loop to do one step */
     146                 :      86082 :     r = addii(v,v1);
     147         [ +  + ]:      86082 :     if (cmpii(r,bmax) > 0) { avma = av; return 0; } /* done, not found */
     148                 :      85545 :     vmax = get_vmax(r, lb, lbb);
     149                 :            :     /* do a Lehmer-Jebelean round */
     150                 :      85545 :     lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, vmax);
     151         [ +  + ]:      85545 :     if (lhmres) /* check progress */
     152                 :            :     { /* apply matrix */
     153 [ +  + ][ +  + ]:      81700 :       if (lhmres == 1 || lhmres == -1)
     154                 :            :       {
     155                 :      20040 :         s = -s;
     156         [ +  + ]:      40080 :         if (xv1 == 1)
     157                 :            :         { /* re-use v+v1 computed above */
     158                 :      10876 :           v = v1; v1 = r;
     159                 :      10876 :           r = subii(d,d1); d = d1; d1 = r;
     160                 :            :         }
     161                 :            :         else
     162                 :            :         {
     163                 :       9164 :           r = subii(d, mului(xv1,d1)); d = d1; d1 = r;
     164                 :       9164 :           r = addii(v, mului(xv1,v1)); v = v1; v1 = r;
     165                 :            :         }
     166                 :            :       }
     167                 :            :       else
     168                 :            :       {
     169                 :      61660 :         r  = subii(muliu(d,xu),  muliu(d1,xv));
     170                 :      61660 :         d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
     171                 :      61660 :         r  = addii(muliu(v,xu),  muliu(v1,xv));
     172                 :      61660 :         v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
     173         [ +  + ]:      61660 :         if (lhmres&1) { togglesign(d); s = -s; } else togglesign(d1);
     174                 :            :       }
     175                 :            :       /* check whether we're done.  Assert v <= bmax here.  Examine v1:
     176                 :            :        * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
     177                 :            :        * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise proceed*/
     178         [ +  + ]:      81700 :       if (cmpii(v1,bmax) > 0)
     179                 :            :       {
     180                 :      23095 :         avma = av;
     181         [ +  + ]:      23095 :         if (cmpii(d,amax) > 0) return 0; /* done, not found */
     182                 :            :         /* done, found */
     183                 :      12066 :         *a = icopy(d); setsigne(*a,-s);
     184                 :      12066 :         *b = icopy(v); return 1;
     185                 :            :       }
     186         [ +  + ]:      58605 :       if (cmpii(d1,amax) <= 0)
     187                 :            :       { /* done, found */
     188                 :      22484 :         avma = av;
     189         [ +  - ]:      22484 :         if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
     190                 :      22484 :         *b = icopy(v1); return 1;
     191                 :            :       }
     192                 :            :     } /* lhmres != 0 */
     193                 :            : 
     194 [ +  + ][ +  - ]:      39966 :     if (lhmres <= 0 && signe(d1))
     195                 :            :     {
     196                 :       3852 :       q = dvmdii(d,d1,&r);
     197                 :            : #ifdef DEBUG_LEHMER
     198                 :            :       err_printf("Full division:\n  q = %Ps\n", q);
     199                 :            : #endif
     200                 :       3852 :       d = d1; d1 = r;
     201                 :       3852 :       r = addii(v, mulii(q,v1));
     202                 :       3852 :       v = v1; v1 = r;
     203                 :       3852 :       s = -s;
     204                 :            :       /* check whether we are done now.  Since we weren't before the div, it
     205                 :            :        * suffices to examine v1 and d1 -- the new d (former d1) cannot cut it */
     206         [ -  + ]:       3852 :       if (cmpii(v1,bmax) > 0) { avma = av; return 0; } /* done, not found */
     207         [ +  + ]:       3852 :       if (cmpii(d1,amax) <= 0) /* done, found */
     208                 :            :       {
     209                 :       1929 :         avma = av;
     210         [ +  - ]:       1929 :         if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
     211                 :       1929 :         *b = icopy(v1); return 1;
     212                 :            :       }
     213                 :            :     }
     214                 :            : 
     215         [ -  + ]:      38037 :     if (low_stack(lim, stack_lim(av,1)))
     216                 :            :     {
     217         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ratlift");
     218                 :          0 :       gerepileall(av1, 4, &d, &d1, &v, &v1);
     219                 :            :     }
     220                 :            :   } /* end while */
     221                 :            : 
     222                 :            :   /* Postprocessing - final sprint.  Since we usually underestimate vmax,
     223                 :            :    * this function needs a loop here instead of a simple conditional.
     224                 :            :    * Note we can only get here when amax fits into one word  (which will
     225                 :            :    * typically not be the case!).  The condition is bogus -- d1 is never
     226                 :            :    * zero at the start of the loop.  There will be at most a few iterations,
     227                 :            :    * so we don't bother collecting garbage */
     228         [ +  - ]:      16052 :   while (signe(d1))
     229                 :            :   {
     230                 :            :     /* Assertions: lgefint(d)==lgefint(d1)==3.
     231                 :            :      * Moreover, we aren't done already, or we would have returned by now.
     232                 :            :      * Recompute vmax */
     233                 :            : #ifdef DEBUG_RATLIFT
     234                 :            :     err_printf("rl-fs: d,d1=%Ps,%Ps\n", d, d1);
     235                 :            :     err_printf("rl-fs: v,v1=%Ps,%Ps\n", v, v1);
     236                 :            : #endif
     237                 :      16052 :     r = addii(v,v1);
     238         [ +  + ]:      16052 :     if (cmpii(r,bmax) > 0) { avma = av; return 0; } /* done, not found */
     239                 :      13624 :     vmax = get_vmax(r, lb, lbb);
     240                 :            :     /* single-word "Lehmer", discarding the gcd or whatever it returns */
     241                 :      13624 :     (void)rgcduu((ulong)*int_MSW(d), (ulong)*int_MSW(d1), vmax, &xu, &xu1, &xv, &xv1, &s0);
     242                 :            : #ifdef DEBUG_RATLIFT
     243                 :            :     err_printf("rl-fs: [%lu,%lu; %lu,%lu] %s\n",
     244                 :            :                xu, xu1, xv, xv1, s0 < 0 ? "-" : "+");
     245                 :            : #endif
     246         [ -  + ]:      13624 :     if (xv1 == 1) /* avoid multiplications */
     247                 :            :     { /* re-use r = v+v1 computed above */
     248                 :          0 :       v = v1; v1 = r;
     249                 :          0 :       r = subii(d,d1); d = d1; d1 = r;
     250                 :          0 :       s = -s;
     251                 :            :     }
     252         [ +  + ]:      13624 :     else if (xu == 0) /* and xv==1, xu1==1, xv1 > 1 */
     253                 :            :     {
     254                 :       2187 :       r = subii(d, mului(xv1,d1)); d = d1; d1 = r;
     255                 :       2187 :       r = addii(v, mului(xv1,v1)); v = v1; v1 = r;
     256                 :       2187 :       s = -s;
     257                 :            :     }
     258                 :            :     else
     259                 :            :     {
     260                 :      11437 :       r  = subii(muliu(d,xu),  muliu(d1,xv));
     261                 :      11437 :       d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
     262                 :      11437 :       r  = addii(muliu(v,xu),  muliu(v1,xv));
     263                 :      11437 :       v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
     264         [ +  + ]:      11437 :       if (s0 < 0) { togglesign(d); s = -s; } else togglesign(d1);
     265                 :            :     }
     266                 :            :     /* check whether we're done, as above.  Assert v <= bmax.
     267                 :            :      * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
     268                 :            :      * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise proceed.
     269                 :            :      */
     270         [ +  + ]:      13624 :     if (cmpii(v1,bmax) > 0)
     271                 :            :     {
     272                 :      12215 :       avma = av;
     273         [ +  + ]:      12215 :       if (cmpii(d,amax) > 0) return 0; /* done, not found */
     274                 :            :       /* done, found */
     275                 :       8229 :       *a = icopy(d); setsigne(*a,-s);
     276                 :       8229 :       *b = icopy(v); return 1;
     277                 :            :     }
     278         [ +  + ]:       1409 :     if (cmpii(d1,amax) <= 0)
     279                 :            :     { /* done, found */
     280                 :        405 :       avma = av;
     281         [ +  + ]:        405 :       if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
     282                 :        405 :       *b = icopy(v1); return 1;
     283                 :            :     }
     284                 :            :   } /* while */
     285                 :            : 
     286                 :            :   /* We have run into d1 == 0 before returning. This cannot happen */
     287                 :          0 :   pari_err_BUG("ratlift failed to catch d1 == 0");
     288                 :     105802 :   return 0; /* NOTREACHED */
     289                 :            : }

Generated by: LCOV version 1.9