[BUG] Numerical integration issues for extreme parameters

[hsloot]

Description

Some extreme parametrization will cause a failure of used numerical integration techniques. This cannot be avoided completely, but sometimes improvements might be possible. This particular case arose when new tests were added and the randomly generated test-parametrizations changed unfavorably.

Reprex

library(rmo)
library(tidyverse)

# parameter
x <- 3.5
alpha <- 0.0902969

# Kendall's tau
(ktau <- 2 - 2 ^ alpha)
#> [1] 0.9354108

# 
try({
    alpha |>
        AlphaStableBernsteinFunction() |>
        valueOf(
            x, 
            method = "stieltjes",
            tolerance = .Machine$double.eps^0.5
        )
})
#> Error in FUN(X[[i]], ...) : 
#>   Numerical integration failed with error: roundoff error is detected in the extrapolation table

integrand_fn <- function(u) {
    (1 - exp(-u * x)) *
        alpha * 1 / gamma(1 - alpha) / u ^ (1 + alpha)
}

ggplot() + 
    xlim(.Machine$double.eps^0.5, .Machine$double.eps^0.1) + 
    geom_function(fun = integrand_fn)

^{Created on 2023-02-07 with reprex v2.0.2}

[hsloot]

Some observations:

• It is unclear whether the roundoff error in the extrapolation table could be ignored if the error estimate is still below the requested tolerance. In any case, I think reporting the issue and letting the user decide is prudent.
• For the primary application, calculating the Bernstein function itself using the integral representation is not required. Is it necessary that this is a test case and potentially causes more issues down the road?
• Is it possible, in general, to scale the integral and improve stability for cases like this?