Skip to the content of the web site.

Project S.3: Bernoulli numbers

The calculation of Bernoulli numbers has a special place in the history of computer programming: the first program written for the Babbage machine by Ada, Countess of Lovelace, was to generate the Bernoulli numbers.

$B_n = {\begin{cases} \phantom{-}1&{\mbox{if }}n=0\\ -\frac{1}{2}&{\mbox{if }}n=1\\ -\sum_{k = 0}^{n - 1} {n \choose k} \frac{B_k}{n - k + 1}&{\mbox{if }}n \geq 2\end{cases}}$.

You will write a program that calculates the Bernoulli numbers from $0$ to $n$ and return a dynamically allocated array storing these numbers.

double *bernoulli( std::size_t n );

Once you take a look at your output, you may notice something peculiar about many of the numbers. After reading the Wikipedia page on Bernoulli numbers, what simplification could you introduce to your function to improve the efficiency and perhaps even the results?

Next, the exact value of $B_{14} = \frac{7}{6} = 1.1666\cdots$. How close is your approximation to this value? Can you think of how you could improve the accuracy of your implementation?

If you are just writing a recursive function, consider using the following main program:

#include <iostream>
int main();
double bernoulli( std::size_t n );

int main() {
    std::size_t const MAX_BERNOULLI{ 20 };

    // Set the output to display 16 significant decimal digits
    //  - this is close to the precision of the double-precision
    //    floating-point representation
    std::cout.precision( 16 );

    for ( std::size_t k{0}; k <= MAX_BERNOULLI; ++k ) {
        // This formatting lines up the output more nicely by ensuring
        // that the equal signs and the leading digits are lined up
        std::cout << "  B_" << k << "\t= "
                  << ((bernoulli( k ) >= 0) ? " " : "")
                  << bernoulli( k ) << std::endl;
    }

    return 0;
}

double bernoulli( std::size_t n ) {
    return 0.0;
}

If, however, you're writing a function that calculates all Bernoulli numbers up to $n$ at once, here is a program you can use:

#include <iostream>
int main();
double *bernoulli( std::size_t n );

int main() {
    std::size_t const MAX_BERNOULLI{ 20 };
    double *bernoulli_numbers{ bernoulli( MAX_BERNOULLI ) };

    // Set the output to display 16 significant decimal digits
    //  - this is close to the precision of the double-precision
    //    floating-point representation
    std::cout.precision( 16 );

    for ( std::size_t k{0}; k <= MAX_BERNOULLI; ++k ) {
        // This formatting lines up the output more nicely by ensuring
        // that the equal signs and the leading digits are lined up
        std::cout << "  B_" << k << "\t= "
                  << ((bernoulli_numbers[k] >= 0) ? " " : "")
                  << bernoulli_numbers[k] << std::endl;
    }

    delete[] bernoulli_numbers;

    return 0;
}

double *bernoulli( std::size_t n ) {
    double result{ new double[n + 1] };

    return result;
}

Finally, how long does it take for your algorithm to print the first 30 Bernoulli numbers? Is it reasonably fast or there a noticable delay of at least a few seconds? What algorithm did you use for calcluating ${n \choose k}$?