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 - basemath - dirichlet.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25819-e703fe1174) Lines: 263 263 100.0 %
Date: 2020-09-18 06:10:04 Functions: 17 17 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2015  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             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**           Dirichlet series through Euler product               **/
      17             : /**                                                                **/
      18             : /********************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : static void
      23          28 : err_direuler(GEN x)
      24          28 : { pari_err_DOMAIN("direuler","constant term","!=", gen_1,x); }
      25             : 
      26             : /* s = t_POL (tolerate t_SER of valuation 0) of constant term = 1
      27             :  * d = minimal such that p^d > X
      28             :  * V indexed by 1..X will contain the a_n
      29             :  * v[1..n] contains the indices nj such that V[nj] != 0 */
      30             : static long
      31       20468 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
      32             : {
      33       20468 :   long i, j, m = n, D = minss(d+2, lg(s));
      34       20468 :   ulong q = 1, X = lg(V)-1;
      35             : 
      36       70245 :   for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
      37             :   {
      38       49777 :     GEN aq = gel(s,i);
      39       49777 :     if (gequal0(aq)) continue;
      40             :     /* j = 1 */
      41       42084 :     gel(V,q) = aq;
      42       42084 :     v[++n] = q;
      43     1766674 :     for (j = 2; j <= m; j++)
      44             :     {
      45     1724590 :       ulong nj = umuluu_le(uel(v,j), q, X);
      46     1724590 :       if (!nj) continue;
      47      133910 :       gel(V,nj) = gmul(aq, gel(V,v[j]));
      48      133910 :       v[++n] = nj;
      49             :     }
      50             :   }
      51       20468 :   return n;
      52             : }
      53             : 
      54             : /* ap != 0 for efficiency, p > sqrt(X) */
      55             : static void
      56      150010 : dirmuleuler_large(GEN V, ulong p, GEN ap)
      57             : {
      58      150010 :   long j, jp, X = lg(V)-1;
      59      150010 :   gel(V,p) = ap;
      60      599879 :   for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
      61      150010 : }
      62             : 
      63             : static ulong
      64        9107 : direulertou(GEN a, GEN fl(GEN))
      65             : {
      66        9107 :   if (typ(a) != t_INT)
      67             :   {
      68          49 :     a = fl(a);
      69          28 :     if (typ(a) != t_INT) pari_err_TYPE("direuler", a);
      70             :   }
      71        9086 :   return signe(a)<=0 ? 0: itou(a);
      72             : }
      73             : 
      74             : static GEN
      75        3556 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
      76             : {
      77        3556 :   long i, l = lg(Sbad);
      78        3556 :   ulong X = lg(V)-1;
      79        3556 :   GEN pbad = gen_1;
      80        9275 :   for (i = 1; i < l; i++)
      81             :   {
      82        5754 :     GEN ai = gel(Sbad,i);
      83             :     ulong q;
      84        5754 :     if (typ(ai) != t_VEC || lg(ai) != 3)
      85          14 :       pari_err_TYPE("direuler [bad primes]",ai);
      86        5740 :     q = gtou(gel(ai,1));
      87        5733 :     if (q <= X)
      88             :     {
      89        4627 :       long d = ulogint(X, q) + 1;
      90        4627 :       GEN s = direuler_factor(gel(ai,2), d);
      91        4613 :       *n = dirmuleuler_small(V, v, *n, q, s, d);
      92        4613 :       pbad = muliu(pbad, q);
      93             :     }
      94             :   }
      95        3521 :   return pbad;
      96             : }
      97             : 
      98             : GEN
      99         504 : direuler_bad(void *E, GEN (*eval)(void *,GEN,long), GEN a,GEN b,GEN c, GEN Sbad)
     100             : {
     101             :   ulong au, bu, X, sqrtX, n, p;
     102         504 :   pari_sp av0 = avma;
     103             :   GEN gp, v, V;
     104             :   forprime_t T;
     105         504 :   au = direulertou(a, gceil);
     106         497 :   bu = direulertou(b, gfloor);
     107         490 :   X = c ? direulertou(c, gfloor): bu;
     108         483 :   if (X == 0) return cgetg(1,t_VEC);
     109         476 :   if (bu > X) bu = X;
     110         476 :   if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
     111         462 :   v = vecsmall_ei(X, 1);
     112         462 :   V = vec_ei(X, 1);
     113         462 :   n = 1;
     114         462 :   if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
     115         427 :   p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
     116        2758 :   while (p <= sqrtX && (p = u_forprime_next(&T)))
     117        2352 :     if (!Sbad || umodiu(Sbad, p))
     118             :     {
     119        2247 :       long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
     120             :       GEN s;
     121        2247 :       gp[2] = p; s = eval(E, gp, d);
     122        2226 :       n = dirmuleuler_small(V, v, n, p, s, d);
     123             :     }
     124       55391 :   while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
     125       54985 :     if (!Sbad || umodiu(Sbad, p))
     126             :     {
     127             :       GEN s;
     128       54978 :       gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
     129       54978 :       if (lg(s) > 3 && !gequal0(gel(s,3)))
     130       23072 :         dirmuleuler_large(V, p, gel(s,3));
     131             :     }
     132         406 :   return gerepilecopy(av0,V);
     133             : }
     134             : 
     135             : /* return a t_SER or a truncated t_POL to precision n */
     136             : GEN
     137       61852 : direuler_factor(GEN s, long n)
     138             : {
     139       61852 :   long t = typ(s);
     140       61852 :   if (is_scalar_t(t))
     141             :   {
     142       31612 :     if (!gequal1(s)) err_direuler(s);
     143       31598 :     return scalarpol_shallow(s,0);
     144             :   }
     145       30240 :   switch(t)
     146             :   {
     147        5467 :     case t_POL: break; /* no need to RgXn_red */
     148       24458 :     case t_RFRAC:
     149             :     {
     150       24458 :       GEN p = gel(s,1), q = gel(s,2);
     151       24458 :       q = RgXn_red_shallow(q,n);
     152       24458 :       s = RgXn_inv(q, n);
     153       24458 :       if (typ(p) == t_POL && varn(p) == varn(q))
     154             :       {
     155          28 :         p = RgXn_red_shallow(p, n);
     156          28 :         s = RgXn_mul(s, p, n);
     157             :       }
     158             :       else
     159       24430 :         if (!gequal1(p)) s = RgX_Rg_mul(s, p);
     160       24458 :       if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
     161       24444 :       break;
     162             :     }
     163         308 :     case t_SER:
     164         308 :       if (!signe(s) || valp(s) || !gequal1(gel(s,2))) err_direuler(s);
     165         308 :       break;
     166           7 :     default: pari_err_TYPE("direuler", s);
     167             :   }
     168       30219 :   return s;
     169             : }
     170             : 
     171             : struct eval_bad
     172             : {
     173             :   void *E;
     174             :   GEN (*eval)(void *, GEN);
     175             : };
     176             : static GEN
     177         420 : eval_bad(void *E, GEN p, long n)
     178             : {
     179         420 :   struct eval_bad *d = (struct eval_bad*) E;
     180         420 :   return direuler_factor(d->eval(d->E, p), n);
     181             : }
     182             : GEN
     183         133 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
     184             : {
     185             :   struct eval_bad d;
     186         133 :   d.E= E; d.eval = eval;
     187         133 :   return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
     188             : }
     189             : 
     190             : static GEN
     191       26614 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
     192             : {
     193       26614 :   GEN P = cgetg(n+1, t_VECSMALL);
     194             :   long i, j;
     195      242816 :   for (i = 1, j = 1; i <= n; i++)
     196             :   {
     197      219863 :     ulong p = u_forprime_next(T);
     198      219863 :     if (!p) { *running = 0; break; }
     199      216202 :     if (Sbad && umodiu(Sbad, p)==0) continue;
     200      211701 :     uel(P,j++) = p;
     201             :   }
     202       26614 :   setlg(P, j);
     203       26614 :   return P;
     204             : }
     205             : 
     206             : GEN
     207        3668 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
     208             : {
     209             :   ulong au, bu, X, sqrtX, n, snX, nX;
     210        3668 :   pari_sp av0 = avma;
     211             :   GEN v, V;
     212             :   forprime_t T;
     213             :   struct pari_mt pt;
     214        3668 :   long running = 1, pending = 0;
     215        3668 :   au = direulertou(a, gceil);
     216        3668 :   bu = direulertou(b, gfloor);
     217        3668 :   X = c ? direulertou(c, gfloor): bu;
     218        3668 :   if (X == 0) return cgetg(1,t_VEC);
     219        3668 :   if (bu > X) bu = X;
     220        3668 :   if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
     221        3661 :   v = vecsmall_ei(X, 1);
     222        3661 :   V = vec_ei(X, 1);
     223        3661 :   n = 1;
     224        3661 :   if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
     225        3661 :   sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
     226        3661 :   if (snX)
     227             :   {
     228        3647 :     GEN P = primelist(&T, Sbad, snX, &running);
     229        3647 :     GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
     230        3647 :     long i, l = lg(P);
     231       17276 :     for (i = 1; i < l; i++)
     232             :     {
     233       13629 :       GEN s = gel(R,i);
     234       13629 :       n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
     235             :     }
     236          14 :   } else snX = 1;
     237        3661 :   mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
     238       29615 :   while (running || pending)
     239             :   {
     240             :     GEN done;
     241       25954 :     GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
     242       25954 :     mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
     243       25954 :     done = mt_queue_get(&pt, NULL, &pending);
     244       25954 :     if (done)
     245             :     {
     246       22967 :       GEN P = gel(done,1), R = gel(done,2);
     247       22967 :       long j, l = lg(P);
     248      221039 :       for (j=1; j<l; j++)
     249             :       {
     250      198072 :         GEN F = gel(R,j);
     251      198072 :         if (degpol(F) && !gequal0(gel(F,3)))
     252      126938 :           dirmuleuler_large(V, uel(P,j), gel(F,3));
     253             :       }
     254             :     }
     255             :   }
     256        3661 :   mt_queue_end(&pt);
     257        3661 :   return gerepilecopy(av0,V);
     258             : }
     259             : 
     260             : /********************************************************************/
     261             : /**                                                                **/
     262             : /**                 DIRPOWERS and DIRPOWERSSUM                     **/
     263             : /**                                                                **/
     264             : /********************************************************************/
     265             : 
     266             : /* [1^B,...,N^B] */
     267             : GEN
     268         259 : vecpowuu(long N, ulong B)
     269             : {
     270             :   GEN v;
     271             :   long p, i;
     272             :   forprime_t T;
     273             : 
     274         259 :   if (B <= 2)
     275             :   {
     276          77 :     if (!B) return const_vec(N,gen_1);
     277          70 :     v = cgetg(N+1, t_VEC); if (N == 0) return v;
     278          70 :     gel(v,1) = gen_1;
     279          70 :     if (B == 1)
     280        3815 :       for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
     281             :     else
     282         644 :       for (i = 2; i <= N; i++) gel(v,i) = sqru(i);
     283          70 :     return v;
     284             :   }
     285         182 :   v = const_vec(N, NULL);
     286         182 :   u_forprime_init(&T, 3, N);
     287        3934 :   while ((p = u_forprime_next(&T)))
     288             :   {
     289             :     long m, pk, oldpk;
     290        3752 :     gel(v,p) = powuu(p, B);
     291        8365 :     for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
     292             :     {
     293        4613 :       if (pk != p) gel(v,pk) = mulii(gel(v,oldpk), gel(v,p));
     294       20797 :       for (m = N/pk; m > 1; m--)
     295       16184 :         if (gel(v,m) && m%p) gel(v, m*pk) = mulii(gel(v,m), gel(v,pk));
     296             :     }
     297             :   }
     298         182 :   gel(v,1) = gen_1;
     299        7602 :   for (i = 2; i <= N; i+=2)
     300             :   {
     301        7420 :     long vi = vals(i);
     302        7420 :     gel(v,i) = shifti(gel(v,i >> vi), B * vi);
     303             :   }
     304         182 :   return v;
     305             : }
     306             : /* [1^B,...,N^B] */
     307             : GEN
     308        2443 : vecpowug(long N, GEN B, long prec)
     309             : {
     310        2443 :   GEN v, pow = NULL;
     311        2443 :   long p, precp = 2, eB, prec0;
     312             :   forprime_t T;
     313        2443 :   if (N == 1) return mkvec(gen_1);
     314        2443 :   eB = gexpo(B);
     315        2443 :   prec0 = eB < 5? prec: prec + nbits2extraprec(eB);
     316        2443 :   u_forprime_init(&T, 2, N);
     317        2443 :   v = const_vec(N, NULL);
     318        2443 :   gel(v,1) = gen_1;
     319      529634 :   while ((p = u_forprime_next(&T)))
     320             :   {
     321             :     long m, pk, oldpk;
     322      527191 :     if (!pow)
     323        2443 :       pow = gpow(utor(p,prec0), B, prec);
     324             :     else
     325             :     {
     326      524748 :       GEN t = gpow(divru(utor(p,prec0), precp), B, prec);
     327      524748 :       pow = gmul(pow, t); /* (p / precp)^B * precp^B */
     328             :     }
     329      527191 :     precp = p;
     330      527191 :     gel(v,p) = pow; /* p^B */
     331      527191 :     if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
     332     1093715 :     for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
     333             :     {
     334      566524 :       if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
     335    13292209 :       for (m = N/pk; m > 1; m--)
     336    12725685 :         if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
     337             :     }
     338             :   }
     339        2443 :   return v;
     340             : }
     341             : 
     342             : GEN
     343         266 : dirpowers(long n, GEN x, long prec)
     344             : {
     345         266 :   pari_sp av = avma;
     346             :   GEN v;
     347         266 :   if (n <= 0) return cgetg(1, t_VEC);
     348         252 :   if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0)
     349           7 :   {
     350          35 :     ulong B = itou(x);
     351          35 :     v = vecpowuu(n, B);
     352          35 :     if (B <= 2) return v;
     353             :   }
     354         217 :   else v = vecpowug(n, x, prec);
     355         224 :   return gerepilecopy(av, v);
     356             : }
     357             : 
     358             : /* P = prime divisors of (squarefree) n */
     359             : static GEN
     360      296824 : smallfact(ulong n, GEN P, ulong sq, GEN V)
     361             : {
     362      296824 :   long i, l = lg(P);
     363             :   ulong p;
     364             :   GEN c;
     365      296824 :   if (l == 1) return gen_1;
     366      288480 :   p = uel(P,l - 1); if (p > sq) return NULL;
     367       49861 :   c = gel(V, p);
     368       66885 :   for (i = l-2; i > 0; i--)
     369             :   {
     370       60088 :     n /= p; if (n <= sq) return gmul(c, gel(V,n));
     371       17024 :     p = uel(P, i); c = gmul(c, gel(V,p));
     372             :   }
     373        6797 :   return c;
     374             : }
     375             : /* sum_{n <= N} n^s. */
     376             : GEN
     377        8372 : dirpowerssum(ulong N, GEN s, long prec)
     378             : {
     379        8372 :   const ulong step = 2048;
     380        8372 :   pari_sp av = avma, av2;
     381             :   GEN P, V, W, F, c2, c3, c6, c12, c123, c1234, tmp, S;
     382             :   forprime_t T;
     383             :   ulong x1, n, sq, p, precp;
     384             :   long prec0;
     385             : 
     386        8372 :   if (!N) return gen_0;
     387        8372 :   if (N < 9UL) return gerepileupto(av, RgV_sum(dirpowers(N, s, prec)));
     388        8344 :   sq = usqrt(N);
     389        8344 :   V = cgetg(sq+1, t_VEC);
     390        8344 :   W = cgetg(sq+1, t_VEC);
     391        8344 :   F = cgetg(sq+1, t_VEC);
     392        8344 :   prec0 = prec + EXTRAPRECWORD;
     393        8344 :   s = gprec_w(s, prec0);
     394        8344 :   gel(V,1) = gel(W,1) = gel(F,1) = gen_1;
     395       46037 :   for (n = 2; n <= sq; n++)
     396             :   {
     397       37693 :     GEN t = divru(utor(n, prec0), n-1);
     398       37693 :     gel(V,n) = gmul(gel(V,n-1), gpow(t, s, prec0));
     399       37693 :     gel(W,n) = gadd(gel(W,n-1), gel(V,n));
     400       37693 :     gel(F,n) = gadd(gel(F,n-1), gsqr(gel(V,n)));
     401             :   }
     402        8344 :   c2 = gel(V,2); c3 = gel(V,3); c6 = gmul(c2, c3);
     403        8344 :   c12 = gaddgs(c2, 1);
     404        8344 :   c123 = gadd(c12, c3);
     405        8344 :   c1234 = gadd(gel(F,2), gadd(c2,c3));
     406        8344 :   precp = 0; tmp = NULL; S = gen_0;
     407        8344 :   u_forprime_init(&T, sq + 1, N);
     408        8344 :   av2 = avma;
     409      139865 :   while ((p = u_forprime_next(&T)))
     410             :   {
     411      131521 :     GEN t = utor(p, prec0);
     412      131521 :     if (precp) t = divru(t, precp);
     413      131521 :     t = gpow(t, s, prec0);
     414      131521 :     tmp = precp? gmul(tmp, t): t; /* p^s */
     415      131521 :     S = gadd(S, gmul(gel(W, N / p), tmp));
     416      131521 :     precp = p;
     417      131521 :     if ((p & 0x1ff) == 1) S = gerepileupto(av2, S);
     418             :   }
     419        8344 :   P = mkvecsmall2(2, 3);
     420        8344 :   av2 = avma;
     421        8344 :   for(x1 = 1;; x1 += step)
     422         329 :   { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
     423        8673 :     ulong j, lv, x2 = (N >= 2*step && N - 2*step >= x1)? x1-1 + step: N;
     424        8673 :     GEN v = vecfactorsquarefreeu_coprime(x1, x2, P);
     425        8673 :     lv = lg(v);
     426      971963 :     for (j = 1; j < lv; j++) if (gel(v,j))
     427             :     {
     428      296824 :       ulong d = x1-1 + j; /* squarefree, coprime to 6 */
     429             :       ulong q;
     430      296824 :       tmp = smallfact(d, gel(v,j), sq, V);
     431      296824 :       if (!tmp) continue;
     432       58205 :       q = N / d;
     433       58205 :       switch(q = N / d)
     434             :       {
     435       15925 :         case 1: break;
     436        6748 :         case 2: tmp = gmul(tmp, c12); break;
     437        3864 :         case 3: tmp = gmul(tmp, c123); break;
     438        5168 :         case 4:
     439        5168 :         case 5: tmp = gmul(tmp, c1234); break;
     440       26500 :         default:
     441             :         {
     442       26500 :           GEN a = gel(F, usqrt(q)), b = gel(F, usqrt(q / 2));
     443       26500 :           GEN c = gel(F, usqrt(q / 3)), d = gel(F, usqrt(q / 6));
     444       26500 :           tmp = gmul(tmp, gadd(gadd(a, gmul(c2, b)),
     445             :                                gadd(gmul(c3, c), gmul(c6, d))));
     446             :         }
     447             :       }
     448       58205 :       S = gadd(S, tmp);
     449             :     }
     450        8673 :     if (x2 == N) break;
     451         329 :     S = gerepileupto(av2, S);
     452             :   }
     453        8344 :   return gerepilecopy(av, S);
     454             : }
     455             : GEN
     456          35 : dirpowerssum0(GEN N, GEN s, long prec)
     457             : {
     458          35 :   if (typ(N) != t_INT) pari_err_TYPE("dirpowerssum", N);
     459          28 :   if (signe(N) <= 0) return gen_0;
     460          14 :   return dirpowerssum(itou(N), s, prec);
     461             : }

Generated by: LCOV version 1.13