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 - basemath - zetamult.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21072-998352a) Lines: 246 250 98.4 %
Date: 2017-09-26 06:25:23 Functions: 21 22 95.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : This file is part of the PARI/GP package.
       3             : 
       4             : PARI/GP is free software; you can redistribute it and/or modify it under the
       5             : terms of the GNU General Public License as published by the Free Software
       6             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       7             : ANY WARRANTY WHATSOEVER.
       8             : 
       9             : Check the License for details. You should have received a copy of it, along
      10             : with the package; see the file 'COPYING'. If not, write to the Free Software
      11             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      12             : 
      13             : /********************************************************************/
      14             : /**                                                                **/
      15             : /**             MULTIPLE ZETA VALUES (AKHILESH ALGORITHM)          **/
      16             : /**                                                                **/
      17             : /********************************************************************/
      18             : #include "pari.h"
      19             : #include "paripriv.h"
      20             : 
      21             : static long
      22        1841 : la(long e, long f) { return (e == f)? 2: (e? 1: 3); }
      23             : static GEN
      24      263669 : lamul(long la, GEN s)
      25             : {
      26      263669 :   switch(la)
      27             :   {
      28      123886 :     case 2: return gmul2n(s,1);
      29       78512 :     case 3: return gmulgs(s,3);
      30       61271 :     default: return s;
      31             :   }
      32             : }
      33             : 
      34             : /* dual of evec[1..l-1] */
      35             : static GEN
      36        1841 : revslice(GEN evec, long l)
      37             : {
      38        1841 :   GEN res = cgetg(l, t_VECSMALL);
      39             :   long i;
      40        1841 :   for (i = 1; i < l; ++i) res[i] = 1 - evec[l-i];
      41        1841 :   return res;
      42             : }
      43             : 
      44             : /* N.B. evec[ne] = 1 */
      45             : static GEN
      46       10871 : etoa(GEN evec)
      47             : {
      48       10871 :   long ct = 0, ctold = 0, i = 1, le = lg(evec);
      49       10871 :   GEN avec = cgetg(le, t_VECSMALL);
      50      138131 :   while (++ct < le)
      51      116389 :     if (evec[ct] == 1) { avec[i++] = ct - ctold; ctold = ct; }
      52       10871 :   setlg(avec, i); return avec;
      53             : }
      54             : 
      55             : static GEN
      56        7315 : atoe(GEN avec)
      57             : {
      58        7315 :   long i, l = lg(avec);
      59        7315 :   GEN evec = cgetg(l, t_VEC);
      60        7315 :   for (i = 1; i < l; i++) { long a = avec[i]; gel(evec,i) = vecsmall_ei(a,a); }
      61        7315 :   return shallowconcat1(evec);
      62             : }
      63             : 
      64             : /* vphi[i] contains phip(j,avec[i..r]) for 1<= j < L
      65             :  * vpow[a][j] = j^-a as a t_INT or t_REAL; j < L */
      66             : static GEN
      67        1099 : phip(GEN avec, GEN vpow)
      68             : {
      69        1099 :   long i, r = lg(avec) - 1;
      70        1099 :   GEN vphi = cgetg(r+1, t_VEC);
      71        1099 :   gel(vphi, r) = gel(vpow, avec[r]);
      72        7077 :   for (i = r-1; i >= 1; i--)
      73             :   {
      74        5978 :     GEN t, u, phi = gel(vphi,i+1), pow = gel(vpow, avec[i]);
      75        5978 :     long j, L = lg(pow);
      76        5978 :     gel(vphi, i) = u = cgetg(L, t_VEC);
      77        5978 :     gel(u,1) = gen_0;
      78        5978 :     gel(u,2) = (i==r-1)? gel(pow,2): gen_0;
      79        5978 :     t = gel(phi,1); /* 0 or 1 */
      80      893704 :     for (j = 3; j < L; j++)
      81             :     {
      82      887726 :       t = mpadd(t, gel(phi,j-1));
      83      887726 :       gel(u,j) = mpmul(t, gel(pow,j)); /* t / j^avec[i] */
      84             :     }
      85             :   }
      86        1099 :   return vphi;
      87             : }
      88             : 
      89             : /* Return 1 if vec2 RHS of vec1, -1 if vec1 RHS of vec2, 0 else */
      90             : static long
      91       14196 : isrhs(GEN v1, GEN v2)
      92             : {
      93       14196 :   long s = 1, i, l1 = lg(v1), l2 = lg(v2);
      94       14196 :   if (l1 < l2) { s = -1; swap(v1,v2); lswap(l1,l2); }
      95       83685 :   for (i = l2-1; i >= 1; --i)
      96       81102 :     if (v2[i] != v1[l1-l2+i]) return 0;
      97        2583 :   return s;
      98             : }
      99             : 
     100             : static long
     101       15295 : istruerhs(GEN v1, GEN v2)
     102             : {
     103       15295 :   long i, l1 = lg(v1), l2 = lg(v2);
     104       15295 :   if (l1 < l2) return 0;
     105       75656 :   for (i = l2-1; i >= 1; --i)
     106       71974 :     if (v2[i] != v1[l1-l2+i]) return 0;
     107        3682 :   return l1-l2+1;
     108             : }
     109             : 
     110             : /* a is a rhs of a unique v[m] */
     111             : static GEN
     112        3682 : isinphi(GEN v, GEN a, GEN vphi)
     113             : {
     114        3682 :   long m, l = lg(v);
     115       15295 :   for (m = 1; m < l; m++)
     116             :   {
     117       15295 :     long s = istruerhs(gel(v,m), a);
     118       15295 :     if (s) return gmael(vphi,m,s);
     119             :   }
     120             :   return NULL; /* LCOV_EXCL_LINE */
     121             : }
     122             : 
     123             : /* If v RHS of LR[i] for some i, return LR. If LR[i] RHS (strict) of v, replace
     124             :  * LR[i] by v. If none, add v to LR. */
     125             : static GEN
     126        3682 : addevec(GEN LR, GEN v)
     127             : {
     128        3682 :   long s, i, l1 = lg(LR);
     129       15295 :   for (i = 1; i < l1; i++)
     130             :   {
     131       14196 :     s = isrhs(gel(LR,i), v);
     132       14196 :     if (s == 1) return LR;
     133       11851 :     if (s ==-1) { gel(LR,i) = v; return LR; }
     134             :   }
     135        1099 :   return vec_append(LR,v);
     136             : }
     137             : 
     138             : /* N > 2 */
     139             : static GEN
     140         140 : get_vbin(long N, long prec)
     141             : {
     142         140 :   GEN v = cgetg(N+1, t_VEC);
     143             :   long n;
     144         140 :   gel(v,1) = gen_0; /* unused */
     145         140 :   gel(v,2) = invr(utor(6,prec));
     146         140 :   for (n = 3; n <= N; n++) gel(v,n) = divru(mulru(gel(v,n-1), n), 4*n-2);
     147         140 :   return v;
     148             : }
     149             : /* m < k */
     150             : static GEN
     151         140 : zetamultinit_i(long k, long m, long bitprec)
     152             : {
     153             :   long i, N, prec;
     154         140 :   GEN vpow = cgetg(m+1, t_VEC);
     155             : 
     156         140 :   bitprec += 64*(1+(k>>5));
     157         140 :   prec = nbits2prec(bitprec);
     158         140 :   N = 5 + bitprec/2;
     159         140 :   gel(vpow,1) = vecpowug(N, gen_m1, prec);
     160         574 :   for (i = 2; i <= m; i++)
     161             :   {
     162         434 :     GEN pow = cgetg(N+1, t_VEC), powm = gel(vpow,i-1);
     163             :     long j;
     164         434 :     gel(pow,1) = gen_1;
     165         434 :     gel(pow,2) = real2n(-i, prec);
     166         434 :     for (j = 3; j <= N; j++) gel(pow,j) = divru(gel(powm,j), j);
     167         434 :     gel(vpow,i) = pow;
     168             :   }
     169         140 :   return mkvec2(vpow, get_vbin(N, prec));
     170             : }
     171             : GEN
     172           7 : zetamultinit(long k, long prec)
     173             : {
     174           7 :   pari_sp av = avma;
     175           7 :   if (k <= 0) pari_err_DOMAIN("zetamultinit", "weight", "<=", gen_0, stoi(k));
     176           7 :   return gerepilecopy(av, zetamultinit_i(k, k-1, prec2nbits(prec)));
     177             : }
     178             : GEN
     179         189 : zetamult0(GEN avec, GEN T, long prec)
     180             : {
     181         189 :   pari_sp ltop = avma;
     182             :   long k, n, i, j, l, lbin;
     183         189 :   GEN vpow, vphi, vbin, S, s, LR, MA, MR, evec = gen_0;
     184             : 
     185         189 :   avec = zetamultconvert(avec, 1);
     186         140 :   if (lg(avec) == 1) return gen_1;
     187         133 :   evec = atoe(avec);
     188         133 :   k = lg(evec)-1; /* weight */
     189         133 :   LR = cgetg(1, t_VEC);
     190         133 :   MA = cgetg(k, t_VEC);
     191         133 :   MR = cgetg(k, t_VEC);
     192        1974 :   for (i = 1; i < k; ++i)
     193             :   {
     194        1841 :     gel(MA,i) = etoa(revslice(evec, i+1));
     195        1841 :     gel(MR,i) = etoa(vecslice(evec, i+1, k));
     196        1841 :     LR = addevec(addevec(LR, gel(MA,i)), gel(MR,i));
     197             :   }
     198         133 :   if (!T)
     199             :   {
     200         133 :     long m = vecvecsmall_max(LR); /* < k */
     201         133 :     T = zetamultinit_i(k, m, prec2nbits(prec));
     202             :   }
     203             :   else
     204             :   {
     205             :     long M;
     206           0 :     if (typ(T) != t_VEC || lg(T) != 3) pari_err_TYPE("zetamult", T);
     207           0 :     M = lg(gel(T,1)); /* need M > m, which is < k */
     208           0 :     if (M < k) pari_err_DOMAIN("zetamult", "weight", ">", utoi(M), utoi(k));
     209             :   }
     210         133 :   vpow = gel(T,1);
     211         133 :   vbin = gel(T,2);
     212         133 :   l = lg(LR); vphi = cgetg(l, t_VEC);
     213         133 :   for (j = 1; j < l; j++) gel(vphi,j) = phip(gel(LR,j), vpow);
     214         133 :   lbin = lg(vbin);
     215         133 :   S = cgetg(lbin, t_VEC);
     216        1974 :   for (i = 1; i < k; i++)
     217             :   {
     218        1841 :     long LA = la(evec[i],evec[i+1]);
     219        1841 :     GEN phi1 = isinphi(LR, gel(MA,i), vphi);
     220        1841 :     GEN phi2 = isinphi(LR, gel(MR,i), vphi);
     221        1841 :     if (i == 1)
     222       17374 :       for (n = 1; n < lbin; n++)
     223       17241 :         gel(S,n) = lamul(LA, mpmul(gel(phi1,n), gel(phi2,n)));
     224             :     else
     225      248136 :       for (n = 1; n < lbin; n++)
     226      246428 :         gel(S,n) = mpadd(gel(S,n), lamul(LA, mpmul(gel(phi1,n), gel(phi2,n))));
     227             :   }
     228         133 :   s = gmul2n(gel(S,1), -1);
     229         133 :   for (n = 2; n < lbin; n++) s = gadd(s, mpmul(gel(S,n), gel(vbin,n)));
     230         133 :   return gerepileuptoleaf(ltop, rtor(s,prec));
     231             : }
     232             : GEN
     233           0 : zetamult(GEN avec, long prec) { return zetamult0(avec, NULL, prec); }
     234             : 
     235             : /**************************************************************/
     236             : /*                         ALL MZV's                          */
     237             : /**************************************************************/
     238             : 
     239             : /* vecsmall to binary */
     240             : static long
     241        7203 : myfd(GEN evec, long ini, long fin)
     242             : {
     243        7203 :   long i, s = 0;
     244        7203 :   for (i = ini; i <= fin; ++i) s = evec[i] | (s << 1);
     245        7203 :   return s;
     246             : }
     247             : 
     248             : /* Given admissible evec w = 0e_2....e_{k-1}1, compute a,b,v such that
     249             :  * w=0{1}_{b-1}v{0}_{a-1}1 with v empty or admissible.
     250             :  * Input: binary vector evec */
     251             : static void
     252          63 : findabv(GEN w, long *pa, long *pb, long *pminit, long *pmmid, long *pmfin)
     253             : {
     254          63 :   long le = lg(w) - 2;
     255          63 :   if (le == 0)
     256             :   {
     257           7 :     *pa = 1;
     258           7 :     *pb = 1;
     259           7 :     *pminit = 2;
     260           7 :     *pmfin = 2;
     261           7 :     *pmmid = 1;
     262             :   }
     263             :   else
     264             :   {
     265             :     long a, b, j, lv;
     266          70 :     for (j = 1; j <= le; ++j)
     267          70 :       if (!w[j+1]) break;
     268          56 :     b = j;
     269         126 :     for (j = le; j >= 1; --j)
     270         105 :       if (w[j+1]) break;
     271          56 :     a = le + 1 - j;
     272          56 :     lv = le + 2 - a - b;
     273          56 :     if (lv > 0)
     274             :     {
     275          21 :       long v = myfd(w, b + 1, le - a + 2);
     276          21 :       *pa = a;
     277          21 :       *pb = b;
     278          21 :       *pminit = (((1 << b) - 1) << (lv - 1)) + (v/2) + 2;
     279          21 :       *pmfin = (((1 << (lv - 1)) + v) << (a - 1)) + 2;
     280          21 :       *pmmid = (1 << (lv - 2)) + (v/2) + 2;
     281             :     }
     282             :     else
     283             :     {
     284          35 :       *pa = a;
     285          35 :       *pb = b;
     286          35 :       *pminit = (1 << (b - 1)) + 1;
     287          35 :       *pmfin = (a == 1) ? 2 : (1 << (a - 2)) + 2;
     288          35 :       *pmmid = 1;
     289             :     }
     290             :   }
     291          63 : }
     292             : 
     293             : /* Returns 'all':
     294             : * all[1] contains zeta(emptyset)_{n-1,n-1},
     295             : * all[2] contains zeta({0})_{n-1,n-1}=zeta({1})_{n-1,n-1} for n >= 2,
     296             : * all[m+2][n] : 1 <= m < 2^{k-2}, 1 <= n <= N + 1
     297             : * contains zeta(w)_{n-1,n-1}, w corresponding to m,n
     298             : * all[m+2] : 2^{k-2} <= m < 2^{k-1} contains zeta(w), w corresponding to m
     299             : (code: w=0y1 iff m=1y). */
     300             : static GEN
     301           7 : fillall(long k, long bitprec)
     302             : {
     303           7 :   long N = 1 + bitprec/2, prec = nbits2prec(bitprec);
     304           7 :   long k1, j, n, m, mbar = 0, K = 1 << (k - 1), K2 = K/2;
     305             :   GEN all, v, p1, p2, r1, pab, S;
     306             : 
     307           7 :   r1 = real_1(prec);
     308           7 :   pab = cgetg(N+1, t_VEC); gel(pab, 1) = gen_0; /* not needed */
     309           7 :   for (n = 2; n <= N; n++) gel(pab, n) = powersr(divru(r1, n), k);
     310             :   /* 1/n^a = gmael(pab, n, a + 1) */
     311           7 :   all = cgetg(K + 2, t_VEC);
     312           7 :   gel(all,1) = v = cgetg(N+1, t_VEC);
     313           7 :   gel(v,1) = gen_0; /* unused */
     314           7 :   gel(v,2) = real2n(-1,prec);
     315           7 :   gel(v,3) = invr(utor(6,prec)); /* cf get_vbin: shifted by 1 :-( */
     316           7 :   for (n = 3; n < N; n++) gel(v,n+1) = divru(mulru(gel(v,n), n), 4*n-2);
     317             : 
     318           7 :   gel(all,2) = p1 = cgetg(N+1, t_VEC);
     319           7 :   gel(p1,1) = gen_0; /* unused */
     320           7 :   for (j = 2; j <= N; j++) gel(p1,j) = divru(gel(v,j), j-1);
     321             : 
     322          56 :   for (m = 1; m < K2; m++)
     323             :   {
     324          49 :     gel(all, m+2) = p1 = cgetg(N+1, t_VEC);
     325          49 :     for (n = 1; n < N; n++) gel(p1, n) = cgetr(prec);
     326          49 :     gel(p1, n) = gen_0;
     327             :   }
     328           7 :   for (m = K2; m < K; m++) gel(all, m+2) = utor(0, prec);
     329          35 :   for (k1 = 2; k1 <= k; k1++)
     330             :   { /* Assume length evec < k1 filled */
     331             :     /* If evec = 0e_2...e_{k_1-1}1 then m = (1e_2...e_{k_1-1})_2 */
     332          28 :     GEN w = cgetg(k1, t_VECSMALL);
     333          28 :     long M = 1 << (k1 - 2);
     334          28 :     pari_sp av = avma;
     335         133 :     for (m = M; m < 2*M; m++)
     336             :     {
     337             :       GEN pinit, pfin, pmid;
     338         105 :       long comp, a, b, minit, mfin, mmid, mc = m, ii = 0;
     339         105 :       p1 = gel(all, m + 2);
     340         343 :       for (j = k1 - 1; j >= 2; --j)
     341             :       {
     342         238 :         w[j] = mc & 1;
     343         238 :         ii = (1 - w[j]) | (ii<<1);
     344         238 :         mc >>= 1;
     345             :       }
     346         105 :       mbar = M + ii;
     347         105 :       comp = mbar - m;
     348         105 :       if (comp < 0) continue;
     349          63 :       p2 = gel(all, mbar + 2);
     350          63 :       findabv(w, &a,&b,&minit,&mmid,&mfin);
     351          63 :       pinit= gel(all, minit);
     352          63 :       pfin = gel(all, mfin);
     353          63 :       pmid = gel(all, mmid);
     354        7056 :       for (n = N-1; n > 1; n--, avma = av)
     355             :       {
     356        6993 :         GEN t = mpmul(gel(pinit,n+1), gmael(pab, n, a+1));
     357        6993 :         GEN u = mpmul(gel(pfin, n+1), gmael(pab, n, b+1));
     358        6993 :         GEN v = mpmul(gel(pmid, n+1), gmael(pab, n, a+b+1));
     359        6993 :         S = mpadd(k1 < k ? gel(p1, n+1) : p1, mpadd(mpadd(t, u), v));
     360        6993 :         if (!signe(S)) S = gen_0;
     361        6993 :         mpaff(S, k1 < k ? gel(p1, n) : p1);
     362        6993 :         if (comp > 0 && k1 < k) mpaff(S, gel(p2, n));
     363             :       }
     364             :       { /* n = 1: same formula simplifies */
     365          63 :         GEN t = gel(pinit,2), u = gel(pfin,2), v = gel(pmid,2);
     366          63 :         S = mpadd(k1 < k ? gel(p1,2) : p1, mpadd(mpadd(t, u), v));
     367          63 :         if (!signe(S)) S = gen_0;
     368          63 :         mpaff(S, k1 < k ? gel(p1,1) : p1);
     369          63 :         if (comp > 0 && k1 < k) mpaff(S, gel(p2, 1));
     370          63 :         avma = av;
     371             :       }
     372          63 :       if (comp > 0 && k1 == k) mpaff(p1, p2);
     373             :     }
     374             :   }
     375         112 :   for (j = 1; j < K; j++)
     376         105 :     gel(all, j) = j < K2 ? gmael(all, j+2, 1) : gel(all, j+2);
     377           7 :   setlg(all, K); return all;
     378             : }
     379             : 
     380             : GEN
     381          21 : zetamultall(long k, long prec)
     382             : {
     383          21 :   pari_sp av = avma;
     384          21 :   if (k < 1) pari_err_DOMAIN("zetamultall", "k", "<", gen_1, stoi(k));
     385          14 :   if (k >= 64) pari_err_OVERFLOW("zetamultall");
     386           7 :   return gerepilecopy(av, fillall(k, prec2nbits(prec) + 32));
     387             : }
     388             : 
     389             : /* m > 0 */
     390             : static GEN
     391        7189 : mtoevec(GEN m)
     392             : {
     393        7189 :   GEN e = vecsmall_append(binary_zv(m), 1);
     394        7189 :   e[1] = 0; return e;
     395             : }
     396             : 
     397             : static GEN
     398        7182 : etoindex(GEN evec)
     399             : {
     400        7182 :   long k = lg(evec) - 1;
     401        7182 :   return utoipos((1 << (k-2)) + myfd(evec, 2, k-1));
     402             : }
     403             : 
     404             : /* Conversions: types are evec, avec, m (if evec=0y1, m=(1y)_2).
     405             :    fl is respectively 0, 1, 2. Type of a is autodetected. */
     406             : GEN
     407       14588 : zetamultconvert(GEN a, long fl)
     408             : {
     409       14588 :   pari_sp av = avma;
     410             :   long i, l;
     411       14588 :   if (fl < 0 || fl > 2) pari_err_FLAG("zetamultconvert");
     412       14588 :   switch(typ(a))
     413             :   {
     414             :     case t_INT:
     415        7196 :       if (signe(a) <= 0) pari_err_TYPE("zetamultconvert",a);
     416        7196 :       switch (fl)
     417             :       {
     418           7 :         case 0: a = mtoevec(a); break;
     419        7182 :         case 1: a = etoa(mtoevec(a)); break;
     420           7 :         case 2: a = icopy(a); break;
     421             :       }
     422        7196 :       break;
     423             :     case t_VEC: case t_COL: case t_VECSMALL:
     424        7385 :       a = gtovecsmall(a);
     425        7385 :       l = lg(a);
     426        7385 :       if (a[1] == 0)
     427             :       {
     428          42 :         if (!a[l-1]) pari_err_TYPE("zetamultconvert", a);
     429         217 :         for (i = 1; i < l; i++)
     430         196 :           if (a[i] & ~1UL) pari_err_TYPE("zetamultconvert", a);
     431          21 :         switch (fl)
     432             :         {
     433           7 :           case 1: a = etoa(a); break;
     434           7 :           case 2: a = etoindex(a);
     435             :         }
     436             :       }
     437             :       else
     438             :       {
     439        7343 :         if (a[1] < 2) pari_err_TYPE("zetamultconvert", a);
     440       36883 :         for (i = 2; i < l; i++)
     441       29561 :           if (a[i] <= 0) pari_err_TYPE("zetamultconvert", a);
     442        7322 :         switch (fl)
     443             :         {
     444           7 :           case 0: a = atoe(a); break;
     445        7175 :           case 2: a = etoindex(atoe(a));
     446             :         }
     447             :       }
     448        7343 :       break;
     449           7 :     default: pari_err_TYPE("zetamultconvert", a);
     450             :   }
     451       14539 :   return gerepileuptoleaf(av, a);
     452             : }

Generated by: LCOV version 1.11