I have forgotten
my Password

Or login with:

  • Facebookhttp://facebook.com/
  • Googlehttps://www.google.com/accounts/o8/id
  • Yahoohttps://me.yahoo.com
get GPL
this unit 1.00
sub units 0.00

Bernoulli B

Calculates array of Bernoulli numbers using an infinite series
Controller: CodeCogs

get GPL add to cart



Bernoulli B

double*dB )
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:

Substitution x=0 gives series expansion for Bernoulli numbers:

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]:

Then Fourier series expansion results in the following:

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:

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

then sequence of new sums

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

or, using

Array dimension should be iMax+1 or greater.


  • 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;
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


iMaxinput maximal index requested
dBoutput pointer on the array of numbers declared in the calling module.


Anatoly Prognimack (Mar 16, 2005)
Developed tested with Borland C++ 3.1 for DOS and Microsoft Visual C++ 5.0, 6.0
Updated by Will Bateman (March 2005)
Source Code

Source code is available when you agree to a GP Licence or buy a Commercial Licence.

Not a member, then Register with CodeCogs. Already a Member, then Login.