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 - ratlift.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.0 lcov report (development 29789-eb252ba62d) Lines: 98 108 90.7 %
Date: 2024-12-03 09:08:48 Functions: 3 3 100.0 %
Legend: Lines: hit not hit

          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; either version 2 of the License, or (at your option) any later
       9             : version. It is distributed in the hope that it will be useful, but WITHOUT
      10             : ANY WARRANTY WHATSOEVER.
      11             : 
      12             : Check the License for details. You should have received a copy of it, along
      13             : with the package; see the file 'COPYING'. If not, write to the Free Software
      14             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      15             : 
      16             : /*==========================================================
      17             :  * Fp_ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax)
      18             :  *==========================================================
      19             :  * Reconstruct rational number from its residue x mod m
      20             :  *    Given t_INT x, m, amax>=0, bmax>0 such that
      21             :  *         0 <= x < m;  2*amax*bmax < m
      22             :  *    attempts to find t_INT a, b such that
      23             :  *         (1) a = b*x (mod m)
      24             :  *         (2) |a| <= amax, 0 < b <= bmax
      25             :  *         (3) gcd(m, b) = gcd(a, b)
      26             :  *    If unsuccessful, it will return 0 and leave a,b unchanged  (and
      27             :  *    caller may deduce no such a,b exist).  If successful, sets a,b
      28             :  *    and returns 1.  If there exist a,b satisfying (1), (2), and
      29             :  *         (3') gcd(m, b) = 1
      30             :  *    then they are uniquely determined subject to (1),(2) and
      31             :  *         (3'') gcd(a, b) = 1,
      32             :  *    and will be returned by the routine.  (The caller may wish to
      33             :  *    check gcd(a,b)==1, either directly or based on known prime
      34             :  *    divisors of m, depending on the application.)
      35             :  * Reference:
      36             :  @article {MR97c:11116,
      37             :      AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.},
      38             :       TITLE = {Efficient rational number reconstruction},
      39             :     JOURNAL = {J. Symbolic Comput.},
      40             :      VOLUME = {20},
      41             :        YEAR = {1995},
      42             :      NUMBER = {3},
      43             :       PAGES = {287--297},
      44             :  }
      45             :  * Preprint available from:
      46             :  * ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz */
      47             : static ulong
      48      109416 : get_vmax(GEN r, long lb, long lbb)
      49             : {
      50      109416 :   long lr = lb - lgefint(r);
      51             :   ulong vmax;
      52      109416 :   if (lr > 1)        /* still more than a word's worth to go */
      53           0 :     vmax = ULONG_MAX;        /* (cannot in fact happen) */
      54             :   else
      55             :   { /* take difference of bit lengths */
      56      109416 :     long lbr = bfffo(*int_MSW(r));
      57      109416 :     lr = lr*BITS_IN_LONG - lbb + lbr;
      58      109416 :     if ((ulong)lr > BITS_IN_LONG)
      59        7678 :       vmax = ULONG_MAX;
      60      101738 :     else if (lr == 0)
      61        8663 :       vmax = 1UL;
      62             :     else
      63       93075 :       vmax = 1UL << (lr-1); /* pessimistic but faster than a division */
      64             :   }
      65      109416 :   return vmax;
      66             : }
      67             : 
      68             : /* assume bmax <= sqrt(m), fast if amax <=sqrt(m) */
      69             : static int
      70     4817081 : Fp_ratlift_hgcd(GEN n, GEN m, GEN amax, GEN bmax, GEN *pa, GEN *pb)
      71             : {
      72     4817081 :   pari_sp av = avma;
      73             :   GEN x, y, a, b;
      74     4817081 :   GEN H = halfgcdii(n, m), M = gel(H,1), V = gel(H,2);
      75     4817434 :   x = gel(V,1); a = gel(V,2); y = gcoeff(M,1,1); b = gcoeff(M,2,1);
      76     5005797 :   while(abscmpii(b, bmax)<=0)
      77             :   {
      78             :     GEN q, r, u;
      79     4689264 :     if (abscmpii(a, amax)<=0)
      80             :     {
      81     4500770 :       if (signe(b)<0)  { a = negi(a); b = negi(b); }
      82     4500503 :       *pa =a; *pb = b;
      83     4500503 :       gerepileall(av, 2, pb, pa); return 1;
      84             :     }
      85      188351 :     q = dvmdii(x, a, &r); x = a; a = r;
      86      188366 :     u = subii(y, mulii(b, q));
      87      188363 :     y = b; b = u;
      88             :   }
      89      316375 :   return gc_bool(av, 0);
      90             : }
      91             : 
      92             : /* Assume x,m,amax >= 0,bmax > 0 are t_INTs, 0 <= x < m, 2 amax * bmax < m */
      93             : int
      94    21135959 : Fp_ratlift(GEN x, GEN m, GEN amax, GEN bmax, GEN *a, GEN *b)
      95             : {
      96             :   GEN d, d1, v, v1, q, r;
      97    21135959 :   pari_sp av = avma, av1;
      98             :   long lb, lbb, s, s0;
      99             :   ulong vmax;
     100             :   ulong xu, xu1, xv, xv1; /* Lehmer stage recurrence matrix */
     101             :   int lhmres;             /* Lehmer stage return value */
     102             : 
     103             :   /* special cases x=0 and/or amax=0 */
     104    21135959 :   if (!signe(x)) { *a = gen_0; *b = gen_1; return 1; }
     105     8877250 :   if (!signe(amax)) return 0;
     106             :   /* assert: m > x > 0, amax > 0 */
     107             : 
     108             :   /* check whether a=x, b=1 is a solution */
     109     8877250 :   if (cmpii(x,amax) <= 0) { *a = icopy(x); *b = gen_1; return 1; }
     110             : 
     111     4902979 :   if (amax == bmax || equalii(amax, bmax))
     112     4816745 :     return Fp_ratlift_hgcd(x, m, amax, bmax, a, b);
     113             : 
     114             :   /* There is no special case for single-word numbers since this is
     115             :    * mainly meant to be used with large moduli. */
     116       86234 :   (void)new_chunk(lgefint(bmax) + lgefint(amax)); /* room for a,b */
     117       86234 :   d = m; d1 = x;
     118       86234 :   v = gen_0; v1 = gen_1;
     119             :   /* assert d1 > amax, v1 <= bmax here */
     120       86234 :   lb = lgefint(bmax);
     121       86234 :   lbb = bfffo(*int_MSW(bmax));
     122       86234 :   s = 1;
     123       86234 :   av1 = avma;
     124             : 
     125             :   /* General case: Euclidean division chain starting with m div x, and
     126             :    * with bounds on the sequence of convergents' denoms v_j.
     127             :    * Just to be different from what invmod and bezout are doing, we work
     128             :    * here with the all-nonnegative matrices [u,u1;v,v1]=prod_j([0,1;1,q_j]).
     129             :    * Loop invariants:
     130             :    * (a) (sign)*[-v,v1]*x = [d,d1] (mod m)  (componentwise)
     131             :    * (sign initially +1, changes with each Euclidean step)
     132             :    * so [a,b] will be obtained in the form [-+d,v] or [+-d1,v1];
     133             :    * this congruence is a consequence of
     134             :    *
     135             :    * (b) [x,m]~ = [u,u1;v,v1]*[d1,d]~,
     136             :    * where u,u1 is the usual numerator sequence starting with 1,0
     137             :    * instead of 0,1  (just multiply the eqn on the left by the inverse
     138             :    * matrix, which is det*[v1,-u1;-v,u], where "det" is the same as the
     139             :    * "(sign)" in (a)).  From m = v*d1 + v1*d and
     140             :    *
     141             :    * (c) d > d1 >= 0, 0 <= v < v1,
     142             :    * we have d >= m/(2*v1), so while v1 remains smaller than m/(2*amax),
     143             :    * the pair [-(sign)*d,v] satisfies (1) but violates (2) (d > amax).
     144             :    * Conversely, v1 > bmax indicates that no further solutions will be
     145             :    * forthcoming;  [-(sign)*d,v] will be the last, and first, candidate.
     146             :    * Thus there's at most one point in the chain division where a solution
     147             :    * can live:  v < bmax, v1 >= m/(2*amax) > bmax,  and this is acceptable
     148             :    * iff in fact d <= amax  (e.g. m=221, x=34 or 35, amax=bmax=10 fail on
     149             :    * this count while x=32,33,36,37 succeed).  However, a division may leave
     150             :    * a zero residue before we ever reach this point  (consider m=210, x=35,
     151             :    * amax=bmax=10),  and our caller may find that gcd(d,v) > 1  (Examples:
     152             :    * keep m=210 and consider any of x=29,31,32,33,34,36,37,38,39,40,41).
     153             :    * Furthermore, at the start of the loop body we have in fact
     154             :    *
     155             :    * (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0,
     156             :    * (and are never done already).
     157             :    *
     158             :    * Main loop is similar to those of invmod() and bezout(), except for
     159             :    * having to determine appropriate vmax bounds, and checking termination
     160             :    * conditions.  The signe(d1) condition is only for paranoia */
     161      109361 :   while (lgefint(d) > 3 && signe(d1))
     162             :   {
     163             :     /* determine vmax for lgcdii so as to ensure v won't overshoot.
     164             :      * If v+v1 > bmax, the next step would take v1 beyond the limit, so
     165             :      * since [+-d1,v1] is not a solution, we give up.  Otherwise if v+v1
     166             :      * is way shorter than bmax, use vmax=ULONG_MAX.  Otherwise, set vmax
     167             :      * to a crude lower approximation of bmax/(v+v1), or to 1, which will
     168             :      * allow the inner loop to do one step */
     169       68721 :     r = addii(v,v1);
     170       68721 :     if (cmpii(r,bmax) > 0) return gc_long(av, 0); /* done, not found */
     171       68228 :     vmax = get_vmax(r, lb, lbb);
     172             :     /* do a Lehmer-Jebelean round */
     173       68228 :     lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, vmax);
     174       68228 :     if (lhmres) /* check progress */
     175             :     { /* apply matrix */
     176       68225 :       if (lhmres == 1 || lhmres == -1)
     177             :       {
     178         767 :         s = -s;
     179         767 :         if (xv1 == 1)
     180             :         { /* re-use v+v1 computed above */
     181           2 :           v = v1; v1 = r;
     182           2 :           r = subii(d,d1); d = d1; d1 = r;
     183             :         }
     184             :         else
     185             :         {
     186         765 :           r = subii(d, mului(xv1,d1)); d = d1; d1 = r;
     187         765 :           r = addii(v, mului(xv1,v1)); v = v1; v1 = r;
     188             :         }
     189             :       }
     190             :       else
     191             :       {
     192       67458 :         r  = subii(muliu(d,xu),  muliu(d1,xv));
     193       67458 :         d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
     194       67458 :         r  = addii(muliu(v,xu),  muliu(v1,xv));
     195       67458 :         v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
     196       67458 :         if (lhmres&1) { togglesign(d); s = -s; } else togglesign(d1);
     197             :       }
     198             :       /* check whether we're done.  Assert v <= bmax here.  Examine v1:
     199             :        * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
     200             :        * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise proceed*/
     201       68225 :       if (cmpii(v1,bmax) > 0)
     202             :       {
     203       31940 :         set_avma(av);
     204       31940 :         if (cmpii(d,amax) > 0) return 0; /* done, not found */
     205             :         /* done, found */
     206       16932 :         *a = icopy(d); setsigne(*a,-s);
     207       16932 :         *b = icopy(v); return 1;
     208             :       }
     209       36285 :       if (cmpii(d1,amax) <= 0)
     210             :       { /* done, found */
     211       13161 :         set_avma(av);
     212       13161 :         if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
     213       13161 :         *b = icopy(v1); return 1;
     214             :       }
     215             :     } /* lhmres != 0 */
     216             : 
     217       23127 :     if (lhmres <= 0 && signe(d1))
     218             :     {
     219           5 :       q = dvmdii(d,d1,&r);
     220           5 :       d = d1; d1 = r;
     221           5 :       r = addii(v, mulii(q,v1));
     222           5 :       v = v1; v1 = r;
     223           5 :       s = -s;
     224             :       /* check whether we are done now.  Since we weren't before the div, it
     225             :        * suffices to examine v1 and d1 -- the new d (former d1) cannot cut it */
     226           5 :       if (cmpii(v1,bmax) > 0) return gc_long(av, 0); /* done, not found */
     227           5 :       if (cmpii(d1,amax) <= 0) /* done, found */
     228             :       {
     229           0 :         set_avma(av);
     230           0 :         if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
     231           0 :         *b = icopy(v1); return 1;
     232             :       }
     233             :     }
     234             : 
     235       23127 :     if (gc_needed(av,1))
     236             :     {
     237           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ratlift");
     238           0 :       gerepileall(av1, 4, &d, &d1, &v, &v1);
     239             :     }
     240             :   } /* end while */
     241             : 
     242             :   /* Postprocessing - final sprint.  Since we usually underestimate vmax,
     243             :    * this function needs a loop here instead of a simple conditional.
     244             :    * Note we can only get here when amax fits into one word  (which will
     245             :    * typically not be the case!).  The condition is bogus -- d1 is never
     246             :    * zero at the start of the loop.  There will be at most a few iterations,
     247             :    * so we don't bother collecting garbage */
     248       45770 :   while (signe(d1))
     249             :   {
     250             :     /* Assertions: lgefint(d)==lgefint(d1)==3.
     251             :      * Moreover, we aren't done already, or we would have returned by now.
     252             :      * Recompute vmax */
     253       45770 :     r = addii(v,v1);
     254       45768 :     if (cmpii(r,bmax) > 0) return gc_long(av, 0); /* done, not found */
     255       41188 :     vmax = get_vmax(r, lb, lbb);
     256             :     /* single-word "Lehmer", discarding the gcd or whatever it returns */
     257       41188 :     (void)rgcduu((ulong)*int_MSW(d), (ulong)*int_MSW(d1), vmax, &xu, &xu1, &xv, &xv1, &s0);
     258       41190 :     if (xv1 == 1) /* avoid multiplications */
     259             :     { /* re-use r = v+v1 computed above */
     260           0 :       v = v1; v1 = r;
     261           0 :       r = subii(d,d1); d = d1; d1 = r;
     262           0 :       s = -s;
     263             :     }
     264       41190 :     else if (xu == 0) /* and xv==1, xu1==1, xv1 > 1 */
     265             :     {
     266        7529 :       r = subii(d, mului(xv1,d1)); d = d1; d1 = r;
     267        7529 :       r = addii(v, mului(xv1,v1)); v = v1; v1 = r;
     268        7528 :       s = -s;
     269             :     }
     270             :     else
     271             :     {
     272       33661 :       r  = subii(muliu(d,xu),  muliu(d1,xv));
     273       33661 :       d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
     274       33661 :       r  = addii(muliu(v,xu),  muliu(v1,xv));
     275       33661 :       v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
     276       33661 :       if (s0 < 0) { togglesign(d); s = -s; } else togglesign(d1);
     277             :     }
     278             :     /* check whether we're done, as above.  Assert v <= bmax.
     279             :      * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
     280             :      * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise proceed.
     281             :      */
     282       41189 :     if (cmpii(v1,bmax) > 0)
     283             :     {
     284       31368 :       set_avma(av);
     285       31368 :       if (cmpii(d,amax) > 0) return 0; /* done, not found */
     286             :       /* done, found */
     287       21430 :       *a = icopy(d); setsigne(*a,-s);
     288       21430 :       *b = icopy(v); return 1;
     289             :     }
     290        9821 :     if (cmpii(d1,amax) <= 0)
     291             :     { /* done, found */
     292        4692 :       set_avma(av);
     293        4692 :       if (signe(d1)) { *a = icopy(d1); setsigne(*a,s); } else *a = gen_0;
     294        4692 :       *b = icopy(v1); return 1;
     295             :     }
     296             :   } /* while */
     297             : 
     298             :   /* We have run into d1 == 0 before returning. This cannot happen */
     299           0 :   pari_err_BUG("ratlift failed to catch d1 == 0");
     300             :   return 0; /* LCOV_EXCL_LINE */
     301             : }

Generated by: LCOV version 1.16