Basically quadrature is never used in optimized special-function implementations, though it’s great to have for comparison and bootstrapping purposes.
Instead, optimized special functions usually employ a combination of polynomial and rational expansions, e.g. continued-fraction expansions for large argument and Taylor series near roots. The tricky part is figuring out what approximation to use where, and how to stitch them together.
Unfortunately, this becomes more and more difficult to implement well as the special functions have more parameters.
Generally, I would recommend against ODE solvers to perform integrals, since that discards a lot of structure from the problem.
If you fix the parameters \alpha, \kappa, p, q and want to compute \gamma(x) and its inverse for some range of x values, then you can do much better. That is, suppose you are computing \gamma^{-1}(y) many times for given \alpha, \kappa, p, q, within a finite range of x and y values, and you can afford some startup cost.
Then, instead of using the DifferentialEquations interpolant, I would just construct a smooth polynomial approximation from a set of quadrature evaluations, e.g. a Chebyshev interpolant using a package like FastChebInterp.jl or ApproxFun.jl. This converges exponentially fast for smooth functions, much more accurate than an ODE interpolant. Given this interpolant for \gamma(x), you can then construct an interpolant for \gamma^{-1}(y) by applying Newton’s method to y values at Chebyshev points in your interval of interest. Thereafter, computing \gamma^{-1}(z) for any z in your interval is just a fast polynomial evaluation.
See also Approximate inverse of a definite integral - #3 by stevengj — as noted there, you can just do Newton’s method directly on the QuadGK results and then compute a Chebyshev interpolant of \gamma^{-1}(y) directly, skipping the intermediate step of constructing an interpolant for interpolant for \gamma(x).