docs: public static double intervalCost(double x1, double y0, double y1)
    {
      if (x1 == 0 || y1 == y0) {
        return 0;
      }
  
      // cost(X, Y) = cost(Y, X), so we swap X and Y to
      // have x_0 <= y_0 and simplify the calculations below
      if (y0 < 0) {
        // swap X and Y
        double tmp = x1;
        x1 = y1 - y0;
        y1 = tmp - y0;
        y0 = -y0;
      }
  
      // since x_0 <= y_0, Y must overlap X if y_0 < x_1
      if (y0 < x1) {
        /**
         * We have two possible cases of overlap:
         *
         * X  = [ A )[ B )[ C )   or  [ A )[ B )
         * Y  =      [   )                 [   )[ C )
         *
         * A is empty if y0 = 0
         * C is empty if y1 = x1
         *
         * cost(X, Y) = cost(A, Y) + cost(B, C) + cost(B, B)
         *
         * cost(A, Y) and cost(B, C) can be calculated using the non-overlapping case,
         * which reduces the overlapping case to computing
         *
         * cost(B, B) = \int_0^{\beta} \int_{0}^{\beta} e^{-|x-y|}dxdy
         *            = 2 \cdot (\beta + e^{-\beta} - 1)
         *
         *            where \beta is the length of interval B
         *
         */
        final double beta;  // b1 - y0, length of interval B
        final double gamma; // c1 - y0, length of interval C
        if (y1 <= x1) {
          beta = y1 - y0;
          gamma = x1 - y0;
        } else {
          beta = x1 - y0;
          gamma = y1 - y0;
        }
        //noinspection SuspiciousNameCombination
        return intervalCost(y0, y0, y1) + // cost(A, Y)
               intervalCost(beta, beta, gamma) + // cost(B, C)
               2 * (beta + FastMath.exp(-beta) - 1); // cost(B, B)
      } else {
        /**
         * In the case where there is no overlap:
         *
         * Given that x_0 <= y_0,
         * then x <= y must be true for all x in [x_0, x_1] and y in [y_0, y_1).
         *
         * therefore,
         *
         * cost(X, Y) = \int_0^{x_1} \int_{y_0}^{y_1} e^{-|x-y|} dxdy
         *            = \int_0^{x_1} \int_{y_0}^{y_1} e^{x-y} dxdy
         *            = (e^{-y_1} - e^{-y_0}) - (e^{x_1-y_1} - e^{x_1-y_0})
         *
         * Note, this expression could be further reduced by factoring out (e^{x_1} - 1),
         * but we prefer to keep the smaller values x_1 - y_0 and x_1 - y_1 in the exponent
         * to avoid numerical overflow caused by calculating e^{x_1}
         */
        final double exy0 = FastMath.exp(x1 - y0);
        final double exy1 = FastMath.exp(x1 - y1);
        final double ey0 = FastMath.exp(0f - y0);
        final double ey1 = FastMath.exp(0f - y1);
  
        return (ey1 - ey0) - (exy1 - exy0);
      }
    }