# Fast implementation of incomplete extended gamma function, QuadGK?

**URL:** <https://discourse.julialang.org/t/fast-implementation-of-incomplete-extended-gamma-function-quadgk/112239>\
**Category:** Numerics\
**Tags:** question, specialfunctions\
**Created:** [March 28, 2024, 4:35pm UTC](https://discourse.julialang.org/t/fast-implementation-of-incomplete-extended-gamma-function-quadgk/112239 "2024-03-28T16:35:17Z")\
**Posts on this page:** 1\
**Showing post:** 2

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 28, 2024, 6:35pm UTC](https://discourse.julialang.org/t/fast-implementation-of-incomplete-extended-gamma-function-quadgk/112239/2 "2024-03-28T18:35:20Z")

</div>

> [@mleprovost](#):
>
> I attach my current implementation based on QuadGK. I would appreciate any feedback to accelerate computations or reduce allocations.

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.

> [@mleprovost](#):
>
> Can it be beneficial to define the function \gamma as the solution of the ODE

Generally, I would recommend against ODE solvers to perform integrals, since that discards a lot of structure from the problem.

> [@mleprovost](#):
>
> Then we can compute its inverse by using a root-finding algorithm combined with the interpolation for free tools from DifferentialEquations

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](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/3) — 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).

---

_[View the full topic](https://discourse.julialang.org/t/fast-implementation-of-incomplete-extended-gamma-function-quadgk/112239)._