Extremely fast generation of gamma variates in Python

Through code profiling, I have identified np.random.gamma as the major bottleneck in my program. My program calls this function on each iteration with a very large shape array as its argument. I have also tried using np.random.Generator with PCG64:

import numpy as np
rg = np.random.Generator(np.random.PCG64())
rg.gamma(shape)

My tests show this is somewhat faster but still nowhere near fast enough. I'm hoping to get at least an order-of-magnitude improvement. Is there a faster way to generate a very large array of gamma-distributed random variates in Python? This is the algorithm being used by the current generator. There's a 2013 paper called Extremely efficient generation of Gamma random variables for α ≥ 1 but I don't know if it has been implemented.