There are two conventions for the Bernoulli numbers, specifically $B_1$ can be defined to be 1/2 or -1/2 (see https://en.wikipedia.org/wiki/Bernoulli_number), depending on whether it is defined using the exponential generating function $\frac{x}{e^x-1}$ or $\frac{xe^x}{e^x - 1}$. SymPy uses the -1/2 definition.

Donald Knuth has recently written that he's been convinced that the -1/2 definition is wrong and historically inaccurate, and he updated Concrete Mathematics accordingly (https://www-cs-faculty.stanford.edu/~knuth/news22.html):

The sequence of Bernoulli numbers plays a prominent role in mathematics, especially in connection with asymptotic expansions related to Euler's summation formula. When I learned about this fascinating sequence, during my undergraduate days, I was taught that $B_1$ is equal to minus one-half. So I duly taught the same to my students, and went on to write books that explained what I thought Bernoulli and Euler had done.

But last year I took a close look at Peter Luschny's Bernoulli manifesto, where he gives more than a dozen good reasons why the value of $B_1$ should really be plus one-half. He explains that some mathematicians of the early 20th century had unilaterally changed the conventions, because some of their formulas came out a bit nicer when the negative value was used. It was their well-intentioned but ultimately poor choice that had led to what I'd been taught in the 1950s.

Luschny's webpage cites, for example, recent treatments of the subject by leading mathematicians such as Terence Tao. And his most compelling argument, from my personal perspective, is the way he unveils the early publications: I learned from him that my own presentation of the story, in The Art of Computer Programming and much more extensively in Concrete Mathematics, was a violation of history! I had put words and thoughts into Bernoulli and Euler's minds that were not theirs at all. This hurt, because I've always tried to present the evolution of ideas faithfully; in this case I'd fooled myself, by trying to conform what they wrote to what I'd learned.

By now, hundreds of books that use the “minus-one-half” convention have unfortunately been written. Even worse, all the major software systems for symbolic mathematics have that 20th-century aberration deeply embedded. Yet Luschny convinced me that we have all been wrong, and that it's high time to change back to the correct definition before the situation gets even worse.

Therefore I changed the definition of $B_1$ in all printings of The Art of Computer Programming during the latter half of 2021. And the new (34th) printing of Concrete Mathematics, released in January 2022, contains the much more extensive changes that are needed to tell a more comprehensive story.

Here is Luschny's page he cites http://luschny.de/math/zeta/The-Bernoulli-Manifesto.html.

So maybe we should reconsider SymPy's definition of bernoulli. At the very least, I think there should be a way to specify which convention is used, $B^+_n$ or $B^-_n$.

As an aside, the docs for bernoulli state that the defining recurrence relation is $$0 = \sum_{k=0}^n \binom{n+1}{k} B_k$$ but for the -1/2 definition it should actually be $$\delta_{m,0} = \sum_{k=0}^n \binom{n+1}{k} B_k$$ (for the +1/2 definition it's $$n + 1 = \sum_{k=0}^n \binom{n+1}{k} B_k,$$ see https://en.wikipedia.org/wiki/Bernoulli_number#Recursive_definition).