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.13.0 lcov report (development 25825-edc109b529) Lines: 194 201 96.5 %
Date: 2020-09-21 06:08:33 Functions: 20 21 95.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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                     BERNOULLI NUMBERS B_2k                     **/
      17             : /**                                                                **/
      18             : /********************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : /* D = divisorsu(n). Return a/b = \sum_{p-1 | 2n: p prime} 1/p
      23             :  * B_2k + a/b in Z [Clausen-von Staudt] */
      24             : static GEN
      25       31013 : fracB2k(GEN D)
      26             : {
      27       31013 :   GEN a = utoipos(5), b = utoipos(6); /* 1/2 + 1/3 */
      28       31289 :   long i, l = lg(D);
      29      166030 :   for (i = 2; i < l; i++) /* skip 1 */
      30             :   {
      31      135281 :     ulong p = 2*D[i] + 1; /* a/b += 1/p */
      32      135281 :     if (uisprime(p)) { a = addii(muliu(a,p), b); b = muliu(b,p); }
      33             :   }
      34       30749 :   return mkfrac(a,b);
      35             : }
      36             : /* precision needed to compute B_k for all k <= N */
      37             : long
      38        4588 : bernbitprec(long N)
      39             : { /* 1.612086 ~ log(8Pi) / 2 */
      40        4588 :   const double log2PI = 1.83787706641;
      41        4588 :   double t = (N + 0.5) * log((double)N) - N*(1 + log2PI) + 1.612086;
      42        4588 :   return (long)ceil(t / M_LN2) + 16;
      43             : }
      44             : static long
      45          21 : bernprec(long N) { return nbits2prec(bernbitprec(N)); }
      46             : /* \sum_{k > M} k^(-n) <= M^(1-n) / (n-1) < 2^-bit_accuracy(prec) */
      47             : static long
      48        4093 : zetamaxpow(long n, long prec)
      49             : {
      50        4093 :   long b = bit_accuracy(prec), M = (long)exp2((double)b/(n-1.0));
      51        4096 :   return M | 1; /* make it odd */
      52             : }
      53             : /* zeta(k) using 'max' precomputed odd powers */
      54             : static GEN
      55       30914 : bern_zeta(long k, GEN pow, long max, long prec)
      56             : {
      57       30914 :   GEN s = gel(pow, max);
      58             :   long j;
      59      208275 :   for (j = max - 2; j >= 3; j -= 2) s = addrr(s, gel(pow,j));
      60       30320 :   return divrr(addrs(s,1), subsr(1, real2n(-k, prec)));
      61             : }
      62             : /* z * j^2 */
      63             : static GEN
      64      209987 : mulru2(GEN z, ulong j)
      65             : { return (j | HIGHMASK)? mulru(mulru(z, j), j)
      66      209987 :                        : mulru(z, j*j); }
      67             : /* 1 <= m <= n, set y[1] = B_{2m}, ... y[n-m+1] = B_{2n} in Q */
      68             : static void
      69         269 : bernset(GEN *y, long m, long n)
      70             : {
      71         269 :   long i, j, k, bit, prec, max, N = n << 1; /* up to B_N */
      72             :   GEN A, C, pow;
      73         269 :   bit = bernbitprec(N);
      74         269 :   prec = nbits2prec(bit);
      75         268 :   A = sqrr(Pi2n(1, prec)); /* (2Pi)^2 */
      76         269 :   C = divrr(mpfactr(N, prec), powru(A, n)); shiftr_inplace(C,1);
      77         268 :   max = zetamaxpow(N, prec);
      78         268 :   pow = cgetg(max+1, t_VEC);
      79        2788 :   for (j = 3; j <= max; j += 2)
      80             :   { /* fixed point, precision decreases with j */
      81        2521 :     long b = bit - N * log2(j), p = b <= 0? LOWDEFAULTPREC: nbits2prec(b);
      82        2521 :     gel(pow,j) = invr(rpowuu(j, N, p));
      83             :   }
      84         267 :   y += n - m;
      85         267 :   for (i = n, k = N;; i--)
      86       30696 :   { /* set B_n, k = 2i */
      87       30963 :     pari_sp av2 = avma;
      88       30963 :     GEN B, z = fracB2k(divisorsu(i)), s = bern_zeta(k, pow, max, prec);
      89             :     long j;
      90             :     /* s = zeta(k), C = 2*k! / (2Pi)^k */
      91       30820 :     B = mulrr(s, C); if (!odd(i)) setsigne(B, -1); /* B ~ B_n */
      92       30950 :     B = roundr(addrr(B, fractor(z,LOWDEFAULTPREC))); /* B - z = B_n */
      93       30983 :     *y-- = gclone(gsub(B, z));
      94       30985 :     if (i == m) break;
      95       30720 :     affrr(divrunu(mulrr(C,A), k-1), C);
      96      239770 :     for (j = max; j >= 3; j -= 2) affrr(mulru2(gel(pow,j), j), gel(pow,j));
      97       30688 :     set_avma(av2);
      98       30689 :     k -= 2;
      99       30689 :     if ((k & 0xf) == 0)
     100             :     { /* reduce precision if possible */
     101        3811 :       long bit2 = bernbitprec(k), prec2 = nbits2prec(bit2), max2;
     102        3827 :       if (prec2 == prec) continue;
     103        3827 :       prec = prec2;
     104        3827 :       max2 = zetamaxpow(k,prec);
     105        3829 :       if (max2 > max) continue;
     106        3813 :       bit = bit2;
     107        3813 :       max = max2;
     108        3813 :       setprec(C, prec);
     109       28755 :       for (j = 3; j <= max; j += 2)
     110             :       {
     111       24953 :         GEN P = gel(pow,j);
     112       24953 :         long b = bit + expo(P), p = b <= 0? LOWDEFAULTPREC: nbits2prec(b);
     113       24946 :         if (realprec(P) > p) setprec(P, p);
     114             :       }
     115             :     }
     116             :   }
     117         265 : }
     118             : /* need B[2..2*nb] as t_INT or t_FRAC */
     119             : void
     120      244363 : constbern(long nb)
     121             : {
     122      244363 :   const pari_sp av = avma;
     123             :   long i, l;
     124             :   GEN B;
     125             :   pari_timer T;
     126             : 
     127      244363 :   l = bernzone? lg(bernzone): 0;
     128      244363 :   if (l > nb) return;
     129             : 
     130         268 :   nb = maxss(nb, l + 127);
     131         269 :   B = cgetg_block(nb+1, t_VEC);
     132         269 :   if (bernzone)
     133        6144 :   { for (i = 1; i < l; i++) gel(B,i) = gel(bernzone,i); }
     134             :   else
     135             :   {
     136         228 :     gel(B,1) = gclone(mkfracss(1,6));
     137         228 :     gel(B,2) = gclone(mkfracss(-1,30));
     138         228 :     gel(B,3) = gclone(mkfracss(1,42));
     139         228 :     gel(B,4) = gclone(mkfracss(-1,30));
     140         228 :     gel(B,5) = gclone(mkfracss(5,66));
     141         227 :     gel(B,6) = gclone(mkfracss(-691,2730));
     142         228 :     gel(B,7) = gclone(mkfracss(7,6));
     143         227 :     gel(B,8) = gclone(mkfracss(-3617,510));
     144         228 :     gel(B,9) = gclone(mkfracss(43867,798));
     145         228 :     gel(B,10)= gclone(mkfracss(-174611,330));
     146         228 :     gel(B,11)= gclone(mkfracss(854513,138));
     147         228 :     gel(B,12)= gclone(mkfracss(-236364091,2730));
     148         227 :     gel(B,13)= gclone(mkfracss(8553103,6)); /* B_26 */
     149         228 :     l = 14;
     150             :   }
     151         269 :   set_avma(av);
     152         269 :   if (DEBUGLEVEL) {
     153           0 :     err_printf("caching Bernoulli numbers 2*%ld to 2*%ld\n", l, nb);
     154           0 :     timer_start(&T);
     155             :   }
     156         269 :   bernset((GEN*)B + l, l, nb);
     157         265 :   if (DEBUGLEVEL) timer_printf(&T, "Bernoulli");
     158         265 :   swap(B, bernzone); guncloneNULL(B);
     159         265 :   set_avma(av);
     160             : }
     161             : /* Obsolete, kept for backward compatibility */
     162             : void
     163           0 : mpbern(long n, long prec) { (void)prec; constbern(n); }
     164             : 
     165             : /* assume n even > 0, if iz != NULL, assume iz = 1/zeta(n) */
     166             : static GEN
     167          21 : bernreal_using_zeta(long n, long prec)
     168             : {
     169          21 :   GEN pi2 = Pi2n(1, prec+EXTRAPRECWORD);
     170          21 :   GEN iz = inv_szeta_euler(n, prec);
     171          21 :   GEN z = divrr(mpfactr(n, prec), mulrr(powru(pi2, n), iz));
     172          21 :   shiftr_inplace(z, 1); /* 2 * n! * zeta(n) / (2Pi)^n */
     173          21 :   if ((n & 3) == 0) setsigne(z, -1);
     174          21 :   return z;
     175             : }
     176             : /* assume n even > 0, B = NULL or good approximation to B_n */
     177             : static GEN
     178          14 : bernfrac_i(long n, GEN B)
     179             : {
     180          14 :   pari_sp av = avma;
     181          14 :   GEN z = fracB2k(divisorsu(n >> 1));
     182          14 :   if (!B) B = bernreal_using_zeta(n, bernprec(n));
     183          14 :   B = roundr( gadd(B, fractor(z,LOWDEFAULTPREC)) );
     184          14 :   return gerepileupto(av, gsub(B, z));
     185             : }
     186             : GEN
     187     9503852 : bernfrac(long n)
     188             : {
     189             :   long k;
     190     9503852 :   if (n <= 1)
     191             :   {
     192        4018 :     if (n < 0) pari_err_DOMAIN("bernfrac", "index", "<", gen_0, stoi(n));
     193        4011 :     return n? mkfrac(gen_m1,gen_2): gen_1;
     194             :   }
     195     9499834 :   if (odd(n)) return gen_0;
     196     9495949 :   k = n >> 1;
     197     9495949 :   if (!bernzone) constbern(0);
     198     9495949 :   if (bernzone && k < lg(bernzone)) return gel(bernzone, k);
     199          19 :   return bernfrac_i(n, NULL);
     200             : }
     201             : GEN
     202         175 : bernvec(long n)
     203             : {
     204             :   long i, l;
     205             :   GEN y;
     206         175 :   if (n < 0) return cgetg(1, t_VEC);
     207         175 :   constbern(n);
     208         175 :   l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
     209         896 :   for (i = 2; i < l; i++) gel(y,i) = gel(bernzone,i-1);
     210         175 :   return y;
     211             : }
     212             : 
     213             : /* x := pol_x(v); B_k(x) = \sum_{i=0}^k binomial(k, i) B_i x^{k-i} */
     214             : static GEN
     215        1988 : bernpol_i(long k, long v)
     216             : {
     217             :   GEN B, C;
     218             :   long i;
     219        1988 :   if (v < 0) v = 0;
     220        1988 :   constbern(k >> 1); /* cache B_2, ..., B_2[k/2] */
     221        1988 :   C = vecbinomial(k);
     222        1988 :   B = cgetg(k + 3, t_POL);
     223       14357 :   for (i = 0; i <= k; ++i) gel(B, k-i+2) = gmul(gel(C,i+1), bernfrac(i));
     224        1988 :   B[1] = evalsigne(1) | evalvarn(v);
     225        1988 :   return B;
     226             : }
     227             : GEN
     228        1918 : bernpol(long k, long v)
     229             : {
     230        1918 :   pari_sp av = avma;
     231        1918 :   if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
     232        1911 :   return gerepileupto(av, bernpol_i(k, v));
     233             : }
     234             : /* x := pol_x(v); return 1^e + ... + x^e = x^e + (B_{e+1}(x) - B_{e+1})/(e+1) */
     235             : static GEN
     236          42 : faulhaber(long e, long v)
     237             : {
     238             :   GEN B;
     239          42 :   if (e == 0) return pol_x(v);
     240          28 :   B = RgX_integ(bernpol_i(e, v)); /* (B_{e+1}(x) - B_{e+1}) / (e+1) */
     241          28 :   gel(B,e+2) = gaddgs(gel(B,e+2), 1); /* add x^e, in place */
     242          28 :   return B;
     243             : }
     244             : /* sum_v T(v), T a polynomial expression in v */
     245             : GEN
     246          49 : sumformal(GEN T, long v)
     247             : {
     248          49 :   pari_sp av = avma, av2;
     249             :   long i, t, d;
     250             :   GEN R;
     251             : 
     252          49 :   T = simplify_shallow(T);
     253          49 :   t = typ(T);
     254          49 :   if (is_scalar_t(t))
     255          14 :     return gerepileupto(av, monomialcopy(T, 1, v < 0? 0: v));
     256          35 :   if (t != t_POL) pari_err_TYPE("sumformal [not a t_POL]", T);
     257          28 :   if (v < 0) v = varn(T);
     258          28 :   av2 = avma;
     259          28 :   R = gen_0;
     260          28 :   d = poldegree(T,v);
     261          91 :   for (i = d; i >= 0; i--)
     262             :   {
     263          63 :     GEN c = polcoef(T, i, v);
     264          63 :     if (gequal0(c)) continue;
     265          42 :     R = gadd(R, gmul(c, faulhaber(i, v)));
     266          42 :     if (gc_needed(av2,3))
     267             :     {
     268           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"sumformal, i = %ld/%ld", i,d);
     269           0 :       R = gerepileupto(av2, R);
     270             :     }
     271             :   }
     272          28 :   return gerepileupto(av, R);
     273             : }
     274             : 
     275             : /* 1/zeta(n) using Euler product. Assume n > 0. */
     276             : GEN
     277        1194 : inv_szeta_euler(long n, long prec)
     278             : {
     279             :   GEN z, res;
     280             :   pari_sp av, av2;
     281             :   double A, D, lba;
     282             :   ulong p, lim;
     283             :   forprime_t S;
     284             : 
     285        1194 :   if (n > prec2nbits(prec)) return real_1(prec);
     286             : 
     287        1187 :   lba = prec2nbits_mul(prec, M_LN2);
     288        1187 :   D = exp((lba - log((double)(n-1))) / (n-1));
     289        1187 :   lim = 1 + (ulong)ceil(D);
     290        1187 :   if (lim < 3) return subir(gen_1,real2n(-n,prec));
     291        1187 :   res = cgetr(prec); incrprec(prec);
     292        1187 :   av = avma;
     293        1187 :   z = subir(gen_1, real2n(-n, prec));
     294             : 
     295        1187 :   (void)u_forprime_init(&S, 3, lim);
     296        1187 :   av2 = avma; A = n / M_LN2;
     297       86198 :   while ((p = u_forprime_next(&S)))
     298             :   {
     299       85011 :     long l = prec2nbits(prec) - (long)floor(A * log((double)p)) - BITS_IN_LONG;
     300             :     GEN h;
     301             : 
     302       85011 :     if (l < BITS_IN_LONG) l = BITS_IN_LONG;
     303       85011 :     l = minss(prec, nbits2prec(l));
     304       85011 :     h = divrr(z, rpowuu(p, (ulong)n, l));
     305       85011 :     z = subrr(z, h);
     306       85011 :     if (gc_needed(av,1))
     307             :     {
     308           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"inv_szeta_euler, p = %lu/%lu", p,lim);
     309           0 :       z = gerepileuptoleaf(av2, z);
     310             :     }
     311             :   }
     312        1187 :   affrr(z, res); set_avma(av); return res;
     313             : }
     314             : 
     315             : /* Return B_n */
     316             : GEN
     317        4641 : bernreal(long n, long prec)
     318             : {
     319             :   GEN B;
     320             :   long p, k;
     321        4641 :   if (n < 0) pari_err_DOMAIN("bernreal", "index", "<", gen_0, stoi(n));
     322        4634 :   if (n == 0) return real_1(prec);
     323        4634 :   if (n == 1) return real_m2n(-1,prec); /*-1/2*/
     324        4634 :   if (odd(n)) return real_0(prec);
     325             : 
     326        4634 :   k = n >> 1;
     327        4634 :   if (!bernzone) constbern(0);
     328        4634 :   if (k < lg(bernzone)) return fractor(gel(bernzone,k), prec);
     329           7 :   p = bernprec(n);
     330           7 :   B = bernreal_using_zeta(n, minss(p, prec));
     331           7 :   if (p < prec) B = fractor(bernfrac_i(n, B), prec);
     332           7 :   return B;
     333             : }
     334             : 
     335             : GEN
     336          56 : eulerpol(long k, long v)
     337             : {
     338          56 :   pari_sp av = avma;
     339             :   GEN B, E;
     340          56 :   if (k < 0) pari_err_DOMAIN("eulerpol", "index", "<", gen_0, stoi(k));
     341          49 :   k++; B = bernpol_i(k, v);
     342          49 :   E = RgX_Rg_mul(RgX_sub(B, RgX_rescale(B, gen_2)), sstoQ(2,k));
     343          49 :   return gerepileupto(av, E);
     344             : }
     345             : GEN
     346           7 : eulervec(long n)
     347             : {
     348             :   pari_sp av;
     349             :   GEN v, E, C;
     350             :   long k;
     351           7 :   if (n < 0) return cgetg(1, t_VEC);
     352           7 :   C = vecbinomial(2*n);
     353           7 :   E = ZX_translate(RgX_rescale(eulerpol(2*n, 0), gen_2), gen_1);
     354           7 :   av = avma; v = cgetg(n + 2, t_VEC); gel(v,1) = gen_1;
     355             :   /* 2^(2n) E_2n(x/2 + 1) = sum_k binomial(2n,k) E_k x^(n-k) */
     356         217 :   for (k = 1; k <= n; k++)
     357         210 :     gel(v,k+1) = diviiexact(gel(E,2*n-2*k+2), gel(C,2*k+1));
     358           7 :   return gerepileupto(av, v);
     359             : }

Generated by: LCOV version 1.13