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 - modules - part.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16937-4bd9b4e) Lines: 149 183 81.4 %
Date: 2014-10-24 Functions: 12 13 92.3 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 98 149 65.8 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2002  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : /* Original code contributed by: Ralf Stephan (ralf@ark.in-berlin.de).
      15                 :            :  * Updated by Bill Allombert (2014) to use Selberg formula for L
      16                 :            :  * following http://dx.doi.org/10.1112/S1461157012001088
      17                 :            :  *
      18                 :            :  * This program is basically the implementation of the script
      19                 :            :  *
      20                 :            :  * Psi(n, q) = my(a=sqrt(2/3)*Pi/q, b=n-1/24, c=sqrt(b));
      21                 :            :  *             (sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c))
      22                 :            :  * L(n,q)=sqrt(k/3)*sum(l=0,2*k-1,
      23                 :            :           if(((3*l^2+l)/2+n)%k==0,(-1)^l*cos((6*l+1)/(6*k)*Pi)))
      24                 :            :  * part(n) = round(sum(q=1,5 + 0.24*sqrt(n),L(n,q)*Psi(n,q)))
      25                 :            :  *
      26                 :            :  * only faster.
      27                 :            :  *
      28                 :            :  * ------------------------------------------------------------------
      29                 :            :  *   The first restriction depends on Pari's maximum precision of floating
      30                 :            :  * point reals, which is 268435454 bits in 2.2.4, since the algorithm needs
      31                 :            :  * high precision exponentials. For that engine, the maximum possible argument
      32                 :            :  * would be in [5*10^15,10^16], the computation of which would need days on
      33                 :            :  * a ~1-GHz computer. */
      34                 :            : 
      35                 :            : #include "pari.h"
      36                 :            : #include "paripriv.h"
      37                 :            : 
      38                 :            : /****************************************************************/
      39                 :            : 
      40                 :            : /* Given c = sqrt(2/3)*Pi*sqrt(N-1/24)
      41                 :            :  * Psi(N, q) = my(a = c/q); sqrt(q) * (a*cosh(a) - sinh(a)) */
      42                 :            : static GEN
      43                 :    1109212 : psi(GEN c, ulong q, long prec)
      44                 :            : {
      45                 :    1109212 :   GEN a = divru(c, q), ea = mpexp(a), invea = invr(ea);
      46                 :    1109212 :   GEN cha = shiftr(addrr(ea, invea), -1);  /* ch(a) */
      47                 :    1109212 :   GEN sha = shiftr(subrr(ea, invea), -1);  /* sh(a) */
      48                 :    1109212 :   return mulrr(sqrtr(stor(q,prec)), subrr(mulrr(a,cha), sha));
      49                 :            : }
      50                 :            : 
      51                 :            : /* L(n,q)=sqrt(k/3)*sum(l=0,2*k-1,
      52                 :            :           if(((3*l^2+l)/2+n)%k==0,(-1)^l*cos((6*l+1)/(6*k)*Pi)))
      53                 :            :  * Never called with q < 3, so ignore this case */
      54                 :            : static GEN
      55                 :    1295452 : L(GEN n, ulong k, long bitprec)
      56                 :            : {
      57                 :            :   ulong r, l, m;
      58                 :    1295452 :   long pr = nbits2prec(bitprec / k + k);
      59                 :    1295452 :   GEN s = stor(0,pr), pi = mppi(pr);
      60                 :    1295452 :   pari_sp av = avma;
      61                 :            : 
      62                 :    1295452 :   r = 2; m = umodiu(n,k);
      63         [ +  + ]:   34034238 :   for (l = 0; l < 2*k; l++)
      64                 :            :   {
      65         [ +  + ]:   32738786 :     if (m == 0)
      66                 :            :     {
      67                 :    2590332 :       GEN c = mpcos(divru(mulru(pi, 6*l+1), 6*k));
      68         [ +  + ]:    2590332 :       if (odd(l)) subrrz(s, c, s); else addrrz(s, c, s);
      69                 :    2590332 :       avma = av;
      70                 :            :     }
      71         [ +  + ]:   32738786 :     m += r; if (m >= k) m -= k;
      72         [ +  + ]:   32738786 :     r += 3; if (r >= k) r -= k;
      73                 :            :   }
      74                 :            :   /* multiply by sqrt(k/3) */
      75         [ +  + ]:    1295452 :   return mulrr(s, sqrtr((k % 3)? rdivss(k,3,pr): utor(k/3,pr)));
      76                 :            : }
      77                 :            : 
      78                 :            : /* Return a low precision estimate of log p(n). */
      79                 :            : static GEN
      80                 :      70000 : estim(GEN n)
      81                 :            : {
      82                 :      70000 :   pari_sp av = avma;
      83                 :      70000 :   GEN p1, pi = mppi (DEFAULTPREC);
      84                 :            : 
      85                 :      70000 :   p1 = divru( itor(shifti(n,1), DEFAULTPREC), 3 );
      86                 :      70000 :   p1 = mpexp( mulrr(pi, sqrtr(p1)) ); /* exp(Pi * sqrt(2N/3)) */
      87                 :      70000 :   p1 = divri (shiftr(p1,-2), n);
      88                 :      70000 :   p1 = divrr(p1, sqrtr( stor(3,DEFAULTPREC) ));
      89                 :      70000 :   return gerepileupto(av, mplog(p1));
      90                 :            : }
      91                 :            : 
      92                 :            : /* c = sqrt(2/3)*Pi*sqrt(n-1/24);  d = 1 / ((2*b)^(3/2) * Pi); */
      93                 :            : static void
      94                 :      70000 : pinit(GEN n, GEN *c, GEN *d, ulong prec)
      95                 :            : {
      96                 :      70000 :   GEN b = divru( itor( subis(muliu(n,24), 1), prec ), 24 ); /* n - 1/24 */
      97                 :      70000 :   GEN sqrtb = sqrtr(b), Pi = mppi(prec), pi2sqrt2, pisqrt2d3;
      98                 :            : 
      99                 :            : 
     100                 :      70000 :   pisqrt2d3 = mulrr(Pi, sqrtr( divru(stor(2, prec), 3) ));
     101                 :      70000 :   pi2sqrt2  = mulrr(Pi, sqrtr( stor(8, prec) ));
     102                 :      70000 :   *c = mulrr(pisqrt2d3, sqrtb);
     103                 :      70000 :   *d = invr( mulrr(pi2sqrt2, mulrr(b,sqrtb)) );
     104                 :      70000 : }
     105                 :            : 
     106                 :            : /* part(n) = round(sum(q=1,5 + 0.24*sqrt(n), L(n,q)*Psi(n,q))) */
     107                 :            : GEN
     108                 :      70021 : numbpart(GEN n)
     109                 :            : {
     110                 :      70021 :   pari_sp ltop = avma, av;
     111                 :            :   GEN sum, est, C, D, p1, p2;
     112                 :            :   long prec, bitprec;
     113                 :            :   ulong q;
     114                 :            : 
     115         [ -  + ]:      70021 :   if (typ(n) != t_INT) pari_err_TYPE("partition function",n);
     116         [ -  + ]:      70021 :   if (signe(n) < 0) return gen_0;
     117         [ +  + ]:      70021 :   if (cmpiu(n, 2) < 0) return gen_1;
     118         [ +  + ]:      70007 :   if (cmpii(n, uu32toi(0x38d7e, 0xa4c68000)) >= 0)
     119                 :          7 :     pari_err_OVERFLOW("numbpart [n < 10^15]");
     120                 :      70000 :   est = estim(n);
     121                 :      70000 :   bitprec = (long)(rtodbl(est)/LOG2) + 32;
     122                 :      70000 :   prec = nbits2prec(bitprec);
     123                 :      70000 :   pinit(n, &C, &D, prec);
     124                 :      70000 :   sum = cgetr (prec); affsr(0, sum);
     125                 :            :   /* Because N < 10^16 and q < sqrt(N), q fits into a long
     126                 :            :    * In fact q < 2 LONG_MAX / 3 */
     127                 :      70000 :   av = avma; togglesign(est);
     128         [ +  + ]:    1365452 :   for (q = (ulong)(sqrt(gtodouble(n))*0.24 + 5); q >= 3; q--, avma=av)
     129                 :            :   {
     130                 :    1295452 :     GEN t = L(n, q, bitprec);
     131         [ +  + ]:    1295452 :     if (absr_cmp(t, mpexp(divru(est,q))) < 0) continue;
     132                 :            : 
     133                 :     969212 :     t = mulrr(t, psi(gprec_w(C, nbits2prec(bitprec / q + 32)), q, prec));
     134                 :     969212 :     affrr(addrr(sum, t), sum);
     135                 :            :   }
     136                 :      70000 :   p1 = addrr(sum, psi(C, 1, prec));
     137                 :      70000 :   p2 = psi(C, 2, prec);
     138         [ +  + ]:      70000 :   affrr(mod2(n)? subrr(p1,p2): addrr(p1,p2), sum);
     139                 :      70014 :   return gerepileuptoint (ltop, roundr(mulrr(D,sum)));
     140                 :            : }
     141                 :            : 
     142                 :            : /* for loop over partitions of integer k.
     143                 :            :  * nbounds can restrict partitions to have length between nmin and nmax
     144                 :            :  * (the length is the number of non zero entries) and
     145                 :            :  * abounds restrict to integers between amin and amax.
     146                 :            :  *
     147                 :            :  * Start from central partition.
     148                 :            :  * By default, remove zero entries on the left.
     149                 :            :  *
     150                 :            :  * Algorithm:
     151                 :            :  *
     152                 :            :  * A partition of k is an increasing sequence v1,... vn with sum(vi)=k
     153                 :            :  * The starting point is the minimal n-partition of k: a,...a,a+1,.. a+1
     154                 :            :  * (a+1 is repeated r times with k = a * n + r).
     155                 :            :  *
     156                 :            :  * The procedure to obtain the next partition:
     157                 :            :  * - find the last index i<n such that v{i-1} != v{i} (that is vi is the start
     158                 :            :  * of the last constant range excluding vn).
     159                 :            :  * - decrease vi by one, and set v{i+1},... v{n} to be a minimal partition (of
     160                 :            :  * the right sum).
     161                 :            :  *
     162                 :            :  * Examples: we highlight the index i
     163                 :            :  * 1 1 2 2 3
     164                 :            :  *     ^
     165                 :            :  * 1 1 1 3 3
     166                 :            :  *       ^
     167                 :            :  * 1 1 1 2 4
     168                 :            :  *       ^
     169                 :            :  * 1 1 1 1 5
     170                 :            :  * ^
     171                 :            :  * 0 2 2 2 3
     172                 :            :  *   ^
     173                 :            :  * This is recursive in nature. Restrictions on upper bounds of the vi or on
     174                 :            :  * the length of the partitions are straightforward to take into account. */
     175                 :            : 
     176                 :            : static void
     177                 :         98 : parse_interval(GEN a, long *amin, long *amax)
     178                 :            : {
     179      [ +  +  - ]:         98 :   switch (typ(a))
     180                 :            :   {
     181                 :            :   case t_INT:
     182                 :         42 :     *amax = itos(a);
     183                 :         42 :     break;
     184                 :            :   case t_VEC:
     185         [ -  + ]:         56 :     if (lg(a) != 3)
     186                 :          0 :       pari_err_TYPE("forpart [expect vector of type [amin,amax]]",a);
     187                 :         56 :     *amin = gtos(gel(a,1));
     188                 :         56 :     *amax = gtos(gel(a,2));
     189 [ +  - ][ +  - ]:         56 :     if (*amin>*amax || *amin<0 || *amax<=0)
                 [ -  + ]
     190                 :          0 :       pari_err_TYPE("forpart [expect 0<=min<=max, 0<max]",a);
     191                 :         56 :     break;
     192                 :            :   default:
     193                 :          0 :     pari_err_TYPE("forpart",a);
     194                 :            :   }
     195                 :         98 : }
     196                 :            : 
     197                 :            : void
     198                 :        126 : forpart_init(forpart_t *T, long k, GEN abound, GEN nbound)
     199                 :            : {
     200                 :            : 
     201                 :            :   /* bound on coefficients */
     202                 :        126 :   T->amin=1;
     203         [ +  + ]:        126 :   if (abound) parse_interval(abound,&T->amin,&T->amax);
     204                 :         70 :   else T->amax = k;
     205                 :            :   /* strip leading zeros ? */
     206                 :        126 :   T->strip = (T->amin > 0) ? 1 : 0;
     207                 :            :   /* bound on number of non-zero coefficients */
     208                 :        126 :   T->nmin=0;
     209         [ +  + ]:        126 :   if (nbound) parse_interval(nbound,&T->nmin,&T->nmax);
     210                 :         84 :   else T->nmax = k;
     211                 :            : 
     212                 :            :   /* non empty if nmin*amin <= k <= amax*nmax */
     213 [ +  - ][ -  + ]:        126 :   if ( T->amin*T->nmin > k || k > T->amax * T->nmax )
     214                 :            :   {
     215                 :          0 :     T->nmin = T->nmax = 0;
     216                 :            :   }
     217                 :            :   else
     218                 :            :   {
     219                 :            :     /* to reach nmin one must have k <= nmin*amax, otherwise increase nmin */
     220         [ +  + ]:        126 :     if ( T->nmin * T->amax < k )
     221                 :         98 :       T->nmin = 1 + (k - 1) / T->amax; /* ceil( k/tmax ) */
     222                 :            :     /* decrease nmax (if strip): k <= amin*nmax */
     223 [ +  + ][ +  + ]:        126 :     if (T->strip && T->nmax > k/T->amin)
     224                 :         14 :       T->nmax = k / T->amin; /* strip implies amin>0 */ /* fixme: take ceil( ) */
     225                 :            :     /* no need to change amin */
     226                 :            :     /* decrease amax if amax + (nmin-1)*amin > k  */
     227         [ +  + ]:        126 :     if ( T->amax + (T->nmin-1)* T->amin > k )
     228                 :         14 :       T->amax = k - (T->nmin-1)* T->amin;
     229                 :            :   }
     230                 :            : 
     231         [ +  + ]:        126 :   if ( T->amax < T->amin )
     232                 :         21 :     T->nmin = T->nmax = 0;
     233                 :            : 
     234                 :        126 :   T->v = zero_zv(T->nmax); /* partitions will be of length <= nmax */
     235                 :        126 :   T->k = k;
     236                 :        126 : }
     237                 :            : 
     238                 :            : GEN
     239                 :    3160759 : forpart_next(forpart_t *T)
     240                 :            : {
     241                 :    3160759 :   GEN v = T->v;
     242                 :    3160759 :   long n = lg(v)-1;
     243                 :            :   long i, s, a, k, vi, vn;
     244                 :            : 
     245 [ +  + ][ +  + ]:    3160759 :   if (n>0 && v[n])
     246                 :            :   {
     247                 :            :     /* find index to increase: i s.t. v[i+1],...v[n] is central a,..a,a+1,..a+1
     248                 :            :        keep s = v[i] + v[i+1] + ... + v[n] */
     249                 :    3160626 :     s = a = v[n];
     250 [ +  + ][ +  + ]:    4636786 :     for(i = n-1; i>0 && v[i]+1 >= a; s += v[i--]);
     251         [ +  + ]:    6321168 :     if (i == 0) {
     252                 :            :       /* v is central [ a, a, .. a, a+1, .. a+1 ] */
     253 [ +  + ][ +  + ]:        798 :       if ((n+1) * T->amin > s || n == T->nmax) return NULL;
     254                 :        714 :       i = 1; n++;
     255                 :        714 :       setlg(v, n+1);
     256                 :        714 :       vi = T->amin;
     257                 :            :     } else {
     258                 :    3159828 :       s += v[i];
     259                 :    3159828 :       vi = v[i]+1;
     260                 :            :     }
     261                 :            :   } else {
     262                 :            :     /* init v */
     263                 :        133 :     s = T->k;
     264         [ +  + ]:        133 :     if (T->amin == 0) T->amin = 1;
     265         [ +  + ]:        133 :     if (T->strip) { n = T->nmin; setlg(T->v, n+1); }
     266         [ +  + ]:        133 :     if (s==0)
     267                 :            :     {
     268 [ +  + ][ +  - ]:         21 :       if (n==0 && T->nmin==0) {T->nmin++; return v;}
     269                 :          7 :       return NULL;
     270                 :            :     }
     271         [ +  + ]:        112 :     if (n==0) return NULL;
     272                 :        105 :     vi = T->amin;
     273         [ +  + ]:        105 :     i = T->strip ? 1 : n + 1 - T->nmin; /* first non-zero index */
     274                 :            :   }
     275                 :            :   /* now fill [ v[i],... v[n] ] with s, start at vi */
     276                 :    3160647 :   vn = s - (n-i)*vi; /* expected value for v[n] */
     277 [ +  - ][ +  + ]:    3160647 :   if (T->amax && vn > T->amax)
     278                 :        140 :   {
     279                 :            :     /* do not exceed amax */
     280                 :            :     long ai, q, r;
     281                 :        140 :     vn -= vi;
     282                 :        140 :     ai = T->amax - vi;
     283                 :        140 :     q = vn / ai; /* number of nmax */
     284                 :        140 :     r = vn % ai; /* value before nmax */
     285                 :            :     /* fill [ v[i],... v[n] ] as [ vi,... vi, vi+r, amax,... amax ] */
     286         [ +  + ]:        469 :     while ( q-- ) v[n--] = T->amax;
     287         [ +  + ]:        140 :     if ( n >= i ) v[n--] = vi + r;
     288         [ +  + ]:        371 :     while ( n >= i ) v[n--] = vi;
     289                 :            :   } else {
     290                 :            :     /* fill as [ v[i], ... v[i], vn ] */
     291         [ +  + ]:    7796019 :     for ( k=i; k<n; v[k++] = vi );
     292                 :    3160507 :     v[n] = vn;
     293                 :            :   }
     294                 :    3160759 :   return v;
     295                 :            : }
     296                 :            : 
     297                 :            : GEN
     298                 :          0 : forpart_prev(forpart_t *T)
     299                 :            : {
     300                 :          0 :   GEN v = T->v;
     301                 :          0 :   long n = lg(v)-1;
     302                 :            :   long j, ni, q, r;
     303                 :            :   long i, s;
     304 [ #  # ][ #  # ]:          0 :   if (n>0 && v[n])
     305                 :            :   {
     306                 :            :     /* find index to decrease: start of last constant sequence, excluding v[n] */
     307                 :          0 :     i = n-1; s = v[n];
     308 [ #  # ][ #  # ]:          0 :     while (i>1 && (v[i-1]==v[i] || v[i+1]==T->amax))
                 [ #  # ]
     309                 :          0 :       s+= v[i--];
     310         [ #  # ]:          0 :     if (!i) return NULL;
     311                 :            :     /* amax condition: cannot decrease i if maximal on the right */
     312         [ #  # ]:          0 :     if ( v[i+1] == T->amax ) return NULL;
     313                 :            :     /* amin condition: stop if below except if strip & try to remove */
     314         [ #  # ]:          0 :     if (v[i] == T->amin) {
     315         [ #  # ]:          0 :       if (!T->strip) return NULL;
     316                 :          0 :       s += v[i]; v[i] = 0;
     317                 :            :     } else {
     318                 :          0 :       v[i]--; s++;
     319                 :            :     }
     320                 :            :     /* zero case... */
     321         [ #  # ]:          0 :     if (v[i] == 0)
     322                 :            :     {
     323         [ #  # ]:          0 :       if (T->nmin > n-i) return NULL; /* need too many non zero coeffs */
     324                 :            :       /* reduce size of v ? */
     325         [ #  # ]:          0 :       if (T->strip) {
     326                 :          0 :         i = 0; n--;
     327                 :          0 :         setlg(v, n+1);
     328                 :            :       }
     329                 :            :     }
     330                 :            :   } else
     331                 :            :   {
     332                 :          0 :     s = T->k;
     333                 :          0 :     i = 0;
     334         [ #  # ]:          0 :     if (s==0)
     335                 :            :     {
     336 [ #  # ][ #  # ]:          0 :       if (n==0 && T->nmin==0) {T->nmin++; return v;}
     337                 :          0 :       return NULL;
     338                 :            :     }
     339 [ #  # ][ #  # ]:          0 :     if (n*T->amax < s || s < T->nmin*T->amin) return NULL;
     340                 :            :   }
     341                 :            :   /* set minimal partition of sum s starting from index i+1 */
     342                 :          0 :   ni = n-i;
     343                 :          0 :   q = s / ni;
     344                 :          0 :   r = s % ni;
     345         [ #  # ]:          0 :   for(j=i+1;   j<=n-r; j++) v[j]=q;
     346         [ #  # ]:          0 :   for(j=n-r+1; j<=n;   j++) v[j]=q + 1;
     347                 :          0 :   return v;
     348                 :            : }
     349                 :            : 
     350                 :            : static long
     351                 :         42 : countpart(long k, GEN abound, GEN nbound)
     352                 :            : {
     353                 :         42 :   pari_sp av = avma;
     354                 :            :   long n;
     355                 :            :   forpart_t T;
     356         [ +  + ]:         42 :   if (k<0) return 0;
     357                 :         35 :   forpart_init(&T, k, abound, nbound);
     358         [ +  + ]:        343 :   for (n=0; forpart_next(&T); n++)
     359                 :        308 :   avma = av;
     360                 :         42 :   return n;
     361                 :            : }
     362                 :            : 
     363                 :            : GEN
     364                 :         42 : partitions(long k, GEN abound, GEN nbound)
     365                 :            : {
     366                 :            :   GEN v;
     367                 :            :   forpart_t T;
     368                 :         42 :   long i, n = countpart(k,abound,nbound);
     369         [ +  + ]:         42 :   if (n==0) return cgetg(1, t_VEC);
     370                 :         28 :   forpart_init(&T, k, abound, nbound);
     371                 :         28 :   v = cgetg(n+1, t_VEC);
     372         [ +  + ]:        336 :   for (i=1; i<=n; i++)
     373                 :        308 :     gel(v,i)=zv_copy(forpart_next(&T));
     374                 :         42 :   return v;
     375                 :            : }
     376                 :            : 
     377                 :            : void
     378                 :         63 : forpart(void *E, long call(void*, GEN), long k, GEN abound, GEN nbound)
     379                 :            : {
     380                 :            :   GEN v;
     381                 :            :   forpart_t T;
     382                 :         63 :   forpart_init(&T, k, abound, nbound);
     383         [ +  + ]:    3160108 :   while ((v=forpart_next(&T)))
     384         [ -  + ]:    3160045 :     if (call(E, v)) break;
     385                 :         63 : }
     386                 :            : 
     387                 :            : void
     388                 :         70 : forpart0(GEN k, GEN code, GEN abound, GEN nbound)
     389                 :            : {
     390                 :         70 :   pari_sp av = avma;
     391         [ -  + ]:         70 :   if (typ(k) != t_INT) pari_err_TYPE("forpart",k);
     392         [ +  + ]:        133 :   if (signe(k)<0) return;
     393                 :         63 :   push_lex(gen_0, code);
     394                 :         63 :   forpart((void*)code, &gp_evalvoid, itos(k), abound, nbound);
     395                 :         63 :   pop_lex(1);
     396                 :         63 :   avma=av;
     397                 :            : }

Generated by: LCOV version 1.9