bernmm is a C++ implementation of the algorithm described in the paper A multimodular algorithm for computing Bernoulli numbers.
It depends on GMP and NTL, and is freely available under a BSD-style license.
It is included in Sage since version 3.1, accessible via the bernoulli function as follows:
sage: bernoulli(100, algorithm="bernmm", num_threads=3) -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330
The code for computing Bk mod p is also available (the following takes about 0.03 seconds on my laptop):
sage: p = 10^8 + 7 sage: bernoulli_mod_p_single(p, 255018) 2