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.16.2 lcov report (development 29395-ef22f77854) Lines: 630 680 92.6 %
Date: 2024-06-15 09:03:40 Functions: 38 40 95.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       28700 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
      33             : {
      34       28700 :   long i, j, m = n, D = minss(d+2, lg(s));
      35       28700 :   ulong q = 1, X = lg(V)-1;
      36             : 
      37       94724 :   for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
      38             :   {
      39       66024 :     GEN aq = gel(s,i);
      40       66024 :     if (gequal0(aq)) continue;
      41             :     /* j = 1 */
      42       53753 :     gel(V,q) = aq;
      43       53753 :     v[++n] = q;
      44     3268013 :     for (j = 2; j <= m; j++)
      45             :     {
      46     3214260 :       ulong nj = umuluu_le(uel(v,j), q, X);
      47     3214260 :       if (!nj) continue;
      48      192017 :       gel(V,nj) = gmul(aq, gel(V,v[j]));
      49      192017 :       v[++n] = nj;
      50             :     }
      51             :   }
      52       28700 :   return n;
      53             : }
      54             : 
      55             : /* ap != 0 for efficiency, p > sqrt(X) */
      56             : static void
      57      308798 : dirmuleuler_large(GEN V, ulong p, GEN ap)
      58             : {
      59      308798 :   long j, jp, X = lg(V)-1;
      60      308798 :   gel(V,p) = ap;
      61     1506547 :   for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
      62      308798 : }
      63             : 
      64             : static ulong
      65       10269 : direulertou(GEN a, GEN fl(GEN))
      66             : {
      67       10269 :   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       10248 :   return signe(a)<=0 ? 0: itou(a);
      73             : }
      74             : 
      75             : static GEN
      76        3724 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
      77             : {
      78        3724 :   long i, l = lg(Sbad);
      79        3724 :   ulong X = lg(V)-1;
      80        3724 :   GEN pbad = gen_1;
      81        9646 :   for (i = 1; i < l; i++)
      82             :   {
      83        5957 :     GEN ai = gel(Sbad,i);
      84             :     ulong q;
      85        5957 :     if (typ(ai) != t_VEC || lg(ai) != 3)
      86          14 :       pari_err_TYPE("direuler [bad primes]",ai);
      87        5943 :     q = gtou(gel(ai,1));
      88        5936 :     if (q <= X)
      89             :     {
      90        4809 :       long d = ulogint(X, q) + 1;
      91        4809 :       GEN s = direuler_factor(gel(ai,2), d);
      92        4795 :       *n = dirmuleuler_small(V, v, *n, q, s, d);
      93        4795 :       pbad = muliu(pbad, q);
      94             :     }
      95             :   }
      96        3689 :   return pbad;
      97             : }
      98             : 
      99             : GEN
     100         672 : 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         672 :   pari_sp av0 = avma;
     104             :   GEN gp, v, V;
     105             :   forprime_t T;
     106         672 :   au = direulertou(a, gceil);
     107         665 :   bu = direulertou(b, gfloor);
     108         658 :   X = c ? direulertou(c, gfloor): bu;
     109         651 :   if (X == 0) return cgetg(1,t_VEC);
     110         644 :   if (bu > X) bu = X;
     111         644 :   if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
     112         630 :   v = vecsmall_ei(X, 1);
     113         630 :   V = vec_ei(X, 1);
     114         630 :   n = 1;
     115         630 :   if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
     116         595 :   p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
     117        8316 :   while (p <= sqrtX && (p = u_forprime_next(&T)))
     118        7742 :     if (!Sbad || umodiu(Sbad, p))
     119             :     {
     120        7637 :       long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
     121             :       GEN s;
     122        7637 :       gp[2] = p; s = eval(E, gp, d);
     123        7616 :       n = dirmuleuler_small(V, v, n, p, s, d);
     124             :     }
     125      740201 :   while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
     126      739627 :     if (!Sbad || umodiu(Sbad, p))
     127             :     {
     128             :       GEN s;
     129      739620 :       gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
     130      739620 :       if (lg(s) > 3 && !gequal0(gel(s,3)))
     131      139153 :         dirmuleuler_large(V, p, gel(s,3));
     132             :     }
     133         574 :   return gerepilecopy(av0,V);
     134             : }
     135             : 
     136             : /* return a t_SER or a truncated t_POL to precision n */
     137             : GEN
     138      752066 : direuler_factor(GEN s, long n)
     139             : {
     140      752066 :   long t = typ(s);
     141      752066 :   if (is_scalar_t(t))
     142             :   {
     143       33194 :     if (!gequal1(s)) err_direuler(s);
     144       33180 :     return scalarpol_shallow(s,0);
     145             :   }
     146      718872 :   switch(t)
     147             :   {
     148        5712 :     case t_POL: break; /* no need to RgXn_red */
     149      712845 :     case t_RFRAC:
     150             :     {
     151      712845 :       GEN p = gel(s,1), q = gel(s,2);
     152      712845 :       q = RgXn_red_shallow(q,n);
     153      712845 :       s = RgXn_inv(q, n);
     154      712845 :       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      712817 :         if (!gequal1(p)) s = RgX_Rg_mul(s, p);
     161      712845 :       if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
     162      712831 :       break;
     163             :     }
     164         308 :     case t_SER:
     165         308 :       if (!signe(s) || valser(s) || !gequal1(gel(s,2))) err_direuler(s);
     166         308 :       break;
     167           7 :     default: pari_err_TYPE("direuler", s);
     168             :   }
     169      718851 :   return s;
     170             : }
     171             : 
     172             : struct eval_bad
     173             : {
     174             :   void *E;
     175             :   GEN (*eval)(void *, GEN);
     176             : };
     177             : static GEN
     178      688303 : eval_bad(void *E, GEN p, long n)
     179             : {
     180      688303 :   struct eval_bad *d = (struct eval_bad*) E;
     181      688303 :   return direuler_factor(d->eval(d->E, p), n);
     182             : }
     183             : GEN
     184         301 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
     185             : {
     186             :   struct eval_bad d;
     187         301 :   d.E= E; d.eval = eval;
     188         301 :   return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
     189             : }
     190             : 
     191             : static GEN
     192       31157 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
     193             : {
     194       31157 :   GEN P = cgetg(n+1, t_VECSMALL);
     195             :   long i, j;
     196      302631 :   for (i = 1, j = 1; i <= n; i++)
     197             :   {
     198      275548 :     ulong p = u_forprime_next(T);
     199      275548 :     if (!p) { *running = 0; break; }
     200      271474 :     if (Sbad && umodiu(Sbad, p)==0) continue;
     201      266791 :     uel(P,j++) = p;
     202             :   }
     203       31157 :   setlg(P, j);
     204       31157 :   return P;
     205             : }
     206             : 
     207             : GEN
     208        4081 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
     209             : {
     210             :   ulong au, bu, X, sqrtX, n, snX, nX;
     211        4081 :   pari_sp av0 = avma;
     212             :   GEN v, V;
     213             :   forprime_t T;
     214             :   struct pari_mt pt;
     215        4081 :   long running = 1, pending = 0;
     216        4081 :   au = direulertou(a, gceil);
     217        4081 :   bu = direulertou(b, gfloor);
     218        4081 :   X = c ? direulertou(c, gfloor): bu;
     219        4081 :   if (X == 0) return cgetg(1,t_VEC);
     220        4081 :   if (bu > X) bu = X;
     221        4081 :   if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
     222        4074 :   v = vecsmall_ei(X, 1);
     223        4074 :   V = vec_ei(X, 1);
     224        4074 :   n = 1;
     225        4074 :   if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
     226        4074 :   sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
     227        4074 :   if (snX)
     228             :   {
     229        4060 :     GEN P = primelist(&T, Sbad, snX, &running);
     230        4060 :     GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
     231        4060 :     long i, l = lg(P);
     232       20349 :     for (i = 1; i < l; i++)
     233             :     {
     234       16289 :       GEN s = gel(R,i);
     235       16289 :       n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
     236             :     }
     237          14 :   } else snX = 1;
     238        4074 :   mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
     239       34212 :   while (running || pending)
     240             :   {
     241             :     GEN done;
     242       30138 :     GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
     243       30138 :     mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
     244       30138 :     done = mt_queue_get(&pt, NULL, &pending);
     245       30138 :     if (done)
     246             :     {
     247       27097 :       GEN P = gel(done,1), R = gel(done,2);
     248       27097 :       long j, l = lg(P);
     249      277599 :       for (j=1; j<l; j++)
     250             :       {
     251      250502 :         GEN F = gel(R,j);
     252      250502 :         if (degpol(F) && !gequal0(gel(F,3)))
     253      169645 :           dirmuleuler_large(V, uel(P,j), gel(F,3));
     254             :       }
     255             :     }
     256             :   }
     257        4074 :   mt_queue_end(&pt);
     258        4074 :   return gerepilecopy(av0,V);
     259             : }
     260             : 
     261             : /********************************************************************/
     262             : /**                                                                **/
     263             : /**                 DIRPOWERS and DIRPOWERSSUM                     **/
     264             : /**                                                                **/
     265             : /********************************************************************/
     266             : 
     267             : /* [1^B,...,N^B] */
     268             : GEN
     269         686 : vecpowuu(long N, ulong B)
     270             : {
     271             :   GEN v;
     272             :   long p, i;
     273             :   forprime_t T;
     274             : 
     275         686 :   if (B <= 8000)
     276             :   {
     277         686 :     if (!B) return const_vec(N,gen_1);
     278         679 :     v = cgetg(N+1, t_VEC); if (N == 0) return v;
     279         679 :     gel(v,1) = gen_1;
     280         679 :     if (B == 1)
     281       92736 :       for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
     282         469 :     else if (B == 2)
     283             :     {
     284             :       ulong o, s;
     285         273 :       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       31073 :         for (i = 2, s = 1, o = 3; i <= N; i++, s += o, o += 2)
     290       30800 :           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         679 :     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       51014 : get_needlog(GEN s)
     331             : {
     332       51014 :   switch(typ(s))
     333             :   {
     334         294 :     case t_REAL: return 2; /* yes but not powcx */
     335       47586 :     case t_COMPLEX: return 1; /* yes using powcx */
     336        3134 :     default: return 0; /* no */
     337             :   }
     338             : }
     339             : /* [1^B,...,N^B] */
     340             : GEN
     341       12031 : vecpowug(long N, GEN B, long prec)
     342             : {
     343       12031 :   GEN v, logp = NULL;
     344       12031 :   long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     345       12031 :   long p, precp = 2, prec0, prec1, needlog;
     346             :   forprime_t T;
     347       12031 :   if (N == 1) return mkvec(gen_1);
     348       12024 :   if (typ(B) == t_INT && lgefint(B) <= 3 && signe(B) >= 0)
     349         168 :     return vecpowuu(N, itou(B));
     350       11856 :   needlog = get_needlog(B);
     351       11856 :   prec1 = prec0 = prec;
     352       11856 :   if (needlog == 1) prec1 = powcx_prec(log2((double)N), B, prec);
     353       11856 :   u_forprime_init(&T, 2, N);
     354       11856 :   v = const_vec(N, NULL);
     355       11856 :   gel(v,1) = gen_1;
     356     1566194 :   while ((p = u_forprime_next(&T)))
     357             :   {
     358             :     long m, pk, oldpk;
     359             :     GEN u;
     360     1554338 :     gp[2] = p;
     361     1554338 :     if (needlog)
     362             :     {
     363       96234 :       if (!logp)
     364       17486 :         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       78748 :         ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
     369       78748 :         GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
     370       78748 :         shiftr_inplace(z, 1); logp = addrr(logp, z);
     371             :       }
     372       92710 :       u = needlog == 1? powcx(gp, logp, B, prec0)
     373       96234 :                       : mpexp(gmul(B, logp));
     374       96234 :       if (p == 2) logp = NULL; /* reset: precp must be odd */
     375             :     }
     376             :     else
     377     1458104 :       u = gpow(gp, B, prec0);
     378     1554338 :     precp = p;
     379     1554338 :     gel(v,p) = u; /* p^B */
     380     1554338 :     if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
     381     3222462 :     for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
     382             :     {
     383     1668124 :       if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
     384    46352335 :       for (m = N/pk; m > 1; m--)
     385    44684211 :         if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
     386             :     }
     387             :   }
     388       11856 :   return v;
     389             : }
     390             : 
     391             : GEN
     392         665 : dirpowers(long n, GEN x, long prec)
     393             : {
     394             :   pari_sp av;
     395             :   GEN v;
     396         665 :   if (n <= 0) return cgetg(1, t_VEC);
     397         651 :   av = avma; v = vecpowug(n, x, prec);
     398         651 :   if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0 && cmpiu(x, 2) <= 0)
     399         133 :     return v;
     400         518 :   return gerepilecopy(av, v);
     401             : }
     402             : 
     403             : static GEN
     404         252 : vecmulsqlv(GEN Q, GEN V)
     405             : {
     406             :   long lq, i;
     407             :   GEN W;
     408         252 :   if (typ(V) != t_VEC) return RgV_Rg_mul(Q, V);
     409          21 :   lq = lg(Q); W = cgetg(lq, t_VEC);
     410         672 :   for (i = 1; i < lq; i++) gel(W, i) = vecmul(gel(Q, i), V);
     411          21 :   return W;
     412             : }
     413             : 
     414             : /* P = prime divisors of (squarefree) n, V[i] = i^s for i <= sq.
     415             :  * Return NULL if n is not sq-smooth, else f(n)n^s */
     416             : static GEN
     417    18206141 : smallfact(ulong n, GEN P, ulong sq, GEN V)
     418             : {
     419             :   long i, l;
     420             :   ulong p, m, o;
     421             :   GEN c;
     422    18206141 :   if (n <= sq) return gel(V,n);
     423    18157204 :   l = lg(P); m = p = uel(P, l-1); if (p > sq) return NULL;
     424     3843789 :   for (i = l-2; i > 1; i--, m = o) { p = uel(P,i); o = m*p; if (o > sq) break; }
     425     3810482 :   c = gel(V,m); n /= m; /* m <= sq, o = m * p > sq */
     426     3810482 :   if (n > sq) { c = vecmul(c, gel(V,p)); n /= p; }
     427     3799465 :   return vecmul(c, gel(V,n));
     428             : }
     429             : 
     430             : static GEN
     431         140 : _Qtor(GEN x, long prec)
     432         140 : { return typ(x) == t_FRAC? fractor(x, prec): x; }
     433             : static GEN
     434         315 : Qtor(GEN x, long prec)
     435             : {
     436         315 :   long tx = typ(x);
     437         357 :   if (tx == t_VEC || tx == t_COL) pari_APPLY_same(_Qtor(gel(x, i), prec));
     438         294 :   return tx == t_FRAC? fractor(x, prec): x;
     439             : }
     440             : 
     441             : /* Here N > 0 is small */
     442             : static GEN
     443         147 : naivedirpowerssum(long N, GEN s, void *E, GEN (*f)(void *, ulong, long),
     444             :                   long prec)
     445             : {
     446         147 :   GEN V = vecpowug(N, s, prec), S;
     447         147 :   if (!f) S = RgV_sum(V);
     448             :   else
     449             :   {
     450             :     long n;
     451          35 :     S = f(E, 1, prec);
     452         308 :     for (n = 2; n <= N; n++) S = gadd(S, gmul(gel(V, n), f(E, n, prec)));
     453             :   }
     454         147 :   return Qtor(S, prec);
     455             : }
     456             : 
     457             : static GEN
     458         126 : smalldirpowerssum(long N, GEN s, void *E, GEN (*f)(void *, ulong, long),
     459             :                   long both, long prec)
     460             : {
     461         126 :   GEN S = naivedirpowerssum(N, s, E, f, prec), SB, sb;
     462         126 :   if (!both) return S;
     463          42 :   sb = gconj(gsubsg(-1, s));
     464          42 :   SB = both==2 && gequal(s,sb)? S: gconj(naivedirpowerssum(N,sb,E,f,prec));
     465          42 :   return mkvec2(S, SB);
     466             : }
     467             : 
     468             : static GEN
     469          63 : dirpowsuminit(GEN s, void *E, GEN (*f)(void *, ulong, long), GEN data,
     470             :               long both, long prec)
     471             : {
     472          63 :   GEN onef = gel(data, 1), zervec = gel(data, 2), sqlpp = gel(data, 3);
     473          63 :   long sq = sqlpp[1], needlog = sqlpp[2], prec0 = sqlpp[3], prec1 = sqlpp[4];
     474          63 :   GEN V = cgetg(sq+1, t_VEC), W = cgetg(sq+1, t_VEC), Q = cgetg(sq+1, t_VEC);
     475          63 :   GEN VB = NULL, WB = NULL, QB = NULL, c2, c2B = NULL;
     476          63 :   GEN Q2, Q3, Q6, Q2B = NULL, Q3B = NULL, Q6B = NULL;
     477          63 :   GEN logp, R, RB = NULL;
     478          63 :   long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     479             :   long n;
     480          63 :   if (both == 1 || (both == 2 && !gequal(real_i(s), gneg(ghalf))))
     481          21 :   { VB = cgetg(sq+1, t_VEC); WB = cgetg(sq+1, t_VEC); QB = cgetg(sq+1, t_VEC);}
     482          63 :   gel(V, 1) = gel(W, 1) = gel(Q, 1) = onef;
     483          63 :   if (VB) { gel(VB, 1) = gel(WB, 1) = gel(QB, 1) = onef; }
     484          63 :   c2 = gpow(gen_2, s, prec0); if (VB) c2B = ginv(gmul2n(gconj(c2), 1));
     485          63 :   if (f)
     486             :   {
     487          42 :     GEN tmp2 = f(E, 2, prec);
     488          42 :     c2 = gmul(c2, tmp2); if (VB) c2B = gmul(c2B, tmp2);
     489             :   }
     490          63 :   gel(V,2) = c2; /* f(2) 2^s */
     491          63 :   gel(W,2) = Qtor(gadd(c2, onef), prec0);
     492          63 :   gel(Q,2) = Qtor(gadd(vecsqr(c2), onef), prec0);
     493          63 :   if (VB)
     494             :   {
     495          21 :     gel(VB, 2) = c2B; gel(WB, 2) = Qtor(gadd(c2B, onef), prec0);
     496          21 :     gel(QB, 2) = Qtor(gadd(vecsqr(c2B), onef), prec0);
     497             :   }
     498          63 :   logp = NULL;
     499        4074 :   for (n = 3; n <= sq; n++)
     500             :   {
     501        4011 :     GEN u = NULL, uB = NULL, ks = f ? f(E, n, prec0) : gen_1;
     502        4011 :     long zks = !gequal0(ks);
     503        4011 :     if (odd(n))
     504             :     {
     505        2023 :       gp[2] = n;
     506        2023 :       if (needlog)
     507             :       {
     508         476 :         if (!logp)
     509          42 :           logp = logr_abs(utor(n, prec1));
     510             :         else
     511             :         { /* log n = log(n-2) + 2 atanh(1 / (n - 1)) */
     512         434 :           GEN z = atanhuu(1, n - 1, prec1);
     513         434 :           shiftr_inplace(z, 1); logp = addrr(logp, z);
     514             :         }
     515         476 :         if (zks)
     516         476 :           u = needlog == 1? powcx(gp, logp, s, prec0) : mpexp(gmul(s, logp));
     517             :       }
     518        1547 :       else if (zks) u = gpow(gp, s, prec0);
     519        2023 :       if (zks)
     520             :       {
     521        2009 :         if (VB) uB = gmul(ginv(gmulsg(n, gconj(u))), ks);
     522        2009 :         u = gmul(u, ks); /* f(n) n^s */
     523             :       }
     524             :     }
     525             :     else
     526             :     {
     527        1988 :       u = vecmul(c2, gel(V, n >> 1));
     528        1988 :       if (VB) uB = vecmul(c2B, gel(VB, n >> 1));
     529             :     }
     530        4011 :     if (zks)
     531             :     { /* V[n]=f(n)n^s, W[n]=sum_{i<=n} f(i)i^s, Q[n]=sum_{i<=n} f(i^2)i^2s */
     532        3983 :       gel(V,n) = u;
     533        3983 :       gel(W,n) = gadd(gel(W, n-1), gel(V,n));
     534        3983 :       gel(Q,n) = gadd(gel(Q, n-1), vecsqr(gel(V,n)));
     535        3983 :       if (VB)
     536             :       {
     537         462 :         gel(VB,n) = uB;
     538         462 :         gel(WB,n) = gadd(gel(WB,n-1), gel(VB,n));
     539         462 :         gel(QB,n) = gadd(gel(QB,n-1), vecsqr(gel(VB,n)));
     540             :       }
     541             :     }
     542             :     else
     543             :     {
     544          28 :       gel(V,n) = zervec; gel(W,n) = gel(W, n-1); gel(Q,n) = gel(Q, n-1);
     545          28 :       if (VB)
     546             :       {
     547           0 :         gel(VB,n) = zervec; gel(WB,n) = gel(WB, n-1);
     548           0 :         gel(QB,n) = gel(QB, n-1);
     549             :       }
     550             :     }
     551             :   }
     552          63 :   Q2 = vecmulsqlv(Q, gel(V,2));
     553          63 :   Q3 = vecmulsqlv(Q, gel(V,3));
     554          63 :   Q6 = vecmulsqlv(Q, gel(V,6));
     555          63 :   if (VB)
     556             :   {
     557          21 :     Q2B = vecmulsqlv(QB, gel(VB,2));
     558          21 :     Q3B = vecmulsqlv(QB, gel(VB,3));
     559          21 :     Q6B = vecmulsqlv(QB, gel(VB,6));
     560             :   }
     561          63 :   R = mkvecn(6, V, W, Q, Q2, Q3, Q6);
     562          63 :   if (VB) RB = mkvecn(6, VB, WB, QB, Q2B, Q3B, Q6B);
     563          63 :   return VB ? mkvec2(R, RB) : mkvec(R);
     564             : }
     565             : 
     566             : static GEN
     567          63 : dirpowsumprimeloop(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
     568             :                    GEN data, GEN W, GEN WB)
     569             : {
     570             :   pari_sp av2;
     571          63 :   GEN zervec = gel(data, 2), S = zervec, SB = zervec, logp = NULL;
     572          63 :   GEN sqlpp = gel(data, 3);
     573             :   forprime_t T;
     574          63 :   long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     575          63 :   long p, precp = 0, sq = sqlpp[1], needlog = sqlpp[2];
     576          63 :   long prec0 = sqlpp[3], prec1 = sqlpp[4];
     577          63 :   u_forprime_init(&T, sq + 1, N);
     578          63 :   av2 = avma;
     579       80969 :   while ((p = u_forprime_next(&T)))
     580             :   {
     581       80906 :     GEN u = NULL, ks = f ? f(E, p, prec1) : gen_1;
     582       80906 :     long zks = !gequal0(ks);
     583       80906 :     gp[2] = p;
     584       80906 :     if (needlog)
     585             :     {
     586        4690 :       if (!logp)
     587          42 :         logp = logr_abs(utor(p, prec1));
     588             :       else
     589             :       { /* log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
     590        4648 :         ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
     591        4648 :         GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
     592        4648 :         shiftr_inplace(z, 1); logp = addrr(logp, z);
     593             :       }
     594        4690 :       if (zks)
     595        4690 :         u = needlog == 1? powcx(gp, logp, s, prec0) : mpexp(gmul(s, logp));
     596             :     }
     597       76216 :     else { if (zks) u = gpow(gp, s, prec0); }
     598       80906 :     if (zks)
     599             :     {
     600       80906 :       S = gadd(S, vecmul(gel(W, N / p), gmul(ks, u)));
     601       80906 :       if (WB)
     602        2345 :         SB = gadd(SB, gdiv(vecmul(ks, gel(WB, N / p)), gmulsg(p, gconj(u))));
     603             :     }
     604       80906 :     precp = p;
     605       80906 :     if ((p & 0x1ff) == 1)
     606             :     {
     607         280 :       if (!logp) gerepileall(av2, SB? 2: 1, &S, &SB);
     608           0 :       else gerepileall(av2, SB? 3: 2, &S, &logp, &SB);
     609             :     }
     610             :   }
     611          63 :   return SB ? mkvec2(S, SB) : mkvec(S);
     612             : }
     613             : 
     614             : static GEN
     615         413 : dirpowsumsqfloop(long N, long x1, long x2, long sq, GEN P, GEN allvecs, GEN S,
     616             :                  GEN allvecsB, GEN SB)
     617             : {
     618         413 :   GEN V = gel(allvecs, 1), Q1 = gel(allvecs, 2), Q2 = gel(allvecs, 3);
     619         413 :   GEN Q3 = gel(allvecs, 4), Q6 = gel(allvecs, 5), Z = gel(allvecs, 6);
     620         413 :   GEN VB = NULL, QB = NULL, Q2B = NULL, Q3B = NULL, Q6B = NULL, ZB = NULL;
     621         413 :   GEN v = vecfactorsquarefreeu_coprime(x1, x2, P);
     622         413 :   long lv = lg(v), j;
     623         413 :   if (allvecsB)
     624             :   {
     625          21 :     VB = gel(allvecsB, 1), QB = gel(allvecsB, 2), Q2B = gel(allvecsB, 3);
     626          21 :     Q3B = gel(allvecsB, 4), Q6B = gel(allvecsB, 5), ZB = gel(allvecsB, 6);
     627             :   }
     628      806813 :   for (j = 1; j < lv; j++)
     629      806400 :     if (gel(v,j))
     630             :     {
     631      245126 :       ulong d = x1 - 1 + j; /* squarefree, coprime to 6 */
     632      245126 :       GEN t = smallfact(d, gel(v,j), sq, V), u;
     633      245126 :       GEN tB = NULL, uB = NULL; /* = d^s */
     634             :       long a, b, c, e, q;
     635      245126 :       if (!t || gequal0(t)) continue;
     636       48748 :       if (VB) tB = vecinv(gmulsg(d, gconj(t)));
     637             :       /* warning: gives 1/conj(f(d)) d^(-1-conj(s)), equal to
     638             :          f(d) d^(-1-conj(s)) only if |f(d)|=1. */
     639             :       /* S += f(d) * d^s * Z[q] */
     640       48748 :       q = N / d;
     641       48748 :       if (q == 1)
     642             :       {
     643       17339 :         S = gadd(S, t); if (VB) SB = gadd(SB, tB);
     644       17339 :         continue;
     645             :       }
     646       31409 :       if (q <= sq) { u = gel(Z, q); if (VB) uB = gel(ZB, q); }
     647             :       else
     648             :       {
     649        1274 :         a = usqrt(q); b = usqrt(q / 2); c = usqrt(q / 3); e = usqrt(q / 6);
     650        1274 :         u = gadd(gadd(gel(Q1,a), gel(Q2,b)), gadd(gel(Q3,c), gel(Q6,e)));
     651        1274 :         if (VB)
     652         161 :           uB = gadd(gadd(gel(QB,a), gel(Q2B,b)), gadd(gel(Q3B,c), gel(Q6B,e)));
     653             :       }
     654       31409 :       S = gadd(S, vecmul(t, u)); if (VB) SB = gadd(SB, vecmul(tB, uB));
     655             :     }
     656         413 :   return VB ? mkvec2(S, SB) : mkvec(S);
     657             : }
     658             : 
     659             : static GEN
     660          63 : dirpowsummakez(GEN V, GEN W, GEN VB, GEN WB, GEN onef, ulong sq)
     661             : {
     662          63 :   GEN Z = cgetg(sq+1, t_VEC), ZB = NULL;
     663             :   ulong a, b, c, e, q;
     664             :   /* a,b,c,e = sqrt(q), sqrt(q/2), sqrt(q/3), sqrt(q/6)
     665             :    * Z[q] = Q[a] + 2^s Q[b] + 3^s Q[c] + 6^s Q[e], with Q[0] = 0 */
     666          63 :   gel(Z, 1) = onef;
     667          63 :   gel(Z, 2) = gel(W, 2);
     668          63 :   gel(Z, 3) = gel(W, 3);
     669          63 :   gel(Z, 4) = gel(Z, 5) = gel(W, 4);
     670          63 :   gel(Z, 6) = gel(Z, 7) = gadd(gel(W, 4), gel(V, 6));
     671          63 :   if (VB)
     672             :   {
     673          21 :     ZB = cgetg(sq+1, t_VEC);
     674          21 :     gel(ZB, 1) = onef;
     675          21 :     gel(ZB, 2) = gel(WB, 2);
     676          21 :     gel(ZB, 3) = gel(WB, 3);
     677          21 :     gel(ZB, 4) = gel(ZB, 5) = gel(WB, 4);
     678          21 :     gel(ZB, 6) = gel(ZB, 7) = gadd(gel(WB, 4), gel(VB, 6));
     679             :   }
     680          63 :   a = 2; b = c = e = 1;
     681        3759 :   for (q = 8; q <= sq; q++)
     682             :   { /* Gray code: at most one of a,b,c,d differs (by 1) from previous value */
     683        3696 :     GEN z = gel(Z, q - 1), zB = NULL;
     684             :     ulong na, nb, nc, ne, na2, nb2, nc2, ne2;
     685        3696 :     if (VB) zB = gel(ZB, q - 1);
     686        3696 :     if ((na = usqrt(q)) != a)
     687         280 :     { a = na; na2 = na * na; z = gadd(z, gel(V, na2));
     688         280 :       if (VB) zB = gadd(zB, gel(VB, na2)); }
     689        3416 :     else if ((nb = usqrt(q / 2)) != b)
     690         203 :     { b = nb; nb2 = 2 * nb * nb; z = gadd(z, gel(V, nb2));
     691         203 :       if (VB) zB = gadd(zB, gel(VB, nb2)); }
     692        3213 :     else if ((nc = usqrt(q / 3)) != c)
     693         161 :     { c = nc; nc2 = 3 * nc * nc; z = gadd(z, gel(V, nc2));
     694         161 :       if (VB) zB = gadd(zB, gel(VB, nc2)); }
     695        3052 :     else if ((ne = usqrt(q / 6)) != e)
     696          98 :     { e = ne; ne2 = 6 * ne * ne; z = gadd(z, gel(V, ne2));
     697          98 :       if (VB) zB = gadd(zB, gel(VB, ne2)); }
     698        3696 :     gel(Z, q) = z; if (VB) gel(ZB, q) = zB;
     699             :   }
     700          63 :   return VB ? mkvec2(Z, ZB) : mkvec(Z);
     701             : }
     702             : 
     703             : static const long step = 2048;
     704             : 
     705             : /* both =
     706             :  * 0: sum_{n<=N}f(n)n^s
     707             :  * 1: sum for (f,s) and (conj(f),-1-s)
     708             :  * 2: sum for (f,s) and (f,-1-s), assuming |f(n)| in {0,1} */
     709             : GEN
     710         203 : dirpowerssumfun(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
     711             :                 long both, long prec)
     712             : {
     713         203 :   pari_sp av = avma, av2;
     714             :   GEN P, V, W, Q, Q2, Q3, Q6, S, Z, onef, zervec;
     715         203 :   GEN VB = NULL, WB = NULL, QB = NULL;
     716         203 :   GEN Q2B = NULL, Q3B = NULL, Q6B = NULL, SB = NULL, ZB = NULL;
     717         203 :   GEN R, RS, data, allvecs, allvecsB = NULL;
     718             :   ulong x1, sq;
     719             :   long prec0, prec1, needlog;
     720             : 
     721         203 :   if ((f && N < 49) || (!f && N < 1000))
     722             :   {
     723         140 :     if (!N)
     724             :     {
     725          14 :       if (!f) return gen_0;
     726           0 :       return gerepileupto(av, gmul(gen_0, f(E, 1, prec)));
     727             :     }
     728         126 :     return gerepilecopy(av, smalldirpowerssum(N, s, E, f, both, prec));
     729             :   }
     730          63 :   onef = f ? f(E, 1, prec) : gen_1;
     731          63 :   zervec = gmul(gen_0, onef);
     732          63 :   sq = usqrt(N);
     733          63 :   prec1 = prec0 = prec + EXTRAPREC64;
     734          63 :   s = gprec_w(s, prec0);
     735          63 :   needlog = get_needlog(s);
     736          63 :   if (needlog == 1) prec1 = powcx_prec(log2((double)N), s, prec);
     737          63 :   data = mkvec3(onef, zervec, mkvecsmall4(sq, needlog, prec0, prec1));
     738          63 :   RS = dirpowsuminit(s, E, f, data, both, prec);
     739          63 :   R = gel(RS, 1); V = gel(R, 1); W = gel(R, 2); Q = gel(R, 3);
     740          63 :   Q2 = gel(R, 4); Q3 = gel(R, 5); Q6 = gel(R, 6);
     741          63 :   if (lg(RS) > 2)
     742             :   {
     743          21 :     GEN RB = gel(RS, 2);
     744          21 :     VB = gel(RB, 1); WB = gel(RB, 2); QB = gel(RB, 3);
     745          21 :     Q2B = gel(RB, 4); Q3B = gel(RB, 5); Q6B = gel(RB, 6);
     746             :   }
     747          63 :   RS = dirpowsumprimeloop(N, s, E, f, data, W, WB);
     748          63 :   S = gel(RS, 1); if (VB) SB = gel(RS, 2);
     749          63 :   RS = dirpowsummakez(V, W, VB, WB, onef, sq);
     750          63 :   Z = gel(RS, 1); if (VB) ZB = gel(RS, 2);
     751          63 :   P = mkvecsmall2(2, 3);
     752          63 :   allvecs = mkvecn(6, V, Q, Q2, Q3, Q6, Z);
     753          63 :   if (VB) allvecsB = mkvecn(6, VB, QB, Q2B, Q3B, Q6B, ZB);
     754          63 :   av2 = avma;
     755          63 :   for(x1 = 1;; x1 += step)
     756         350 :   { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
     757         413 :     ulong x2 = (N >= 2*step && N - 2*step >= x1)? x1-1 + step: N;
     758         413 :     RS = dirpowsumsqfloop(N, x1, x2, sq, P, allvecs, S, allvecsB, SB);
     759         413 :     S = gel(RS, 1); if (VB) SB = gel(RS, 2);
     760         413 :     if (x2 == N) break;
     761         350 :     gerepileall(av2, SB? 2: 1, &S, &SB);
     762             :   }
     763          63 :   if (both == 0) return gerepileupto(av, S);
     764          42 :   return gerepilecopy(av, mkvec2(S, gconj(VB? SB: S)));
     765             : }
     766             : 
     767             : GEN
     768         133 : dirpowerssum(ulong N, GEN s, long both, long prec)
     769         133 : { return dirpowerssumfun(N, s, NULL, NULL, both, prec); }
     770             : 
     771             : static GEN
     772       13818 : gp_callUp(void *E, ulong x, long prec)
     773             : {
     774       13818 :   long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     775       13818 :   court[2] = x; return gp_callprec(E, court, prec);
     776             : }
     777             : 
     778             : GEN
     779         154 : dirpowerssum0(GEN N, GEN s, GEN f, long both, long prec)
     780             : {
     781         154 :   if (typ(N) != t_INT) pari_err_TYPE("dirpowerssum", N);
     782         147 :   if (signe(N) <= 0) N = gen_0;
     783         147 :   if (!f) return dirpowerssum(itou(N), s, both, prec);
     784          70 :   if (typ(f) != t_CLOSURE) pari_err_TYPE("dirpowerssum", f);
     785          70 :   return dirpowerssumfun(itou(N), s, (void*)f, gp_callUp, both, prec);
     786             : }
     787             : 
     788             : /*******************************************************************/
     789             : /*                     Parallel dirpowerssumfun                    */
     790             : /*******************************************************************/
     791             : /* f is a totally multiplicative function of modulus 0 or 1
     792             :  * (essentially a Dirichlet character). Compute simultaneously
     793             :  * sum_{0 < n <= N} f(n)n^s and sum_{0 < n <= N} f(n)n^{-1-conj(s)}
     794             :  * Warning: s is conjugated, but not f. Main application for Riemann-Siegel,
     795             :  * where we need R(chi,s) and conj(R(chi,1-conj(s))). */
     796             : 
     797             : static GEN
     798      405385 : mycallvec(GEN f, ulong n, long prec)
     799             : {
     800      405385 :   if (typ(f) == t_CLOSURE) return gp_callUp((void*)f, n, prec);
     801      405385 :   return gel(f, (n-1) % (lg(f)-1) + 1);
     802             : }
     803             : 
     804             : static GEN
     805          70 : _smalldirpowerssum(long N, GEN s, long fl, long prec)
     806             : {
     807          70 :   GEN vg = vecpowug(N, s, prec), S1 = _Qtor(RgV_sum(vg), prec);
     808             :   long n;
     809          70 :   if (!fl) return mkvec2(S1, S1);
     810         616 :   for (n = 1; n <= N; n++) gel(vg, n) = ginv(gmulsg(n, gconj(gel(vg, n))));
     811          28 :   return mkvec2(S1, _Qtor(RgV_sum(vg), prec));
     812             : }
     813             : 
     814             : static GEN
     815        2919 : gmulvecsqlv(GEN Q, GEN V)
     816             : {
     817             :   long lq, i;
     818             :   GEN W;
     819        2919 :   if (typ(V) != t_VEC) return RgV_Rg_mul(Q, V);
     820        1155 :   lq = lg(Q); W = cgetg(lq, t_VEC);
     821       61005 :   for (i = 1; i < lq; i++) gel(W, i) = vecmul(gel(Q, i), V);
     822        1155 :   return W;
     823             : }
     824             : 
     825             : GEN
     826       29132 : parsqfboth_worker(GEN gk, GEN vZ, GEN vVQ, GEN vV, GEN P, GEN Nsq)
     827             : {
     828       29132 :   pari_sp av2 = avma;
     829       29132 :   GEN Z1 = gel(vZ, 1), Z2 = gel(vZ, 2), V1 = gel(vV, 1), V2 = gel(vV, 2);
     830       29132 :   GEN VQ1 = gel(vVQ, 1), VQN, v, S, SB;
     831       29132 :   GEN Q1 = gel(VQ1, 1), Q2 = gel(VQ1, 2), Q3 = gel(VQ1, 3), Q6 = gel(VQ1, 4);
     832       29132 :   GEN QB = NULL, Q2B = NULL, Q3B = NULL, Q6B = NULL;
     833       29132 :   long k = itos(gk), N = Nsq[1];
     834       29132 :   long x1 = 1 + step * k, x2, j, lv, fl = !gcmp0(V2);
     835       29132 :   ulong sq = Nsq[2];
     836       29132 :   if (typ(gel(V1, 1)) == t_VEC)
     837        2016 :   { long lvv = lg(gel(V1, 1)) - 1; S = zerovec(lvv); SB = zerovec(lvv); }
     838       27116 :   else { S = gen_0; SB = gen_0; }
     839       29132 :   if (fl)
     840           0 :   { VQN = gel(vVQ, 2); QB = gel(VQN, 1); Q2B = gel(VQN, 2);
     841           0 :     Q3B = gel(VQN, 3); Q6B = gel(VQN, 4); }
     842             :   /* beware overflow, fuse last two bins (avoid a tiny remainder) */
     843       29132 :   x2 = (N >= 2*step && N - 2*step >= x1)? x1 - 1 + step: N;
     844             : 
     845       29132 :   v = vecfactorsquarefreeu_coprime(x1, x2, P);
     846       29133 :   lv = lg(v);
     847    59133259 :   for (j = 1; j < lv; j++)
     848    59105257 :     if (gel(v,j))
     849             :     {
     850    17963288 :       ulong d = x1 - 1 + j; /* squarefree, coprime to 6 */
     851    17963288 :       GEN t = smallfact(d, gel(v,j), sq, V1), u;
     852    17927171 :       GEN tB = NULL, uB = NULL; /* = f(d) d^s */
     853             :       ulong a, b, c, e, q;
     854    17927171 :       if (!t || gequal0(t)) continue;
     855     3731392 :       if (fl) tB = vecinv(gmulsg(d, gconj(t)));
     856             :       /* warning: gives 1/conj(f(d)) d^(-1-conj(s)), equal to
     857             :          f(d) d^(-1-conj(s)) only if |f(d)|=1. */
     858             :       /* S += f(d) d^s * Z[q] */
     859     3808519 :       q = N / d;
     860     3808519 :       if (q == 1)
     861             :       {
     862     1453224 :         S = gadd(S, t); if (fl) SB = gadd(SB, tB);
     863     1457986 :         continue;
     864             :       }
     865     2355295 :       if (q <= sq) { u = gel(Z1, q); if (fl) uB = gel(Z2, q); }
     866             :       else
     867             :       {
     868       46904 :         a = usqrt(q); b = usqrt(q / 2); c = usqrt(q / 3); e = usqrt(q / 6);
     869       46977 :         u = gadd(gadd(gel(Q1,a), gel(Q2,b)), gadd(gel(Q3,c), gel(Q6,e)));
     870       46977 :         if (fl)
     871           0 :           uB = gadd(gadd(gel(QB,a), gel(Q2B,b)), gadd(gel(Q3B,c), gel(Q6B,e)));
     872             :       }
     873     2355368 :       S = gadd(S, vecmul(t, u)); if (fl) SB = gadd(SB, vecmul(tB, uB));
     874             :     }
     875       28002 :   return gerepilecopy(av2, mkvec2(S, SB));
     876             : }
     877             : 
     878             : GEN
     879       38128 : parsumprimeWfunboth_worker(GEN gk, GEN s, GEN W1, GEN W2, GEN f, GEN Nsqprec)
     880             : {
     881             :   pari_sp av2;
     882       38128 :   GEN S1, S2 = gen_0, logp;
     883             :   forprime_t T;
     884       38128 :   long k = itou(gk), N = Nsqprec[1], sq = Nsqprec[2], precp;
     885       38123 :   long STEP = Nsqprec[3], prec0 = Nsqprec[4], prec1 = Nsqprec[5], p;
     886       38123 :   long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     887       38123 :   long needlog = get_needlog(s), fl = !gequal0(W2), lv;
     888             : 
     889       38113 :   if (f && gequal0(f)) f = NULL;
     890       38090 :   if (!f) lv = 0;
     891             :   else
     892             :   {
     893       15966 :     GEN tmp = mycallvec(f, 1, prec1);
     894       15966 :     lv = typ(tmp) == t_VEC ? lg(tmp) - 1 : 0;
     895             :   }
     896       38090 :   precp = 0; logp = NULL;
     897       38090 :   if (lv) { S1 = const_vec(lv, real_0(prec1)); if (fl) S2 = const_vec(lv, real_0(prec1)); }
     898       22124 :   else { S1 = real_0(prec1); if (fl) S2 = real_0(prec1); }
     899       38087 :   u_forprime_init(&T, k * STEP + sq + 1, minss(N, (k + 1) * STEP + sq));
     900       38070 :   av2 = avma;
     901     5200823 :   while ((p = u_forprime_next(&T)))
     902             :   {
     903     5164931 :     GEN u = gen_0, ks = f ? mycallvec(f, p, prec1) : gen_1;
     904     5164925 :     long zks = !gequal0(ks);
     905     5164859 :     gp[2] = p;
     906     5164859 :     if (needlog)
     907             :     {
     908     5164988 :       if (!logp)
     909       30292 :         logp = logr_abs(utor(p, prec1));
     910             :       else
     911             :       { /* log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
     912     5134696 :         ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
     913     5134696 :         GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
     914     5131642 :         shiftr_inplace(z, 1); logp = addrr(logp, z);
     915             :       }
     916     5162351 :       if (zks)
     917     5162404 :         u = needlog == 1? powcx(gp, logp, s, prec0) : mpexp(gmul(s, logp));
     918             :     }
     919           0 :     else { if (zks) u = gpow(gp, s, prec0); }
     920     5164738 :     if (zks)
     921             :     {
     922     5164921 :       S1 = gadd(S1, vecmul(gel(W1, N / p), gmul(ks, u)));
     923     5162936 :       if (fl)
     924           0 :         S2 = gadd(S2, gdiv(vecmul(ks, gel(W2, N / p)), gmulsg(p, gconj(u))));
     925             :     }
     926     5162753 :     precp = p;
     927     5162753 :     if ((p & 0x1ff) == 1)
     928             :     {
     929       19285 :       if (!logp) gerepileall(av2, 2, &S1, &S2);
     930       19285 :       else gerepileall(av2, 3, &S1, &S2, &logp);
     931             :     }
     932             :   }
     933       37937 :   return gcopy(mkvec2(S1, S2));
     934             : }
     935             : 
     936             : static GEN
     937          14 : smalldirpowerssumfunvec_i(GEN f, ulong N, GEN s, long prec)
     938             : {
     939          14 :   GEN S = real_0(prec);
     940             :   ulong n;
     941          14 :   if (f && gequal0(f)) f = NULL;
     942          14 :   if (f)
     943             :   {
     944          14 :     GEN tmp = mycallvec(f, 1, prec);
     945          14 :     if (typ(tmp) == t_VEC) S = const_vec(lg(tmp) - 1, real_0(prec));
     946             :   }
     947          28 :   for (n = 1; n <= N; n++)
     948             :   {
     949          14 :     GEN ks = mycallvec(f, n, prec);
     950          14 :     if (!gequal0(ks)) S = gadd(S, gmul(ks, gpow(utoipos(n), s, prec)));
     951             :   }
     952          14 :   return S;
     953             : }
     954             : 
     955             : static GEN
     956          14 : smalldirpowerssumfunvec(GEN f, ulong N, GEN s, long fl, long prec)
     957             : {
     958          14 :   GEN tmp = smalldirpowerssumfunvec_i(f, N, s, prec);
     959          14 :   if (fl) return mkvec2(tmp, smalldirpowerssumfunvec_i(f, N, gsubgs(gneg(gconj(s)), 1), prec));
     960          14 :   else return mkvec2(tmp, tmp);
     961             : }
     962             : 
     963             : static GEN
     964        2163 : getscal(GEN V, long isscal)
     965             : {
     966        2163 :   if (isscal != 1 || typ(V) != t_VEC) return V;
     967         588 :   switch (lg(V))
     968             :   {
     969           0 :     case 2: return gel(V,1);
     970         588 :     case 3: return mkvec2(getscal(gel(V,1), 1), getscal(gel(V,2), 1));
     971             :   }
     972           0 :   pari_err_TYPE("getscal", V);
     973             :   return NULL; /*LCOV_EXCL_LINE*/
     974             : }
     975             : 
     976             : GEN
     977        1085 : pardirpowerssumfun(GEN f, ulong N, GEN s, long both, long prec)
     978             : {
     979        1085 :   pari_sp av = avma;
     980        1085 :   GEN P, V1, W1, Q1, V2 = gen_0, W2 = gen_0, QN = gen_0, c2, c2N;
     981        1085 :   GEN Q2, Q2N = gen_0, Q3, Q3N = gen_0, Q6, Q6N = gen_0;
     982        1085 :   GEN S1, S2, RES, Z1, Z2 = gen_0, logp;
     983        1085 :   long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     984             :   ulong a, b, c, e, q, n, sq, fl;
     985        1085 :   long prec0, prec1, needlog, lv = 0, isscal = 1;
     986        1085 :   GEN unvec = gen_1, zervec = gen_0, re0, re1, tmp2 = NULL;
     987             : 
     988        1085 :   if (f)
     989             :   {
     990         399 :     if (typ(f) == t_CLOSURE) isscal = -1; else isscal = 0;
     991         399 :     tmp2 = mycallvec(f, 2, prec + EXTRAPRECWORD);
     992         399 :     if (typ(tmp2) == t_VEC)
     993             :     {
     994         399 :       lv = lg(tmp2) - 1;
     995         399 :       unvec = const_vec(lv, gen_1);
     996         399 :       zervec = const_vec(lv, gen_0);
     997             :     }
     998             :   }
     999        1085 :   if (N <= 0)
    1000          28 :   { GEN v = mkvec2(gen_0, gen_0); return f ? const_vec(lv, v) : v; }
    1001        1057 :   fl = both && !gequal(real_i(s), gneg(ghalf));
    1002        1057 :   if (f && N < 49)
    1003          14 :     return gerepilecopy(av, getscal(smalldirpowerssumfunvec(f, N, s, fl, prec), isscal));
    1004        1043 :   if (!f && N < 10000UL)
    1005          70 :     return gerepilecopy(av, _smalldirpowerssum(N, s, fl, prec));
    1006         973 :   sq = usqrt(N);
    1007         973 :   V1 = cgetg(sq+1, t_VEC); W1 = cgetg(sq+1, t_VEC); Q1 = cgetg(sq+1, t_VEC);
    1008         973 :   if (fl)
    1009             :   {
    1010           0 :     V2 = cgetg(sq+1, t_VEC); W2 = cgetg(sq+1, t_VEC); QN = cgetg(sq+1, t_VEC);
    1011             :   }
    1012         973 :   prec1 = prec0 = prec + EXTRAPRECWORD;
    1013         973 :   s = gprec_w(s, prec0);
    1014         973 :   needlog = get_needlog(s);
    1015         973 :   if (needlog == 1) prec1 = powcx_prec(log2((double)N), s, prec);
    1016         973 :   gel(V1,1) = gel(W1,1) = gel(Q1,1) = unvec;
    1017         973 :   if (fl) { gel(V2,1) = gel(W2,1) = gel(QN,1) = unvec; }
    1018         973 :   c2 = gpow(gen_2, s, prec0); c2N = ginv(gmul2n(gconj(c2), 1));
    1019         973 :   re0 = real_0(prec0); re1 = real_1(prec0);
    1020         973 :   if (f) { c2 = gmul(c2, tmp2); c2N = gmul(c2N, tmp2); }
    1021         973 :   gel(V1,2) = c2; /* f(2) 2^s */
    1022         973 :   gel(W1,2) = gmul(re1, gadd(c2, unvec));
    1023         973 :   gel(Q1,2) = gmul(re1, gadd(vecsqr(c2), unvec));
    1024         973 :   if (fl)
    1025             :   {
    1026           0 :     gel(V2,2) = c2N; gel(W2,2) = gmul(re1, gadd(c2N, unvec));
    1027           0 :     gel(QN,2) = gmul(re1, gadd(vecsqr(c2N), unvec));
    1028             :   }
    1029         973 :   logp = NULL;
    1030      154623 :   for (n = 3; n <= sq; n++)
    1031             :   {
    1032      153650 :     GEN u1 = zervec, u2 = zervec, ks = f ? mycallvec(f, n, prec) : gen_1;
    1033      153650 :     long zks = !gequal0(ks);
    1034      153650 :     if (odd(n))
    1035             :     {
    1036       77287 :       gp[2] = n;
    1037       77287 :       if (needlog)
    1038             :       {
    1039       77287 :         if (!logp)
    1040         973 :           logp = logr_abs(utor(n, prec1));
    1041             :         else
    1042             :         { /* log n = log(n-2) + 2 atanh(1 / (n - 1)) */
    1043       76314 :           GEN z = atanhuu(1, n - 1, prec1);
    1044       76314 :           shiftr_inplace(z, 1); logp = addrr(logp, z);
    1045             :         }
    1046       77287 :         if (zks)
    1047       75712 :           u1 = needlog == 1? powcx(gp, logp, s, prec0) : mpexp(gmul(s, logp));
    1048             :       }
    1049           0 :       else { if (zks) u1 = gpow(gp, s, prec0); }
    1050       77287 :       if (zks)
    1051             :       {
    1052       75712 :         if(fl) u2 = gmul(ks, ginv(gmulsg(n, gconj(u1))));
    1053       75712 :         u1 = gmul(ks, u1); /* f(n) n^s */
    1054             :       }
    1055             :     }
    1056             :     else
    1057             :     {
    1058       76363 :       u1 = vecmul(c2, gel(V1, n >> 1));
    1059       76363 :       if (fl) u2 = vecmul(c2N, gel(V2, n >> 1));
    1060             :     }
    1061      153650 :     gel(V1,n) = u1; /* f(n) n^s */
    1062      153650 :     gel(W1,n) = gadd(gel(W1,n-1), gel(V1,n));       /* = sum_{i<=n} f(i)i^s */
    1063      153650 :     gel(Q1,n) = gadd(gel(Q1,n-1), vecsqr(gel(V1,n))); /* = sum_{i<=n} f(i^2)i^2s */
    1064      153650 :     if (fl)
    1065             :     {
    1066           0 :       gel(V2,n) = u2;
    1067           0 :       gel(W2,n) = gadd(gel(W2,n-1), gel(V2,n));
    1068           0 :       gel(QN,n) = gadd(gel(QN,n-1), vecsqr(gel(V2,n)));
    1069             :     }
    1070             :   }
    1071         973 :   Q2 = gmulvecsqlv(Q1, gel(V1,2));
    1072         973 :   Q3 = gmulvecsqlv(Q1, gel(V1,3));
    1073         973 :   Q6 = gmulvecsqlv(Q1, gel(V1,6));
    1074         973 :   if (fl)
    1075             :   {
    1076           0 :     Q2N = gmulvecsqlv(QN, gel(V2,2));
    1077           0 :     Q3N = gmulvecsqlv(QN, gel(V2,3));
    1078           0 :     Q6N = gmulvecsqlv(QN, gel(V2,6));
    1079             :   }
    1080         973 :   S1 = S2 = typ(zervec) == t_VEC? const_vec(lv, re0): re0;
    1081         973 :   RES = mkvec2(S1, S2);
    1082             :   {
    1083         973 :     GEN fspec = f ? f : gen_0;
    1084         973 :     long m = mt_nbthreads();
    1085         973 :     long STEP = maxss(N / (m * m), 1);
    1086         973 :     GEN VS = mkvecsmalln(5, N, sq, STEP, prec0, prec1);
    1087         973 :     GEN FUN = snm_closure(is_entry("_parsumprimeWfunboth_worker"),
    1088             :                           mkvec5(s, W1, W2, fspec, VS));
    1089         973 :     RES = gadd(RES, parsum(gen_0, utoipos((N - 1) / STEP), FUN));
    1090             :   }
    1091         973 :   P = mkvecsmall2(2, 3);
    1092         973 :   Z1 = cgetg(sq+1, t_VEC);
    1093             :   /* a,b,c,e = sqrt(q), sqrt(q/2), sqrt(q/3), sqrt(q/6)
    1094             :    * Z[q] = Q[a] + 2^s Q[b] + 3^s Q[c] + 6^s Q[e], with Q[0] = 0 */
    1095         973 :   gel(Z1, 1) = unvec;
    1096         973 :   gel(Z1, 2) = gel(W1, 2);
    1097         973 :   gel(Z1, 3) = gel(W1, 3);
    1098         973 :   gel(Z1, 4) = gel(Z1, 5) = gel(W1, 4);
    1099         973 :   gel(Z1, 6) = gel(Z1, 7) = gadd(gel(W1, 4), gel(V1, 6));
    1100         973 :   if (fl)
    1101             :   {
    1102           0 :     Z2 = cgetg(sq+1, t_VEC);
    1103           0 :     gel(Z2, 1) = unvec;
    1104           0 :     gel(Z2, 2) = gel(W2, 2);
    1105           0 :     gel(Z2, 3) = gel(W2, 3);
    1106           0 :     gel(Z2, 4) = gel(Z2, 5) = gel(W2, 4);
    1107           0 :     gel(Z2, 6) = gel(Z2, 7) = gadd(gel(W2, 4), gel(V2, 6));
    1108             :   }
    1109         973 :   a = 2; b = c = e = 1;
    1110      149758 :   for (q = 8; q <= sq; q++)
    1111             :   { /* Gray code: at most one of a,b,c,d differs (by 1) from previous value */
    1112      148785 :     GEN z1 = gel(Z1, q - 1), z2 = gen_0;
    1113             :     ulong na, nb, nc, ne, na2, nb2, nc2, ne2;
    1114      148785 :     if (fl) z2 = gel(Z2, q - 1);
    1115      148785 :     if ((na = usqrt(q)) != a)
    1116        8911 :     { a = na; na2 = na * na; z1 = gadd(z1, gel(V1, na2)); if (fl) z2 = gadd(z2, gel(V2, na2)); }
    1117      139874 :     else if ((nb = usqrt(q / 2)) != b)
    1118        5999 :     { b = nb; nb2 = 2 * nb * nb; z1 = gadd(z1, gel(V1, nb2)); if (fl) z2 = gadd(z2, gel(V2, nb2)); }
    1119      133875 :     else if ((nc = usqrt(q / 3)) != c)
    1120        5257 :     { c = nc; nc2 = 3 * nc * nc; z1 = gadd(z1, gel(V1, nc2)); if (fl) z2 = gadd(z2, gel(V2, nc2)); }
    1121      128618 :     else if ((ne = usqrt(q / 6)) != e)
    1122        3010 :     { e = ne; ne2 = 6 * ne * ne; z1 = gadd(z1, gel(V1, ne2)); if (fl) z2 = gadd(z2, gel(V2, ne2)); }
    1123      148785 :     gel(Z1,q) = z1; if (fl) gel(Z2,q) = z2;
    1124             :   }
    1125             :   {
    1126         973 :     GEN vQ1 = mkvec4(Q1, Q2, Q3, Q6);
    1127         973 :     GEN vQ2 = fl? mkvec4(QN, Q2N, Q3N, Q6N): gen_0;
    1128         973 :     GEN worker = snm_closure(is_entry("_parsqfboth_worker"),
    1129             :                    mkvec5(mkvec2(Z1, Z2), mkvec2(vQ1, vQ2), mkvec2(V1, V2),
    1130             :                    P, mkvecsmall2(N, sq)));
    1131         973 :     RES = gadd(RES, parsum(gen_0, utoipos(maxss((N-1) / step - 1, 0)), worker));
    1132             :   }
    1133         973 :   if (!fl) gel(RES, 2) = gel(RES, 1);
    1134         973 :   return gerepilecopy(av, getscal(RES, isscal));
    1135             : }
    1136             : 
    1137             : GEN
    1138           0 : pardirpowerssum0(GEN N, GEN s, GEN f, long both, long prec)
    1139             : {
    1140           0 :   pari_sp av = avma;
    1141             :   GEN R;
    1142           0 :   if (typ(N) != t_INT) pari_err_TYPE("pardirpowerssum", N);
    1143           0 :   R = pardirpowerssumfun(f, itou(N), s, both, prec);
    1144           0 :   return both ? R : gerepilecopy(av, gel(R, 1));
    1145             : }
    1146             : 
    1147             : GEN
    1148           0 : pardirpowerssum(ulong N, GEN s, long prec)
    1149             : {
    1150           0 :   pari_sp av = avma;
    1151           0 :   return gerepilecopy(av, gel(pardirpowerssumfun(NULL, N, s, 0, prec), 1));
    1152             : }

Generated by: LCOV version 1.16