CodeCogs

Calculates array of Bernoulli numbers using an infinite series

You're viewing an older version of this page (#82). View the current version.

Interface

#include <codecogs/maths/discrete/number_theory/bernoulli_b.h>

using namespace Maths::Discrete::Number_Theory;

for( int j = 3;
const double eps = DBL_EPSILON;
const double TwoPi = 6.2831853071795860;
for( int n = 1;

Overview

Bernoulli numbers $B_n$ are particular values of Bernoulli polynomials $B_n(x): B_n=B_n(0)$ . Polynomials $B_n(x)$ are expandable in the Fourier series:

$B_{2n}(x)=2{{(-1)}^{n+1}}(2n)!\sum _{r=1}^{\infty }{{(2 \pi r)}^{-2n}}cos(2 \pi r x)$

(1)

$\begin{array}{l} {B_{2n+1}}(x)= 2{{(-1)}^{n+1}}(2n+1)!\sum _{r=1}^{\infty }{{(2 \pi r)}^{-2n-1}} sin(2 {\pi rx}) \end{array}$

(2)

Substitution $x=0$ gives series expansion for Bernoulli numbers:

$\begin{array}{l} {B_{2n}}=2{{(-1)}^{n+1}}(2n)! \sum _{r=1}^{\infty }{{(2 \pi r)}^{-2n}} \\ n=1,2,3,... \\ {B_{2n+1}}=0,\\ n=1,2,3,...\\ \end{array}$

(3)

However, series like this is not quite suitable for numerical evaluation, especially at small n, in view of slow convergence and uneasy accuracy estimation. One needs to use another relation [2,23.10.21]:

${B_{n(}} \big ( \frac{1}{2}\big)=-(1-{2^{1-n}}){B_n}$

(4)

Then Fourier series expansion results in the following:

$\begin{array}{l} {B_{2n}}=\frac{{{(-1)}^{n+1}}2(2n)!}{(1-{2^{1-2n}}){{(2\pi )}^{2n}}}\sum _{r=1}^{\infty }\frac{{{(-1)}^{r-1}}}{{r^{2n}}}, \\ n=1,2,3,....\\ \end{array}$

(5)

Alternating series converges faster and gives possibility for simple accuracy estimation: error arising due to series truncation is less than first term neglected. Then N, the amount of series terms to obtain accuracy $(\epsilon )$ may be estimated as follows:

$N={{\epsilon }^{-1/2n}}$

(6)

In the worst case n=2 reasonable accuracy $(\epsilon ={10}^{-12} )$ results in $({10}^6)$ series terms to sum. To significantly reduce a mount of summands and speed up convergence it is suitable to use Aitken transformation of partial series sums. Let

${S_n}=\sum _{k=1}^{n}{f_k}$

(7)

then sequence of new sums

${T_n}= \frac{S_{n+1}S_{n-1} - S^2_n}{ S_{n+1} + S_{n-1} - 2{S_n}}$

(8)

may converge faster compared with $(S_n)$ . To avoid subtractions of near values and losses of significant digits expression should be rewritten as follows:

${T_n}= {S_n}-\frac{f_{n+1}f_n}{f_{n+1}-f_n}$

(9)

or, using

${f_n}= {(-1)^n}/{g_n}$

(10)

${T_n}= S_n + \frac{{(-1)}^n}{g_{n+1}} + g_n$

(11)

Array dimension should be iMax+1 or greater.

References

Higher Transcendental Functions, vol.1, (1.13) by H.Bateman and A.Erdelyi (Bateman Manuscript Project), 1953
M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions, 1964 chapt.23
Yu.Luke, Mathematical functions and their approximations, 1975 chapt.14.2

Example 1

#include <stdio.h>
#include <codecogs/maths/discrete/number_theory/bernoulli_b.h>

#define MAX_INDEX 16

int main()
{
  double dBernoulli[MAX_INDEX+1];

  printf( "%8s%2c%20s\n", " ", 'n', "Bn" );
  printf( "%8s", " " );
  for(int i = 0; i < 22; i++ )
    printf( "%c", '-' );
  printf( "\n" );

  Maths::NumberTheory::bernoulli_B( MAX_INDEX, dBernoulli );

  printf( "%10d%20.12f\n", 0, dBernoulli[0] );
  printf( "%10d%20.12f\n", 1, dBernoulli[1] );
  for(int i = 2; i <= MAX_INDEX; i += 2 )
    printf( "%10d%20.12f\n", i, dBernoulli[i] );
  return 0;
}

Output:

n                  Bn
----------------------
0      1.000000000000
1     -0.500000000000
2      0.166666666667
4     -0.033333333333
6      0.023809523810
8     -0.033333333333
10     0.075757575758
12    -0.253113553114
14     1.166666666667
16    -7.092156862745

Parameters

iMax

input maximal index requested

output pointer on the array of numbers declared in the calling module.

FUNCTION

for

odd index elements

DECLARATION

eps

computer dependent value ultimate calculation accuracy

DECLARATION

TwoPi

minimal eps such that 1.0 + eps != 1.0