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 - bern.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 340 358 95.0 %
Date: 2024-03-29 08:06:26 Functions: 35 36 97.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2018  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             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_bern
      19             : 
      20             : /********************************************************************/
      21             : /**                                                                **/
      22             : /**                     BERNOULLI NUMBERS B_2k                     **/
      23             : /**                                                                **/
      24             : /********************************************************************/
      25             : 
      26             : /* D = divisorsu(n). Return a/b = \sum_{p-1 | 2n: p prime} 1/p
      27             :  * B_2k + a/b in Z [Clausen-von Staudt] */
      28             : static GEN
      29       81534 : fracB2k(GEN D)
      30             : {
      31       81534 :   GEN a = utoipos(5), b = utoipos(6); /* 1/2 + 1/3 */
      32       81552 :   long i, l = lg(D);
      33      618998 :   for (i = 2; i < l; i++) /* skip 1 */
      34             :   {
      35      537502 :     ulong p = 2*D[i] + 1; /* a/b += 1/p */
      36      537502 :     if (uisprime(p)) { a = addii(muliu(a,p), b); b = muliu(b,p); }
      37             :   }
      38       81496 :   return mkfrac(a,b);
      39             : }
      40             : /* precision needed to compute B_k for all k <= N */
      41             : long
      42        1870 : bernbitprec(long N)
      43             : { /* 1.612086 ~ log(8Pi) / 2 */
      44        1870 :   const double log2PI = 1.83787706641;
      45        1870 :   double logN = log((double)N);
      46        1870 :   double t = (N + 4) * logN - N*(1 + log2PI) + 1.612086;
      47        1870 :   return (long)ceil(t / M_LN2) + 10;
      48             : }
      49             : static long
      50          21 : bernprec(long N) { return nbits2prec(bernbitprec(N)); }
      51             : /* \sum_{k > M} k^(-n) <= M^(1-n) / (n-1) < 2^-b */
      52             : static long
      53        1373 : zetamaxpow(long n)
      54             : {
      55        1373 :   long M = (long)ceil(n / (2 * M_PI * M_E));
      56        1373 :   return M | 1; /* make it odd */
      57             : }
      58             : /* v * zeta(k) using r precomputed odd powers */
      59             : static GEN
      60       81510 : bern_zeta(GEN v, long k, GEN pow, long r, long p)
      61             : {
      62       81510 :   GEN z, s = gel(pow, r);
      63             :   long j;
      64    10360793 :   for (j = r - 2; j >= 3; j -= 2) s = addii(s, gel(pow,j));
      65       81418 :   z = s = itor(s, nbits2prec(p));
      66       81508 :   shiftr_inplace(s, -p); /* zeta(k)(1 - 2^(-k)) - 1*/
      67       81508 :   s = addrs(s, 1); shiftr_inplace(s, -k);
      68             :   /* divide by 1 - 2^(-k): s + s/2^k + s/2^(2k) + ... */
      69      361153 :   for (; k < p; k <<= 1) s = addrr(s, shiftr(s, -k));
      70       81518 :   return addrr(v, mulrr(v, addrr(z, s)));
      71             : }
      72             : /* z * j^2 */
      73             : static GEN
      74    50915811 : muliu2(GEN z, ulong j)
      75    50915811 : { return (j | HIGHMASK)? mulii(z, sqru(j)): muliu(z, j*j); }
      76             : /* 1 <= m <= n, set y[1] = B_{2m}, ... y[n-m+1] = B_{2n} in Q */
      77             : static void
      78         295 : bernset(GEN *y, long m, long n)
      79             : {
      80         295 :   long i, j, k, p, prec, r, N = n << 1; /* up to B_N */
      81             :   GEN u, b, v, t;
      82         295 :   p = bernbitprec(N); prec = nbits2prec(p);
      83         295 :   u = sqrr(Pi2n(1, prec)); /* (2Pi)^2 */
      84         295 :   v = divrr(mpfactr(N, prec), powru(u, n)); shiftr_inplace(v,1);
      85         295 :   r = zetamaxpow(N);
      86         295 :   t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
      87        5558 :   for (j = 3; j <= r; j += 2)
      88             :   {
      89        5263 :     GEN z = cgeti(nbits2lg(p));
      90        5263 :     pari_sp av2 = avma;
      91        5263 :     affii(divii(b, powuu(j, N)), z);
      92        5263 :     gel(t,j) = z; set_avma(av2);
      93             :   }
      94         295 :   y += n - m;
      95         295 :   for (i = n, k = N;; i--)
      96       81208 :   { /* set B_n, k = 2i */
      97       81503 :     pari_sp av2 = avma;
      98       81503 :     GEN z = fracB2k(divisorsu(i)), B = bern_zeta(v, k, t, r, p);
      99             :     long j;
     100             :     /* B = v * zeta(k), v = 2*k! / (2Pi)^k */
     101       81509 :     if (!odd(i)) setsigne(B, -1); /* B ~ B_n */
     102       81509 :     B = roundr(addrr(B, fractor(z,LOWDEFAULTPREC))); /* B - z = B_n */
     103       81510 :     *y-- = gclone(gsub(B, z));
     104       81518 :     if (i == m) break;
     105       81223 :     affrr(divrunextu(mulrr(v,u), k-1), v);
     106    10439670 :     for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
     107       81202 :     set_avma(av2); k -= 2;
     108       81207 :     if (((N - k) & 0x7f) == 0x7e)
     109             :     { /* reduce precision if possible */
     110        1078 :       long p2 = p, prec2 = prec;
     111        1078 :       p = bernbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
     112        1078 :       setprec(v, prec); r = zetamaxpow(k);
     113      158372 :       for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
     114        1078 :       set_avma(av2);
     115             :     }
     116             :   }
     117         295 : }
     118             : /* need B[2..2*nb] as t_INT or t_FRAC */
     119             : void
     120      237804 : constbern(long nb)
     121             : {
     122      237804 :   const pari_sp av = avma;
     123             :   long i, l;
     124             :   GEN B;
     125             :   pari_timer T;
     126             : 
     127      237804 :   l = bernzone? lg(bernzone): 0;
     128      237804 :   if (l > nb) return;
     129             : 
     130         295 :   nb = maxss(nb, l + 127);
     131         295 :   B = cgetg_block(nb+1, t_VEC);
     132         295 :   if (bernzone)
     133        7040 :   { for (i = 1; i < l; i++) gel(B,i) = gel(bernzone,i); }
     134             :   else
     135             :   {
     136         247 :     gel(B,1) = gclone(mkfracss(1,6));
     137         247 :     gel(B,2) = gclone(mkfracss(-1,30));
     138         247 :     gel(B,3) = gclone(mkfracss(1,42));
     139         247 :     gel(B,4) = gclone(mkfracss(-1,30));
     140         247 :     gel(B,5) = gclone(mkfracss(5,66));
     141         247 :     gel(B,6) = gclone(mkfracss(-691,2730));
     142         247 :     gel(B,7) = gclone(mkfracss(7,6));
     143         247 :     gel(B,8) = gclone(mkfracss(-3617,510));
     144         247 :     gel(B,9) = gclone(mkfracss(43867,798));
     145         247 :     gel(B,10)= gclone(mkfracss(-174611,330));
     146         247 :     gel(B,11)= gclone(mkfracss(854513,138));
     147         247 :     gel(B,12)= gclone(mkfracss(-236364091,2730));
     148         247 :     gel(B,13)= gclone(mkfracss(8553103,6)); /* B_26 */
     149         247 :     l = 14;
     150             :   }
     151         295 :   set_avma(av);
     152         295 :   if (DEBUGLEVEL) {
     153           0 :     err_printf("caching Bernoulli numbers 2*%ld to 2*%ld\n", l, nb);
     154           0 :     timer_start(&T);
     155             :   }
     156         295 :   bernset((GEN*)B + l, l, nb);
     157         295 :   if (DEBUGLEVEL) timer_printf(&T, "Bernoulli");
     158         295 :   swap(B, bernzone); guncloneNULL(B);
     159         295 :   set_avma(av);
     160             : #if 0
     161             :   if (nb > 200000)
     162             : #endif
     163             :   {
     164         295 :     const ulong p = 4294967291UL;
     165         295 :     long n = 2 * nb + 2;
     166         295 :     GEN t = const_vecsmall(n+1, 1);
     167         295 :     t[1] = evalvarn(0); t[2] = 0;
     168         295 :     t = Flx_shift(Flx_invLaplace(t, p), -1); /* t = (exp(x)-1)/x */
     169         295 :     t = Flx_Laplace(Flxn_inv(t, n, p), p);
     170       92022 :     for (i = 1; i <= nb; i++)
     171       91727 :       if (Rg_to_Fl(bernfrac(2*i), p) != uel(t,2*i+2))
     172             :       {
     173           0 :         gunclone(bernzone); bernzone = NULL;
     174           0 :         pari_err_BUG(stack_sprintf("B_{2*%ld}", i));
     175             :       }
     176         295 :     set_avma(av);
     177             :   }
     178             : }
     179             : /* Obsolete, kept for backward compatibility */
     180             : void
     181           0 : mpbern(long n, long prec) { (void)prec; constbern(n); }
     182             : 
     183             : /* assume n even > 0, if iz != NULL, assume iz = 1/zeta(n) */
     184             : static GEN
     185          21 : bernreal_using_zeta(long n, long prec)
     186             : {
     187          21 :   GEN pi2 = Pi2n(1, prec+EXTRAPREC64);
     188          21 :   GEN iz = inv_szeta_euler(n, prec);
     189          21 :   GEN z = divrr(mpfactr(n, prec), mulrr(powru(pi2, n), iz));
     190          21 :   shiftr_inplace(z, 1); /* 2 * n! * zeta(n) / (2Pi)^n */
     191          21 :   if ((n & 3) == 0) setsigne(z, -1);
     192          21 :   return z;
     193             : }
     194             : /* assume n even > 0, B = NULL or good approximation to B_n */
     195             : static GEN
     196          14 : bernfrac_i(long n, GEN B)
     197             : {
     198          14 :   GEN z = fracB2k(divisorsu(n >> 1));
     199          14 :   if (!B) B = bernreal_using_zeta(n, bernprec(n));
     200          14 :   B = roundr( gadd(B, fractor(z,LOWDEFAULTPREC)) );
     201          14 :   return gsub(B, z);
     202             : }
     203             : GEN
     204     9233039 : bernfrac(long n)
     205             : {
     206             :   pari_sp av;
     207             :   long k;
     208     9233039 :   if (n <= 1)
     209             :   {
     210        4431 :     if (n < 0) pari_err_DOMAIN("bernfrac", "index", "<", gen_0, stoi(n));
     211        4424 :     return n? mkfrac(gen_m1,gen_2): gen_1;
     212             :   }
     213     9228608 :   if (odd(n)) return gen_0;
     214     9224684 :   k = n >> 1;
     215     9224684 :   if (!bernzone) constbern(0);
     216     9224684 :   if (bernzone && k < lg(bernzone)) return gel(bernzone, k);
     217          16 :   av = avma;
     218          16 :   return gerepileupto(av, bernfrac_i(n, NULL));
     219             : }
     220             : GEN
     221          77 : bernvec(long n)
     222             : {
     223             :   long i, l;
     224             :   GEN y;
     225          77 :   if (n < 0) return cgetg(1, t_VEC);
     226          77 :   constbern(n);
     227          77 :   l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
     228       53536 :   for (i = 2; i < l; i++) gel(y,i) = gel(bernzone,i-1);
     229          77 :   return y;
     230             : }
     231             : 
     232             : /* x := pol_x(v); B_k(x) = \sum_{i=0}^k binomial(k, i) B_i x^{k-i} */
     233             : static GEN
     234        2198 : bernpol_i(long k, long v)
     235             : {
     236             :   GEN B, C;
     237             :   long i;
     238        2198 :   if (v < 0) v = 0;
     239        2198 :   constbern(k >> 1); /* cache B_2, ..., B_2[k/2] */
     240        2198 :   C = vecbinomial(k);
     241        2198 :   B = cgetg(k + 3, t_POL);
     242       14966 :   for (i = 0; i <= k; ++i) gel(B, k-i+2) = gmul(gel(C,i+1), bernfrac(i));
     243        2198 :   B[1] = evalsigne(1) | evalvarn(v);
     244        2198 :   return B;
     245             : }
     246             : GEN
     247        2086 : bernpol(long k, long v)
     248             : {
     249        2086 :   pari_sp av = avma;
     250        2086 :   if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
     251        2079 :   return gerepileupto(av, bernpol_i(k, v));
     252             : }
     253             : GEN
     254          98 : bernpol_eval(long k, GEN x)
     255             : {
     256          98 :   pari_sp av = avma;
     257             :   GEN B;
     258          98 :   if (!x) return bernpol(k, 0);
     259          49 :   if (gequalX(x)) return bernpol(k, varn(x));
     260          49 :   if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
     261          42 :   B = bernpol_i(k, fetch_var_higher());
     262          42 :   return gerepileupto(av, poleval(B, x));
     263             : }
     264             : 
     265             : /* x := pol_x(v); return 1^e + ... + x^e = x^e + (B_{e+1}(x) - B_{e+1})/(e+1) */
     266             : static GEN
     267          49 : faulhaber(long e, long v)
     268             : {
     269             :   GEN B;
     270          49 :   if (e == 0) return pol_x(v);
     271          35 :   B = RgX_integ(bernpol_i(e, v)); /* (B_{e+1}(x) - B_{e+1}) / (e+1) */
     272          35 :   gel(B,e+2) = gaddgs(gel(B,e+2), 1); /* add x^e, in place */
     273          35 :   return B;
     274             : }
     275             : /* sum_v T(v), T a polynomial expression in v */
     276             : GEN
     277          49 : sumformal(GEN T, long v)
     278             : {
     279          49 :   pari_sp av = avma, av2;
     280             :   long i, t, d;
     281             :   GEN R;
     282             : 
     283          49 :   T = simplify_shallow(T);
     284          49 :   t = typ(T);
     285          49 :   if (is_scalar_t(t))
     286          14 :     return gerepileupto(av, monomialcopy(T, 1, v < 0? 0: v));
     287          35 :   if (t != t_POL) pari_err_TYPE("sumformal [not a t_POL]", T);
     288          28 :   if (v < 0) v = varn(T);
     289          28 :   av2 = avma;
     290          28 :   R = gen_0;
     291          28 :   d = poldegree(T,v);
     292          91 :   for (i = d; i >= 0; i--)
     293             :   {
     294          63 :     GEN c = polcoef_i(T, i, v);
     295          63 :     if (gequal0(c)) continue;
     296          42 :     R = gadd(R, gmul(c, faulhaber(i, v)));
     297          42 :     if (gc_needed(av2,3))
     298             :     {
     299           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"sumformal, i = %ld/%ld", i,d);
     300           0 :       R = gerepileupto(av2, R);
     301             :     }
     302             :   }
     303          28 :   return gerepileupto(av, R);
     304             : }
     305             : 
     306             : /* 1/zeta(n) using Euler product. Assume n > 0. */
     307             : GEN
     308        1194 : inv_szeta_euler(long n, long prec)
     309             : {
     310        1194 :   long bit = prec2nbits(prec);
     311             :   GEN z, res;
     312             :   pari_sp av, av2;
     313             :   double A, D, lba;
     314             :   ulong p, lim;
     315             :   forprime_t S;
     316             : 
     317        1194 :   if (n > bit) return real_1(prec);
     318             : 
     319        1187 :   lba = prec2nbits_mul(prec, M_LN2);
     320        1187 :   D = exp((lba - log((double)(n-1))) / (n-1));
     321        1187 :   lim = 1 + (ulong)ceil(D);
     322        1187 :   if (lim < 3) return subir(gen_1,real2n(-n,prec));
     323        1187 :   res = cgetr(prec); av = avma; incrprec(prec);
     324             : 
     325        1187 :   (void)u_forprime_init(&S, 3, lim);
     326        1187 :   av2 = avma; A = n / M_LN2; z = subir(gen_1, real2n(-n, prec));
     327       86205 :   while ((p = u_forprime_next(&S)))
     328             :   {
     329       85018 :     long l = bit - (long)floor(A * log((double)p));
     330             :     GEN h;
     331             : 
     332       85018 :     if (l < BITS_IN_LONG) l = BITS_IN_LONG;
     333       85018 :     l = minss(prec, nbits2prec(l));
     334       85018 :     h = divrr(z, rpowuu(p, (ulong)n, l));
     335       85018 :     z = subrr(z, h);
     336       85018 :     if (gc_needed(av,1))
     337             :     {
     338           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"inv_szeta_euler, p = %lu/%lu", p,lim);
     339           0 :       z = gerepileuptoleaf(av2, z);
     340             :     }
     341             :   }
     342        1187 :   affrr(z, res); set_avma(av); return res;
     343             : }
     344             : 
     345             : /* Return B_n */
     346             : GEN
     347        7917 : bernreal(long n, long prec)
     348             : {
     349             :   pari_sp av;
     350             :   GEN B;
     351             :   long p, k;
     352        7917 :   if (n < 0) pari_err_DOMAIN("bernreal", "index", "<", gen_0, stoi(n));
     353        7910 :   if (n == 0) return real_1(prec);
     354        7910 :   if (n == 1) return real_m2n(-1,prec); /*-1/2*/
     355        7910 :   if (odd(n)) return real_0(prec);
     356             : 
     357        7910 :   k = n >> 1;
     358        7910 :   if (!bernzone) constbern(0);
     359        7910 :   if (k < lg(bernzone)) return fractor(gel(bernzone,k), prec);
     360           7 :   p = bernprec(n); av = avma;
     361           7 :   B = bernreal_using_zeta(n, minss(p, prec));
     362           7 :   if (p < prec) B = fractor(bernfrac_i(n, B), prec);
     363           7 :   return gerepileuptoleaf(av, B);
     364             : }
     365             : 
     366             : GEN
     367          49 : eulerpol(long k, long v)
     368             : {
     369          49 :   pari_sp av = avma;
     370             :   GEN B, E;
     371          49 :   if (k < 0) pari_err_DOMAIN("eulerpol", "index", "<", gen_0, stoi(k));
     372          42 :   k++; B = bernpol_i(k, v);
     373          42 :   E = RgX_Rg_mul(RgX_sub(B, RgX_rescale(B, gen_2)), uutoQ(2,k));
     374          42 :   return gerepileupto(av, E);
     375             : }
     376             : 
     377             : /*******************************************************************/
     378             : /**                      HARMONIC NUMBERS                         **/
     379             : /*******************************************************************/
     380             : /* 1/a + ... + 1/(b-1); a < b <= 2^(BIL-1) */
     381             : static GEN
     382        4676 : hrec(ulong a, ulong b)
     383             : {
     384             :   ulong m;
     385        4676 :   switch(b - a)
     386             :   {
     387        4501 :     case 1: retmkfrac(gen_1, utoipos(a));
     388         133 :     case 2: if (a < 65536) retmkfrac(utoipos(2*a + 1), utoipos(a * a + a));
     389           0 :       retmkfrac(utoipos(2*a + 1), muluu(a, a+1));
     390             :   }
     391          42 :   m = (a + b) >> 1;
     392          42 :   return gadd(hrec(a, m), hrec(m, b));
     393             : }
     394             : /* exact Harmonic number H_n, n < 2^(BIL-1).
     395             :  * Could use H_n = sum_k 2^(-k) H^odd_{n \ 2^k} */
     396             : GEN
     397        4592 : harmonic(ulong n)
     398             : {
     399        4592 :   pari_sp av = avma;
     400        4592 :   return n? gerepileupto(av, hrec(1, n+1)): gen_0;
     401             : }
     402             : 
     403             : /* 1/a^k + ... + 1/(b-1)^k; a < b */
     404             : static GEN
     405          77 : hreck(ulong a, ulong b, ulong k)
     406             : {
     407             :   ulong m;
     408          77 :   switch(b - a)
     409             :   {
     410             :     GEN x, y;
     411          14 :     case 1: retmkfrac(gen_1, powuu(a, k));
     412          28 :     case 2:
     413          28 :       x = powuu(a, k); y = powuu(a + 1, k);
     414          28 :       retmkfrac(addii(x, y), mulii(x, y));
     415             :   }
     416          35 :   m = (a + b) >> 1;
     417          35 :   return gadd(hreck(a, m, k), hreck(m, b, k));
     418             : }
     419             : GEN
     420          56 : harmonic0(ulong n, GEN k)
     421             : {
     422          56 :   pari_sp av = avma;
     423             :   ulong r;
     424          56 :   if (!n) return gen_0;
     425          49 :   if (n & HIGHBIT) pari_err_OVERFLOW("harmonic");
     426          49 :   if (!k) return harmonic(n);
     427          21 :   if (typ(k) != t_INT) pari_err_TYPE("harmonic", k);
     428          14 :   if (signe(k) < 0)
     429             :   {
     430           7 :     GEN H = poleval(faulhaber(-itos(k), 0), utoipos(n));
     431           7 :     return gerepileuptoint(av, H);
     432             :   }
     433           7 :   r = itou(k);
     434           7 :   if (!r) return utoipos(n);
     435           7 :   if (r == 1) return harmonic(n);
     436           7 :   return gerepileupto(av, hreck(1, n+1, r));
     437             : }
     438             : 
     439             : 
     440             : /**************************************************************/
     441             : /*                      Euler numbers                         */
     442             : /**************************************************************/
     443             : 
     444             : /* precision needed to compute E_k for all k <= N */
     445             : static long
     446         798 : eulerbitprec(long N)
     447             : { /* 1.1605 ~ log(32/Pi) / 2 */
     448         798 :   const double logPIS2 = 0.4515827;
     449         798 :   double t = (N + 1) * log((double)N) - N*(1 + logPIS2) + 1.1605;
     450         798 :   return (long)ceil(t / M_LN2) + 10;
     451             : }
     452             : static long
     453          14 : eulerprec(long N) { return nbits2prec(eulerbitprec(N)); }
     454             : 
     455             : /* sum_{k > M, k odd} (-1)^((k-1)/2)k^(-n) < M^(-n) < 2^-b */
     456             : static long
     457         784 : lfun4maxpow(long n)
     458             : {
     459         784 :   long M = (long)ceil(2 * n / (M_E * M_PI));
     460         784 :   return M | 1; /* make it odd */
     461             : }
     462             : 
     463             : /* lfun4(k) using r precomputed odd powers */
     464             : static GEN
     465       49791 : euler_lfun4(GEN v, GEN pow, long r, long p)
     466             : {
     467       49791 :   GEN s = ((r & 3L) == 1)? gel(pow, r): negi(gel(pow, r));
     468             :   long j;
     469    40556894 :   for (j = r - 2; j >= 3; j -= 2)
     470    40507103 :     s = ((j & 3L) == 1)? addii(s, gel(pow,j)): subii(s, gel(pow,j));
     471       49791 :   s = mulri(v, s); shiftr_inplace(s, -p);
     472       49791 :   return addrr(v, s);
     473             : }
     474             : 
     475             : /* 1 <= m <= n, set y[1] = E_{2m}, ... y[n-m+1] = E_{2n} in Z */
     476             : static void
     477          21 : eulerset(GEN *y, long m, long n)
     478             : {
     479          21 :   long i, j, k, p, prec, r, N = n << 1, N1 = N + 1; /* up to E_N */
     480             :   GEN b, u, v, t;
     481          21 :   p = eulerbitprec(N); prec = nbits2prec(p);
     482          21 :   u = sqrr(Pi2n(-1, prec)); /* (Pi/2)^2 */
     483          21 :   v = divrr(mpfactr(N, prec), mulrr(powru(u, n), Pi2n(-2,prec)));
     484          21 :   r = lfun4maxpow(N1);
     485          21 :   t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
     486       11921 :   for (j = 3; j <= r; j += 2)
     487             :   {
     488       11900 :     GEN z = cgeti(nbits2lg(p));
     489       11900 :     pari_sp av2 = avma;
     490       11900 :     affii(divii(b, powuu(j, N+1)), z);
     491       11900 :     gel(t,j) = z; set_avma(av2);
     492             :   }
     493          21 :   y += n - m;
     494          21 :   for (i = n, k = N1;; i--)
     495       49770 :   { /* set E_n, k = 2i + 1 */
     496       49791 :     pari_sp av2 = avma;
     497       49791 :     GEN E = euler_lfun4(v, t, r, p);
     498             :     long j;
     499             :     /* E = v * lfun4(k), v = (4/Pi)*k! / (Pi/2)^k */
     500       49791 :     E = roundr(E); if (odd(i)) setsigne(E, -1); /* E ~ E_n */
     501       49791 :     *y-- = gclone(E);
     502       49791 :     if (i == m) break;
     503       49770 :     affrr(divrunextu(mulrr(v,u), k-2), v);
     504    40606202 :     for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
     505       49770 :     set_avma(av2); k -= 2;
     506       49770 :     if (((N1 - k) & 0x7f) == 0x7e)
     507             :     { /* reduce precision if possible */
     508         763 :       long p2 = p, prec2 = prec;
     509         763 :       p = eulerbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
     510         763 :       setprec(v, prec); r = lfun4maxpow(k);
     511      622923 :       for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
     512         763 :       set_avma(av2);
     513             :     }
     514             :   }
     515          21 : }
     516             : 
     517             : /* need E[2..2*nb] as t_INT */
     518             : static void
     519          28 : constreuler(long nb)
     520             : {
     521          28 :   const pari_sp av = avma;
     522             :   long i, l;
     523             :   GEN E;
     524             :   pari_timer T;
     525             : 
     526          28 :   l = eulerzone? lg(eulerzone): 0;
     527          28 :   if (l > nb) return;
     528             : 
     529          21 :   nb = maxss(nb, l + 127);
     530          21 :   E = cgetg_block(nb+1, t_VEC);
     531          21 :   if (eulerzone)
     532         896 :   { for (i = 1; i < l; i++) gel(E,i) = gel(eulerzone,i); }
     533             :   else
     534             :   {
     535          14 :     gel(E,1) = gclone(stoi(-1));
     536          14 :     gel(E,2) = gclone(stoi(5));
     537          14 :     gel(E,3) = gclone(stoi(-61));
     538          14 :     gel(E,4) = gclone(stoi(1385));
     539          14 :     gel(E,5) = gclone(stoi(-50521));
     540          14 :     gel(E,6) = gclone(stoi(2702765));
     541          14 :     gel(E,7) = gclone(stoi(-199360981));
     542          14 :     l = 8;
     543             :   }
     544          21 :   set_avma(av);
     545          21 :   if (DEBUGLEVEL) {
     546           0 :     err_printf("caching Euler numbers 2*%ld to 2*%ld\n", l, nb);
     547           0 :     timer_start(&T);
     548             :   }
     549          21 :   eulerset((GEN*)E + l, l, nb);
     550          21 :   if (DEBUGLEVEL) timer_printf(&T, "Euler");
     551          21 :   swap(E, eulerzone); guncloneNULL(E);
     552          21 :   set_avma(av);
     553             : }
     554             : 
     555             : /* 1/lfun(-4,n) using Euler product. Assume n > 0. */
     556             : static GEN
     557          14 : inv_lfun4(long n, long prec)
     558             : {
     559          14 :   long bit = prec2nbits(prec);
     560             :   GEN z, res;
     561             :   pari_sp av, av2;
     562             :   double A;
     563             :   ulong p, lim;
     564             :   forprime_t S;
     565             : 
     566          14 :   if (n > bit) return real_1(prec);
     567             : 
     568          14 :   lim = 1 + (ulong)ceil(exp2((double)bit / n));
     569          14 :   res = cgetr(prec); av = avma; incrprec(prec);
     570             : 
     571          14 :   (void)u_forprime_init(&S, 3, lim);
     572          14 :   av2 = avma; A = n / M_LN2; z = real_1(prec);
     573         369 :   while ((p = u_forprime_next(&S)))
     574             :   {
     575         355 :     long l = bit - (long)floor(A * log((double)p));
     576             :     GEN h;
     577             : 
     578         355 :     if (l < BITS_IN_LONG) l = BITS_IN_LONG;
     579         355 :     l = minss(prec, nbits2prec(l));
     580         355 :     h = rpowuu(p, (ulong)n, l); if ((p & 3UL) == 1) setsigne(h, -1);
     581         355 :     z = addrr(z, divrr(z, h)); /* z *= 1 - chi_{-4}(p) / p^n */
     582         355 :     if (gc_needed(av,1))
     583             :     {
     584           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"inv_lfun4, p = %lu/%lu", p,lim);
     585           0 :       z = gerepileuptoleaf(av2, z);
     586             :     }
     587             :   }
     588          14 :   affrr(z, res); set_avma(av); return res;
     589             : }
     590             : /* assume n even > 0, E_n = (-1)^(n/2) (4/Pi) n! lfun4(n+1) / (Pi/2)^n */
     591             : static GEN
     592          14 : eulerreal_using_lfun4(long n, long prec)
     593             : {
     594          14 :   GEN pisur2 = Pi2n(-1, prec+EXTRAPREC64);
     595          14 :   GEN iz = inv_lfun4(n+1, prec);
     596          14 :   GEN z = divrr(mpfactr(n, prec), mulrr(powru(pisur2, n+1), iz));
     597          14 :   if ((n & 3L) == 2) setsigne(z, -1);
     598          14 :   shiftr_inplace(z, 1); return z;
     599             : }
     600             : /* Euler numbers: 1, 0, -1, 0, 5, 0, -61,... */
     601             : GEN
     602         217 : eulerfrac(long n)
     603             : {
     604             :   pari_sp av;
     605             :   long k;
     606             :   GEN E;
     607         217 :   if (n <= 0)
     608             :   {
     609          21 :     if (n < 0) pari_err_DOMAIN("eulerfrac", "index", "<", gen_0, stoi(n));
     610          14 :     return gen_1;
     611             :   }
     612         196 :   if (odd(n)) return gen_0;
     613         189 :   k = n >> 1;
     614         189 :   if (!eulerzone) constreuler(0);
     615         189 :   if (eulerzone && k < lg(eulerzone)) return gel(eulerzone, k);
     616          14 :   av = avma; E = eulerreal_using_lfun4(n, eulerprec(n));
     617          14 :   return gerepileuptoleaf(av, roundr(E));
     618             : }
     619             : GEN
     620          14 : eulervec(long n)
     621             : {
     622             :   long i, l;
     623             :   GEN y;
     624          14 :   if (n < 0) return cgetg(1, t_VEC);
     625          14 :   constreuler(n);
     626          14 :   l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
     627       49224 :   for (i = 2; i < l; i++) gel(y,i) = gel(eulerzone,i-1);
     628          14 :   return y;
     629             : }
     630             : 
     631             : /* Return E_n */
     632             : GEN
     633          21 : eulerreal(long n, long prec)
     634             : {
     635          21 :   pari_sp av = avma;
     636             :   GEN B;
     637             :   long p, k;
     638          21 :   if (n < 0) pari_err_DOMAIN("eulerreal", "index", "<", gen_0, stoi(n));
     639          14 :   if (n == 0) return real_1(prec);
     640          14 :   if (odd(n)) return real_0(prec);
     641             : 
     642          14 :   k = n >> 1;
     643          14 :   if (!eulerzone) constreuler(0);
     644          14 :   if (k < lg(eulerzone)) return itor(gel(eulerzone,k), prec);
     645           0 :   p = eulerprec(n);
     646           0 :   B = eulerreal_using_lfun4(n, minss(p, prec));
     647           0 :   if (p < prec) B = itor(roundr(B), prec);
     648           0 :   return gerepileuptoleaf(av, B);
     649             : }

Generated by: LCOV version 1.14