How does scipy compute random variates for the distributions given in scipy.stats?

I am currently working on a small probability distributions package (primarily as an exercise) that I am basing on the probability distribution classes in SciPy. Specifically, scipy.stats and its continuous and discrete probability distribution classes.

Each of these classes have the "rvs" method that generates random variates for a given distribution. I scoured the SciPy GitHub repo to try and track down how this is done, but I am getting lost in the vast source code of SciPy. Here is my trail so far (using the normal distribution as an example):

scipy.stats._continuous_distns.norm_gen._rvs

scipy.stats._distn_infrastructure.rv_frozen.rvs

and this is where my reasoning ends. I assume I could use Newton's method or the bisection method with inverse transform sampling, but I sense there could be better optimized algorithms for particular distributions.

tl;dr: What is the best way to implement random variates in Python?

Is the best solution to use these numerical methods, or is there an existing catalog for specific random variate algorithms per type of distribution?