Cumulative Distribution Function of the skewed generalized error distribution I´m trying to calculate the cumulative distribution function of the skewed generalized error distribution with the probability density function:

where u = y - m.
From Theodossiou in (https://www.mfsociety.org/modules/modDashboard/uploadFiles/journals/MJ~0~p1a4fjq38m1k2p45t6481fob7rp4.pdf). I tried some calculations but they are not worth mentioning.

Applying the general recommendations for recognizing distributions set out at https://stats.stackexchange.com/a/332844/919, we see (without any calculation) that up to a change of scale (differing according to whether $u$ is positive or negative), these densities are proportional to $\exp(-|u|^k).$
The absolute value and the presence of the signum ("sign of") function suggest breaking the integration of the density into positive and negative regions (and reversing the sign of $u$ on the negative region).  This reduces the question to the following:

Evaluate $\int_0^u \exp(-t^k)\mathrm dt$ for $u \gt 0.$

The most direct and obvious route to simplification is to make $z = t^k$ the variable of integration.  This is expedited by using the logarithmic differential $\mathrm dz/z = \mathrm d(t^k)/t^k = k\,\mathrm dt/t$ (see https://stats.stackexchange.com/a/185709/919):
$$\exp\left(-t^k\right)\mathrm dt = \exp\left(-t^k\right)\ t \frac{\mathrm dt}{t} = \frac{1}{k}\exp(-z)\ z^{1/k}\frac{\mathrm dz}{z}.$$
Its integral (without the factor of $1/k$) defines the lower incomplete Gamma function with shape parameter $1/k,$ often written $\gamma.$  Thus, since the upper limit of the integral is $z^{1/k} = t = u,$ we must evaluate this integral from $z = 0$ to $z = u^k:$

$$\int_0^u \exp\left(-t^k\right)\mathrm dt = \frac{1}{k}\gamma\left(\frac{1}{k}, u^k\right).$$

(Evidently $k$ must be positive.)  The rest is straightforward algebra.
As far as implementing $\gamma$ goes, it's available in most (if not all) statistical computing platforms because when $X$ has a $\Gamma(1/k)$ distribution and $x \gt 0,$
$$\Gamma(1/k)\Pr(X\le x) = \int_0^x z^{1/k}\exp(-z)\frac{\mathrm dz}{z} = \gamma\left(\frac{1}{k}, x\right).$$
Compute the factors on the left to obtain the value of $\gamma.$

As a check and a working example, in R we may numerically integrate $\exp(-t^k)$
f <- Vectorize(function(u, k) {integrate(\(t) exp(-t^k), 0, u)$value}, "u")

and compare that to this solution where $\Gamma$ is implemented with lgamma (the logarithm of the Gamma function, used to avoid numeric overflows for small $k$) and $\Pr(X\le x)$ is implemented with pgamma:
g <- function(u, k) exp(-log(k) + lgamma(1/k) + pgamma(u^k, 1/k, log.p = TRUE))

Graphs of the functions f and g coincide (at least for non-extreme values of $k;$ for very small values the numeric integration becomes problematic).
