docs: def gpinv(p, k, sigma):
      """Inverse Generalised Pareto distribution function."""
      x = np.empty(p.shape)
      x.fill(np.nan)
      if sigma <= 0:
          return x
      ok = (p > 0) & (p < 1)
      if np.all(ok):
          if np.abs(k) < np.finfo(float).eps:
              np.negative(p, out=x)
              np.log1p(x, out=x)
              np.negative(x, out=x)
          else:
              np.negative(p, out=x)
              np.log1p(x, out=x)
              x *= -k
              np.expm1(x, out=x)
              x /= k
          x *= sigma
      else:
          if np.abs(k) < np.finfo(float).eps:
              # x[ok] = - np.log1p(-p[ok])
              temp = p[ok]
              np.negative(temp, out=temp)
              np.log1p(temp, out=temp)
              np.negative(temp, out=temp)
              x[ok] = temp
          else:
              # x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
              temp = p[ok]
              np.negative(temp, out=temp)
              np.log1p(temp, out=temp)
              temp *= -k
              np.expm1(temp, out=temp)
              temp /= k
              x[ok] = temp
          x *= sigma
          x[p == 0] = 0
          if k >= 0:
              x[p == 1] = np.inf
          else:
              x[p == 1] = -sigma / k
      return x