Inconsistent results for Inf integrands

I am looking at some SciML problems where, depending on the stability of the ODE, the integrand used w/ HCubature.jl may return Inf. I am noticing inconsistent behavior across `hquadrature` and `hcubature` for single and multi-dim domains.

```julia
using HCubature

g(x) = x[1] > 0.0 ? Inf : x[1]

hquadrature(g, -1.0, 1.0) #Inf
hcubature(g, -ones(1), ones(1)) #Inf
hcubature(g, -ones(2), ones(2)) #NaN
```

I see similar behavior in Cubature.jl as well, but `Cubature.hcubature` returns Inf up to dim = 4 and then NaN, while `Cubature.pcubature` returns Inf up dim = 17 (haven't test passed that).

Is this the expected behavior? I assume this is due to some Inf*0 that is propagating down.

Thanks.

Start by running the provided hquadrature and hcubature examples, then compare their Inf and NaN results across single- and multi-dimensional domains. Read the hquadrature and hcubature entry points to trace how these results are handled; done means the expected behavior is established and the inconsistency has a reproducible resolution.

Written by the indexing model from the issue text.