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.14.0 lcov report (development 27775-aca467eab2) Lines: 313 342 91.5 %
Date: 2022-07-03 07:33:15 Functions: 21 21 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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**           Dirichlet series through Euler product               **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : static void
      24          28 : err_direuler(GEN x)
      25          28 : { pari_err_DOMAIN("direuler","constant term","!=", gen_1,x); }
      26             : 
      27             : /* s = t_POL (tolerate t_SER of valuation 0) of constant term = 1
      28             :  * d = minimal such that p^d > X
      29             :  * V indexed by 1..X will contain the a_n
      30             :  * v[1..n] contains the indices nj such that V[nj] != 0 */
      31             : static long
      32       23247 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
      33             : {
      34       23247 :   long i, j, m = n, D = minss(d+2, lg(s));
      35       23247 :   ulong q = 1, X = lg(V)-1;
      36             : 
      37       79982 :   for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
      38             :   {
      39       56735 :     GEN aq = gel(s,i);
      40       56735 :     if (gequal0(aq)) continue;
      41             :     /* j = 1 */
      42       48559 :     gel(V,q) = aq;
      43       48559 :     v[++n] = q;
      44     2283547 :     for (j = 2; j <= m; j++)
      45             :     {
      46     2234988 :       ulong nj = umuluu_le(uel(v,j), q, X);
      47     2234988 :       if (!nj) continue;
      48      167566 :       gel(V,nj) = gmul(aq, gel(V,v[j]));
      49      167566 :       v[++n] = nj;
      50             :     }
      51             :   }
      52       23247 :   return n;
      53             : }
      54             : 
      55             : /* ap != 0 for efficiency, p > sqrt(X) */
      56             : static void
      57      193249 : dirmuleuler_large(GEN V, ulong p, GEN ap)
      58             : {
      59      193249 :   long j, jp, X = lg(V)-1;
      60      193249 :   gel(V,p) = ap;
      61      789404 :   for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
      62      193249 : }
      63             : 
      64             : static ulong
      65        9877 : direulertou(GEN a, GEN fl(GEN))
      66             : {
      67        9877 :   if (typ(a) != t_INT)
      68             :   {
      69          49 :     a = fl(a);
      70          28 :     if (typ(a) != t_INT) pari_err_TYPE("direuler", a);
      71             :   }
      72        9856 :   return signe(a)<=0 ? 0: itou(a);
      73             : }
      74             : 
      75             : static GEN
      76        3696 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
      77             : {
      78        3696 :   long i, l = lg(Sbad);
      79        3696 :   ulong X = lg(V)-1;
      80        3696 :   GEN pbad = gen_1;
      81        9562 :   for (i = 1; i < l; i++)
      82             :   {
      83        5901 :     GEN ai = gel(Sbad,i);
      84             :     ulong q;
      85        5901 :     if (typ(ai) != t_VEC || lg(ai) != 3)
      86          14 :       pari_err_TYPE("direuler [bad primes]",ai);
      87        5887 :     q = gtou(gel(ai,1));
      88        5880 :     if (q <= X)
      89             :     {
      90        4774 :       long d = ulogint(X, q) + 1;
      91        4774 :       GEN s = direuler_factor(gel(ai,2), d);
      92        4760 :       *n = dirmuleuler_small(V, v, *n, q, s, d);
      93        4760 :       pbad = muliu(pbad, q);
      94             :     }
      95             :   }
      96        3661 :   return pbad;
      97             : }
      98             : 
      99             : GEN
     100         504 : direuler_bad(void *E, GEN (*eval)(void *,GEN,long), GEN a,GEN b,GEN c, GEN Sbad)
     101             : {
     102             :   ulong au, bu, X, sqrtX, n, p;
     103         504 :   pari_sp av0 = avma;
     104             :   GEN gp, v, V;
     105             :   forprime_t T;
     106         504 :   au = direulertou(a, gceil);
     107         497 :   bu = direulertou(b, gfloor);
     108         490 :   X = c ? direulertou(c, gfloor): bu;
     109         483 :   if (X == 0) return cgetg(1,t_VEC);
     110         476 :   if (bu > X) bu = X;
     111         476 :   if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
     112         462 :   v = vecsmall_ei(X, 1);
     113         462 :   V = vec_ei(X, 1);
     114         462 :   n = 1;
     115         462 :   if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
     116         427 :   p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
     117        2765 :   while (p <= sqrtX && (p = u_forprime_next(&T)))
     118        2359 :     if (!Sbad || umodiu(Sbad, p))
     119             :     {
     120        2254 :       long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
     121             :       GEN s;
     122        2254 :       gp[2] = p; s = eval(E, gp, d);
     123        2233 :       n = dirmuleuler_small(V, v, n, p, s, d);
     124             :     }
     125       57533 :   while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
     126       57127 :     if (!Sbad || umodiu(Sbad, p))
     127             :     {
     128             :       GEN s;
     129       57120 :       gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
     130       57120 :       if (lg(s) > 3 && !gequal0(gel(s,3)))
     131       23611 :         dirmuleuler_large(V, p, gel(s,3));
     132             :     }
     133         406 :   return gerepilecopy(av0,V);
     134             : }
     135             : 
     136             : /* return a t_SER or a truncated t_POL to precision n */
     137             : GEN
     138       64148 : direuler_factor(GEN s, long n)
     139             : {
     140       64148 :   long t = typ(s);
     141       64148 :   if (is_scalar_t(t))
     142             :   {
     143       33194 :     if (!gequal1(s)) err_direuler(s);
     144       33180 :     return scalarpol_shallow(s,0);
     145             :   }
     146       30954 :   switch(t)
     147             :   {
     148        5705 :     case t_POL: break; /* no need to RgXn_red */
     149       24934 :     case t_RFRAC:
     150             :     {
     151       24934 :       GEN p = gel(s,1), q = gel(s,2);
     152       24934 :       q = RgXn_red_shallow(q,n);
     153       24934 :       s = RgXn_inv(q, n);
     154       24934 :       if (typ(p) == t_POL && varn(p) == varn(q))
     155             :       {
     156          28 :         p = RgXn_red_shallow(p, n);
     157          28 :         s = RgXn_mul(s, p, n);
     158             :       }
     159             :       else
     160       24906 :         if (!gequal1(p)) s = RgX_Rg_mul(s, p);
     161       24934 :       if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
     162       24920 :       break;
     163             :     }
     164         308 :     case t_SER:
     165         308 :       if (!signe(s) || valp(s) || !gequal1(gel(s,2))) err_direuler(s);
     166         308 :       break;
     167           7 :     default: pari_err_TYPE("direuler", s);
     168             :   }
     169       30933 :   return s;
     170             : }
     171             : 
     172             : struct eval_bad
     173             : {
     174             :   void *E;
     175             :   GEN (*eval)(void *, GEN);
     176             : };
     177             : static GEN
     178         420 : eval_bad(void *E, GEN p, long n)
     179             : {
     180         420 :   struct eval_bad *d = (struct eval_bad*) E;
     181         420 :   return direuler_factor(d->eval(d->E, p), n);
     182             : }
     183             : GEN
     184         133 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
     185             : {
     186             :   struct eval_bad d;
     187         133 :   d.E= E; d.eval = eval;
     188         133 :   return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
     189             : }
     190             : 
     191             : static GEN
     192       31066 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
     193             : {
     194       31066 :   GEN P = cgetg(n+1, t_VECSMALL);
     195             :   long i, j;
     196      302386 :   for (i = 1, j = 1; i <= n; i++)
     197             :   {
     198      275366 :     ulong p = u_forprime_next(T);
     199      275366 :     if (!p) { *running = 0; break; }
     200      271320 :     if (Sbad && umodiu(Sbad, p)==0) continue;
     201      266672 :     uel(P,j++) = p;
     202             :   }
     203       31066 :   setlg(P, j);
     204       31066 :   return P;
     205             : }
     206             : 
     207             : GEN
     208        4053 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
     209             : {
     210             :   ulong au, bu, X, sqrtX, n, snX, nX;
     211        4053 :   pari_sp av0 = avma;
     212             :   GEN v, V;
     213             :   forprime_t T;
     214             :   struct pari_mt pt;
     215        4053 :   long running = 1, pending = 0;
     216        4053 :   au = direulertou(a, gceil);
     217        4053 :   bu = direulertou(b, gfloor);
     218        4053 :   X = c ? direulertou(c, gfloor): bu;
     219        4053 :   if (X == 0) return cgetg(1,t_VEC);
     220        4053 :   if (bu > X) bu = X;
     221        4053 :   if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
     222        4046 :   v = vecsmall_ei(X, 1);
     223        4046 :   V = vec_ei(X, 1);
     224        4046 :   n = 1;
     225        4046 :   if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
     226        4046 :   sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
     227        4046 :   if (snX)
     228             :   {
     229        4032 :     GEN P = primelist(&T, Sbad, snX, &running);
     230        4032 :     GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
     231        4032 :     long i, l = lg(P);
     232       20286 :     for (i = 1; i < l; i++)
     233             :     {
     234       16254 :       GEN s = gel(R,i);
     235       16254 :       n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
     236             :     }
     237          14 :   } else snX = 1;
     238        4046 :   mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
     239       34614 :   while (running || pending)
     240             :   {
     241             :     GEN done;
     242       30568 :     GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
     243       30568 :     mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
     244       30568 :     done = mt_queue_get(&pt, NULL, &pending);
     245       30568 :     if (done)
     246             :     {
     247       27034 :       GEN P = gel(done,1), R = gel(done,2);
     248       27034 :       long j, l = lg(P);
     249      277452 :       for (j=1; j<l; j++)
     250             :       {
     251      250418 :         GEN F = gel(R,j);
     252      250418 :         if (degpol(F) && !gequal0(gel(F,3)))
     253      169638 :           dirmuleuler_large(V, uel(P,j), gel(F,3));
     254             :       }
     255             :     }
     256             :   }
     257        4046 :   mt_queue_end(&pt);
     258        4046 :   return gerepilecopy(av0,V);
     259             : }
     260             : 
     261             : /********************************************************************/
     262             : /**                                                                **/
     263             : /**                 DIRPOWERS and DIRPOWERSSUM                     **/
     264             : /**                                                                **/
     265             : /********************************************************************/
     266             : 
     267             : /* [1^B,...,N^B] */
     268             : GEN
     269         644 : vecpowuu(long N, ulong B)
     270             : {
     271             :   GEN v;
     272             :   long p, i;
     273             :   forprime_t T;
     274             : 
     275         644 :   if (B <= 8000)
     276             :   {
     277         644 :     if (!B) return const_vec(N,gen_1);
     278         637 :     v = cgetg(N+1, t_VEC); if (N == 0) return v;
     279         637 :     gel(v,1) = gen_1;
     280         637 :     if (B == 1)
     281       92106 :       for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
     282         462 :     else if (B == 2)
     283             :     {
     284             :       ulong o, s;
     285         266 :       if (N & HIGHMASK)
     286           0 :         for (i = 2, o = 3; i <= N; i++, o += 2)
     287           0 :           gel(v,i) = addiu(gel(v,i-1), o);
     288             :       else
     289       32145 :         for (i = 2, s = 1, o = 3; i <= N; i++, s += o, o += 2)
     290       31879 :           gel(v,i) = utoipos(s + o);
     291             :     }
     292         196 :     else if (B == 3)
     293         840 :       for (i = 2; i <= N; i++) gel(v,i) = powuu(i, B);
     294             :     else
     295             :     {
     296         182 :       long k, Bk, e = expu(N);
     297        7553 :       for (i = 3; i <= N; i += 2) gel(v,i) = powuu(i, B);
     298        1239 :       for (k = 1; k <= e; k++)
     299             :       {
     300        1057 :         N >>= 1; Bk = B * k;
     301        8498 :         for (i = 1; i <= N; i += 2) gel(v, i << k) = shifti(gel(v, i), Bk);
     302             :       }
     303             :     }
     304         637 :     return v;
     305             :   }
     306           0 :   v = const_vec(N, NULL);
     307           0 :   u_forprime_init(&T, 3, N);
     308           0 :   while ((p = u_forprime_next(&T)))
     309             :   {
     310             :     long m, pk, oldpk;
     311           0 :     gel(v,p) = powuu(p, B);
     312           0 :     for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
     313             :     {
     314           0 :       if (pk != p) gel(v,pk) = mulii(gel(v,oldpk), gel(v,p));
     315           0 :       for (m = N/pk; m > 1; m--)
     316           0 :         if (gel(v,m) && m%p) gel(v, m*pk) = mulii(gel(v,m), gel(v,pk));
     317             :     }
     318             :   }
     319           0 :   gel(v,1) = gen_1;
     320           0 :   for (i = 2; i <= N; i+=2)
     321             :   {
     322           0 :     long vi = vals(i);
     323           0 :     gel(v,i) = shifti(gel(v,i >> vi), B * vi);
     324             :   }
     325           0 :   return v;
     326             : }
     327             : 
     328             : /* does n^s require log(x) ? */
     329             : static long
     330       11625 : get_needlog(GEN s)
     331             : {
     332       11625 :   switch(typ(s))
     333             :   {
     334         294 :     case t_REAL: return 2; /* yes but not powcx */
     335        8316 :     case t_COMPLEX: return 1; /* yes using powcx */
     336        3015 :     default: return 0; /* no */
     337             :   }
     338             : }
     339             : /* [1^B,...,N^B] */
     340             : GEN
     341       11751 : vecpowug(long N, GEN B, long prec)
     342             : {
     343       11751 :   GEN v, logp = NULL;
     344       11751 :   long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     345       11751 :   long p, precp = 2, prec0, prec1, needlog;
     346             :   forprime_t T;
     347       11751 :   if (N == 1) return mkvec(gen_1);
     348       11744 :   if (typ(B) == t_INT && lgefint(B) <= 3 && signe(B) >= 0)
     349         133 :     return vecpowuu(N, itou(B));
     350       11611 :   needlog = get_needlog(B);
     351       11611 :   prec1 = prec0 = prec;
     352       11611 :   if (needlog == 1) prec1 = powcx_prec(log2((double)N), B, prec);
     353       11611 :   u_forprime_init(&T, 2, N);
     354       11611 :   v = const_vec(N, NULL);
     355       11611 :   gel(v,1) = gen_1;
     356      758716 :   while ((p = u_forprime_next(&T)))
     357             :   {
     358             :     long m, pk, oldpk;
     359             :     GEN u;
     360      747105 :     gp[2] = p;
     361      747105 :     if (needlog)
     362             :     {
     363       90613 :       if (!logp)
     364       17220 :         logp = logr_abs(utor(p, prec1));
     365             :       else
     366             :       { /* Assuming p and precp are odd,
     367             :          * log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
     368       73393 :         ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
     369       73393 :         GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
     370       73393 :         shiftr_inplace(z, 1); logp = addrr(logp, z);
     371             :       }
     372       87089 :       u = needlog == 1? powcx(gp, logp, B, prec0)
     373       90613 :                       : mpexp(gmul(B, logp));
     374       90613 :       if (p == 2) logp = NULL; /* reset: precp must be odd */
     375             :     }
     376             :     else
     377      656492 :       u = gpow(gp, B, prec0);
     378      747105 :     precp = p;
     379      747105 :     gel(v,p) = u; /* p^B */
     380      747105 :     if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
     381     1598707 :     for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
     382             :     {
     383      851602 :       if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
     384    16711815 :       for (m = N/pk; m > 1; m--)
     385    15860213 :         if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
     386             :     }
     387             :   }
     388       11611 :   return v;
     389             : }
     390             : 
     391             : GEN
     392         602 : dirpowers(long n, GEN x, long prec)
     393             : {
     394             :   pari_sp av;
     395             :   GEN v;
     396         602 :   if (n <= 0) return cgetg(1, t_VEC);
     397         588 :   av = avma; v = vecpowug(n, x, prec);
     398         588 :   if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0 && cmpiu(x, 2) <= 0)
     399         105 :     return v;
     400         483 :   return gerepilecopy(av, v);
     401             : }
     402             : 
     403             : /* P = prime divisors of (squarefree) n, V[i] = i^s for i <= sq.
     404             :  * Return NULL if n is not sq-smooth, else f(n)n^s */
     405             : static GEN
     406      234087 : smallfact(ulong n, GEN P, ulong sq, GEN V)
     407             : {
     408             :   long i, l;
     409             :   ulong p, m, o;
     410             :   GEN c;
     411      234087 :   if (n <= sq) return gel(V,n);
     412      233198 :   l = lg(P); m = p = uel(P, l-1); if (p > sq) return NULL;
     413       46739 :   for (i = l-2; i > 1; i--, m = o) { p = uel(P,i); o = m*p; if (o > sq) break; }
     414       46431 :   c = gel(V,m); n /= m; /* m <= sq, o = m * p > sq */
     415       46431 :   if (n > sq) { c = gmul(c, gel(V,p)); n /= p; }
     416       46431 :   return gmul(c, gel(V,n));
     417             : }
     418             : static GEN
     419         105 : Qtor(GEN x, long prec)
     420         105 : { return typ(x) == t_FRAC? fractor(x, prec): x; }
     421             : /* sum_{n <= N} f(n)n^s. */
     422             : GEN
     423          91 : dirpowerssumfun(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
     424             :                 long prec)
     425             : {
     426          91 :   const ulong step = 2048;
     427          91 :   pari_sp av = avma, av2;
     428             :   GEN P, V, W, Q, c2, Q2, Q3, Q6, S, Z, logp;
     429             :   forprime_t T;
     430          91 :   long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     431             :   ulong a, b, c, e, q, x1, n, sq, p, precp;
     432             :   long prec0, prec1, needlog;
     433             : 
     434          91 :   if (!N) return gen_0;
     435          91 :   if (f)
     436             :   {
     437          14 :     if (N < 49)
     438             :     {
     439           7 :       V = vecpowug(N, s, prec);
     440           7 :       S = gen_1;
     441          28 :       for (n = 2; n <= N; n++)
     442          21 :         S = gadd(S, gmul(gel(V,n), f(E, n, prec)));
     443           7 :       return gerepileupto(av, Qtor(S, prec));
     444             :     }
     445             :   }
     446          77 :   else if (N < 1000UL)
     447          70 :     return gerepileupto(av, Qtor(RgV_sum(vecpowug(N, s, prec)), prec));
     448          14 :   sq = usqrt(N);
     449          14 :   V = cgetg(sq+1, t_VEC);
     450          14 :   W = cgetg(sq+1, t_VEC);
     451          14 :   Q = cgetg(sq+1, t_VEC);
     452          14 :   prec1 = prec0 = prec + EXTRAPRECWORD;
     453          14 :   s = gprec_w(s, prec0);
     454          14 :   needlog = get_needlog(s);
     455          14 :   if (needlog == 1) prec1 = powcx_prec(log2((double)N), s, prec);
     456          14 :   gel(V,1) = gel(W,1) = gel(Q,1) = gen_1;
     457          14 :   c2 = gpow(gen_2, s, prec0);
     458          14 :   if (f) c2 = gmul(c2, f(E, 2, prec));
     459          14 :   gel(V,2) = c2; /* f(2) 2^s */
     460          14 :   gel(W,2) = Qtor(gaddgs(c2, 1), prec0);
     461          14 :   gel(Q,2) = Qtor(gaddgs(gsqr(c2), 1), prec0);
     462          14 :   logp = NULL;
     463        2898 :   for (n = 3; n <= sq; n++)
     464             :   {
     465             :     GEN u;
     466        2884 :     if (odd(n))
     467             :     {
     468        1442 :       gp[2] = n;
     469        1442 :       if (needlog)
     470             :       {
     471           0 :         if (!logp)
     472           0 :           logp = logr_abs(utor(n, prec1));
     473             :         else
     474             :         { /* log n = log(n-2) + 2 atanh(1 / (n - 1)) */
     475           0 :           GEN z = atanhuu(1, n - 1, prec1);
     476           0 :           shiftr_inplace(z, 1); logp = addrr(logp, z);
     477             :         }
     478           0 :         u = needlog == 1? powcx(gp, logp, s, prec0)
     479           0 :                         : mpexp(gmul(s, logp));
     480             : 
     481             :       }
     482             :       else
     483        1442 :         u = gpow(gp, s, prec0);
     484        1442 :       if (f) u = gmul(u, f(E, n, prec0)); /* f(n) n^s */
     485             :     }
     486             :     else
     487        1442 :       u = gmul(c2, gel(V, n>> 1));
     488        2884 :     gel(V,n) = u; /* f(n) n^s */
     489        2884 :     gel(W,n) = gadd(gel(W,n-1), gel(V,n));       /* = sum_{i<=n} f(i)i^s */
     490        2884 :     gel(Q,n) = gadd(gel(Q,n-1), gsqr(gel(V,n))); /* = sum_{i<=n} f(i^2)i^2s */
     491             :   }
     492          14 :   Q2 = RgV_Rg_mul(Q, gel(V,2));
     493          14 :   Q3 = RgV_Rg_mul(Q, gel(V,3));
     494          14 :   Q6 = RgV_Rg_mul(Q, gel(V,6));
     495          14 :   precp = 0; logp = NULL; S = gen_0;
     496          14 :   u_forprime_init(&T, sq + 1, N);
     497          14 :   av2 = avma;
     498       75131 :   while ((p = u_forprime_next(&T)))
     499             :   {
     500             :     GEN u;
     501       75117 :     gp[2] = p;
     502       75117 :     if (needlog)
     503             :     {
     504           0 :       if (!logp)
     505           0 :         logp = logr_abs(utor(p, prec1));
     506             :       else
     507             :       { /* log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
     508           0 :         ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
     509           0 :         GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
     510           0 :         shiftr_inplace(z, 1); logp = addrr(logp, z);
     511             :       }
     512           0 :       u = needlog == 1? powcx(gp, logp, s, prec0)
     513           0 :                       : mpexp(gmul(s, logp));
     514             :     }
     515             :     else
     516       75117 :       u = gpow(gp, s, prec0);
     517       75117 :     if (f) u = gmul(u, f(E, p, prec0)); /* f(p) p^s */
     518       75117 :     S = gadd(S, gmul(gel(W, N / p), u));
     519       75117 :     precp = p;
     520       75117 :     if ((p & 0x1ff) == 1)
     521             :     {
     522         280 :       if (!logp) S = gerepileupto(av2,S);
     523           0 :       else gerepileall(av2, 2, &S, &logp);
     524             :     }
     525             :   }
     526          14 :   P = mkvecsmall2(2, 3);
     527          14 :   Z = cgetg(sq+1, t_VEC);
     528             :   /* a,b,c,e = sqrt(q), sqrt(q/2), sqrt(q/3), sqrt(q/6)
     529             :    * Z[q] = Q[a] + 2^s Q[b] + 3^s Q[c] + 6^s Q[e], with Q[0] = 0 */
     530          14 :   gel(Z, 1) = gen_1;
     531          14 :   gel(Z, 2) = gel(W, 2);
     532          14 :   gel(Z, 3) = gel(W, 3);
     533          14 :   gel(Z, 4) = gel(Z, 5) = gel(W,4);
     534          14 :   gel(Z, 6) = gel(Z, 7) = gadd(gel(W,4), gel(V,6));
     535          14 :   a = 2; b = c = e = 1;
     536        2828 :   for (q = 8; q <= sq; q++)
     537             :   { /* Gray code: at most one of a,b,c,d differs (by 1) from previous value */
     538        2814 :     GEN z = gel(Z, q - 1);
     539             :     ulong na, nb, nc, ne;
     540        2814 :     if ((na = usqrt(q)) != a)
     541         161 :     { a = na; z = gadd(z, gel(V, na * na)); }
     542        2653 :     else if ((nb = usqrt(q / 2)) != b)
     543         119 :     { b = nb; z = gadd(z, gel(V, 2 * nb * nb)); }
     544        2534 :     else if ((nc = usqrt(q / 3)) != c)
     545          91 :     { c = nc; z = gadd(z, gel(V, 3 * nc * nc)); }
     546        2443 :     else if ((ne = usqrt(q / 6)) != e)
     547          63 :     { e = ne; z = gadd(z, gel(V, 6 * ne * ne)); }
     548        2814 :     gel(Z,q) = z;
     549             :   }
     550          14 :   av2 = avma;
     551          14 :   for(x1 = 1;; x1 += step)
     552         350 :   { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
     553         364 :     ulong j, lv, x2 = (N >= 2*step && N - 2*step >= x1)? x1-1 + step: N;
     554         364 :     GEN v = vecfactorsquarefreeu_coprime(x1, x2, P);
     555         364 :     lv = lg(v);
     556      770364 :     for (j = 1; j < lv; j++) if (gel(v,j))
     557             :     {
     558      234087 :       ulong d = x1-1 + j; /* squarefree, coprime to 6 */
     559      234087 :       GEN t = smallfact(d, gel(v,j), sq, V), u; /* = d^s */
     560      234087 :       if (!t) continue;
     561             :       /* S += d^s * Z[q] */
     562       47320 :       q = N / d;
     563       47320 :       if (q == 1) { S = gadd(S, t); continue; }
     564       30184 :       if (q <= sq) u = gel(Z, q);
     565             :       else
     566             :       {
     567         889 :         a = usqrt(q); b = usqrt(q / 2); c = usqrt(q / 3); e = usqrt(q / 6);
     568         889 :         u = gadd(gadd(gel(Q,a), gel(Q2,b)), gadd(gel(Q3,c), gel(Q6,e)));
     569             :       }
     570       30184 :       S = gadd(S, gmul(t, u));
     571             :     }
     572         364 :     if (x2 == N) break;
     573         350 :     S = gerepileupto(av2, S);
     574             :   }
     575          14 :   return gerepileupto(av, S);
     576             : }
     577             : GEN
     578          77 : dirpowerssum(ulong N, GEN s, long prec)
     579          77 : { return dirpowerssumfun(N, s, NULL, NULL, prec); }
     580             : static GEN
     581        8799 : gp_callUp(void *E, ulong x, long prec)
     582             : {
     583        8799 :   long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     584        8799 :   court[2] = x; return gp_callprec(E, court, prec);
     585             : }
     586             : GEN
     587          56 : dirpowerssum0(GEN N, GEN s, GEN f, long prec)
     588             : {
     589          56 :   if (typ(N) != t_INT) pari_err_TYPE("dirpowerssum", N);
     590          49 :   if (signe(N) <= 0) return gen_0;
     591          35 :   if (!f) return dirpowerssum(itou(N), s, prec);
     592          14 :   if (typ(f) != t_CLOSURE) pari_err_TYPE("dirpowerssum", f);
     593          14 :   return dirpowerssumfun(itou(N), s, (void*)f, gp_callUp, prec);
     594             : }

Generated by: LCOV version 1.13