Hello—

Using AD with the Gauss 2F1 function has come up in some work I am doing. I am incredibly impressed by the performance and functionality of HypergeometricFunctions.jl, but I’m having a strange issue where only for only a somewhat narrow band of values for the second argument I’m getting a stack overflow error. I’m having a lot of trouble narrowing down where in the source I should be looking and I’m hoping somebody can provide advice.

Here is a MWE minus the unicode, which I can’t get to work here:

```
using HypergeometricFunctions, ForwardDiff
fun(arg, v) = _\_2 F \_1 (0.5, 0.5*(v+1), 1.5, -(arg^2)/v)
ForwardDiff.derivative(v->fun(2.4, v), 1.825) # error
ForwardDiff.derivative(v->fun(2.4, v), 1.8) # works fine
```

The error is a stack overflow on `unsafe_gamma`

on a line which defines this:

```
unsafe_gamma(x::Real) = unsafe_gamma(float(x))
```

and then this function gets called from this one:

```
G(z::Number, ϵ::Number) = ϵ == 0 ? digamma(z)/unsafe_gamma(z) : (inv(unsafe_gamma(z))-inv(unsafe_gamma(z+ϵ)))/ϵ
```

From some print debugging, it is the call to `inv_unsafe(z+\epsilon)`

that gets caught in an infinite loop.

So here’s where I can’t figure out where to go: `unsafe_gamma`

is defined for duals, using `DualNumbers`

. And when I print `\epsilon`

, I get a dual. But that branch of the dispatch doesn’t seem to be getting called, and `isreal(\epsilon)`

even seems to return `true`

, which is what I think is causing the issue. I’m pretty stumped about what could be causing that. Can somebody help me?

Thank you so much in advance for reading.