Bernoulli numbers with Boost are different from Mathematica

In the latest Boost, there is a function to compute the Bernoulli number, but I miss what it does exactly.

For example, Mathematica, Python mpmath and www.bernoulli.org say that:

BernoulliB[1] == -1/2

but the boost version

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/math/special_functions/bernoulli.hpp>

boost::multiprecision::cpp_dec_float_100 x =  bernoulli_b2n<boost::multiprecision::cpp_dec_float_100>(1);

returns 0.166667

Why this difference? Am I missing something?

All odd Bernoulli numbers are zero, apart of B1, which you know is -1/2. So, boost::math::bernoulli_b2n returns the only even (2nth) Bernoulli numbers.

For example, to get B4 you need to actually pass 2:

std::cout
    << std::setprecision(std::numeric_limits<double>::digits10)
    << boost::math::bernoulli_b2n<double>(2) << std::endl;

and if you pass 1, you get B2.

See docs: http://www.boost.org/doc/libs/1_56_0/libs/math/doc/html/math_toolkit/number_series/bernoulli_numbers.html

Of course, you can make a simple wrapper, to imitate preferred syntax¹:

double bernoulli(int n)
{
    if (n == 1) return -1.0 / 2.0; // one
    if (n % 2) return 0; // odd 
    return boost::math::bernoulli_b2n<double>(n / 2);
}

int main()
{
    std::cout << std::setprecision(std::numeric_limits<double>::digits10);
    for (int i = 0; i < 10; ++i)
    {
        std::cout << "B" << i << "\t" << bernoulli(i) << "\n";
    }
}

or even a class with overloaded operator[] (for demanding persons ;) ):

class Bernoulli
{
public:
    double operator[](int n)
    {
        return bernoulli(n);
    }
};

or even make use of template magic and do all this checks at compile time (I will left it as an exercise for a reader ;) ).

_{¹Please note, that this exact function body is not well verified and can contains mistakes. But I hope you've got the idea of how you can make a wrapper.}