\section{Introduction}\label{s:intro}
Stable L\'evy processes and densities appear in various fields of natural sciences, engineering and finance.
 See, e.g.,
\cite{Chandrasekhar,Chavanis09,Mand1,Mand2,Zolotarev,SamorodnitskyTaqqu94,Weron96,MitRachev,Nolan97,Nolan98,Nolan03,NikiasShao95,SignalProc10,Saenko16,AmentONeil18} and the bibliographies therein. Efficient calculation of (cumulative) probability distribution functions  (cpdf and pdf) and expectations in stable L\'evy models are non-trivial even in the one-dimensional case (1D case).  The slow decay of the transition kernels at infinity and very large derivatives of the kernels near the peak make accurate calculations in the state space difficult. If the Fourier transform technique is used, standard numerical quadratures face serious difficulties because of  the oscillation and slow decay of the integrands at infinity. In particular,
standard tools such as the fast Fourier transform produce very large errors unless the index $\alpha$ of the stable L\'evy process is  close to 2, hence, the process is close to the Brownian motion (BM).
The methods in op.cit. are very difficult to generalize to subordinated stable L\'evy processes and mixtures of stable distributions,
which arise in applications to signal processing  \cite{SignalProc10}. See \cite{ConfAccelerationStable} for details.
 Generalizations to distributions of extrema of a L\'evy processes
and joint distributions, which we consider
in the present paper, is essentially impossible. It is necessary to evaluate
  double and triple integrals, hence, the detailed analysis of the integrands used in the 1D case does not seem feasible.

In the present paper, as in \cite{ConfAccelerationStable} and  \cite{EfficientLevyExtremum}, we use the conformal deformation technique.
The conformal deformations in  \cite{EfficientLevyExtremum} are applicable to L\'evy processes with the characteristic exponents admitting analytic continuation to a strip; in  \cite{ConfAccelerationStable},  variations of the conformal deformation technique suitable  for evaluation pdf and cpdf of stable
L\'evy distribution on ${\mathbb R}$ are developed.  The new elements of the present
paper are efficient methods for evaluation of the Wiener-Hopf factors and 2D-3D integrals in the formula for
expectations  $V(f,T;x_1,x_2)={\mathbb E}[f(x_1+X_T, \max\{x_2, x_1+{\bar X}_T\})]$ of functions of a stable L\'evy process  $X$  and its
supremum process ${\bar X}$ both starting at zero. The generalization to the case when the infimum process replaces supremum process is by symmetry. The formula is derived in \cite{EfficientLevyExtremum} for any L\'evy process; numerical realizations have to be modified for the case of stable L\'evy processes.
In the nutshell, the main idea is to deform the contours  of integration using  appropriate conformal maps, which must be in a certain agreement. The corresponding changes of variables and simplified trapezoid rule allow us to satisfy a small error tolerance, at a small CPU cost.   The efficiency of the simplified
trapezoid rule  stems from the fact that the discretization error of the  infinite trapezoid rule is  $O(\exp[-2\pi d/\zeta])$, where $d$ is the half-width of the strip of analyticity around the line of integration and $\zeta$ is the step size. As in \cite{ConfAccelerationStable}, certain
preliminary regularizations of the integrals are needed. In many cases, the method is faster and more accurate than the saddle point method and methods based on the reduction to an appropriate cut in the complex plane (see, e.g., \cite{Fedoryuk,NG-MBS,IAC}).
The schemes that we design use  the general analytical properties of the characteristic exponent and can be easily modified to more involved cases of mixtures of stable distributions and subordinated
stable L\'evy processes, and  to solution of several types of boundary problems for fractional
differential equations; generalizations to higher dimensions are possible.

Special cases which we consider in more detail in the paper, are: (1) $f$ depends on the first argument only; (2) $f(x_1,x_2)={\bf 1}_{(\infty,a]}(x_2)$;
(3) $f(x_1, x_2)={\bf 1}_{(-\infty, a_1]}(x_1){\bf 1}_{(-\infty, a_2]}(x_2)$, where $x_1\le x_2\le a_2, a_1\le a_2$; (4)
$f(x_1,x_2)=(\beta x_1-x_2)_+e^{-\lambda x_2}$, where $\beta>1$ and $\lambda>0$.
Case (1), which we consider in Section \ref{s:1D}, is a straightforward generalization of the case
of the cpdf of a stable L\'evy process considered in \cite{ConfAccelerationStable}. Numerical schemes outlined for this case are used as  blocks in
more complicated situations.
First, as in \cite{ConfAccelerationStable}, we calculate
the expectation ${\mathbb E}[f(x+X_T)]$ using the Fourier transform technique.
We reduce the integral in the Fourier inversion formula to an integral over
over a half-axis, rotate the half-axis to $e^{i\omega}{\mathbb R}_+$, where the choice of $\omega\in (-\pi/2,\pi/2)$ is determined by  the domains of analyticity and oscillation of the characteristic function of the process and Fourier transform ${\hat f}$. Then we make the change of variables $\xi=e^{i\omega+y}$ and apply the simplified trapezoid rule. The exponential rate of decay of the error of the infinite trapezoid rule as a function of $1/\zeta$ allows one to easily satisfy a small error tolerance unless $d$, the half-width of the strip of analyticity, is very small. The $d$ is very small in the case of strongly
 asymmetric L\'evy processes of index $\alpha\in (0,2), \alpha\neq 1$; for asymmetric L\'evy processes of index 1, the rotation is impossible for
  $x$ in a certain semi-infinite interval unless high precision arithmetic is used. To alleviate these difficulties, in \cite{ConfAccelerationStable}, we used two additional families of conformal deformations; the CPU time significantly increased. In the present paper, we calculate of
${\mathbb E}[f(x+X_T)]$ using the Laplace transform w.r.t. $T$ and Fourier transform w.r.t. $x$. The resulting 2D integral is calculated using either the Gaver-Wynn Rho algorithm (GWR algorithm) for the Laplace inversion or the conformal deformationof the line of integration in the Bromwich integral  used in \cite{Contrarian,EfficientLevyExtremum}.  The deformation of the same type, namely, the sinh-acceleration, was introduced in \cite{Sinh} in the context of evaluation of special functions and used in \cite{SINHregular} to price European options in wide classes of L\'evy and affine models. The inner integral (the inverse Fourier transform) can be calculated using the rotation of the half-axis of integration and the exponential change of variable $\xi=e^{i\omega+y}$. The half-width of the strip of analyticity $d=\pi/8$ in the new coordinate $y$ is not small, hence, the number of terms in the simplified trapezoid rule for the inner integral is small,
which compensates for the drawback of calculations for several dozen of points in the Laplace inversion formula.

Case (2) is the probability distribution of the supremum process. We calculate the Laplace transform w.r.t. $T$ using
the Wiener-Hopf factorization technique, and evaluate the Bromwich integral. The result is a 2D integral, the integrand
being expressed in terms of the Wiener-Hopf factors. An efficient numerical procedure uses the conformal deformations
in the formula for the Wiener-Hopf factor and in the 2D integral; all deformations must be in a certain agreement.
A simpler but less accurate version uses GWR algorithm to evaluate the Bromwich integral.
The elements of the Wiener-Hopf factorization technique used in the paper are in Section \ref{ss:WHF_gen}.
In particular, we introduce a decomposition of the Wiener-Hopf factors, which is necessary for  efficient realizations of the
general formula for $V(f,T;x_1,x_2)$. The formula is derived in \cite{EfficientLevyExtremum}, for any L\'evy process, and
an efficient numerical realization is developed for L\'evy processes with exponentially decaying tails of the L\'evy density.
The decomposition allows us to repeat the constructions in \cite{EfficientLevyExtremum} almost verbatim.
The main difference is that in \cite{EfficientLevyExtremum}, the sinh-changes of variables on ${\mathbb R}$, of the form
$\xi=i\omega_1+b\sinh(i\omega+y)$,
where $\omega_1\in {\mathbb R}$, $b>0$ and $\omega\in (0,\pi/2)$ are used; in the present paper,
the exponential changes of variables on rays $e^{i\omega}{\mathbb R}_+$ are used: $\xi=e^{i\omega +y}$.
As it was demonstrated in  \cite{ConfAccelerationStable}, typically, the sinh-change of variables and exponential change of variables are the most efficient ones but, in some cases, less efficient families are more efficient or even indispensable.

Case (3) is the joint cpdf of the stable L\'evy process and its extremum. If GWR algorithm is used, then, for each $q>0$ in GWR, we need
to evaluate 2D integral; if the sinh-acceleration is applied to the Browmwich integral, then 3D integral. The latter version allows us to
achieve much higher precision than the former but it is possible only if the index $\alpha>1$ or $\alpha=1$ and the jump part is symmetric or $\alpha\in (0,1)$ and $\mu=0$. Note that in Case (3), $V(f,T;x_1,x_2)$ is independent of $x_2$.
In Case (4), $V(f,T;x_1,x_2)$ depends on both $x_1$ and $x_2$. Furthermore, Case (4) is used to illustrate how the general theorem formulated and proved for bounded measurable functions can be used
when $f$ is unbounded: introduce the dampening factor, derive the integral representation and prove that it is possible to pass to the limit $\lambda\to 0$. In Case (4), the limit is finite iff $\alpha>1$.

The integral representations are derived in Section \ref{int_repr}. In Section \ref{numer_real}, efficient numerical realizations are
discussed, the detailed algorithms for Case (3) are formulated, and  numerical examples are  produced.
Section \ref{concl} concludes; technical details and tables with numerical examples are relegated to appendices \ref{s:techn} and
\ref{s:tables}, respectively.

\section{Evaluation of expectations $E[f(x+X_T)]$. Conformal deformations technique}\label{s:1D}

