efficient sampling from beta-binomial distribution in python

for a stochastic simulation I need to draw a lot of random numbers which are beta binomial distributed.

At the moment I implemented it this way (using python):

import scipy as scp
from scipy.stats import rv_discrete

class beta_binomial(rv_discrete):
       """
       creating betabinomial distribution by defining its pmf
       """
       def _pmf(self, k, a, b, n):
          return scp.special.binom(n,k)*scp.special.beta(k+a,n-k+b)/scp.special.beta(a,b)

so sampling a random number x can be done by:

betabinomial = beta_binomial(name="betabinomial")
x = betabinomial.rvs(0.5,0.5,3) # with some parameter 

The problem is, that sampling one random number takes ca. 0.5ms, which is in my case dominating the whole simulation speed. The limiting element is the evaluation of the beta functions (or gamma functions within these).

Does anyone has a great idea how to speed up the sampling?