(Why) is the std::binomial_distribution biased for large probabilities p and slow for small n?

I want to generate binomially distributed random numbers in c++. Speed is a major concern. Not knowing a lot about random number generators, I use the standard libraries' tools. My code looks like something below:

#include <random>

static std::random_device random_dev;
static std::mt19937 random_generator{random_dev()};
std::binomial_distribution<int> binomial_generator;

void RandomInit(int s) {
    //I create the generator object here to save time. Does this make sense?
    binomial_generator = std::binomial_distribution<int>(1, 0.5);
    random_generator.seed(s);
}

int binomrand(int n, double p) {
    binomial_generator.param(std::binomial_distribution<int>::param_type(n, p));
    return binomial_generator(random_generator);
}

To test my implementation, I have built a cython wrapper and then executed and timed the function from within python. For reference I have also implemented a "stupid" binomial distribution, which just returns the sum of Bernoulli trials.

int binomrand2(int n, double p) {
    int result = 0;
    for (int i = 0; i<n; i++) {
        if (_Random() < p) //_Random is a thoroughly tested custom random number generator on U[0,1)
            result++;
    }

    return result;
}

Timing showed that the latter implementation is about 50% faster than the former if n < 25. Furthermore, for p = 0.95, the former yielded significantly biased results (the mean over 1000000 trials for n = 40 was 38.23037; standard deviation is 0.0014; the result was reproducable with different seeds).

Is this a (known) issue with the standard library's functions or is my implementation wrong? What could I do to achieve my goal of obtaining accurate results with high efficiency?

The parameter n will mostly be below 100 and smaller values will occur more frequently.

I am open to suggestions outside the realm of the standard library, but I may not be able to use external software libraries.

I am using the VC 2019 compiler on 64bit Windows.

Edit

I have also tested the bias without using python:

double binomrandTest(int n, double p, long long N) {
    long long result = 0;
    for (long long i = 0; i<N; i++) {
        result += binomrand(n, p);
    }
    return ((double) result) / ((double) N);
}

The result remained biased (38.228045 for the parameters above, where something like 38.000507 would be expected).