  For L\'evy processes with exponentially decaying tails of the L\'evy density,
we derive integral representations for the joint cpdf $V$ of $(X_T, \bar
X_T,\tau_T)$ (the process, its supremum evaluated at $T<+\infty$, and the first
time at which $X$ attains its supremum). The first representation is a
Riemann-Stieltjes integral in terms of the (cumulative) probability
distribution of the supremum process and joint probability distribution
function of the process and its supremum process. The integral is evaluated
using a combination an analog of the trapezoid rule. The second representation
is amenable to more accurate albeit slower calculations. We calculate
explicitly the Laplace-Fourier transform of $V$ w.r.t. all arguments, apply the
inverse transforms, and reduce the problem to evaluation of the sum of 5D
integrals. The integrals can be evaluated using the summation by parts in the
infinite trapezoid rule and simplified trapezoid rule; the inverse Laplace
transforms can be calculated using the Gaver-Wynn-Rho algorithm. Under
additional conditions on the domain of analyticity of the characteristic
exponent, the speed of calculations is greatly increased using the conformal
deformation technique. For processes of infinite variation, the program in
Matlab running on a Mac with moderate characteristics achieves the precision
better than E-05 in a fraction of a second; the precision better than E-10 is
achievable in dozens of seconds. As the order of the process (the analog of the
Blumenthal-Getoor index) decreases, the CPU time increases, and the best
accuracy achievable with double precision arithmetic decreases.
