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

Generated by: LCOV version 1.14