\subsection{Stable L\'evy processes: main notation}\label{ss:stable_notation}
Let $X$ be a one-dimensional stable L\'evy process on the filtered probability space $({\Omega}, {\mathcal F}, \{{\mathcal F}_t\}_{t\ge 0}, {\mathbb P})$
satisfying the usual conditions, and let ${\mathbb E}$ be the expectation operator under ${\mathbb P}$.
We use the same parametrization of stable L\'evy processes on ${\mathbb R}$ and representations of the characteristic exponent as in \cite{ConfAccelerationStable}.
Let $c_\pm\ge 0$, $c_++c_->0$, $\alpha\in (0,2)$, and let $X$ be a stable L\'evy process of index $\alpha$, with the L\'evy density
\[
F(dx)=c_+ x^{-\alpha-1}{\bf 1}_{(0,+\infty)}(x)dx+c_-|x|^{-\alpha-1}{\bf 1}_{(-\infty,0)}(x)dx.
\]
If $\alpha\neq 1$,  the characteristic exponent is of the form
\begin{equation}\label{stdr}
\psi_{st}(\xi)=-i\mu\xi+\psi^0_{st}(\xi),
\end{equation}
where $\mu\in{\mathbb R}$,
\begin{eqnarray}\label{SLnuneq10}
\psi^0_{st}(\xi)&=&{C_+} |\xi|^\alpha{\bf 1}_{(0,+\infty)}(\xi)+ {C_-}|\xi|^\alpha{\bf 1}_{(-\infty,0)}(\xi),\ \xi\in {\mathbb R},
\\\label{Cp}
{C_+}&=&{C_+}(\alpha,{c_+},{c_-})=-{c_+}\Gamma(-\alpha)e^{-i\pi\alpha/2}-{c_-}\Gamma(-\alpha)e^{i\pi\alpha/2},
\end{eqnarray}
and ${C_-}=\overline{{C_+}}$.
If $\alpha=1$, then the formula for $\psi^0_{st}$ is more involved:
  \begin{equation}\label{SLnu1}
\psi^0_{st}(\xi)=\sigma|\xi|(1+i(2\beta/\pi)\,\mathrm{sign}\,\xi\ln|\xi|),\ \xi\in {\mathbb R},
\end{equation}
where $\sigma=({c_+}+{c_-})\pi/2$, $\beta=({c_+}-{c_-})/({c_+}+{c_-})$. This is a version of Zolotarev's  parametrizations \cite{Zolotarev}
for stable processes of index 1.  In the symmetric case $c=c_+=c_-$, \eq{SLnu1} reduces to
  \eq{SLnuneq10} with ${C_+}={C_-}=c\pi$ and $\alpha=1$.
\begin{prop}\label{prop:anal_psi_st}
\begin{enumerate}[(a)]
\item
$\psi^0_{st}(\xi)$ admits analytic continuation from $(0,+\infty)$ to the right half-plane (and to an appropriate Riemann surface).
\item
$\psi^0_{st}(\xi)$ admits analytic continuation from $(-\infty,0)$ to the left half-plane (and to an appropriate Riemann surface).
\end{enumerate}
\end{prop}
\begin{proof} It suffices to note that  $\xi^\alpha=\exp[\alpha\ln \xi]$ for $\xi\in{\mathbb C}\setminus (-\infty,0]$.
\end{proof}
We use  the notation
 ${\mathcal C}_{\gamma_-,\gamma_+}=\{e^{i\varphi}\rho\ |\ \rho> 0, \varphi\in (\gamma_-,\gamma_+)\cup (\pi-\gamma_+,\pi-\gamma_-)\}$,
 ${\mathcal C}_{\gamma}=\{e^{i\varphi}\rho\ |\ \rho> 0, \varphi\in (-\gamma,\gamma)\}$. The proof of the following lemma is by inspection.
\begin{lem}\label{lem:asymp_psi_st}  For any $\gamma\in (0,\pi/2)$, as $({\mathcal C}_\gamma\ni)\xi\to \infty$,
\begin{enumerate}[(a)]
\item
if $\alpha\in (1,2)$ or $\alpha\in(0,1)$ and $\mu=0$,
\begin{equation}\label{aspsimain}
\psi_{st}(\xi)={C_+} \xi^\alpha
(1+|\mu|O(|\xi|^{1-\alpha})),
\end{equation}
and
\begin{equation}\label{aspsimainRe}
\Re\psi_{st}(\xi)=|{C_+}|\cos(\varphi_0+\varphi) |\xi|^\alpha
(1+|\mu|O(|\xi|^{1-\alpha})),
\end{equation}
where $\varphi_0=\mathrm{arg}\,{C_+}$ and $\varphi=\mathrm{arg}\,\xi$;
\item
if $\alpha\in (0,1)$ and $\mu\neq 0$,
\begin{equation}\label{aspsi01}
\psi_{st}(\xi)=-i\mu \xi(1+O(|\xi|^{\alpha-1})),
\end{equation}
and
\begin{equation}\label{aspsi01Re}
\Re\psi_{st}(\xi)=\mu \sin(\varphi)|\xi|(1+O(|\xi|^{\alpha-1}));
\end{equation}
\item
if $\alpha=1$ and $c_+=c_-=c$, then
\begin{equation}\label{aspsi1eq}
\psi_{st}(\xi)=(c\pi-i\mu) \xi,
\end{equation}
and
\begin{equation}\label{aspsi1eqRe}
\Re\psi_{st}(\xi)=\sqrt{(c\pi)^2+\mu^2}\cos(-\arctan(\mu/(c\pi))+\varphi) |\xi|;
\end{equation}
\item
if $\alpha=1$ and $c_+\neq c_-$, then
\begin{equation}\label{aspsi1neq}
\psi_{st}(\xi)=((c_++c_-)\pi/2-i\mu)\xi+i(c_+-c_-)\xi\ln\xi,
\end{equation}
and
\begin{equation}\label{aspsi1neqRe}
\Re\psi_{st}(\xi)=(c_--c_+)\sin(\varphi)|\xi|\ln|\xi|+O(|\xi|).
\end{equation}
Furthermore, if $c_+>c_-$, then there exists $R>0$ s.t. for   $\xi\in \{\xi \in {\mathcal C}_{-\pi/2,0}\ | \ |\xi|\ge R\}$,
\begin{equation}\label{aspsi1neqRe1p}
\Re\psi_{st}(\xi)\ge \sqrt{\sigma^2+\mu^2}\cos(\varphi-\arctan(\mu/\sigma)|\xi|,
\end{equation}
and if $c_+<c_-$, then \eq{aspsi1neqRe1p} holds for   $\xi\in \{\xi \in{\mathcal C}_{0,\pi/2}\ |\ |\xi|\ge R\}$.
\end{enumerate}
\end{lem}

\begin{lem}\label{SLroots}
For
any $\xi\in {\mathbb C}\setminus i{\mathbb R}$ and $q>0$, $q+\psi_{st}(\xi)\neq 0$.
\end{lem}
\begin{proof}
This property is proved in \cite{EfficientAmenable} for a wide class of L\'evy processes
called {\em Stieltjes-L\'evy processes} (SL-processes). It is also proved that stable L\'evy processes are SL-processes.
\end{proof}

If $\alpha\ge 1$ or $\alpha\in (0,1)$ and $\mu=0$, set $\bar\alpha=\alpha$. If $\alpha\in (0,1)$ and $\mu\neq 0$, set $\bar\alpha=1$.

The following lemma is immediate from Lemma \ref{lem:asymp_psi_st}.
\begin{lem}\label{lem:psi_st_bounds}
\begin{enumerate}[(a)]
\item
For any $\gamma\in (0,\pi/2)$ and any $q>0$, there exists $C>0$ such that
\begin{equation}\label{bound_psi_inv}
q/|q+\psi_{st}(\xi)|\le C(1+|\xi|)^{-\bar\alpha},\ \xi\in {\mathcal C}_\gamma.
\end{equation}
\item
There exist $-\pi/2<\gamma_-\le 0\le \gamma_+<\pi/2$ and $c, R>0$ such that $\gamma_-<\gamma_+$ and
\begin{equation}\label{lim_bound_Re_psi}
\Re\psi^0_{st}(\xi)\ge c|\xi|^{\alpha},\ \xi\in {\mathcal C}_{\gamma_-,\gamma_+}, |\xi|\ge R.
\end{equation}
\item
\begin{enumerate}[(i)]
\item
If $\alpha\in (0,1)\cup (1,2)$ or $\alpha=1$ and $c_+=c_-$, then $\gamma_-<0<\gamma_+$.
\item
If $\alpha=1$ and $c_+>c_-$ (resp., $c_+<c_-$), then $\gamma_+=0$ (resp., $\gamma_-=0$).
\end{enumerate}
\end{enumerate}
\end{lem}
\begin{rem}\label{rem:Re vs upper bound}{\rm
The line of integration in the Fourier inversion formula in the next subsection can be deformed into a cone where \eq{lim_bound_Re_psi}
holds, and the lines of integration in the formulas in the following sections - into a cone where
\eq{bound_psi_inv} holds. This fact implies that if GWR algorithm is used to evaluate the Bromwich integral, the restrictions on the type
of deformations in the formulas for the Laplace transform ${\tilde V}(f;q,x_1,x_2)$ are milder than when the Fourier inversion is used
to evaluate ${\mathbb E}[f(x+X_T)]$. One can use the Laplace-Fourier inversion to evaluate ${\mathbb E}[f(x+X_T)]$ which means that
we  replace a 1D integral
with 2D. We will analyze in a separate publication the advantages of this  approach in applications to evaluation of pdf and cpdf of strongly asymmetric stable L\'evy processes, of index $\alpha=1$ or close to 1 especially.
}
\end{rem}

\subsection{Evaluation of $V(T,x)={\mathbb E}[f(x+X_T)]$}\label{ss:1D} Representing $f:{\mathbb R}\to {\mathbb C}$ as a linear combination of positive functions supported on either  $(-\infty,0]$ or  $[0,+\infty)$, we may restrict ourselves to non-negative $f$ supported on $[0,+\infty)$ only.

\vskip0.1cm
\noindent
{\sc Assumption $f$.} The following conditions are satisfied:
\begin{enumerate}[(I)]
\item
${\mathrm{supp}}\, f\in [0,+\infty)$, and,
 for any $\lambda>0$, $f(\lambda,\cdot):=f(\cdot)e^{-\lambda\cdot}\in L_1({\mathbb R};
 {\mathbb R}_+)$;
 \item  there exist $a\in {\mathbb R}$ and $\gamma_0\in (\pi/2,\pi)$ such that the Fourier transform ${\hat f}$ is of the form ${\hat f}(\xi)=e^{-ia\xi}{\hat f}_0(\xi)$,
where ${\hat f}_0(\xi)$ admits analytic continuation to  $-i{\mathcal C}_{\gamma_{0}}$;

\item there exist  $C>0$, $\alpha_0<\alpha$ and $\alpha_\infty\in {\mathbb R}$, such that
\begin{equation}\label{boundfD1}
|{\hat f}_0(\xi)|\le C(|\xi|^{-\alpha_0-1}{\bf 1}_{(0,1]}(|\xi|)+|\xi|^{\alpha_\infty-1}{\bf 1}_{[1,+\infty)}(|\xi|)),\ \xi\in -i{\mathcal C}_{\gamma_{0}}.
\end{equation}
\end{enumerate}
Set $V(\lambda;T,x)={\mathbb E}[f(x+X_T)e^{-\lambda (x+X_T)}]$.
We have $V(T,x)=\lim_{\lambda\downarrow 0}V(\lambda;T,x)$. Let $Y$ be the BM of variance $2|{C_+}|$, with drift $\mu$, independent of $X$. Since ${\hat f}_0(\lambda;\xi)$
admits analytic continuation to $\lambda-i{\mathcal C}_{\gamma_{0}}$ and is polynomially bounded at infinity, we can use the Fourier transform technique.
Set $x'=x+\mu T-a$. We have
\begin{eqnarray*}
V(\lambda;T,x)&=&\frac{1}{2\pi}\int_{\mathbb R} e^{ix'\xi-T\psi^0_{st}(\xi)}{\hat f}_0(\lambda;\xi)d\xi\\
&=&\frac{1}{2\pi}\int_{\mathbb R} e^{ix'\xi}\left(e^{-T\psi^0_{st}(\xi)}-e^{-T|{C_+}|\xi^2}\right){\hat f}_0(\lambda;\xi)d\xi +{\mathbb E}[f(\lambda;x+Y_T)].
\end{eqnarray*}
As $\xi\to 0$, $e^{-T\psi_{st}(\xi)}-e^{-T|{C_+}|\xi^2}=O(|\xi|^\alpha)$, and ${\hat f}(\lambda,\xi)=O(|\xi|^{-\alpha_0-1})$, uniformly in $\lambda$. Hence, we can pass to the limit $\lambda\to 0$, and obtain $V(T,x)=V_1(T,x)+{\mathbb E}[f(x+Y_T)]$, where
\begin{equation}\label{eq:VT}
V_1(T,x)=\frac{1}{2\pi}\int_{\mathbb R} e^{ix'\xi}\left(e^{-T\psi^0_{st}(\xi)}-e^{-T|{C_+}|\xi^2}\right){\hat f}_0(\xi)d\xi.
\end{equation}
For the BM model, efficient numerical procedures are well-known (see, e.g., \cite{SINHregular} and the bibliography therein),
hence, it suffices to evaluate  $V_1(T,x)$.
Since $\overline{e^{ix'\xi}}=e^{-ix'\bar\xi}$, $\overline{\psi^0_{st}(\xi)}=\psi^0_{st}(-\bar\xi)$ and $\overline{{\hat f}_0(\xi)}={\hat f}_0(-\bar\xi)$,
we have
\begin{equation}\label{V1Tx}
V_1(T,x)=\frac{1}{\pi}\Re\int_0^{+\infty} e^{ix'\xi}\left(e^{-T\psi^0_{st}(\xi)}-e^{-T|{C_+}|\xi^2}\right){\hat f}_0(\xi)d\xi.
\end{equation}
Assume that either $\alpha\neq 1$ or $\alpha=1$ and $c_+=c_-$. Then there exist $\gamma_-<0<\gamma_+$ such that \eq{lim_bound_Re_psi} holds.
We redefine $\gamma_+:=\min\{\gamma_+,\gamma_0-\pi/2\}$, and, if $x'>0$,  take $\omega\in (0,\gamma_+/2)$, deform the line of integration into $e^{i\omega}{\mathbb R}_+$, and change the variable $\xi=e^{i\omega+y}$:
\begin{equation}\label{V1Tx2}
V_1(T,x)=\frac{1}{\pi}\Re\int_{-\infty}^{+\infty} e^{ix'\xi(y)}\left(e^{-T\psi^0_{st}(\xi(y))}-e^{-T|{C_+}|\xi(y)^2}\right){\hat f}_0(\xi(y))e^{i\omega+y}dy.
\end{equation}
If $x'<0$ (resp., $x'=0$), we use \eq{V1Tx2} with $\omega\in (\gamma_-/2,0)$ (resp., $\omega=(\gamma_++\gamma_-)/2$).  We evaluate
the integral on the RHS of \eq{V1Tx2} using the infinite trapezoid rule.

The integrand on the RHS of \eq{V1Tx2}, denote it  $g$, is analytic in the strip
$S_{(-d,d)}:=\{\xi\ | \Im\xi\in (-d,d)\}$, for any $d\in (0,\gamma_+/2)$ (if $x'>0$, $d\in (0,-\gamma_-/2)$; if $x'=0$,
$d\in (0, (\gamma_+-\gamma_-)/2)$) and decays at infinity sufficiently fast so that
$\lim_{A\to \pm\infty}\int_{-d}^d |g(i a+A)|da=0,$
and
\begin{equation}\label{Hnorm}
H(g,d):=\|g\|_{H^1(S_{(-d,d)})}:=\lim_{a\downarrow -d}\int_{\mathbb R}|g(i a+ y)|dy+\lim_{a\uparrow d}\int_{\mathbb R}|g(i a+y)|dy<\infty
\end{equation}
is finite. We write $g\in H^1(S_{(-d,d)})$. The integral
$I=\int_{\mathbb R} g(\xi)d\xi$
can be evaluated using the infinite trapezoid rule
\begin{equation}\label{inftrap}
I\approx \zeta\sum_{j\in {\mathbb Z}} g(j\zeta),
\end{equation}
where $\zeta>0$.
The following key lemma is proved in \cite{stenger-book} using the heavy machinery of sinc-functions. A simple proof can be found in
\cite{paraHeston}.
\begin{lem}[\cite{stenger-book}, Thm.3.2.1] Let $g\in H^1(S_{(-d,d)})$.
The error of the infinite trapezoid rule \eq{inftrap} admits an upper bound
\begin{equation}\label{Err_inf_trap}
{\rm Err}_{\rm disc}\le H(g,d)\frac{\exp[-2\pi d/\zeta]}{1-\exp[-2\pi d/\zeta]}.
\end{equation}
\end{lem}
Once
an approximately bound for $H(g,d)$ is derived, it becomes possible to satisfy the desired error tolerance
with a good accuracy.  The infinite sum is truncated
\begin{equation}\label{inftrap_trunc}
I\approx \zeta\sum_{j=-N_-}^{N_+} g(j\zeta).
\end{equation}
As $j\to +\infty$, $g(j\zeta)$ exhibits a super-exponential decay, hence, a fairly small $N_+$ can be chosen to satisfy even a very small error tolerance $\epsilon$. In a neighborhood of  $-\infty$, $g(j\zeta)$ admits a bound via $C\exp[(\alpha-\alpha_0)j\zeta]$, where an approximate upper bound for $C$ can be easily derived. Hence, we choose $N_-=\mathrm{ceil}\,(\ln(C/\epsilon)/(\alpha-\alpha_0))$. If ${\hat f}_0$ admits an asymptotic expansion
as $\xi\to 0$, the asymptotic expansions of $e^{ix'\xi}, e^{-T\psi^0_{st}(\xi)}$ and $e^{-T|{C_+}|\xi^2}$ as $\xi\to 0$ can be used
to calculate an asymptotic expansion
\[
g(j\zeta)=\sum_{k=1}^{M-1} c_k e^{\beta_k j\zeta}+ O(e^{\beta_M j\zeta}),
\]
where $c_k$ and $\beta_k>0$ depend on ${\hat f}_0, \alpha, c_\pm$ and $T$.  After that we calculate explicitly the sums
$\sum_{j=-\infty}^{-N_--1}e^{\beta_k j\zeta}$, and choose significantly smaller $N_-=\mathrm{ceil}\,(\ln(C/\epsilon)/\beta_M)$.
See \cite{ConfAccelerationStable} for details.
\begin{rem}\label{rem:other families}{\rm The complexity of the scheme based on the exponential change of variables admits a bound via
 $C(\ln(1/\epsilon))^2$. As we demonstrated in \cite{ConfAccelerationStable}, if either $\gamma_+$ or $-\gamma_-$ are very small,
then the half-width of the strip of analyticity in the $y$-coordinate is very small, $\zeta$ is very small, and the constant $C$ is extremely
large. In the result, even for a fairly small error tolerance, e.g., E-15, it may be
advantageous to use a polynomial change of variables. The complexity of the corresponding scheme is of the order of
$C\epsilon^{-m}\ln(1/\epsilon)$, where $m>0$. Finally, if $\alpha=1$ and $c_+\neq c_-$, then, depending on the sign of $x'(c_+-c_-)$,
neither the exponential change of variables nor polynomial one can be used, and only a sub-polynomial increase of the rate of the decay of
the integrand can be achieved. The number of terms and CPU time are much larger.

}
\end{rem}

\subsection{An alternative evaluation of $V(T,x)={\mathbb E}[f(x+X_T)]$ and GWR algorithm}\label{ss:1D_alt} The trick in this subsection allows us to use the exponential change of variables in all cases, at the cost of an additional integration (Laplace inversion).
Let $f$  satisfy the same conditions as above, with an additional restriction on $\alpha_\infty$ in \eq{boundfD1}. Namely,
we assume that $\alpha_\infty<\bar\alpha$, where
$\bar\alpha$ is as in Lemma \ref{lem:psi_st_bounds}. Somewhat counterintuitively, we represent the 1D integral \eq{eq:VT}
as a 2D integral
\begin{equation}\label{eq:V1TBr}
V(T,x)=\frac{1}{2\pi i}\int_{\Re q=\sigma}dq\, e^{qT}\frac{1}{2\pi}\int_{\mathbb R} e^{ix'\xi}\left(\frac{1}{q+\psi^0_{st}(\xi)}-\frac{1}{q+|{C_+}|\xi^2}\right){\hat f}_0(\xi)d\xi.
\end{equation}
If the outer integral (Bromwich integral) is evaluated using GWR algorithm, it suffices to evaluate the inner integral for a moderate number of values of $q>0$. Similarly to \eq{V1Tx}, for each $q$, we simplify the inner integral
\begin{equation}\label{tVq1D}
{\tilde V}(q,x) =\frac{1}{\pi}\Re \int_0^{+\infty}e^{ix'\xi}\left(\frac{1}{q+\psi^0_{st}(\xi)}-\frac{1}{q+|{C_+}|\xi^2}\right){\hat f}_0(\xi)d\xi.
\end{equation}
On the strength \eq{bound_psi_inv}, we can rotate the ray of integration to $e^{i\omega}{\mathbb R}_+$ using any $\omega\in (0,\pi/4)$ (resp., $\omega\in (-\pi/4,0)$;
$\omega\in (-\pi/4,\pi/4)$) if $x'>0$ (resp., $x'<0$; $x'=0$), and the choice $\omega=\pi/8$ (resp., $\omega=-\pi/8$; $\omega=0$) is optimal.
After the change of variables $\xi=e^{i\omega+y}$, we apply the simplified trapezoid rule. The half-width of the strip of analyticity
of the integrand in the $y$-coordinate, which can be used to derive the recommendation for the choice of $\zeta$, is
an arbitrary $d\in (0,\pi/8)$, if $x'\neq 0$, and $d\in (0,\pi/4)$, if $x'=0$. Thus, in all cases, we can use $\zeta$ of the order of
$(\pi^2/4)/\ln(1/\epsilon)$, whereas in Section \ref{ss:1D}, if either $\gamma_+$ or $-\gamma_-$ is very small, $\zeta$ must be very small as well.
Thus,   the number of terms in the simplified trapezoid rule for the evaluation of \eq{tVq1D} is many times smaller than
the number of terms in Section \ref{ss:1D}. Although we need to evaluate ${\tilde V}(q,x)$ for all $q$'s used in the algorithm chosen for the Laplace inversion, the number of $q$'s needed to satisfy even a very small error tolerance is measured in dozens, and
the calculations can be easily parallelized. Set
\[
g(y)=e^{ix'\xi(y)}\left(\frac{1}{q+\psi^0_{st}(\xi(y))}-\frac{1}{q+|{C_+}|\xi(y)^2}\right){\hat f}_0(\xi(y))e^{i\omega+y}.
\]
The truncation parameters $N_\pm$ in the simplified trapezoid rule
\begin{equation}\label{tVqsimp}
{\tilde V}(q,x)\approx \frac{\zeta}{\pi}\Re\sum_{j=-N_-}^{N_+}g(j\zeta)
\end{equation}
are chosen taking into account the exponential asymptotics of $g(y)$ as $y\to\pm \infty$. Contrary to Section \ref{ss:1D}, we can explicitly calculate several terms of of the truncated parts of the infinite sum in both neighborhoods of $j=-\infty$ and $j=+\infty$, and
choose  smaller $N_\pm$.

\subsection{ Sinh-acceleration in the Bromwich integral}\label{ss:1D_alt_Br}
Instead of GWR-algorithm, which typically requires high precision arithmetic to achieve the precision of the order better than E-7-E-8,
we can use the sinh-acceleration in the Bromwich integral (the outer integral on the RHS of \eq{eq:V1TBr}) and the exponential acceleration in the inner integral. This is possible under an additional conditions on the parameters of the process.
\begin{lem}\label{q_xi_cone}
\begin{enumerate}[(a)]
\item
Let either $\alpha\in (1,2)$ or $\alpha=1$ and $c_+=c_-$ or $\alpha\in (0,1)$ and $\mu=0$.

Then
there exist $\sigma>0, \gamma_0\in (0,\pi/2)$ and $\gamma_-<0<\gamma_+$ such that
\begin{equation}\label{asym_1}
q+\psi_{st}(\xi)\not\in (-\infty,0],\ \forall\ q\in [\sigma,+\infty)+({\mathcal C}_{\pi/2+\gamma_0}\cup\{0\}),\  \xi\in {\mathcal C}_{\gamma_-,\gamma_+}.
\end{equation}
\item
If $\alpha\neq 1$ or $\alpha=1$ and $c_+=c_-$, then there exist $\sigma>0, \gamma_0\in (0,\pi/2)$ and $\gamma_-<0<\gamma_+$ such that $\psi^0_{st}$
satisfies \eq{asym_1}.
\item
If $\alpha=1$ and $c_+\neq c_-$ or $\alpha\in (0,1)$ and $\mu\neq 0$, then a collection  of $\sigma>0, \gamma_0\in (\pi/2,\pi)$ and $\gamma_-<0<\gamma_+$ such that
\eq{asym_1} holds does not exist.

\end{enumerate}
\end{lem}
\begin{proof} (a) It follows from Lemma \ref{lem:asymp_psi_st} that for any $\varphi\in (-\pi/2,\pi/2)$,
\[
\psi_{st}(\rho e^{i\varphi})\sim C e^{i(\varphi_1+\alpha\varphi)}\rho^\alpha,\ \rho\to +\infty,
\]
where $C>0$ and $\varphi_1\in (-\pi/2,\pi/2)$. Hence, for any $\delta>0$, there exist $\gamma_-<0<\gamma_+$ and $\sigma>0$ such that
for all $\xi\in {\mathcal C}_{\gamma_-,\gamma_+}$, $\mathrm{arg}\, \psi_{st}(\xi)\in {\mathcal C}_{\pi/2-\delta}$. Then \eq{asym_1} holds for any $\gamma_0\in (0,\delta)$.

(b) is immediate from (a).

(c) If $\alpha=1$ and $c_+\neq c_-$, then, for any $\omega\in (-\pi/2,\pi/2)$, there exists $c(\omega)\in {\mathbb R}\setminus\{0\}$
such that
  $
  \psi_{sr}(\rho e^{i\omega})\sim c(\omega) e^{i(\pi/2+\omega)}\rho \ln\rho$ as $ \rho\to+\infty.
  $
  If $\alpha\in (0,1)$ and $\mu\neq 0$, then  the  asymptotic formula holds without the factor $\ln\rho$.
  It follows that if $a>0$ is sufficiently large, then there exists $\rho>0$ s.t. either $\sigma+ia+\psi_{sr}(\rho e^{i\omega})\in (-\infty,0]$
  or $\sigma-ia+\psi_{sr}(\rho e^{i\omega})\in (-\infty,0]$.
   \end{proof}
Assume that either $\alpha\neq 1$ or $\alpha=1$ and $c_+=c_-$. Then $\psi^0_{st}$ satisfies the condition in (a).
 We find $\sigma, \gamma_0, -\pi/4<\gamma_-<0<\gamma_+<\pi/4$ such that \eq{asym_1} holds for $\psi^0_{st}$
and $\psi^0_{BM}(\xi^2)$. Then choose $\sigma_\ell>\sigma$,
$\omega_\ell\in (0,\gamma_0)$, next,  $b_\ell>0$ such that  $\sigma_\ell-b_\ell\sin\omega_\ell>0$, and then define
 \begin{equation}\label{eq:sinhLapl}
\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}(y)=\sigma_\ell +i b_\ell\sinh(i\omega_\ell+y).
\end{equation}
 Deform
 the line of integration in the Bromwich integral to
${\mathcal L}^{(L)}=\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}(R)$. For $q\in {\mathcal L}^{(L)}$, we can calculate
${\tilde V}(q,x)$ choosing $\omega\in (\gamma_-,\gamma_+)$ (the choice depends on the sign of $x'$ as in the case $q>0$) and deforming
the line of integration on the RHS of \eq{eq:V1TBr} into the angle $e^{i\omega}{\mathbb R}_+\cup e^{i(\pi-\omega)}{\mathbb R}_+$:
\[
V (T,x)=\frac{1}{2\pi i}\int_{{\mathcal L}^{(L)}}dq\, e^{qT}\frac{1}{2\pi}\left(\int_{e^{i(\pi-\omega)}{\mathbb R}_+}+\int_{e^{i\omega}{\mathbb R}_+}\right)e^{ix'\xi}\left(\frac{1}{q+\psi^0_{st}(\xi)}-\frac{1}{q+|{C_+}|\xi^2}\right){\hat f}_0(\xi)d\xi.
\]
We make the change of variables $\xi=e^{i(\pi-\omega)+y'}$ and $\xi=e^{i\omega+y'}$ in the inner integrals, and
$q=\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}(y)$ in the outer integral. Each of the resulting integrals is calculated using the simplified trapezoid rule.

\section{Wiener-Hopf factorization for L\'evy processes and main theorem}\label{ss:WHF_gen}
The definitions and results in this Section are valid for any L\'evy process $X$ on ${\mathbb R}$. Let  ${\bar X}_t=\sup_{0\le s\le t}X_s$ and ${\underline X}_t=\inf_{0\le s\le t}X_s$ be the supremum
and infimum processes (defined path-wise, a.s.), and $X_0={\bar X}_0={\underline X}_0=0$.
Let $q>0$ and let $T_q$
be an exponentially distributed random variable of mean $1/q$, independent of $X$. Set ${\phi^+_q}(\xi)={\mathbb E}[e^{i\xi{\bar X}_{T_q}}]$,
${\phi^-_q}(\xi)={\mathbb E}[e^{i\xi{\underline X}_{T_q}}]$, and define the (normalized) expected present value operators ${\mathcal E_q}$ and ${\mathcal E}^\pm_q$ by
${\mathcal E_q} u(x)={\mathbb E}[u(x+X_{T_q})]$, ${\mathcal E^+_q} u(x)={\mathbb E}[u(x+{\bar X}_{T_q})]$, ${\mathcal E^-_q} u(x)={\mathbb E}[u(x+{\underline X}_{T_q})]$. Evidently,
${\mathcal E_q}, {\mathcal E}^\pm_q: L_\infty({\mathbb R})\to L_\infty({\mathbb R})$ are bounded operators and, for $\xi\in{\mathbb R}$,
\[
{\mathcal E_q} e^{ix\xi}=\frac{q}{q+\psi(\xi)}e^{ix\xi},\ {\mathcal E}^\pm_q e^{ix\xi}=\phi^\pm_q(\xi)e^{ix\xi}.
\]
The probabilisitic form
 \begin{equation}\label{whf0}
 \frac{q}{q+\psi(\xi)}={\phi^+_q}(\xi){\phi^-_q}(\xi)
 \end{equation}
 of the Wiener-Hopf factorization identity is a special case of the Wiener-Hopf factorization of functions of wider classes than
 the ones arising in probability. Likewise, the operator form
  \begin{equation}\label{eq:operWHF}
 {\mathcal E_q}={\mathcal E^+_q}{\mathcal E^-_q}={\mathcal E^-_q}{\mathcal E^+_q}
 \end{equation}
 is a special case of the Wiener-Hopf factorziation used in the general theory of boundary problems
 for pseudo-differential operators (pdo), where more general classes of functions and operators appear
(see, e.g., \cite{eskin}) but additional regularity conditions on functions are imposed. The probabilistic version \eq{eq:operWHF} was proved in   \cite{BLSIAM02,NG-MBS,barrier-RLPE,EPV,IDUU} under additional regularity conditions on the process, and in \cite{single}, for any L\'evy process.
Both forms \eq{whf0} and \eq{eq:operWHF} are immediate from the following lemma; however, \eq{whf0} was derived earlier. See \cite{sato} for the bibliography.
\begin{lem}\label{l:deep} (\cite[Lemma 2.1]{greenwood-pitman}, and
\cite[p.81]{RW})
Let $X$ and $T_q$ be as above.  Then
\begin{enumerate}[(a)]
\item the random variables ${\bar X}_{T_q}$ and
$X_{T_q}-{\bar X}_{T_q}$ are independent; and

\item the random variables ${\underline X}_{T_q}$ and
$X_{T_q}-{\bar X}_{T_q}$ are identical in law.
\end{enumerate}
\end{lem}
(By symmetry, the statements (a), (b) are valid with ${\bar X}$ and ${\underline X}$ interchanged).
 By definition, part (a) amounts to the statement that the
probability distribution of the ${\mathbb R}^2$-valued random variable
$({\bar X}_{T_q}, X_{T_q}-{\bar X}_{T_q})$ is equal to the product (in the sense
of ``product measure'') of the distribution of ${\bar X}_{T_q}$ and the
distribution of $X_{T_q}-{\bar X}_{T_q}$. Using this observation, in \cite{EfficientLevyExtremum}, we proved the following theorem.
In the formulation of the theorem, we use the notation
$U_+=\{(x\in {\mathbb R}^2\ |\ x_1\le x_2\}$;  $f_+$ is the extension of $f:U\to {\mathbb C}$ by 0 to a function on ${\mathbb R}^2$, $I$ denotes the identity operator, and $\Delta$ is the diagonal map: $\Delta(x)=(x,x)$.
\begin{thm}\label{thm:X_barX_exp}
 Let $X$ be a L\'evy process on ${\mathbb R}$, $q>0$, and let $f:U_+\to {\mathbb R}$ be a measurable and uniformly bounded  function
s.t.   $(({\mathcal E^-_q}\otimes I)f)\circ \Delta:{\mathbb R}\to{\mathbb R}$ is measurable.
 Then
 \begin{enumerate}
  \item
 for any $x_1\le x_2$,
 \begin{eqnarray}\label{tVq0}
 q{\tilde V}(f,q;x_1,x_2)&=&(({\mathcal E_q}\otimes I)f_+)(x_1,x_2)+({\mathcal E^+_q} w(f;q,\cdot, x_2))(x_1),
 \end{eqnarray}
 where
 \begin{equation}\label{eq:wqVtq}
 w(f;q,y,x_2)={\bf 1}_{[x_2,+\infty)}(y)((({\mathcal E^-_q}\otimes I)f_+)(y,y)-(({\mathcal E^-_q}\otimes I)f_+)(y,x_2));
 \end{equation}
  \item
   the RHS' of \eq{tVq0} and \eq{eq:wqVtq} admit analytic continuation w.r.t. $q$ to the right half-plane.
   \end{enumerate}
 \end{thm}
 If an efficient numerical procedure for ${\tilde V}(f,q;x_1,x_2)$ is developed, we can calculate the expectation $V(f;T,x_1,x_2)$ as the Bromwich integral
\begin{equation}\label{tVBrom}
V(f;T;x_1,x_2)=\frac{1}{2\pi i}\int_{\Re q=\sigma}e^{qT}{\tilde V}(f;q;x_1,x_2)\,dq,
\end{equation}
where $\sigma>0$ is arbitrary.

 \begin{rem}\label{rem:LapltVq}{\rm
   In \cite{EfficientLevyExtremum}, for the case of L\'evy processes with characteristic exponent admitting analytical continuation to a strip,
we proved Theorem \ref{thm:X_barX_exp} for exponentially increasing $f$. In the case of stable L\'evy processes,
the restrictions on the rate of increase of $f$ at infinity are more stringent. If $f$ is unbounded, we recommend to calculate the expectation of the exponentially damped $f$, apply Theorem \ref{thm:X_barX_exp}, and calculate the limit of the resulting analytic expression as the dampening parameter tends to 0. See and
Example \ref{ex:exchange_stable}.}
\end{rem}
Note the following special cases of \eq{tVq0}.
\begin{enumerate}[I.]
\item
European option on the underlying. Let $f(x_1,x_2)=g(x_1)$. Then
\begin{equation}\label{perp_opt_X}
{\tilde V}(f,q,x_1,+\infty)=q^{-1}({\mathcal E_q} g)(x_1).
\end{equation}
We considered this case in Section \ref{s:1D}.
\item
European option on the supremum. Let $f(x_1,x_2)=g(x_2)$. Then, for $x_1\le x_2$,
\begin{equation}\label{perp_opt_barX}
{\tilde V}(f,q,x_1,x_2)=q^{-1}({\mathcal E^+_q} g)(x_2).
\end{equation}
\item
No-touch option with the upper barrier $a_2$ and payoff $g(X_T)$ at maturity. Then $f(x_1,x_2)=g(x_1){\bf 1}_{(-\infty,a_2]}(x_2)$,
and  $q{\tilde V}(f,q,x_1,x_2)$ is given by the RHS of \eq{tVq0} or, equivalently, ${\tilde V}(f,q,x_1,x_2)$ is given by either of the two formulas
\begin{eqnarray}\label{perp_opt_X_barr2}
{\tilde V}(f,q,x_1,x_2)&=&q^{-1}({\mathcal E_q} g)(x_1)-q^{-1}({\mathcal E^+_q}{\bf 1}_{(a_2,+\infty)}{\mathcal E^-_q}  g)(x_1),
\end{eqnarray}
and
\begin{eqnarray}
\label{perp_opt_X_barr}
{\tilde V}(f,q,x_1,x_2)&=&q^{-1}({\mathcal E^+_q}{\bf 1}_{(-\infty,a_2]}{\mathcal E^-_q}  g)(x_1).
\end{eqnarray}
 Eq. \eq{perp_opt_X_barr} was derived
in \cite{KoBoL,NG-MBS,BLSIAM02,barrier-RLPE} for wide classes of L\'evy processes and generalized for all L\'evy processes in \cite{single}.
\end{enumerate}

We finish the section with the outline of a scheme of a numerical calculation of the first term on the RHS of \eq{tVq0}; the numerical tools for the calculations of the second term are analyzed in the rest of the main body of the paper.
 Assume that $\widehat{(f_+)}(\xi_1,x_2)$, the Fourier transform of $f_+(x_1,x_2)$ w.r.t. $x_1$, is a sum of products of oscillating exponents multiplied by functions
analytic in a cone ${\mathcal C}_\gamma\cup (-{\mathcal C}_\gamma), \gamma\in (0,\pi/2),$ around ${\mathbb R}\setminus 0$ and admitting sufficiently good bounds in neighborhoods of 0 and $\infty$. Explicitly,
\begin{equation}\label{reprfplus}
\widehat{(f_+)}(\xi_1,x_2)=\sum_{k=1}^N e^{-ia_k\xi_1}{\hat g}_k(\xi_1,x_2),
\end{equation}
and, for all $x_2\in{\mathbb R}$ and $\xi\in {\mathcal C}_\gamma\cup (-{\mathcal C}_\gamma)$,
\begin{equation}\label{boundg_k}
|{\hat g}_k(\xi_1,x_2)|\le C_k(x_2)(|\xi_1|^{\alpha_0-1}{\bf 1}_{|\xi|\le 1}+|\xi_1|^{-1+\delta}{\bf 1}_{|\xi|\ge 1}),
\end{equation}
where  $\alpha_0>-\alpha$, $ \delta<\alpha$, are independent of $x_2$, $k$ and $\xi$, and $C_k(x_2)$ depend on $k$ and $x_2$ only.

We represent $(({\mathcal E_q}\otimes I)f_+)(x_1,x_2)$ as a sum according to \eq{reprfplus}, and evaluate each term
as in Sections \ref {ss:1D_alt}-\ref{ss:1D_alt_Br} ($x_2$ is regarded as a parameter).  Then we apply
GWR algorithm. If $\alpha\in (1,2)$ or $\alpha=1$ and $c_+=c_-$ or $\alpha\in(0,1)$ and $\mu=0$, then the sinh-acceleration can be applied.

\section{Wiener-Hopf factorization for stable L\'evy processes}\label{s:WHF-StLevy}

\subsection{Case of real $q$}\label{ss:form_WHF_real}
In \cite{EfficientAmenable}, we defined a class of Stieltjes-L\'evy processes (SL-processes), and proved that
stable L\'evy processes are SL-processes. We also proved that the characteristic exponent $\psi$ of an SL-process
admits analytic continuation to ${\mathbb C}\setminus i{\mathbb R}$ and
enjoys the following important property:
\begin{equation}\label{eq:eqSL}
\forall\ q>0\ \mathrm{and}\ \xi\in {\mathbb C}\setminus i{\mathbb R},\ q+\psi(\xi)\neq 0.
\end{equation}
Using \eq{eq:eqSL} and Proposition \ref{prop:anal_psi_st}, we derive
\begin{thm}\label{tm:anal_psi_st} Let $q>0$. Then
\begin{enumerate}[(a)]
\item
 ${\phi^+_q}(\xi)$ admits analytic continuation w.r.t. $\xi$ from the upper half-plane to the complex plane with the cut $i(-\infty,0]$ (and an appropriate Riemann surface)
 by
 \begin{equation}\label{anal_cont_phipq}
 {\phi^+_q}(\xi)=q/((q+\psi_{st}(\xi)){\phi^-_q}(\xi));
 \end{equation}
 \item
${\phi^-_q}(\xi)$ admits analytic continuation w.r.t. $\xi$ from the lower half-plane to the complex plane with the cut $i[0,+\infty)$ (and an appropriate Riemann surface)
 by
  \begin{equation}\label{anal_cont_phimq}
  {\phi^-_q}(\xi)=q/((q+\psi_{st}(\xi)){\phi^+_q}(\xi)).
  \end{equation}
 \end{enumerate}
 \end{thm}
For $\om_+\in (0,\pi/2)$ and $\om_-\in (-\pi/2,0)$, define the contours ${\mathcal L}^+_{\om_+}=e^{i(\pi-\om_+)}{\mathbb R}_+\cup e^{i\om_+}{\mathbb R}_+$ and
${\mathcal L}^-_{\om_-}=e^{i(-\pi-\om_-)}{\mathbb R}_+\cup e^{i\om_-}{\mathbb R}_+$.
 The direction on each contour is from the left to the right. For $q>0$ and $\xi\in {\mathbb R}\setminus\{0\}$, define
\begin{eqnarray}\label{phip1}
\phi^{+,'}_q(\xi)&=&\exp\left[\frac{1}{2\pi i}\int_{{\mathcal L}^-_{\om_-}}\frac{\xi\ln(1+\psi_{st}(\eta)/q)}{\eta(\xi-\eta)}d\eta\right],\\
\label{phim1}
\phi^{-,'}_q(\xi)&=&\exp\left[-\frac{1}{2\pi i}\int_{{\mathcal L}^+_{\om_+}}\frac{\xi\ln(1+\psi_{st}(\eta)/q)}{\eta(\xi-\eta)}d\eta\right],
\end{eqnarray}
and set $\phi^{\pm,'}_q(0)=1$. On the strength of Lemma \ref{SLroots},
 for any $\omega\in(0,\pi/2)$, the integrands on the RHS' of \eq{phip1} and  \eq{phim1} are well-defined.
Since $\ln (1+\psi_{st}(\xi)/q)=O(|\eta|^\alpha)$ as $\eta\to 0$, and $=O(|\eta|^\epsilon)$ as $\eta\to \infty$, for any $\epsilon>0$,
the integrals are finite.
\begin{lem}\label{lem:stableWHF}
For any $q>0$ and $\xi\in {\mathbb R}$, $\phi^\pm_q(\xi)=\phi^{\pm,'}_q(\xi)$.
\end{lem}
\begin{proof} It suffices to consider $\xi\in {\mathbb R}\setminus\{0\}$. If $\xi>0$, we note that both integrands are analytic in
${\mathcal C}_{\om_-,\om_+}$, decay at infinity faster than $|\xi|^{-3/2}$ and continuous up to the boundary of ${\mathcal C}_{\om_-,\om_+}$.
Hence, the  residue theorem is applicable, and
we obtain
\[
\frac{1}{2\pi i}\int_{{\mathcal L}^-_{\om_-}}\frac{\xi\ln(1+\psi_{st}(\eta)/q)}{\eta(\xi-\eta)}d\eta-\frac{1}{2\pi i}\int_{{\mathcal L}^+_{\om_+}}\frac{\xi\ln(1+\psi_{st}(\eta)/q)}{\eta(\xi-\eta)}d\eta=\ln\frac{q}{q+\psi_{st}(\xi)},
\]
therefore, \eq{whf0} holds with the product $\phi^{+,'}_q(\xi)\phi^{-,'}_q(\xi)$ on the RHS. It is known that the Wiener-Hopf factors $\phi^\pm_q(\xi)$
and their reciprocals are polynomially bounded as $\xi\to\infty$ in the corresponding half-plane.  The proof of Lemma \ref{asympWHF} below implies that  $\phi^{+,'}_q(\xi)$ and $\phi^{-,'}_q(\xi)$ are polynomially bounded, and a modification of same proof
can be used to prove that the reciprocals are polynomially bounded. The standard argument based on Moreira's theorem gives
$\phi^\pm_q(\xi)=\phi^{\pm,'}_q(\xi)$.
\end{proof}

  Define
$\varphi_0= \varphi_0(\alpha,{c_+},{c_-})=\mathrm{arg}\,{C_+}(\alpha,{c_+},{c_-})$, and
\begin{enumerate}[i.]
\item
if $\alpha\in (1,2)$ or $\alpha\in (0,1)$ and $\mu=0$, set $\alpha^\pm=\alpha/2\pm \varphi_0/\pi$;
\item
if $\alpha\in (0,1)$ and $\mu>0$, set $\alpha_+=1, \alpha_-=0$;
\item
if $\alpha\in (0,1)$ and $\mu<0$, set $\alpha_+=0, \alpha_-=1$;
\item
if $\alpha=1$ and $c_+=c_-$, set $\alpha_+=\alpha_-=1/2$.
\end{enumerate}
 The following lemma is proved  in Section \ref{proof:lem_asympWHF}; the asymptotic formulas are valid uniformly in $\xi\in {\mathcal C}_\gamma$,
 for any $\gamma\in (0,\pi/2)$.
 \begin{lem}\label{asympWHF}
  Let $q>0$ and $\gamma\in (0,\pi/2)$. Then,
 \begin{enumerate}[(1)]
 \item
 as $({\mathcal C}_\gamma\ni)\xi\to 0$,
 \begin{equation}\label{asphim0}
 \phi^\pm(\xi)=1+O(|\xi|^{\min\{1,\alpha\}});
 \end{equation}
 \item
if either
$\alpha\in (1,2)$ and both $c_\pm>0$ or $\alpha\in(0,1)$ and $\mu=0$ or $\alpha=1$ and $c_+=c_-$, then both
$\alpha_\pm>0$ and,
 as $({\mathcal C}_\gamma\ni)\xi\to \infty$,
  \begin{equation}\label{asphipmqmain}
\phi^\pm_q(\xi)=O(|\xi|^{-\alpha_\pm});
\end{equation}
\item
\begin{enumerate}[(a)]
\item
if $\alpha\in (1,2)$ and $c_-=0$, then $\alpha_+=\alpha, \alpha_-=0$, and  $\exists\ \delta>0$ s.t. as $({\mathcal C}_\gamma\ni)\xi\to \infty$,
\begin{eqnarray}\label{asphipqp12}
{\phi^+_q}(\xi)&=&O(|\xi|^{-\alpha}),\\\label{asphimqp12}
{\phi^-_q}(\xi)&=&a^-_q + O(|\xi|^{-\delta}),
\end{eqnarray}
where
\begin{equation}\label{amq12}
a^-_q=\exp\left[-\frac{1}{2\pi i}\int_{{\mathcal L}^+_{\om_+}}\frac{\ln(1-i\mu\eta/(q+\psi^0_{st}(\eta))}{\eta}d\eta\right],
\end{equation}
and $\om_+\in (0,\pi/2)$ is arbitrary;
\item
if $\alpha\in (1,2)$ and $c_+=0$, then $\alpha_+=0, \alpha_-=\alpha$, and $\exists\ \delta>0$ s.t. as $({\mathcal C}_\gamma\ni)\xi\to \infty$,

\begin{eqnarray}\label{asphipqm12}
{\phi^+_q}(\xi)&=&a^+_q + O(|\xi|^{-\delta}),\\\label{asphimqm12}
{\phi^-_q}(\xi)&=&O(|\xi|^{-\alpha}),
\end{eqnarray}
where
\begin{equation}\label{apq12}
a^+_q=\exp\left[\frac{1}{2\pi i}\int_{{\mathcal L}^-_{\om_-}}\frac{\ln(1-i\mu\eta/(q+\psi^0_{st}(\eta))}{\eta}d\eta\right],
\end{equation}
and $\om_-\in (-\pi/2,0)$ is arbitrary;
\item
if $\alpha\in (0,1)$ and $\mu> 0$, then there exists $\delta>0$ such that as $({\mathcal C}_\gamma\ni)\xi\to \infty$,
\begin{eqnarray}\label{asphipqp01}
{\phi^+_q}(\xi)&=&O(|\xi|^{-1}),\\\label{asphimqp01}
{\phi^-_q}(\xi)&=&a^-_q + O(|\xi|^{-\delta}),
\end{eqnarray}
where
\begin{equation}\label{amq}
a^-_q=\exp\left[-\frac{1}{2\pi i}\int_{{\mathcal L}^+_{\om_+}}\frac{\ln((1+\psi^0(\eta)/q)/(1-i\mu\eta/q))}{\eta}d\eta\right],
\end{equation}
and $\om_+\in (0,\pi/2)$ is arbitrary;

\item
if $\alpha\in (0,1)$ and $\mu< 0$, then then there exists $\delta>0$ such that as $({\mathcal C}_\gamma\ni)\xi\to \infty$,
\begin{eqnarray}\label{asphipqm01}
{\phi^+_q}(\xi)&=&a^+_q+ O(|\xi|^{-\delta}),\\\label{asphimqm01}
{\phi^-_q}(\xi)&=&O(|\xi|^{-1}),
\end{eqnarray}
where
\begin{equation}\label{apq}
a^+_q=\exp\left[\frac{1}{2\pi i}\int_{{\mathcal L}^-_{\om_-}}\frac{\ln((1+\psi^0(\eta)/q)/(1-i\mu\eta/q))}{\eta}d\eta\right],
\end{equation}
and $\om_-\in (-\pi/2,0)$ is arbitrary;

\item
if $\alpha=1$ and $c_+>c_-$, then, for any $\epsilon>0$, as $({\mathcal C}_\gamma\ni)\xi\to \infty$,
\begin{eqnarray}\label{asphipq1p}
{\phi^+_q}(\xi)&=&O(|\xi|^{-1+\epsilon}),\\\label{asphimq1p}
{\phi^-_q}(\xi)&=& O(|\xi|^{\epsilon})
\end{eqnarray}
(and ${\phi^-_q}(\xi)$ is uniformly bounded on the lower half-plane);
\item
if $\alpha=1$ and $c_+<c_-$, then, for any $\epsilon>0$,
\begin{eqnarray}\label{asphipq1m}
{\phi^+_q}(\xi)&=&O(|\xi|^{\epsilon}),\\\label{asphimq1m}
{\phi^-_q}(\xi)&=& O(|\xi|^{-1-\epsilon})
\end{eqnarray}
(and ${\phi^+_q}(\xi)$ is uniformly bounded on the upper half-plane).
\end{enumerate}
\end{enumerate}
\end{lem}

 \begin{cor}\label{cor:asphipq}
 Let $X$ be either a stable L\'evy process of index $\alpha\neq 1$ or index $\alpha=1$ and $c_+=c_-$.
 Then, for any $q>0$, there exist $a^\pm_q\ge 0$ and $\delta_\pm>0$ such that, for any $\gamma\in (-\pi/2,\pi/2)$,
 \begin{equation}\label{eq:asqpmq}
 \phi^\pm_q(\xi)=a^\pm_q+O(|\xi|^{-\delta_\pm}), \ ({\mathcal C}_\gamma\ni)\xi\to \infty.
 \end{equation}
\end{cor}
\begin{rem}\label{rem:errorFT}{\rm We use \eq{asphim0} and \eq{eq:asqpmq} to regularize the action of the EPV operators ${\mathcal E}^\pm_q$ as PDO: ${\mathcal E}^\pm_q u(x)={\mathcal F}^{-1}_{\xi\to x}\phi^\pm_q(\xi){\mathcal F}_{x\to\xi} u(x)$. Without the regularization, the integrands in the integral representation  that we derive decay too slowly at infinity, and accurate numerical realizations require unnecessarily long grids.
 }
 \end{rem}
 For $\beta >0$, introduce   the convolution operators  $I^\pm_{\beta}$ by
\begin{eqnarray}\label{def:Ip}
I^+_\beta u(x)&=&\beta\int_{0}^{\infty} e^{-\beta y}u(x+y)dy\\\label{def:Im}
I^-_{\beta} u(x)&=&\beta \int_{-\infty}^{0} e^{\beta y}u(x+y)dy.
\end{eqnarray}
The symbol of $I^\pm_{\beta}$ is $(1\mp i\xi/\beta)^{-1}$. The following proposition is evident.
\begin{lem}\label{lem_supp} Let $a_1, a_2\in {\mathbb C}$, $\beta>0$ and let there exist $\beta_1<\beta$ such that $e^{-\beta_1|\cdot|}u\in L_1$. Then
\begin{enumerate}[(a)]
\item
if $u(x)=0$ for all $x< h$ (resp., $x\le h$), then $((a_1I+a_2I^-_{\beta})u)(x)=0, x<h$ (resp., $x\le h$);
\item
if $u(x)=0$ for all $x> h$ (resp., $x\ge h$), then $((a_1I+a_2I^+_{\beta})u)(x)=0, x>h$ (resp., $x\ge h$).

\end{enumerate}
\end{lem}
We calculate the integrals on the RHS' of the formulas \eq{phip1}-\eq{phim1}  for the Wiener-Hopf factors and \eq{amq12}, \eq{apq12}, \eq{amq},
\eq{apq} for the constants $a^\pm_q$ representing each integral as a sum of integrals over two rays of the form
$e^{i\gamma}{\mathbb R}_+$. In the integral over a ray of the form $e^{i\gamma}{\mathbb R}_+$, we change the variable $\eta=e^{i\gamma+y}$ and apply the simplified trapezoid rule. In the process of deformation, the expressions under the logarithm sign must be well-defined and may not assume values in $(-\infty,0]$.

\subsection{Case of complex $q$}\label{ss:form_WHF_complex}
If the conformal deformation technique is applied to the Bromwich integral, then the Wiener-Hopf factors need to be calculated
for $q$ on the (deformed) contour of integration, and additional conditions  on the parameters of the process (see  Lemma \ref{q_xi_cone}(a)) and the deformed  contour have to be imposed. The restrictions on the parameters of the process and the rotation parameter $\omega$ in the formulas for the Wiener-Hopf factors and formulas for ${\tilde V}(f;q;x_1,x_2)$ are similar to but different from to the restrictions in Section \ref{ss:1D_alt_Br}.
The condition in Lemma \ref{q_xi_cone} (a)  must hold for $\psi_{st}$ because \eq{asym_1} must hold for $\psi_{st}$ whereas in Section \ref{ss:1D_alt_Br},
we need \eq{asym_1} to hold for $\psi^0_{st}$ and $\psi_{BM}$.

\begin{lem}\label{lem:anal_cont_WHF}
Let the condition in Lemma \ref{q_xi_cone} (a)  hold. Then there exist $\sigma>0, \gamma_0\in (0,\pi/2)$ and $\gamma_-<0<\gamma_+$ such that
\begin{enumerate}[(a)]
\item
${\phi^+_q}(\xi)$ admits analytic continuation to $\{(q,\xi)\ |\ q\in \sigma+{\mathcal C}_{\pi/2+\gamma_0},
\xi\in i{\mathcal C}_{\pi/2-\gamma_-}\}$. For any $\om_-\in (\gamma_-,0)$ and $\gamma^{--}\in (\om_-,0)$, the restriction of ${\phi^+_q}(\xi)$ on $\{(q,\xi)\ |\ q\in \sigma+{\mathcal C}_{\pi/2+\gamma_0},
\xi\in i{\mathcal C}_{\pi/2-\gamma^{--}}\}$ is defined by
 \eq{phip1};
\item
${\phi^-_q}(\xi)$ admits analytic continuation to $\{(q,\xi)\ |\ q\in \sigma+{\mathcal C}_{\pi/2+\gamma_0},
\xi\in -i{\mathcal C}_{\pi/2+\gamma_+}\}$. For any $\om_+\in (0,\gamma_+)$ and $\gamma^{++}\in (0,\om_+)$, the restriction of ${\phi^-_q}(\xi)$ on $\{(q,\xi)\ |\ q\in \sigma+{\mathcal C}_{\pi/2+\gamma_0},
\xi\in -i{\mathcal C}_{\pi/2+\gamma_{++}}\}$ is defined by
 \eq{phim1};
\item
if $\alpha\in (1,2)$ and $c_-=0$, then $a^-_q$ defined by \eq{amq12} is an analytic function on  $\sigma+{\mathcal C}_{\pi/2+\gamma_0}$;
\item
if $\alpha\in (1,2)$ and $c_+=0$, then $a^+_q$ defined by \eq{apq12} is an analytic function on  $\sigma+{\mathcal C}_{\pi/2+\gamma_0}$.
\end{enumerate}
\end{lem}
\begin{proof} (a) By Lemma \ref{q_xi_cone}, the condition in Lemma \ref{q_xi_cone} (a)  implies that there exist $\sigma, \gamma_0, -\pi/2<\gamma_-<0<\gamma_+<\pi/2$ such that \eq{asym_1} holds. Take $\om_-\in (\gamma_-,0)$, and, for $q\in \sigma+{\mathcal C}_{\pi/2+\gamma_0}$ and $\xi$ above ${\mathcal L}^-_{\om_-}$, define ${\phi^+_q}(\xi)$ by \eq{phip1}. It follows from \eq{asym_1} that ${\phi^+_q}(\xi)$ is analytic on $\{q\in \sigma+{\mathcal C}_{\pi/2+\gamma_0},
\xi\in i{\mathcal C}_{\pi/2-\om_-}\}$. It remains to note that $\om_-\in (\gamma_-,0)$ is arbitrary.

(b) The proof is by symmetry.

(c,d) If $\psi_{st}$ satisfies the condition in Lemma \ref{q_xi_cone} (a), then  $\psi^0_{st}$ satisfies this condition  as well. Hence, there exist $\sigma, \gamma_0, -\pi/2<\gamma_-<0<\gamma_+<\pi/2$ such that   \eq{asym_1} holds for $\psi_{st}$ and $\psi^0_{st}$.

\end{proof}

\section{Integral representations}\label{int_repr
In the remainder of the main body of the paper, we assume that either  $\alpha\neq 1$ or $\alpha=1$ and $c_+=c_-$, and we use the representation \eq{eq:asqpmq}. We exclude the case of asymmetric stable L\'evy processes of order $\alpha=1$ because
of the irregular behavior of the Wiener-Hopf factors.

\subsection{Calculation of the cpdf of ${\bar X}_T$}\label{ss:cpdf_barXT}
Let $a\ge x$, and $T>0$. We calculate
\begin{equation}\label{cpdfsup0}
V(T,a;x)={\mathbb E}[{\bf 1}_{x+{\bar X}_T\le a}]=1-{\mathbb E}[{\bf 1}_{x+{\bar X}_T> a}]=1-{\mathbb E}[{\bf 1}_{{\bar X}_T> a-x}].
\end{equation}
Take $\sigma>0$ and use the Bromwich integral:
\begin{eqnarray}\label{VcpdfbarX}
{\mathbb E}[{\bf 1}_{x+{\bar X}_T> a}]&=&\frac{1}{2\pi i}\int_{\Re q=\sigma}q^{-1}({\mathcal E^+_q} {\bf 1}_{(a,+\infty)})(x).
\end{eqnarray}
On the strength of Lemma \ref{lem_supp} (a), we can replace ${\mathcal E^+_q}$ with ${\mathcal E^+_q}-(a^+_qI+(1-a^+_qI^-_1)$. The symbol of the modified operator
$\phi^+_{q,mod}(\xi_1):={\phi^+_q}(\xi_1)-(a^+_q+(1-a^+_q)/(1+i\xi_1))$ is $O(|\xi_1|^{\min\{1,\alpha\}})$ as $\xi_1\to 0$, and $O(|\xi_1|^{-\delta_+})$ as $\xi_1\to\infty$, where $\delta_+>0$ is from \eq{eq:asqpmq},
and $\widehat{{\bf 1}_{(a,+\infty)}}(\xi_1)=e^{-ia\xi_1}/(i\xi_1)$.
Therefore, for $x\le a$,
\begin{equation}\label{cpdfsup}
({\mathcal E^+_q} {\bf 1}_{(a,+\infty)})(x)=\frac{1}{2\pi}\int_{\mathbb R} e^{i(x-a)\xi_1}\frac{{\phi^+_q}(\xi_1)-(a^+_q+(1-a^+_q)(1+i\xi_1)^{-1})}{i\xi_1}d\xi_1,
\end{equation} where for any $\gamma\in (0,\pi/2)$, the integrand
is uniformly bounded on $\{\mathrm{arg}\, \xi_1\in [-\pi, -\pi+\gamma]\cup [-\gamma,0]\}$ and
 decays as $|\xi_1|^{-1-\delta_+}$ as
$\xi_1\to \infty$ in $\{\xi_1\ |\ \mathrm{arg}\, \xi_1\in [-\pi, -\pi+\gamma]\cup [-\gamma,0]\}$.

\subsubsection{Numerical realization based on GWR algorithm}\label{sss:cpdf_barXT_numer_GWR}
For each $q>0$ used in GWR algorithm, we calculate $({\mathcal E^+_q} {\bf 1}_{(a,+\infty)})(x), x\le a$.  Since $x-a\le 0$, we can take any $\om_-\in (-\pi/2,0)$, and deform
the contour of integration into ${\mathcal L}^-_{\om_-}=e^{i\om_-}{\mathbb R}_+\cup e^{i(\pi-\om_-)}{\mathbb R}_+$:
\begin{equation}\label{cpdfsup2}
({\mathcal E^+_q} {\bf 1}_{(a,+\infty)})(x)=\frac{1}{2\pi}\int_{{\mathcal L}^-_{\om_-}} e^{i(x-a)\xi_1}\frac{{\phi^+_q}(\xi_1)-(a^+_q+(1-a^+_q)(1+i\xi_1)^{-1})}{i\xi_1}d\xi_1.
\end{equation}
We use  $\overline{{\phi^+_q}(\xi)}={\phi^+_q}(-\bar\xi)$ and $\overline{e^{i(x-a)\xi}}=e^{-i(x-a)\bar\xi}$ to obtain
\begin{equation}\label{cpdfsup3}
({\mathcal E^+_q} {\bf 1}_{(a,+\infty)})(x)=\frac{1}{\pi}\Re\int_{e^{i\om_-}{\mathbb R}_+} e^{i(x-a)\xi}\frac{{\phi^+_q}(\xi)-(a^+_q+(1-a^+_q)(1+i\xi)^{-1})}{i\xi}d\xi.
\end{equation}
Then we change the variable $\xi=\xi(y)=e^{i\om_-+y}$
\begin{equation}\label{cpdfsup4}
({\mathcal E^+_q} {\bf 1}_{(a,+\infty)})(x)=\frac{1}{\pi}\Im\int_{-\infty}^{+\infty} e^{i(x-a)\xi(y)}
({\phi^+_q}(\xi(y))-a^+_q-(1-a^+_q)/(1+i\xi(y)))dy,
\end{equation}
and evaluate the integral on the RHS of \eq{cpdfsup4} using the simplified trapezoid rule. We apply the bounds and constructions analogous to the ones
in \cite{ConfAccelerationStable}, where the pdf and cpdf of stable distributions were efficiently evaluated. In order to find the number of terms in the simplified trapezoid rule sufficient to achieve the target precision, we use the bound for $e^{i(x-a)\xi(y)}\phi^+_{q,mod}(\xi(y))$:
\begin{eqnarray}\label{eq:asymp0cpdfminf}
e^{i(x-a)\xi(y)}\phi^+_{q,mod}(\xi(y))&=&O(e^{\min\{\alpha,1\}y}), \ y\to-\infty,\\\label{eq:asymp0cpdfpinf}
e^{i(x-a)\xi(y)}\phi^+_{q,mod}(\xi(y))&=&O(e^{-\min\{\delta,1\}y}\cdot e^{\sin(\om_-)(a-x)e^y}), \ y\to+\infty,
 \end{eqnarray}
 where  $\delta>0$ is the same as in in Lemma \ref{asympWHF}. In particular,
 if \eq{asphipmqmain} holds, $\delta=\alpha_+$.

 \subsubsection{Numerical realization based on the sinh-acceleration applied to the Bromwich integral}\label{sss:cpdf_barXT_numer_sinh}
 The sinh-acceleration can be applied to the Bromwich integral only if analytic continuation of ${\phi^+_q}(\xi)$ w.r.t. $q$ is possible, and if we deform the contour in the formula for $({\mathcal E^+_q} {\bf 1}_{(a,+\infty)})(x), x\le a$, then analytic continuation w.r.t $(q,\xi)$ must be possible.
 Assume that the condition in Lemma \ref{q_xi_cone} (a) holds. Using Lemma \ref{lem:anal_cont_WHF}, we can
find $\sigma, \gamma_0, -\pi/2<\gamma_-<0<\gamma_+<\pi/2$ such that \eq{asym_1} holds. Then choose $\sigma_\ell>\sigma$,
$\omega_\ell\in (0,\gamma_0)$, next, choose $b_\ell>0$ such that  $\sigma_\ell-b_\ell\sin\omega_\ell>0$ and define the  map
$\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}$ by \eq{eq:sinhLapl}. Choose $\om_-\in (\gamma_-,0)$. For $q=q(y')=\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}(y')$
 on the contour ${\mathcal L}^{(L)}=\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}(R)$, and $\xi\in {\mathcal L}^-_{\om_-}$,
we calculate
 $\phi^\pm_q(\xi)$ and $a^\pm_q$ as in Lemma \ref{lem:anal_cont_WHF}, and then
 \begin{eqnarray}\label{VcpdfbarXsinh}
{\mathbb E}[{\bf 1}_{x+{\bar X}_T> a}]&=&\frac{b_\ell}{2\pi }\int_{{\mathcal L}^{(L)}}dy'\,e^{q(y')T}\frac{\cosh(i\omega_\ell +y')}{q(y')}({\mathcal E}^+_{q(y')} {\bf 1}_{(a,+\infty)})(x),
\end{eqnarray}
where $({\mathcal E}^+_{q(y')} {\bf 1}_{(a,+\infty)})(x)$ is given by \eq{cpdfsup2}. We apply the simplified trapezoid rule to the integral on the RHS of
\eq{VcpdfbarXsinh}, and, for each $y'$ on the grid, evaluate
\[
({\mathcal E}^+_{q(y')} {\bf 1}_{(a,+\infty)})(x)=\sum_{\gamma=\om_-, \pi-\om_-}I_\gamma(y',x),
\]
where
\begin{equation}\label{cpdfsup5}
I_\gamma(y',x)=\frac{1}{2\pi}\int_{e^{i\gamma}{\mathbb R}_+} e^{i(x-a)\xi_1}
\frac{\phi^+_{q(y')}(\xi_1)-(a^+_{q(y')}+(1-a^+_{q(y')})(1+i\xi_1)^{-1})}{i\xi_1}d\xi_1
\end{equation}
(for each $\gamma$, the direction is from the left to the right), making the change of variables $\xi=e^{i\gamma +y}$ on the RHS of \eq{cpdfsup5}, and then applying the simplified trapezoid rule.

 \subsection{Calculation of the joint cpdf of $X_T$ and ${\bar X}_T$}\label{ss:joint_stable}
Let $a_2>0, a_1\le a_2$, $x_1\le x_2\le a_2$,  and $T>0$. Consider
\[V(T,a_1,a_2;x_1,x_2)=
{\mathbb E}[{\bf 1}_{(-\infty,a_1]}(x_1+X_T){\bf 1}_{(-\infty,a_2]}(\max\{x_2, x_1+{\bar X}_T\})].
\] If $a_1=a_2$, we have the cpdf of ${\bar X}_T$,
which we considered in Section \ref{ss:cpdf_barXT}. Hence, we assume that $a_1<a_2$. We have
$
V(T,a_1,a_2;x_1,x_2)={\mathbb E}[{\bf 1}_{(-\infty,a_1]}(x_1+X_T)]-V_1(T,a_1,a_2;x_1,x_2)$,
where \[V_1(T,a_1,a_2;x_1,x_2)={\mathbb E}[{\bf 1}_{(-\infty,a_1]}(x_1+X_T){\bf 1}_{(a_2,+\infty)}(\max\{x_2, x_1+{\bar X}_T\})].\]
Function $a_1\mapsto {\mathbb E}[{\bf 1}_{(-\infty,a_1]}(x_1+X_T)]$ is the cpdf of the stable L\'evy process starting at $x_1$, for which we developed explicit efficient algorithms in \cite{ConfAccelerationStable} and Section \ref{s:1D}.
Hence, it suffices to consider $V_1(T,a_1,a_2;x_1,x_2)$. Using \eq{perp_opt_X_barr2}, we obtain
\[
{\tilde V}_1(q,a_1,a_2;x_1,x_2)=q^{-1}({\mathcal E^+_q}{\bf 1}_{(a_2,+\infty)}{\mathcal E^-_q}  {\bf 1}_{(-\infty,a_1]})(x_1).
\]
Since $a_2>a_1$ and $x_1\le a_2$, Lemma \ref{lem_supp} allows us to replace
${\bf 1}_{(a_2,+\infty)}{\mathcal E^-_q}  {\bf 1}_{(-\infty,a_1]}$ and ${\mathcal E^+_q}{\bf 1}_{(a_2,+\infty)}$
with ${\bf 1}_{(a_2,+\infty)}{\mathcal E}^{-}_{q,mod}{\bf 1}_{(-\infty,a_1]}$ and
${\mathcal E}^{+}_{q,mod}{\bf 1}_{(a_2,+\infty)}$, respectively, where
\[
{\mathcal E}^{-}_{q,mod}={\mathcal E^-_q}-(a^-_qI+(1-a^-_q)I^+_1), \
{\mathcal E}^{+}_{q,mod}={\mathcal E^+_q}-(a^+_qI+(1-a^+_q)I^-_1).
\]
We have $\widehat{{\bf 1}_{(-\infty,a_1]})}(\xi_1)=e^{-ia_1\xi_1}/(-i\xi_1)$, and the symbols of the modified EPV operators
$\phi^+_{q,mod}(\xi_1):={\phi^+_q}(\xi_1)-(a^+_q+(1-a^+_q)/(1+i\xi_1))$ and
$\phi^-_{q,mod}(\xi_1):={\phi^-_q}(\xi_1)-(a^-_q+(1-a^-_q)/(1-i\xi_1))$
are $O(|\xi_1|^{\min\{1,\alpha\}})$ as $\xi_1\to 0$, and $O(|\xi_1|^{-\delta_\pm})$ as $\xi_1\to\infty$, where $\delta_\pm>0$ are from \eq{eq:asqpmq}.
Therefore, we can use the pdo representation of the modified EPV operators
to calculate, for $y\ge a_2$,
\[
({\mathcal E}^{-}_{q,mod}{\bf 1}_{(-\infty,a_1]})(y)=\frac{1}{2\pi}\int_{{\mathbb R}}e^{i(y-a_1)\xi_1}\frac{\phi^-_{q,mod}(\xi_1)}{-i\xi_1}d\xi_1.
\]
Let $q>0$. Since $y-a_1>0$, we can deform the contour of integration: for any $\om_+\in (0,\pi/2)$,
\[
({\mathcal E}^{-}_{q,mod}{\bf 1}_{(-\infty,a_1]})(y)=\frac{1}{2\pi}\int_{{\mathcal L}^+}e^{i(y-a_1)\xi_1}\frac{\phi^-_{q,mod}(\xi_1)}{-i\xi_1}d\xi_1.
\]
where ${\mathcal L}^+={\mathcal L}^+_{\om_+}=e^{i\om_+}{\mathbb R}_+\cup e^{i(\pi-\om_+)}{\mathbb R}_+$. (In the case of $q\not\in (0,+\infty)$, we need to impose the condition in Lemma \ref{q_xi_cone} (a), and choose
$\om_+$ so that $q+\psi(\xi_1)\neq 0$ for all $q$ and $\xi_1$ arising in the process of deformation).
Then we calculate the Fourier transform $\hat w(\eta)$ of $w:={\bf 1}_{(a_2,+\infty)}{\mathcal E}^{-}_{q,mod}{\bf 1}_{(-\infty,a_1]}$ for
$\{\eta\ |\ \Im\eta\le 0, \eta\neq 0\}$;  since the integrand below decays exponentially as a function of $y$ as $y\to+\infty$, and as $|\xi|^{-1-\delta_-}$ as $\xi_1\to\infty$, where $\delta_->0$, Fubini's theorem is applicable, and we obtain
\begin{eqnarray*}
\hat w(\eta)&=&\int_{a_2}^{+\infty} dy e^{-iy\eta}\frac{1}{2\pi}\int_{{\mathcal L}^+}e^{i(y-a_1)\xi_1}\frac{\phi^-_{q,mod}(\xi_1)}{-i\xi_1}d\xi_1\\
&=&e^{-ia_2\eta}\frac{1}{2\pi}\int_{{\mathcal L}^+}e^{i(a_2-a_1)\xi_1}\frac{\phi^-_{q,mod}(\xi_1)}{(-i\xi_1)i(\eta-\xi_1)}d\xi_1.
\end{eqnarray*}
Let $q>0$. Fix $\om_-\in (-\pi/2,0)$, and introduce
${\Omega}=\{(\xi_1,\eta)\ |\ \xi_1\in ({\mathcal L}^+\setminus\{0\}), \mathrm{arg}\, \eta\in [-\pi,-\pi-\om_-]\cup [\om_-,0]\}$
and ${\mathcal L}^-={\mathcal L}^-_{\om_-}=e^{i\om_-}{\mathbb R}_+\cup e^{i(\pi-\om_-)}{\mathbb R}_+$.
Since $x_1\le a_2$, the absolute value of the integrand of the double integral
\[{\tilde V}_1(q,a_1,a_2;x_1,x_2)
=\frac{1}{(2\pi)^2 q}\int_{\mathbb R} d\eta\, e^{i(x_1-a_2)\eta}\phi^+_{q,mod}(\eta)\int_{{\mathcal L}^+}e^{i(a_2-a_1)\xi_1}\frac{\phi^-_{q,mod}(\xi_1)}{\xi_1(\eta-\xi_1)}d\xi_1
\]
 is bounded by
$C(q,\om_-,\om_+)g(|\xi_1|,|\eta|)$, where $C(q,\om_-,\om_+)$ is independent of  $(\xi_1,\eta)\in {\Omega}$ and
\[
g(|\xi_1|,|\eta|)=(|\eta|^{-\delta_+}{\bf 1}_{|\eta|\ge 1}+|\eta|^{\min\{1,\alpha\}}{\bf 1}_{|\eta|\le 1})(|\eta|+|\xi|)^{-1}
(|\xi_1|^{-1-\delta_-}{\bf 1}_{|\xi_1|\ge 1}+|\xi_1|^{\min\{0,\alpha-1\}}{\bf 1}_{|\xi_1|\le 1}).
\]
Considering separately the subregions defined by the pairs of inequalities (1) $|\xi_1|\le1, |\eta|\le 1$; (2)
$|\xi_1|\ge1, |\eta|\le 1$; (3) $|\xi_1|\le1, |\eta|\ge 1$; (4) $|\xi_1|\ge1, |\eta|\ge 1$; and taking into account that
$\delta_\pm>0$ and $\alpha>0$, we conclude that
$g\in L_1({\mathbb R}^2)$. Hence, we can change the order of integration and obtain
\begin{equation}\label{tV1q} {\tilde V}_1(q,a_1,a_2;x_1,x_2)=\frac{1}{(2\pi)^2 q}\int_{{\mathcal L}^-} d\eta\, e^{i(x_1-a_2)\eta}\phi^+_{q,mod}(\eta)
\int_{{\mathcal L}^+}\frac{e^{i(a_2-a_1)\xi_1}\phi^-_{q,mod}(\xi_1)}{\xi(\eta-\xi_1)}d\xi_1.
\end{equation}

\subsection{General case and example}\label{int_repr_gen} We evaluate $V(f;T;x_1,x_2)={\mathbb E}[f_+(x_1+X_T, \max\{x_2, x_1+{\bar X}_T\})]$.
The first term on the RHS of \eq{tVq0} is calculated at the end of Section \ref{ss:WHF_gen}.  Below, we consider in detail the evaluation of the second term on the RHS of \eq{tVq0}.
The proof of the following theorem is a straightforward modification of the proof for the joint cpdf. We assume that either $q>0$ or the condition in Lemma \ref{q_xi_cone} (a) holds,
 and the parameters of the deformations are chosen so that $q+\psi_{st}(\xi_1), q+\psi_{st}(\eta)\not\in (-\infty,0]$ for all $(q,\xi_1,\eta)$ arising in the process of deformation of all contours of integration.
\begin{thm}\label{thm:gen_stable}
Let $x_1\le x_2$, and let the following conditions hold
\begin{enumerate}[(i)]
\item
$\exists$\ $\alpha_0>-\alpha$, $ \delta<\delta_-$, $\gamma_+\in (0,\pi/2)$ such that $\forall$\
$\xi_1\in\{\xi_1\ |\  \mathrm{arg}\,\xi_1\in [0,\gamma_+]\cup[\pi,\pi-\gamma_+]\}$,
\begin{equation}\label{boundf1}
|e^{ix_2\xi_1}\widehat{(f_+)}(\xi_1,x_2)|\le C_1(x_2)(|\xi_1|^{\alpha_0-1}{\bf 1}_{|\xi_1|\le 1}+|\xi_1|^{-1+\delta}{\bf 1}_{|\xi_1|\ge 1}),
\end{equation}
where $C_1(x_2)$ depends on $x_2$ but not on $\xi_1$;
\item
$\exists$\ $\alpha_0>-\alpha$, $ \delta<\delta_+$, $\gamma_-\in (-\pi/2,0)$ such that $\forall$\
$\eta\in\{\eta\ |\ \mathrm{arg}\,\eta\in [\gamma_-,0]\cup[-\pi,-\pi-\gamma_-]\}$,
\begin{equation}\label{boundw1}
\left|e^{ix_2\eta}\int_{x_2}^{+\infty} dy\, e^{iy\eta}(({\mathcal E}^-_{q,mod}\otimes I)f_+)(y,y)\right|\le C_2(x_2)(|\eta|^{\alpha_0-1}{\bf 1}_{|\eta|\le 1}+|\eta|^{-1+\delta}{\bf 1}_{|\eta|\ge 1}),
\end{equation}
where $C_2(x_2)$ depends on $x_2$ but not on $\eta$.
\end{enumerate}
Then, for any $x_1\le x_2$, $\om_-\in (\gamma_-,0)$ and $\om_+\in (0,\gamma_+)$, the second term on the RHS of \eq{tVq0}
admits the representation
\begin{eqnarray}\label{eq:gen_stable}
&&({\mathcal E}^+_{q,mod}{\bf 1}_{(x_2,+\infty)}(y)((({\mathcal E}^-_{q,mod}\otimes I)f_+)(y,y)-
(({\mathcal E}^-_{q,mod}\otimes I)f_+)(y,x_2)))(x_1)\\\nonumber
&=&\frac{1}{2\pi}\int_{{\mathcal L}^-_{\om_-}} d\eta\, e^{i(x_1-x_2)\eta}\phi^+_{q,mod}(\eta)
\int_{x_2}^{+\infty} dy\, e^{i(x_2-y)\eta}\\\nonumber
&&\cdot\frac{1}{2\pi}\int_{{\mathcal L}^+_{\om_+}} e^{iy\xi_1}\phi^-_{q,mod}(\xi_1)
(\widehat{(f_+)}_1(\xi_1,y)-\widehat{(f_+)}_1(\xi_1,x_2)) d\xi_1.
\end{eqnarray}
\end{thm}

\begin{example}\label{ex:exchange_stable}{\rm For $\beta>1$, $\lambda>0$, set
$f(x_1,x_2):=f(\lambda,\beta;x_1,x_2)=(\beta x_1-x_2)_+ {\bf 1}_{x_1\le x_2}e^{-\lambda x_2}$. We introduce the dampening factor $e^{-\lambda x_2}$ so that
$V(f(\lambda,\beta;\cdot,\cdot);T;x_1,x_2)$ be finite for any $\alpha\in (0,2)$. By inspection of the integral representation for ${\tilde V}(f,q,x_1,x_2)$
that we derive, it is seen that if $\alpha>1$, the integrands are uniformly bounded by an absolutely integrable
function independent of $\lambda$, and converge point-wise as $\lambda\to0$. Hence,
the limits ${\tilde V}(f(0+,\beta;\cdot,\cdot); q,x_1,x_2)$ and $V(f(0+,\beta;\cdot,\cdot); T,x_1,x_2)$  exist and and are given by
the same integrals with $\lambda=0$.

The function
\begin{eqnarray*}
\widehat{(f_+)}_1(\xi_1,x_2)&=&e^{-\lambda x_2}\int_{x_2/\beta}^{x_2} e^{-ix_1\xi_1}(\beta x_1-x_2)dx_1\\
&=&\frac{e^{-\lambda x_2}}{-i\xi_1}\left[(\beta-1)x_2e^{-ix_2\xi_1}-\beta\int_{x_2/\beta}^{x_2} e^{-ix_1\xi_1}d\xi_1\right]\\
&=&e^{-\lambda x_2}\left([(-i\xi_1)^{-1}(\beta-1)x_2+\xi_1^{-2}\beta] e^{-ix_2\xi_1}-\beta\xi_1^{-2}e^{-ix_2\xi_1/\beta}\right)
\end{eqnarray*}
satisfies \eq{boundf1} with $\alpha_0=1, \delta=0$.
For $y>x_2$, $e^{i(y-x_2)\xi_1}$ is uniformly bounded in
the upper half-plane, and exponentially decays as $\xi_1\to \infty$ along any ray $e^{i\om_+}{\mathbb R}_+$, $\om_+\in (0,\pi)$. Since
$1/\beta<1$, $e^{i(y-x_2/\beta)\xi_1}$ enjoys the same property since $x_2\ge 0$.
Let $q>0$. We take $\om_+\in (0,\pi/2)$, and calculate
\begin{eqnarray*}
&&{\mathcal F}_{y\to\eta}{\bf 1}_{(x_2,+\infty)}(y)(({\mathcal E}^-_{q,mod}\otimes I)f_+)(y,x_2)\\
&=&\int_{x_2}^{\infty}dy\, e^{-(\lambda+i\eta)y}\frac{1}{2\pi}\int_{\mathbb R} \phi^-_{q,mod}(\xi_1)\left[e^{i(y-x_2)\xi_1}
\left(\frac{(\beta-1)x_2}{-i\xi_1}+\frac{\beta}{\xi_1^2}\right)-e^{i(y-x_2/\beta)\xi_1}\frac{\beta}{\xi_1^2}\right] d\xi_1\\
&=&\int_{x_2}^{\infty}dy\, e^{-(\lambda+i\eta)y}\frac{1}{2\pi}\int_{{\mathcal L}^+_{\om_+}} \phi^-_{q,mod}(\xi_1)\left[e^{i(y-x_2)\xi_1}
\left(\frac{(\beta-1)x_2}{-i\xi_1}+\frac{\beta}{\xi_1^2}\right)-e^{i(y-x_2/\beta)\xi_1}\frac{\beta}{\xi_1^2}\right] d\xi_1
\end{eqnarray*}
For $\eta\in \{\Im\eta\le 0, \eta\neq 0\}$, the integrand decays exponentially w.r.t. $\xi_1$ and $y$, hence, the double integral is absolutely convergent. Applying Fubini's theorem, and calculating the integral w.r.t. $y$, we continue
\begin{eqnarray*}
&=&\frac{e^{-(\lambda+i\eta)x_2}}{2\pi}\int_{{\mathcal L}^+_{\om_+}} \frac{\phi^-_{q,mod}(\xi_1)}{\lambda+i(\eta-\xi_1)}\left[
\left(\frac{(\beta-1)x_2}{-i\xi_1}+\frac{\beta}{\xi_1^2}\right)-e^{ix_2(1-1/\beta)\xi_1}\frac{\beta}{\xi_1^2}\right] d\xi_1.
\end{eqnarray*}
Similarly,
\begin{eqnarray*}
&&{\mathcal F}_{y\to\eta}{\bf 1}_{(x_2,+\infty)}(y)(({\mathcal E}^-_{q,mod}\otimes I)f_+)(y,y) \\
&=&\int_{x_2}^{\infty}dy\, e^{-(\lambda+i\eta)y}\frac{1}{2\pi}\int_{{\mathcal L}^+_{\om_+}} \phi^-_{q,mod}(\xi_1)\left[
\left(\frac{(\beta-1)y}{-i\xi_1}+\frac{\beta}{\xi_1^2}\right)-e^{iy(1-1/\beta)\xi_1}\frac{\beta}{\xi_1^2}\right] d\xi_1\\
&=&\frac{e^{-(\lambda+i\eta) x_2}}{2\pi}\int_{{\mathcal L}^+_{\om_+}} \phi^-_{q,mod}(\xi_1)\left[
\left(\left(\frac{1}{\lambda+i\eta}-\frac{1}{(\lambda+i\eta)^2}\right)\frac{(\beta-1)}{-i\xi_1}+\frac{1}{\lambda+i\eta}\frac{\beta}{\xi_1^2}\right)
\right.\\&& \hskip2cm\left.
-\frac{\beta e^{ix_2(1-1/\beta)\xi_1}}{\xi_1^2(\lambda+i(\eta-(1-1/\beta)\xi_1))}\right] d\xi_1.
\end{eqnarray*}
Substituting into \eq{eq:gen_stable}, we obtain
\begin{eqnarray}\label{eq:gen_stable_example}
&&({\mathcal E}^+_{q,mod}{\bf 1}_{(x_2,+\infty)}(y)((({\mathcal E}^-_{q,mod}\otimes I)f_+)(y,y)-
(({\mathcal E}^-_{q,mod}\otimes I)f_+)(y,x_2)))(x_1)\\\nonumber
&=&\frac{e^{-\lambda x_2}}{(2\pi)^2}\int_{{\mathcal L}^-_{\om_-}} d\eta\, e^{i(x_1-x_2)\eta}\phi^+_{q,mod}(\eta)
\int_{{\mathcal L}^+_{\om_+}}d\xi_1\,\phi^-_{q,mod}(\xi_1)\\\nonumber
&&\cdot\left[\left(\frac{1}{\lambda+i\eta}-\frac{1}{(\lambda+i\eta)^2}\right)\frac{(\beta-1)}{-i\xi_1}-\frac{(\beta-1)x_2}{-i\xi_1(\lambda+i(\eta-\xi_1)}
+\frac{\beta}{i\xi_1(\lambda+i\eta)(\lambda+i(\eta-\xi_1))}\right.\\\nonumber
&&\left. +\frac{\beta e^{ix_2(1-1/\beta)\xi_1}(1-1/\beta)}{-i\xi_1(\lambda+i(\eta-(1-1/\beta)\xi_1)(\lambda+i\eta)}\right]\\\nonumber
&=&\frac{e^{-\lambda x_2}}{(2\pi)^2}\int_{{\mathcal L}^-_{\om_-}} d\eta\, e^{i(x_1-x_2)\eta}\phi^+_{q,mod}(\eta)
\int_{{\mathcal L}^+_{\om_+}}d\xi_1\,\frac{\phi^-_{q,mod}(\xi_1)}{-i\xi_1}\\\nonumber
&&\cdot\left[\left(\frac{1}{\lambda+i\eta}-\frac{1}{(\lambda+i\eta)^2}\right)(\beta-1)-\frac{(\beta-1)x_2}{\lambda+i(\eta-\xi_1)}
-\frac{\beta}{(\lambda+i\eta)(\lambda+i(\eta-\xi_1))}\right.\\\nonumber
&&\left. +\frac{\beta e^{ix_2(1-1/\beta)\xi_1}(1-1/\beta)}{(\lambda+i(\eta-(1-1/\beta)\xi_1)(\lambda+i\eta)}\right]\
\end{eqnarray}

}
\end{example}

\section{Numerical evaluation of the cpdf of the supremum process and joint cpdf}\label{numer_real}

\subsection{Cpdf of the supremum process}\label{ss:supremum_numer:main_blocks}
Evaluation of  ${\mathbb E}[{\bf 1}_{x+{\bar X}_T\le a}]$ is reducible to evaluation of \eq{cpdfsup4} (see \eq{cpdfsup0}). The integral on
the RHS of  \eq{cpdfsup4} is evaluated  using the simplified trapezoid rule exactly as the outer integral in the case of the joint cpdf is evaluated in Sect. \ref{ss:joint_stable_numer}.  The Laplace inversion is performed exactly as in the case of the joint cpdf.

 \subsection{Joint cpdf of $X_T$ and ${\bar X}_T$: numerical realization based on GWR algorithm}\label{ss:joint_stable_numer}
 We use $V(T,a_1,a_2;x_1,x_2)={\mathbb E}[{\bf 1}_{(-\infty,a_1]}(x_1+X_T)]-V_1(T,a_1,a_2;x_1,x_2)$. Detailed algorithms for and numerous examples of numerical evaluation of ${\mathbb E}[{\bf 1}_{(-\infty,a_1]}(x_1+X_T)]$ can be found in
 \cite{ConfAccelerationStable}, therefore, it suffices to provide an algorithm for numerical evaluation of $V_1$; the Laplace transform
 $ {\tilde V}_1(q,a_1,a_2;x_1,x_2)$ is given by \eq{tV1q}. We proved \eq{tV1q} assuming that
  either $\alpha\neq 1$ or $\alpha=1$ and $c_-=c_+$. If $\alpha=1$ and $c_+\neq c_-$, then one can use \eq{tV1q} with $\phi^\pm_q$ instead
  of the modified $\phi^\pm_{q,mod}$ but then the convergence of the integrals and infinite sums in the simplified trapezoid rule
  are worse, longer grids are needed, and errors are more difficult to control.

 Set $\om_+=\pi/8, \om_-=-\pi/8$, construct grids
  \begin{eqnarray}\label{plus:grids}
 \vec{y^{+}}&=&\zeta^+*(-N^+_-, N^+_+),\ \vec{\xi^{+-}}=\exp[i*(\pi-\om_+)+ \vec{y^{+}}],\ \vec{\xi^{++}}=\exp[i*\om_++ \vec{y^{+}}],
 \\\label{minus:grids}
 \vec{y^{-}}&=&\zeta^-*(-N^-_-, N^-_+),\ \vec{\xi^{--}}=\exp[i*(-\pi-\om_-)+ \vec{y^{-}}],\ \vec{\xi^{-+}}=\exp[i*\om_-+ \vec{y^{-}}];
 \end{eqnarray}
and
 calculate the arrays $\psi_{st}(\vec{\xi^{+-}}), \psi_{st}(\vec{\xi^{++}}), \psi_{st}(\vec{\xi^{--}}), \psi_{st}(\vec{\xi^{-+}})$.

 Then,
 for each $q>0$ used in GWR algorithm, calculate
 \begin{enumerate}[(I)]
  \item
 the Wiener-Hopf factors ${\phi^-_q}(\vec{\xi^{--}})$ and ${\phi^-_q}(\vec{\xi^{-+}})$ and $a^-_q\ge 0$ using the formulas in Lemma \ref{asympWHF}; the integrals over the rays $e^{i\om_+}{\mathbb R}_+$ and $e^{i(\pi-\om_+)}{\mathbb R}_+$ are calculated using the exponential changes
 of variables and simplified trapezoid rule, the grids \eq{plus:grids} being used;
 \item
 the Wiener-Hopf factors ${\phi^+_q}(\vec{\xi^{+-}})$ and ${\phi^+_q}(\vec{\xi^{++}})$ and $a^+_q\ge 0$ using the formulas in Lemma \ref{asympWHF}; the integrals over the rays $e^{i\om_-}{\mathbb R}_+$ and $e^{i(-\pi-\om_-)}{\mathbb R}_+$ are calculated using the exponential changes
 of variables and simplified trapezoid rule, the grids \eq{minus:grids} being used;
 \item
 arrays ${\phi^+_q}(\vec{\xi^{--}})$, ${\phi^+_q}(\vec{\xi^{-+}})$ and ${\phi^-_q}(\vec{\xi^{+-}})$, ${\phi^-_q}(\vec{\xi^{++}})$ using the arrays
 in (I) and (II) and the Wiener-Hopf identity;
 \item
 arrays  $\phi^+_{q,mod}(\vec{\xi^{--}})$, $\phi^+_{q, mod}(\vec{\xi^{-+}})$,
 $\phi^-_{q,mod}(\vec{\xi^{+-}})$ and $\phi^-_{q, mod}(\vec{\xi^{++}})$;
 \item
 making the exponential changes of variables on each of four arrays comprising the contours in the double integral on the right-most side  of
 \eq{tV1q} and applying the simplified trapezoid rule to each of 4 integrals, calculate ${\tilde V}_1(q,a_1,a_2;x_1,x_2)$.
  \end{enumerate}
  Finally, apply GWR algorithm to calculate $V_1(T,a_1,a_2;x_1,x_2)$.

 \subsection{Sinh-acceleration in the Bromwich integral}\label{ss:SINH-Bromwich}
Let the condition in Lemma \ref{q_xi_cone} (a) hold. Then we find $\sigma>0, \gamma_0\in (0,\pi/2)$ and $\gamma_-<0<\gamma_+$ as indicated in Lemma \ref{lem:anal_cont_WHF}.
Next, we choose  $\sigma_\ell>\sigma$,
$\omega_\ell\in (0,\gamma_0)$, and $b_\ell>0$ such that  $\sigma_\ell-b_\ell\sin\omega_\ell>0$ and define the  map
$\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}$ by \eq{eq:sinhLapl}. In the Bromwich integral for $V_1(T, a_1,a_2;x_1,x_2)$, we make the change
of variables $q=q(y')=\chi_{L; \sigma_\ell,b_\ell,\omega_\ell}(y')$, and choose the grid $\vec{y'}=(j\zeta)_{j=-N}^N$ for the simplified trapezoid rule.
We calculate the factors ${\tilde V}_1(q(j\zeta),a_1,a_2;x_1,x_2)$ in the approximation
\begin{equation}\label{sinhBrCPDF}
V_1(T,a_1,a_2;x_1,x_2)\approx \frac{b_\ell\zeta}{\pi}\Re \sum_{j=0}^N e^{q(j\zeta)T}\cosh(i\omega_\ell+j\zeta))
{\tilde V}_1(q(j\zeta),a_1,a_2;x_1,x_2)(1-\delta_{j0}/2)
\end{equation}
as for $q>0$ above using $\om_+=\gamma_+/2$ and $\om_-=\gamma_-/2$ and Lemma \ref{lem:anal_cont_WHF}
instead of Lemma \ref{asympWHF}.

 \subsection{Numerical examples}\label{ss:numer_ex}
 The calculations in the paper
were performed in MATLAB 2017b-academic use, on a MacPro Chip Apple M1 Max Pro chip
with 10-core CPU, 24-core GPU, 16-core Neural Engine 32GB unified memory,
1TB SSD storage.

The CPU times shown in the tables in Sect. \ref{s:tables} can be significantly improved
using
parallelization, especially in the blocks where calculations are performed for each $q$ in the Laplace inversion procedure, and
 asymptotic expansions and summation by parts (similarly to the ones used to evaluate cpdf of $X_T$ in Sect. \ref{s:1D} and \cite{ConfAccelerationStable}) in the blocks where the Wiener-Hopf factors. Note that when the cpdf of ${\bar X}_T$ is evaluated, the CPU time for the evaluation of the Wiener-Hopf factors is 80\% of the total CPU time, approximately.
As one would expect, the length of arrays and CPU time depends on parameters of the process, $T$, the barrier $a_2$ and distance to barrier $a_2-a_1$. Hence, the comprehensive list of examples would be extremely long; we choose examples to show that if
the parameters of the process, $T$ and $a_2$, $a_2-a_1$ are in a not very favorable region of the parameter space, the method
of the paper can achieve good precision fairly fast. If the sinh-acceleration is applicable to the Bromwich integral,  we achieve the precision
better than E-15 in seconds in the case of the joint distribution, and in a fraction of a second for cpdf of the supremum process.
The precision of the order E-10 is achievable in hundreds and dozens of milliseconds, respectively. When GWR acceleration is applied, we fix the order $2M=16$ and choose other parameters of the numerical scheme so that the order of the errors does not improve if we choose longer and finer grids. Thus, the errors shown are, essentially, the errors of GWR acceleration itself. If $\alpha<1$ and $\mu\neq 0$, then we can apply only GWR acceleration, and we show the differences between the results of two schemes with the integration over different curves.

As in the case of the evaluation of cpdf of $X_T$, the length of grids and CPU time needed to satisfy a chosen error tolerance greatly
increases if $\beta\neq 0$ and $\alpha$  is close to 0 or 1; the effect depends on the sign of $\beta$. In the examples below, $\alpha=0.2$ and
$\alpha=1.2$ are
close but not very close to 0 and 1, and $\beta=-0.2$ and $\mu=-0.02$, $\mu=0$ and $\mu=0.02$ are not large in absolute value.
In the majority of tables, $T=0.25$ is moderate, and the upper boundary $a$ for the supremum process is small, in the range $[0.0125, 0.075]$.
In Tables 4-6, $T=0.004, 1, 10$; in Tables 5-6, $a\in [100,600]$ and $a\in [0.5,3]$, respectively. The results for the joint cpdf (Tables 6 and 7) show that
$V_1$ can be calculated more accurately than cpdf of the supremum process using the same grids although the CPU time in the former case (triple integral) is much larger than in the latter case (double integral). The reason is that $V_1$ is smoother than cpdf of the supremum process.

 \section{Conclusion}\label{concl}
 In Sect. \ref{s:1D}, we generalized the  efficient numerical method
 for evaluation of pdfs and cpdfs of stable distributions developed in \cite{ConfAccelerationStable}
 to the case of expectations ${\mathbb E}[f(x+X_T)$. We also suggested
 a new method, which can be used for if the parameters of the distribution are in the regions where the method of  \cite{ConfAccelerationStable} requires extremely  grids.
In the following sections, we derived integral representations for expectations of functions $f(X_T, {\bar X}_T)$ of a stable L\'evy process and the corresponding supremum process and design efficient numerical procedures for the evaluation of the resulting double and triple integrals.
 The main block is the integral representation for the Laplace transform ${\tilde V}(q)$ of the expectation w.r.t. $T$;
 the form of the representation depends on the properties of $f$.
 The integral representations are  analogs of the representations derived in \cite{EfficientLevyExtremum}
 for L\'evy processes with exponentially decaying tails of the density of jumps.
As in  \cite{EfficientLevyExtremum}, the resulting formulas are in terms
 of the Wiener-Hopf factors $\phi^\pm_q(\xi)$. Since the characteristic exponent $\psi$ of a stable L\'evy process is non-smooth at 0,
 the justification of the integral representations in the present paper are more involved, and a novel regularization of the Wiener-Hopf factors
 in the resulting formulas is needed. The regularization relies on a sufficiently regular behavior of $\phi^\pm_q(\xi)$ as $\xi\to\infty$.
 In the case of asymmetric processes of index 1, we were unable to derive the asymptotics of the Wiener-Hopf factors, hence,
 we have neither a rigorous proof of the integral representations in terms of $\phi^\pm_q(\xi)$ nor  representations
 in terms of the modified Wiener-Hopf factors. The general theory of pseudo-differential operators (see, e.g., \cite{eskin,NG-MBS})
 can be used to justify the integral representations in the sense of generalized functions, and study the regularity of the integral representations. The numerical schemes that we develop can be used in this case as well but the convergence of the integrals and infinite sums can be poor without the regularization.

 Numerical realizations of the formulas for the Wiener-Hopf factors and integral representations are based on appropriate rotations of
 rays of integration, equivalently, on changes of variables of the form $\xi=e^{i\omega+y}$; in \cite{EfficientLevyExtremum}, conformal deformations of the lines of integration and the sinh-changes of variables of the form $x\i=i\omega_1+b\sinh(i\omega+y)$ were used.
In repeated integrals, the changes of variables must be in a certain agreement.
 As in  \cite{EfficientLevyExtremum}, after the changes of the variables are made, the simplified trapezoid rule is applied.
In  \cite{EfficientLevyExtremum} and the present paper, the complexity of the scheme of the evaluation of the Wiener-Hopf factors is of the order of $O((\ln(1/\epsilon))^2)$, where $\epsilon$ is the error tolerance, and  the complexity of the scheme of the evaluation
 of ${\tilde V}(q)$ is of the order $O((\ln(1/\epsilon))^{2m})$, where $m$ is the number of repeated integrations.

 As in \cite{EfficientLevyExtremum}, we consider two types of the Laplace inversion: the Gaver-Wynn Rho acceleration and
 the sinh-acceleration applied to the Bromwich integral. The latter can be applied in the index $\alpha>1$ or $\alpha=1$ in the symmetric case or $\alpha<1$ and zero drift. Using the sinh-acceleration, it is possible to achieve the precision of the order E-15 in a fraction of a second if
 $V$ is represented as a double integral, and in several seconds if $V$ is represented as a triple integral; the precision of the order of $E-10$ can be achieved in dozens of milliseconds and a fraction of a second, respectively. The scheme based on the Gaver-Wynn Rho acceleration is not so accurate but simpler to implement because for $q>0$, the deformations of the contours of deformations
 enjoy better properties, and the parameters of the deformations are much easier to choose. The reason is the following useful property
 established in \cite{EfficientAmenable} for a wide class of L\'evy processes called Stieltjes-L\'evy processes (SL processes). Namely,
 if $X$ is a SL process, then $q+\psi(\xi)\neq 0$ for all $\xi\not\in i{\mathbb R}$. It is proved in  \cite{EfficientAmenable} that essentially all popular classes of L\'evy processes, stable ones including, are SL processes.

 The numerical examples in the paper are provided for the cpdf of the supremum process and joint cpdf of the stable L\'evy process and supremum process. As it is well-known (see, e.g., \cite{ConfAccelerationStable} for an extensive list of examples), for certain regions in the parameter space (relatively small ones), the evaluation even of cpdf of the stable L\'evy process is very difficult and requires additional tricks. Similar difficulties are expected in the cases considered in the paper. However, as numerical examples in the paper demonstrate,
 if the set of parameters is not too close to the bad region in the parameter space, the method of the paper works well.