For Lévy processes with exponentially decaying tails of the Lévy density, we derive integral representations for the joint cpdf $V$ of $(X_T, \bar X_T,τ_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.