This is an electronic reprint of the original article published in American Journal of Computational Mathematics, 2023, Vol. 13, No. 4, 476–504, 10.4236/ajcm.2023.134026. This reprint differs from the original in reference style, pagination, and typographic details.

§ INTRODUCTION

Let $\Phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^x e^{-\frac{1}{2}t^2}\,\mathrm{d}t$ be the standard normal cumulative distribution function (cdf) and $\varphi(x)=\frac{1}{\sqrt{2\pi}}\;e^{-\frac{1}{2}x^2}$ the standard normal probability density function (pdf). The bivariate standard normal cdf

\begin{equation}\label{Phirxy1}
\Phi_\varrho^2(x,y)=\tfrac{1}{2\pi\,\sqrt{1-\varrho^2}}\;\int_{-\infty}^x \int_{-\infty}^y e^{-\frac{s^2-2\varrho st+t^2}{2\,(1-\varrho^2)}}\,\mathrm{d}s\;\mathrm{d}t
\end{equation}

with a correlation coefficient $\varrho\in(-1,1)$ has the bivariate standard normal pdf

\begin{equation}\label{phirxy}
\varphi_\varrho^2(x,y)=\tfrac{\partial^2\Phi_\varrho^2(x,y)}{\partial x\,\partial y}= \tfrac{1}{2\pi\,\sqrt{1-\varrho^2}}\;e^{-\frac{x^2-2\varrho xy+y^2}{2\,(1-\varrho^2)}}.
\end{equation}

If $\varrho=0$, then $\Phi_0^2(x,y)=\Phi(x)\,\Phi(y)$, while for $\abs{\varrho}=1$ we have the degenerated cases

\begin{equation}\label{Phi1xy}
\Phi_{\pm 1}^2(x,y)=
\begin{cases}
\min\{\Phi(x),\Phi(y)\},     & \mbox{if }\varrho=+1,\\
\max\{\Phi(x)+\Phi(y)-1,0\}, & \mbox{if }\varrho=-1.
\end{cases}
\end{equation}

Additionally, for $\varrho\in(-1,1)$ define

\begin{equation*}
L(x,y,\varrho)=\tfrac{1}{2\pi\,\sqrt{1-\varrho^2}}\,\int_x^{\infty}\int_y^{\infty} e^{-\frac{s^2-2\varrho st+t^2}{2\,(1-\varrho^2)}}\,\mathrm{d}s\;\mathrm{d}t,
\end{equation*}

$L(x,y,-1)=\max\{1-\Phi(x)-\Phi(y),0\}$ and $L(x,y,1)=1-\max\{\Phi(x),\Phi(y)\}$. The functions $\Phi_\varrho^2(x,y)$ and $L(x,y,\varrho)$ are symmetric regarding $x$ and $y$, and because they are complementary, it suffices to analyze only one of them. In this respect, $L$ is often preferred, but not in this paper where some equations for $L$ are also rewritten using $\Phi_\varrho^2(x,y)=L(-x,-y,\varrho)$ and some other auxiliary equations.

Noting that (<ref>) implies $\Phi_1^2(x,x)=\Phi(x)=2\,\Phi_0^2(x,0)$ and $\Phi_{-1}^2(x,-x)=0=\Phi_0^2(x,0)+\Phi_0^2(-x,0)-\frac{1}{2}$, the equations from <cit.>, which do not cover these cases if $x\ne 0$, as we get an indefinite expression $\frac{0}{0}$ for the auxiliary correlation coefficients $\varrho_x$ and $\varrho_y$, can be transformed into

\begin{equation}\label{Phirxy2}
\Phi_\varrho^2(x,y)=
\begin{cases}
\tfrac{1}{4}+\tfrac{1}{2\pi}\arcsin\varrho,            & \mbox{if }x=y=0,\\
\Phi_{\varrho_x}^2(x,0)+\Phi_{\varrho_y}^2(y,0)-\beta, & \mbox{if }x\ne 0\mbox{ or }y\ne 0,
\end{cases}
\end{equation}

where $\beta=0$ if $xy>0$ or $xy=0$ and $x+y\ge 0$, and $\beta=\frac{1}{2}$ otherwise,

\begin{equation}\label{Phirxy3}
\varrho_x=\left\{
\begin{array}{lcr}
0 & \mbox{if }\abs{\varrho}=1\mbox{ and }x\ne 0\mbox{ and }x-y\,\sgn(\varrho)=0 & 0\\
\tfrac{(\varrho x-y)\,\sgn(x)}{\sqrt{x^2-2\varrho xy+y^2}} & \mbox{otherwise} &
\tfrac{(\varrho y-x)\,\sgn(y)}{\sqrt{x^2-2\varrho xy+y^2}}
\end{array}
\right\}=\varrho_y,
\end{equation}

where $\sgn(x)=1$, if $x\ge 0$ and $\sgn(x)=-1$, if $x<0$. Note that although $\abs{\varrho}$ is not close to $1$, $\abs{\varrho_x}$ or $\abs{\varrho_y}$ can be, and that $\abs{\varrho}=1$ and the lower equations (<ref>) imply $\abs{\varrho_x}=1$ and $\abs{\varrho_y}=1$.

Since <cit.> imply

\begin{equation}\label{Phirx0}
\Phi_\varrho^2(-x,0)=\Phi_\varrho^2(x,0)-\Phi(x)+\tfrac{1}{2}\quad\text{and}\quad
\Phi_{-\varrho}^2(x,0)=\Phi(x)-\Phi_\varrho^2(x,0),
\end{equation}

and $\Phi_{-\varrho}^2(-x,0)=\frac{1}{2}-\Phi_\varrho^2(x,0)$, the problem how to compute $\Phi_\varrho^2(x,y)$ is reduced to the computation of $\Phi_\varrho^2(x,0)$ for $x\ge 0$ and $\varrho\in[0,1)$.

Extensive older bibliography exists about how $\Phi_\varrho^2(x,y)$ can be computed <cit.>. Many old methods, e.g. [Donnelly, 1973]'s ([Donnelly, 1973]) method based on [Owen, 1956]'s ([Owen, 1956]) results, [Divgi, 1979]'s ([Divgi, 1979]) method and the method of Drezner and Wesolowsky, 1990, are sufficiently effective and accurate for the majority of present needs, but new methods still appear, e.g. [Patefield and Tandy, 2000, Genz, 2004, Meyer, 2013, Fayed and Atiya, 2014]. More details and comparisons of some methods can be found in [Terza and Welland, 1991, Gai, 2001, Aǧca and Chance, 2003].

The Owen's $T$ function is defined by <cit.>

\begin{equation}\label{Owen1}
\begin{split}
T(h,a)&=\tfrac{1}{2\pi}\int_0^a\tfrac{e^{-\frac{1}{2}h^2(1+x^2)}}{1+x^2}\, \mathrm{d}x\\
&=\tfrac{\arctan a}{2\pi}-\tfrac{1}{2\pi}\sum_{k=0}^\infty\tfrac{(-1)^k a^{2k+1}}{2k+1}\left(1-e^{-\frac{1}{2}h^2}\sum_{i=0}^k\tfrac{h^{2i}}{2^i i!}\right)
\end{split}
\end{equation}

for $a,h\in\mathbb{R}$. It appears frequently in the extensive Owen's table of integrals involving Gaussian functions <cit.> and can be used for computing $\Phi_\varrho^2(x,y)$, as will be explained in the next section. The series (<ref>) is used in [Donnelly, 1973]'s ([Donnelly, 1973]) algorithm and, based on it, is also realized in pbinorm() function from the VGAM package [Yee, 2023]. Since it alternates, there are problems with numerical stability, e.g. using double precision arithmetic, it may provide only single precision results. More about Owen's $T$ function and older ways of calculating it are given by Patefield and Tandy, 2000, who developed a hybrid method to compute it. To achieve double precision accuracy, they use the series (<ref>) or one of five other methods. Which method is chosen and which its parameters are used depends on which of the $18$ regions $(h,a)$ belongs to <cit.>. Their method is also implemented in OwenT() function from the OwenQ package [Laurent, 2023], hereafter OwenQ::OwenT().

The method presented in this paper extends the theoretical knowledge and can be easily implemented. Its great advantage over the already mentioned methods, which enable $\Phi_\varrho^2(x,y)$ computation via Owen's $T$ function, is its simplicity, faster convergence and ensured numerical stability. Regarding other methods that have proven to be good or even excellent in terms of $\Phi_\varrho^2(x,y)$ computation speed and accuracy, e.g. the one implemented in the pmvnorm() function from the mvtnorm package [Genz and Bretz, 2009], the main advantages are simplicity and independence from various constants. Both enable easy portability between different environments and using high-precision arithmetic. There is no need to adjust the number and precision of various constants, such as weights and knots in quadrature formulas, which are used in many $\Phi_\varrho^2(x,y)$ computational methods.

In Section <ref>, the background results are presented, which are the starting point for the development of a new method, as well as for the acceleration of the tetrachoric series, which was used for the benchmark values computation. In Section <ref>, auxiliary results are derived. In core Section <ref>, new series for $\Phi_\varrho^2(x,0)$ are derived in Theorems <ref> and <ref>, and in Corollary <ref> they are combined in a unified series for $\Phi_\varrho^2(x,y)$. However, in all cases, the series must be supplemented by the standard normal cdf and some of them also by the arctangent function computation. Four series for computing the Owen's $T$ function are presented in Theorem <ref>. Among them (<ref>) can be viewed as a modified Euler's arctangent series. In Sections <ref> and <ref>, algorithmic and numeric issues are discussed, and testing results are presented.

Using double precision arithmetic and testing on a random sample of one million triplets $(x,y,\varrho)$ with uniformly distributed components, the maximum absolute error was $3.45\cdot 10^{-16}$, which is approximately $3.11$-times the machine epsilon $2^{-53}$. On the same random sample, but with $\varrho$ transformed so that more than half of $\abs{\varrho}$ were greater than $0.9999$, and by equivalent computing of two function values instead of $\Phi_\varrho^2(x,y)$ in rare cases with $\varphi_\varrho^2(x,y)>1$, the maximum absolute error was even slightly smaller.

All new series are suitable for vector and high-precision computation. The fastest one asymptotically converges as $\sum_k\frac{q^{k+1}e^{-q}}{(k+1)!}\,\frac{\min\{a^{2k},1\}}{\sqrt{k+1}\;(1+a^2)^k}$, where $q=\frac{1}{2}(1+a^2)h^2$. Using the Rmpfr package [Maechler, 2023] and $1024$-bit precision computation, $199$ iterations were needed to compute $\Phi_{\sqrt{2}/2}^2(2.1,0)$, achieving an absolute error smaller than $2^{-1024}\approx 5.56\cdot 10^{-309}$.

§ BACKGROUND RESULTS

The bivariate standard normal cdf can be expanded into tetrachoric series

\begin{equation}\label{TetrachA}
\Phi_\varrho^2(x,y)=\Phi(x)\,\Phi(y)+\varphi(x)\,\varphi(y)
\sum_{k=0}^\infty H\!e_k(x)\,H\!e_k(y)\,\tfrac{\varrho^{k+1}}{(k+1)!},
\end{equation}

where $H\!e_k(x)$ are the “probabilist's” Hermite polynomials defined by

\begin{equation*}
\end{equation*}

The tetrachoric series (<ref>) converges a bit faster than a geometric series with a quotient $\varrho$ and slowly converges for $\abs{\varrho}$ close to or equal to $1$.

Let $\mathbb{N}_{-1}=\mathbb{N}\cup\{-1,0\}$. Extending the usual double factorial function definition for $k\in\mathbb{N}_0$ to $k\in\mathbb{N}_{-1}$ with $(-1)!!=1$, we have

\begin{equation*}
\begin{cases}
0,                       & \mbox{if $k$ is odd},\\
(-1)^\frac{k}{2}(k-1)!!, & \mbox{if $k$ is even}.
\end{cases}
\end{equation*}

Setting $y=0$ in (<ref>), transforming the summation index $\frac{k}{2}\mapsto k$, and using $(2k-1)!!=\frac{(2k)!}{2^k k!}$ for $k\in\mathbb{N}_0$, we get a slightly faster converging series

\begin{equation}\label{TetrachB}
\Phi_\varrho^2(x,0)=\tfrac{1}{2}\,\Phi(x)+\tfrac{\varphi(x)}{\sqrt{2\pi}}
\sum_{k=0}^\infty H\!e_{2k}(x)\,\tfrac{(-1)^k\varrho^{2k+1}}{2^k k!\,(2k+1)}.
\end{equation}

<cit.> found an alternative series for the computation of $\Phi_\varrho^2(x,0)$, which is a good alternative to the series (<ref>) if $\varrho^2>\frac{1}{2}$.

Integrating the well known equation <cit.>

\begin{equation}\label{Phirxy4}
\tfrac{\partial}{\partial\varrho}\,\Phi_\varrho^2(x,y)=
\tfrac{\partial^2\Phi_\varrho^2(x,y)}{\partial x\,\partial y}=\varphi_\varrho^2(x,y)
\end{equation}

with respect to $\varrho$, using equations $\Phi_0^2(x,y)=\Phi(x)\,\Phi(y)$ and (<ref>), we get

\begin{equation}\label{Phirxy5}
\begin{split}
\Phi_\varrho^2(x,y)&=\Phi(x)\,\Phi(y)+\int_0^\varrho\varphi_s^2(x,y)\, \mathrm{d}s\\
\end{split}
\end{equation}

Due to historical reasons, the variates $h$ and $k$ are often used instead of $x$ and $y$. We will use $h$ instead of $x$ but will not need $k$ in the role of $y$.

In the sequel, let $\varrho\in(-1,1)$ and $r\in\mathbb{R}$ be permanently connected by the one-to-one correspondence $r=\frac{\varrho}{\sqrt{1-\varrho^2}}$ and $\varrho=\frac{r}{\sqrt{1+r^2}}$. Notice that (<ref>) and $\arcsin\varrho=\arctan r$ imply $\Phi_\varrho^2(0,0)=\frac{1}{4}+\frac{1}{2\pi}\arctan r$. Setting $x=h$, $y=0$, and writing $e^{-\frac{1}{2}h^2}\int_a^b e^{\frac{1}{2}h^2}\varphi_s^2(h,0)\,\mathrm{d}s$ instead of $\int_a^b\varphi_s^2(h,0)\,\mathrm{d}s$, the equations (<ref>) can be rewritten to

\begin{equation}\label{Phirh0a}
\begin{split}
\Phi_\varrho^2(h,0)
\end{split}
\end{equation}

where, using a substitution $s=\frac{x}{\sqrt{1+x^2}}$,

\begin{equation*}
e^{\frac{ws^2}{1-s^2}}\,\mathrm{d}s=\tfrac{1}{\sqrt{2\pi}} \int_0^r\tfrac{e^{wx^2}}{1+x^2}\,\mathrm{d}x\quad(w\in\mathbb{R}).
\end{equation*}

Analogously, $\widetilde{T}_r(w)$ and $\overline{T}_r(w)$ are defined for $w\le 0$ by integrals with limits $(\varrho,1)$ and $(-1,\varrho)$, and $(r,\infty)$ and $(-\infty,r)$ respectively, which do not exist for $w>0$. The computation of $\Phi_\varrho^2(h,0)$ is reduced to the computation of one of the integrals $T_r(-\frac{1}{2}h^2)$, $\widetilde{T}_r(-\frac{1}{2}h^2)$ and $\overline{T}_r(-\frac{1}{2}h^2)$.

The special cases (<ref>) of the equations (<ref>) in combination with (<ref>) are a starting point for some computation methods for $\Phi_\varrho^2(x,y)$ which are based on the integration. Owen's is one of them – note that $\varphi(h)\,T_a(-\frac{1}{2}h^2)=T(h,a)$. Using $r$ instead of $a$ in the sequel and resolving the problem with $\frac{0}{0}$ as before, Owen's equation <cit.> can be rewritten as

\begin{equation}\label{Owen2}
\Phi_\varrho^2(x,y)=
\begin{cases}
\tfrac{1}{4}+\tfrac{1}{2\pi}\arcsin\varrho,            & \mbox{if }x=y=0,\\
\tfrac{1}{2}(\Phi(x)+\Phi(y))-T(x,r_x)-T(y,r_y)-\beta, & \mbox{if }x\ne 0\mbox{ or }y\ne 0,
\end{cases}
\end{equation}

where $\beta=0$ if $xy>0$ or $xy=0$ and $x+y\ge 0$, and $\beta=\frac{1}{2}$ otherwise, and

\begin{equation}\label{Owen3}
\begin{array}{lcr}
0 & \mbox{if }\abs{\varrho}=1\mbox{ and }x\ne 0\mbox{ and }x-y\,\sgn(\varrho)=0 & 0\\
\tfrac{y-\varrho x}{x\,\sqrt{1-\varrho^2}} & \mbox{otherwise} &
\frac{x-\varrho y}{y\,\sqrt{1-\varrho^2}}
\end{array}
\right\}=r_y.
\end{equation}

Another Owen's important equation is <cit.>

\begin{equation}\label{Owen4}
\end{equation}

where $\beta=0$ if $r\ge 0$ and $\beta=\frac{1}{2}$ if $r<0$. It implies $\varphi(h)\,T_1(-\frac{1}{2}h^2)=T(h,1)=\frac{1}{2}\Phi(h)\,(1-\Phi(h))$.

If $\varrho=\frac{\sqrt{2}}{2}$, then $r=1$. Using (<ref>) and the right equation (<ref>), we get

\begin{equation}\label{Phirh0b}
\Phi_{\sqrt{2}/2}^2(h,0)=\Phi(h)(1-\tfrac{1}{2}\,\Phi(h))\quad\text{and}\quad \Phi_{-\sqrt{2}/2}^2(h,0)=\tfrac{1}{2}(\Phi(h))^2.
\end{equation}

Let $\varrho\in(-1,1)$, $\overline{\varrho}=\sgn(\varrho)\sqrt{1-\varrho^2}$,

$r=\frac{\varrho}{\sqrt{1-\varrho^2}}$, $\overline{r}=\frac{\overline{\varrho}}{\sqrt{1-\overline{\varrho}^2}}=\frac{1}{r}$ and $h\in\mathbb{R}$. By adding $\frac{1}{2}(\Phi(h)+\Phi(rh))$ on both sides of (<ref>) and using (<ref>) and (<ref>), it follows

\begin{equation}\label{Owen5}
\Phi_\varrho^2(h,0)+\,\Phi_{\overline{\varrho}}^2(rh,0)=\Phi(h)+\Phi(rh)- \Phi(h)\,\Phi(rh)-\beta,
\end{equation}

where $\beta=0$ if $r\ge 0$ and $\beta=\frac{1}{2}$ if $r<0$.

If $\abs{\varrho}>\frac{\sqrt{2}}{2}$, then $\abs{r}>1$, $\abs{\overline{r}}<1$ and $\abs{\overline{\varrho}}<\frac{\sqrt{2}}{2}$. $\Phi_\varrho^2(h,0)$ can be computed by computing $\Phi_{\overline{\varrho}}^2(rh,0)$ and using (<ref>). Consequently, the right equation (<ref>) and the transformation of the parameters $rh\mapsto\overline{h}$ and $\frac{1}{r}\mapsto\overline{r}$, if needed, reduce the condition $\varrho\in(-1,1)$ to $\varrho\in[0,\frac{\sqrt{2}}{2}]$ with a corresponding $r\in[0,1]$. The problem how to compute $\Phi_\varrho^2(x,y)$ is reduced to the computation of $\Phi_\varrho^2(h,0)$ for $h\ge 0$ and $\varrho\in[0,\frac{\sqrt{2}}{2}]$.

The theorems in the sequel do not depend on the restriction to $\varrho^2\le\frac{1}{2}$ with a corresponding $\abs{r}\le 1$ found in [Patefield and Tandy, 2000]. However, if $\abs{r}>1$, the transformation of $rh\mapsto\overline{h}$ and $\frac{1}{r}\mapsto\overline{r}$ speeds up convergence significantly. This is, for example, also the case for the Vasicek's series and the tetrachoric series (<ref>), which transforms into

\begin{equation}\label{TetrachC}
\Phi_\varrho^2(h,0)=\Phi(h)+\tfrac{1}{2}\Phi(rh)-\Phi(h)\Phi(rh)
-\tfrac{\varphi(rh)}{\sqrt{2\pi}}\sum_{k=0}^\infty H\!e_{2k}(rh)\, \tfrac{(-1)^k\overline{\varrho}^{2k+1}}{2^k k!\,(2k+1)}-\beta,
\end{equation}

where $\overline{\varrho}=\sgn(\varrho)\sqrt{1-\varrho^2}$, $\beta=0$ if $r\ge 0$ and $\beta=\frac{1}{2}$ if $r<0$. Since the results calculated with the Vasicek's series in Section <ref> only supplement benchmark values calculated with the tetrachoric series, details about it are not given here.

§ AUXILIARY RESULTS

This paper pays attention to the estimation of the method error due to truncation of infinite series, while the analysis of rounding errors is out of its scope. The following lemma is intended for this. In it and all the theorems below, the phrase “truncated after the $n^{th}$ term” also means omitting the entire series. In these cases, set $n=-1$ and assume the empty sum is $0$.

If the Euler's series

\begin{equation}\label{Euler1}
\arctan r=\tfrac{r}{1+r^2}\sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!} \left(\tfrac{r^2}{1+r^2}\right)^k\quad(r\in\mathbb{R})
\end{equation}

is truncated after the $n^{th}$ term, then the remainder

\begin{equation*}
R_n=\tfrac{r}{1+r^2}\sum_{k=n+1}^\infty\tfrac{(2k)!!}{(2k+1)!!} \left(\tfrac{r^2}{1+r^2}\right)^k \quad(n\in\mathbb{N}_{-1})
\end{equation*}

has the same sign as $r$ and is bounded in absolute value by

\begin{equation}\label{Euler2}
B_n\le\left\lvert R_n\right\rvert\le(1+r^2)\,B_n,\quad\text{where}\quad B_n=\tfrac{(2n+2)!!}{(2n+3)!!}\,\tfrac{\abs{r}}{1+r^2} \left(\tfrac{r^2}{1+r^2}\right)^{n+1}.
\end{equation}

Let $\alpha\in(0,1]$. Noting that $k\ge n+1\ge 0$ and $\int_0^1(1-x^2)^k\,\mathrm{d}x= \int_0^\frac{\pi}{2}\cos^{2k+1}t\,\mathrm{d}t=\frac{(2k)!!}{(2k+1)!!}$ for $k\in\mathbb{N}_0$ <cit.>, and using Tonelli's theorem, it follows

\begin{equation}\label{Euler3}
\begin{split}
S(r,\alpha)&=\sum_{k=n+1}^\infty\tfrac{(2k)!!}{(2k+1)!!} \left(\tfrac{\alpha\,r^2}{1+r^2}\right)^k=
\sum_{k=n+1}^\infty\left(\tfrac{\alpha\,r^2}{1+r^2}\right)^k \int_0^1(1-x^2)^k\,\mathrm{d}x\\
&=\int_0^1\sum_{k=n+1}^\infty\left(\tfrac{\alpha\,r^2\,(1-x^2)}{1+r^2} \right)^k\mathrm{d}x
=(1+r^2)\left(\tfrac{\alpha\,r^2}{1+r^2}\right)^{n+1} \int_0^1\tfrac{(1-x^2)^{n+1}\,\mathrm{d}x} {1+(1-\alpha)\,r^2+\alpha\,r^2x^2}.
\end{split}
\end{equation}

The denominator in the last integrand is an increasing function of $x$. Inserting $1$ and $0$ instead of $x$ implies

\begin{equation}\label{Euler4}
\tfrac{(2n+2)!!}{(2n+3)!!}\,\tfrac{\alpha\,\abs{r}}{1+r^2} \left(\tfrac{\alpha\,r^2}{1+r^2}\right)^{n+1}\le \left\lvert\tfrac{\alpha\,r}{1+r^2}\,S(r,\alpha)\right\rvert\le \tfrac{(2n+2)!!}{(2n+3)!!}\,\tfrac{\alpha\,\abs{r}}{1+(1-\alpha)\,r^2} \left(\tfrac{\alpha\,r^2}{1+r^2}\right)^{n+1}
\end{equation}

and, together with $B_n=\frac{r}{1+r^2}\,S(r,1)$, (<ref>).

If $r\ne 0$, $\arctan r=\sgn(r)\,\frac{\pi}{2}-\arctan\frac{1}{r}$, (<ref>) with $\frac{1}{r}$ instead of $r$, and $r=\sgn(r)\abs{r}$ imply

\begin{equation*}
\arctan r=\sgn(r)\left(\tfrac{\pi}{2}-\tfrac{\abs{r}}{1+r^2} \sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!}\,\tfrac{1}{(1+r^2)^k}\right)
\end{equation*}

and (<ref>) with $B_n=\frac{(2n+2)!!}{(2n+3)!!}\,\frac{\abs{r}}{(1+r^2)^{n+2}}$. This series for $\abs{r}>1$ converges faster than (<ref>), so we can already sense some kind of similarity with (<ref>).

In the sequel, the lower and upper regularized gamma functions $P(a,x)=\frac{\gamma(a,x)}{\Gamma(a)}$ and $Q(a,x)=\frac{\Gamma(a,x)}{\Gamma(a)}$ respectively, where $\gamma(a,x)=\int_0^x t^{a-1}e^{-t}\,\mathrm{d}t$ and $\Gamma(a,x)=\int_x^\infty t^{a-1}e^{-t}\,\mathrm{d}t$ are the lower and upper incomplete gamma functions respectively, both defined for $a>0$ and $x\ge 0$, will be very important, as well as their lower and upper bounds.

In the Alzer's inequality

\begin{equation}\label{Alzer1}
(1-e^{-s_{k+1}x})^a\le P(a,x)\le(1-e^{-x})^a\quad(a\ge 1,x\ge 0),
\end{equation}

where $s_{k+1}=((k+1)!)^{-\frac{1}{k+1}}$ , allowing $a=1$ resulted in strict inequalities becoming non-strict. It implies

\begin{equation}\label{Alzer2}
1-(1-e^{-x})^a\le Q(a,x)\le 1-(1-e^{-s_{k+1}x})^a.
\end{equation}

For $n\in\mathbb{N}_0$ and $x\ge 0$, we have $I_n(x)=e^x\,P(n+1,x)$ <cit.>,

where $I_n(x)=e^x-\sum_{k=0}^n\frac{x^k}{k!}$ is the remainder of the Taylor series of the exponential function. The Lagrangian form of the remainder $I_n(x)=\frac{x^{n+1}\,e^{\xi_n(x)}}{(n+1)!}$, where $0<\xi_n(x)<x$, implies $I_n(x)\ge \frac{x^{n+1}}{(n+1)!}$ for $x\ge 0$. Since $I_n(x)<\frac{n+2}{n+2-x}\,\frac{x^{n+1}}{(n+1)!}$ for $x\in(0,n+2)$ <cit.>, it follows

\begin{equation}\label{Mitrinovic}
\tfrac{x^{k+1}\,e^{-x}}{(k+1)!}\le P(k+1,x)\le\tfrac{2\,x^{k+1}\,e^{-x}}{(k+1)!} \quad\left(k\in\mathbb{N}_0,0\le x\le\tfrac{k+2}{2}\right),
\end{equation}

where the equality holds only for $x=0$.

Let $\hat{f}(x)=\sum_{k=0}^\infty\hat{c}_k\,x^{2k}$ and $\tilde{f}(x)=\sum_{k=0}^\infty\tilde{c}_k\,x^{2k+1}$ be absolutely convergent series for $x\in\mathbb{R}$. Then

\begin{equation}\label{Eq:Theo1a}
\hat{F}(x)=\int_0^x\varphi(\beta t)\,\hat{f}(t)\,\mathrm{d}t=\tfrac{1}{2}\sum_{k=0}^\infty\tfrac{(2k-1)!!\, \hat{c}_k}{\abs{\beta}^{2k+1}}\,P(k+\tfrac{1}{2},\tfrac{1}{2}\beta^2x^2)
\end{equation}

\begin{equation}\label{Eq:Theo1b}
\tilde{F}(x)=\int_0^x\varphi(\beta t)\,\tilde{f}(t)\,\mathrm{d}t=\tfrac{1}{\sqrt{2\pi}}\sum_{k=0}^\infty \tfrac{(2k)!!\,\tilde{c}_k}{\beta^{2k+2}}\,P(k+1,\tfrac{1}{2}\beta^2x^2)
\end{equation}

for $x\ge 0$ and $\beta\ne 0$.

Let $f_k(t)=\hat{c}_k\,\varphi(\beta t)\,t^{2k}$ and $F_k(x)=\int_0^x f_k(t)\,\mathrm{d}t$ for $x\ge 0$. Using $(2k-1)!!=\frac{2^k}{\sqrt{\pi}}\,\Gamma(k+\frac{1}{2})$ and a substitution $\frac{1}{2}\beta^2t^2=u$, we get

\begin{equation*}
\begin{split}
F_k(x)&=\tfrac{\hat{c}_k}{\sqrt{2\pi}}\int_0^x e^{-\frac{1}{2}\beta^2t^2}t^{2k}\,\mathrm{d}t=\tfrac{2^{k-1} \hat{c}_k}{\sqrt{\pi}\;\abs{\beta}^{2k+1}}\int_0^{\frac{1}{2}\beta^2x^2} u^{k-\frac{1}{2}}e^{-u}\,\mathrm{d}u\\
&=\tfrac{(2k-1)!!\,\hat{c}_k}{2\,\abs{\beta}^{2k+1}}\, \tfrac{\gamma(k+\frac{1}{2},\frac{1}{2}\beta^2x^2)}{\Gamma(k+\frac{1}{2})}= \tfrac{(2k-1)!!\,\hat{c}_k}{2\,\abs{\beta}^{2k+1}}\, P(k+\tfrac{1}{2},\tfrac{1}{2}\beta^2x^2).
\end{split}
\end{equation*}

Setting $j=\lceil\beta^2x^2\rceil+1$ implies $0\le u<\frac{j+2}{2}$ for $u\in[0,\frac{1}{2}\beta^2x^2]$. Then $\sum_{k=0}^\infty\abs{F_k(x)}= \sum_{k=0}^j\abs{F_k(x)}+\sum_{k=j+1}^\infty \abs{F_k(x)}$, where the first sum is finite. Using the fact that $P(a,x)$ is a decreasing function of $a$ <cit.>, (<ref>) and the absolute convergence assumption of $\hat{f}(x)$, it follows

\begin{equation*}
P(k+\tfrac{1}{2},\tfrac{1}{2}\beta^2x^2)\le P(k,\tfrac{1}{2}\beta^2x^2)\le \tfrac{2\,\beta^{2k}\,x^{2k}\,e^{-\frac{1}{2}\beta^2x^2}}{2^k k!}\quad(k\ge j+1)
\end{equation*}

\begin{equation*}
\sum_{k=j+1}^\infty\abs{F_k(x)}=\sum_{k=j+1}^\infty\tfrac{(2k-1)!!\,
\abs{\hat{c}_k}}{2\,\abs{\beta}^{2k+1}}\,P(k+\tfrac{1}{2},\tfrac{1}{2}\beta^2x^2)
\le\tfrac{e^{-\frac{1}{2}\beta^2x^2}}{\abs{\beta}} \sum_{k=j+1}^\infty\tfrac{(2k-1)!!}{(2k)!!}\;\abs{\hat{c}_k}\,x^{2k}<\infty.
\end{equation*}

Noting that $\int_0^x\abs{f_k(t)}\,\mathrm{d}t=\abs{F_k(x)}$, we proved that $\sum_{k=0}^\infty\int_0^x\abs{f_k(t)}\,\mathrm{d}t<\infty$. By Fubini's theorem it follows that $\int_0^x \varphi(\beta t)\,\hat{f}(t)\,\mathrm{d}t=\sum_{k=0}^\infty F_k(x)$, hence (<ref>).

Let $f_k(t)=\tilde{c}_k\,\varphi(\beta t)\,t^{2k+1}$ and $F_k(x)=\int_0^x f_k(t)\,\mathrm{d}t$ for $x\ge 0$. Analogously to above, and using $(2k)!!=2^k k!=2^k\,\Gamma(k+1)$, we get $F_k(x)=\frac{(2k)!!\,\tilde{c}_k}{\sqrt{2\pi}\; \beta^{2k+2}}\,P(k+1,\frac{1}{2}\beta^2x^2)$,

\begin{equation*}
\sum_{k=j+1}^\infty\abs{F_k(x)}\le\tfrac{\abs{x}\,e^{-\frac{1}{2}\beta^2x^2}} {\sqrt{2\pi}}\sum_{k=j+1}^\infty\tfrac{\abs{\tilde{c}_k}\, \abs{x}^{2k+1}}{k+1}<\infty
\end{equation*}

and (<ref>).

Note that in the same way as the equations (<ref>) and (<ref>) were derived, analogous equations can be derived if the integration limits $0$ and $x$ are replaced by $x$ and $\infty$ respectively, and lower regularized gamma function is replaced by the upper one. However, for each case it should be checked whether the order of integration and summation can be reversed.

The usefulness of Theorem <ref> is that if the sequences $\{\hat{c}_k\}$ and $\{\tilde{c}_k\}$ can be computed recursively, then the terms of the series (<ref>) and (<ref>) can be too. Let $q=\frac{1}{2}\beta^2x^2$. Crucial is the calculation of $\hat{d}_k=Q(k+\frac{1}{2},q)$ and $\tilde{d}_k=P(k+1,q)$, $k\in\mathbb{N}_0$, the rest is trivial. Since $P(k+1,q)=1-Q(k+1,q)$, we will calculate $d_k=Q(k+1,q)$ instead of $\tilde{d}_k=1-d_k$.

Using $P(\frac{1}{2},q)=2\,\Phi(\sqrt{2q}\,)-1$ <cit.>, we get $\hat{d}_0=1-P(\frac{1}{2},q)=2\,(1-\Phi(\sqrt{2q}\,))$, and $d_0=e^{-q}$. Since $P(a+1,x)=P(a,x)-\frac{x^a\,e^{-x}}{\Gamma(a+1)}$ <cit.> implies $Q(a+1,x)=Q(a,x)+\frac{x^a\,e^{-x}}{\Gamma(a+1)}$, using auxiliary variables $\hat{b}_k=\frac{q^{k-\frac{1}{2}}\,e^{-q}}{\Gamma(k+\frac{1}{2})}$ and $b_k=\frac{q^k\,e^{-q}}{\Gamma(k+1)}$, it follows

\begin{equation*}
\hat{b}_0=\tfrac{e^{-q}}{\sqrt{q\,\pi}},\quad\hat{d}_0=2\,(1-\Phi(\sqrt{2\,q}\,)) \quad\text{and}\quad
\hat{b}_{k+1}=\tfrac{q\,\hat{b}_k}{k+\frac{1}{2}},\quad\hat{d}_{k+1}= \hat{d}_k+\hat{b}_{k+1}
\end{equation*}

for $q\ne 0$ and $k=0,1,2,\dotsc$, and

\begin{equation}\label{Eq:Theo1c}
b_0=e^{-q},\quad d_0=b_0\quad\text{and}\quad b_{k+1}=\tfrac{q\,b_k}{k+1},\quad d_{k+1}=d_k+b_{k+1}
\end{equation}

for $k=0,1,2,\dotsc$

§ MAIN RESULTS

The core of this article are Theorems <ref>, <ref> and <ref>. Since the first two do not specify variants for the transformed parameters, the last one, together with (<ref>) and (<ref>), is more practically useful. It also justifies the title of the article.

Let $h\in\mathbb{R}$, $\varrho\in(-1,1)$, $r=\frac{\varrho}{\sqrt{1-\varrho^2}}$ and $q=\frac{1}{2}(1+r^2)h^2=\frac{h^2}{2\,(1-\varrho^2)}$. Then

\begin{equation}\label{Eq:Theo2a}
\begin{split}
\Phi_\varrho^2(h,0)&=\tfrac{1}{2}\,\Phi(h)+\tfrac{\arctan r}{2\pi}-\tfrac{r}{2\pi\,(1+r^2)}\sum_{k=0}^\infty c_k\left(\tfrac{r^2}{1+r^2}\right)^k\\
&=\tfrac{1}{2}\,\Phi(h)+\tfrac{\arcsin\varrho}{2\pi}-\tfrac{\sqrt{1-\varrho^2}} {2\pi}\, \sum_{k=0}^\infty c_k\,\varrho^{2k+1},
\end{split}
\end{equation}

where $c_k=\frac{(2k)!!}{(2k+1)!!}\,P(k+1,q)$. The series $S=\sum_{k=0}^\infty c_k\left(\frac{r^2}{1+r^2}\right)^k$ converges for every $r,h\in\mathbb{R}$. If it is truncated after the $n^{th}$ term, then

\begin{equation}\label{Eq:Theo2b}
\Phi_\varrho^2(h,0)=\tfrac{1}{2}\,\Phi(h)+\tfrac{\arctan r}{2\pi}-\tfrac{r\,S_n}{2\pi\,(1+r^2)}-\sgn(r)R_n\quad(n\in\mathbb{N}_{-1}),
\end{equation}

where $S_n=\sum_{k=0}^n c_k\left(\frac{r^2}{1+r^2}\right)^k$ and $R_n=\frac{\abs{r}}{2\pi\,(1+r^2)} \sum_{k=n+1}^\infty c_k\left(\frac{r^2}{1+r^2}\right)^k$ is bounded by

\begin{equation}\label{Eq:Theo2c}
0\le R_n\le\tfrac{B_n}{1+e^{-q}r^2}\le B_n,
\end{equation}

where $B_n=\frac{(2n+2)!!}{(2n+3)!!}\,\frac{\abs{r}}{2\pi}\,(1-e^{-q})^{n+2} \left(\frac{r^2}{1+r^2}\right)^{n+1}$.

From (<ref>) and $\Phi(x)-\frac{1}{2}=\varphi(x)\sum_{k=0}^\infty\frac{x^{2k+1}}{(2k+1)!!}$ for $x\in\mathbb{R}$ <cit.>, it follows

\begin{equation}\label{Eq:Theo2d}
\tfrac{\partial T(h,r)}{\partial h}=-\,\tfrac{h}{2\pi}\int_0^r e^{-\frac{1}{2}h^2(1+x^2)}\,\mathrm{d}x=-\varphi(h)\left(\Phi(rh)-\tfrac{1}{2}\right)= \varphi(\beta h)\sum_{k=0}^\infty\tilde{c}_k\,h^{2k+1},
\end{equation}

$r,h\in\mathbb{R}$, where $\beta=\sqrt{1+r^2}$ and $\tilde{c}_k=-\,\frac{r^{2k+1}}{\sqrt{2\pi}\;(2k+1)!!}$. Since

\begin{equation*}
\sum_{k=0}^\infty\tilde{c}_k\,h^{2k+1}=-\,\tfrac{\varphi(h)}{\varphi(\beta h)}\left(\Phi(rh)-\tfrac{1}{2}\right)= -\,e^{\frac{1}{2}r^2h^2}\left(\Phi(rh)-\tfrac{1}{2}\right),
\end{equation*}

the series on the left side converges for every $r,h\in\mathbb{R}$ and is an odd function of $rh$. If $rh<0$, then all its terms are positive, implying that it converges absolutely.

All requirements of Theorem <ref> are fulfilled. Assuming $r\ge 0$ and using $T(0,r)=\frac{\arctan r}{2\pi}$, (<ref>) and (<ref>), it follows

\begin{equation*}
T(h,r)=\tfrac{\arctan r}{2\pi}-\tfrac{1}{2\pi}\sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!}\, \tfrac{r^{2k+1}}{\beta^{2k+2}}\,P(k+1,\tfrac{1}{2}\beta^2h^2),
\end{equation*}

which, together with (<ref>), (<ref>), $\frac{r^2}{\beta^2}=\frac{r^2}{1+r^2}=\varrho^2$, $\frac{r}{1+r^2}=\varrho\,\sqrt{1-\varrho^2}$ and $\frac{1}{2}\beta^2h^2=q$, imply (<ref>) for $r\ge 0$. Since (<ref>) implies $T(h,-r)=-T(h,r)$, (<ref>) is also valid for $r<0$. Let $\alpha=1-e^{-q}$. Using (<ref>), (<ref>) and (<ref>), it follows

\begin{equation*}
0\le R_n\le\tfrac{\alpha\,\abs{r}}{2\pi\,(1+r^2)}\sum_{k=n+1}^\infty \tfrac{(2k)!!}{(2k+1)!!}\left(\tfrac{\alpha\,r^2}{1+r^2}\right)^k= \tfrac{\alpha\,\abs{r}}{2\pi\,(1+r^2)}\,S(r,\alpha)
\end{equation*}

and (<ref>).

Note that if $rh\ne 0$, then $c_k>0$ and $\frac{c_{k+1}}{c_k}=\frac{2k+2}{2k+3}\,\frac{P(k+2,q)}{P(k+1,q)}<1$. Hence the sequence $\left\{c_k\left(\frac{r^2}{1+r^2}\right)^k\right\}$ is decreasing.

For every $h\in\mathbb{R}$ we have

\begin{equation}\label{Eq:Corr1a}
\Phi(h)=\tfrac{1}{2}+\sgn(h)\sqrt{S},\quad\text{where}\quad S=\tfrac{1}{2\pi}\sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!}\, \tfrac{P(k+1,h^2)}{2^k}.
\end{equation}

Using (<ref>) with $\varrho=-\,\frac{\sqrt{2}}{2}$, hence $r=-1$, $q=h^2$ and the right equation (<ref>), it follows

\begin{equation*}
\tfrac{1}{2}\,\Phi(h)-\tfrac{1}{8}+\tfrac{1}{4\pi}\sum_{k=0}^\infty \tfrac{(2k)!!}{(2k+1)!!}\,\tfrac{P(k+1,h^2)}{2^k}=\tfrac{1}{2}(\Phi(h))^2
\end{equation*}

and a quadratic equation $(\Phi(h))^2-\Phi(h)+\frac{1}{4}-S=0$. If $h=0$, then $S=0$ and (<ref>) is valid. If $h\ne 0$, since $\Phi(h)>\frac{1}{2}$ if $h>0$ and $\Phi(h)<\frac{1}{2}$ if $h<0$, only $\frac{1}{2}+\sgn(h)\sqrt{S}$ is a proper root of the two. Note that sending $h\to\infty$ implies $\sum_{k=0}^\infty\frac{k!}{(2k+1)!!}=\frac{\pi}{2}$.

Let $h\in\mathbb{R}$, $\varrho\in(-1,1)$, $r=\frac{\varrho}{\sqrt{1-\varrho^2}}$ and $q=\frac{1}{2}(1+r^2)h^2=\frac{h^2}{2\,(1-\varrho^2)}$. Then

\begin{equation}\label{Eq:Theo3a}
\begin{split}
\Phi_\varrho^2(h,0)&=\tfrac{1}{2}\,\Phi(h)+\tfrac{r}{2\pi\,(1+r^2)} \sum_{k=0}^\infty c_k\left(\tfrac{r^2}{1+r^2}\right)^k\\
&=\tfrac{1}{2}\,\Phi(h)+\tfrac{\sqrt{1-\varrho^2}}{2\pi}\,\sum_{k=0}^\infty c_k\,\varrho^{2k+1},
\end{split}
\end{equation}

where $c_k=\frac{(2k)!!}{(2k+1)!!}\,Q(k+1,q)$. The series $S=\sum_{k=0}^\infty c_k\left(\frac{r^2}{1+r^2}\right)^k$ converges for every $r,h\in\mathbb{R}$. If it is truncated after the $n^{th}$ term, then

\begin{equation}\label{Eq:Theo3b}
\Phi_\varrho^2(h,0)=\tfrac{1}{2}\,\Phi(h)+\tfrac{r\,S_n}{2\pi\,(1+r^2)}+ \sgn(r)R_n\quad(n\in\mathbb{N}_{-1}),
\end{equation}

where $S_n=\sum_{k=0}^n c_k\left(\frac{r^2}{1+r^2}\right)^k$ and $R_n=\frac{\abs{r}}{2\pi\,(1+r^2)} \sum_{k=n+1}^\infty c_k\left(\frac{r^2}{1+r^2}\right)^k$ is bounded by

\begin{equation}\label{Eq:Theo3c}
0\le R_n\le B_n,\quad\text{where}\quad B_n=\tfrac{(2n+2)!!}{(2n+3)!!}\,\tfrac{\abs{r}}{2\pi} \left(\tfrac{r^2}{1+r^2}\right)^{n+1}.
\end{equation}

Using (<ref>), (<ref>), $P(k+1,x)=1-Q(k+1,x)$, $\frac{r^2}{1+r^2}=\varrho^2$ and $\frac{r}{1+r^2}=\varrho\,\sqrt{1-\varrho^2}$, it follows (<ref>). Noting that $R_n=\frac{\abs{r}}{2\pi\,(1+r^2)}\,S(r,1)$ and (<ref>) imply (<ref>).

Note that the upper bound (<ref>) is not suitable for a simple estimation of $R_n$ here because $\lim_{k\to\infty}s_{k+1}=0$, hence $Q(k+1,q)\le 1$ was implicitly used instead. $Q(a,x)$ is an increasing function of $a$. The ratio $\frac{c_{k+1}}{c_k}= \frac{2k+2}{2k+3}\,\frac{Q(k+2,q)}{Q(k+1,q)}$ can be greater than $1$, e.g. if $h=3$ and $r=1$, we get $\frac{c_2}{c_1}=\frac{4\,Q(3,9)}{5\,Q(2,9)}\approx 4.04$. In Theorem <ref>, the sequence $\left\{c_k\left(\frac{r^2}{1+r^2}\right)^k\right\}$ is not always decreasing, but, in any case, the well-behaved $\arctan\,\abs{r}$ series (<ref>) divided by $2\pi$ is a majorant series for $\frac{\abs{r}\,S}{2\pi\,(1+r^2)}$ from Theorems <ref> and <ref>.

Let $x,y\in\mathbb{R}$, $x^2+y^2>0$, $\varrho\in(-1,1)$, $\varrho_x=\frac{(\varrho x-y)\,\sgn(x)}{\sqrt{x^2-2\varrho xy+y^2}}$, $\varrho_y=\frac{(\varrho y-x)\,\sgn(y)} {\sqrt{x^2-2\varrho xy+y^2}}$ and $q=\frac{x^2-2\varrho xy+y^2}{2\,(1-\varrho^2)}$. Let $c_k=\frac{(2k)!!}{(2k+1)!!}\,P(k+1,q)$ and $\tilde{c}_k=\frac{(2k)!!}{(2k+1)!!}\,Q(k+1,q)$ for $k\in\mathbb{N}_0$. Then

\begin{equation}\label{Eq:Corr2a}
\begin{split}
\Phi_\varrho^2(x,y)&=\tfrac{1}{2}\,(\Phi(x)+\Phi(y))+\tfrac{\arcsin\varrho_x+ \arcsin\varrho_y}{2\pi}-S-\beta\\
\end{split}
\end{equation}

\begin{equation*}
S=\tfrac{1}{2\pi\,\sqrt{2q}}\sum_{k=0}^\infty c_k\left(\abs{x}\varrho_x^{2k+1}+\abs{y}\varrho_y^{2k+1}\right),
\end{equation*}

$\tilde{S}$ is $S$ with $\tilde{c}_k$ instead of $c_k$, and $\beta=0$ if $xy>0$ or $xy=0$ and $x+y\ge 0$, and $\beta=\frac{1}{2}$ otherwise.

Since $q_x=\frac{x^2}{2\,(1-\varrho_x^2)}=\frac{x^2-2\varrho xy+y^2}{2\,(1-\varrho^2)}=q$ and $q_y=\frac{y^2}{2\,(1-\varrho_y^2)}=\frac{x^2-2\varrho xy+y^2}{2\,(1-\varrho^2)}=q$, $\sqrt{1-\varrho_{x\phantom{y}}^2}=\frac{\abs{x}\,\sqrt{1-\varrho^2}} {\sqrt{x^2-2\varrho xy+y^2}}=\frac{\abs{x}}{\sqrt{2q}}$ and $\sqrt{1-\varrho_y^2}=\frac{\abs{y}}{\sqrt{2q}}$, (<ref>), (<ref>), (<ref>) and (<ref>) imply (<ref>).

Note that $\frac{r}{1+r^2}=\frac{\overline{r}}{1+\overline{r}^2}$ and $q=\frac{1}{2}(1+r^2)h^2=\frac{1}{2}(1+\overline{r}^2)\overline{h}^2=\overline{q}$, therefore they are invariant on the transformation of $rh\mapsto\overline{h}$ and $\frac{1}{r}\mapsto\overline{r}$, but $S_n$, $R_n$ and $B_n$ from Theorems <ref> and <ref> are not because $\varrho^2= \frac{r^2}{1+r^2}\mapsto\frac{1}{1+r^2}=\frac{\overline{r}^2}{1+\overline{r}^2}= \overline{\varrho}^2$ and $\varrho^2\ne\overline{\varrho}^2$ if $\abs{r}\ne 1$. They transform to $\overline{S}_n$, $\overline{R}_n$ and $\overline{B}_n$. By transforming the parameters and using (<ref>), the equations (<ref>) and (<ref>) can be rewritten to

\begin{equation}\label{Eq:Theo2e}
\Phi_\varrho^2(h,0)=\tfrac{1}{2}\,\Phi(rh)+\Phi(h)(1-\Phi(rh))- \tfrac{\arctan\overline{r}}{2\pi}+\tfrac{r\,\overline{S}_n}{2\pi\,(1+r^2)}+ \sgn(r)\overline{R}_n-\beta
\end{equation}

\begin{equation}\label{Eq:Theo3d}
\Phi_\varrho^2(h,0)=\tfrac{1}{2}\,\Phi(rh)+\Phi(h)(1-\Phi(rh))- \tfrac{r\,\overline{S}_n}{2\pi\,(1+r^2)}-\sgn(r)\overline{R}_n-\beta
\end{equation}

respectively, where $\beta=0$ if $r\ge 0$ and $\beta=\frac{1}{2}$ if $r<0$.

The transformation makes sense if $\abs{r}>1$, implying $\overline{\varrho}^2<\frac{1}{2}$ and faster convergence. For example, the inequality $0\le R_n\le B_n$ from Theorem <ref> transforms into $0\le\overline{R}_n\le\overline{B}_n$, where the upper bound $\overline{B}_n$ is $r^{2n+4}$-times smaller than a corresponding upper bound $B_n$. The noteworthy cost for faster convergence is the additional computation of $\Phi(rh)$.

If $\abs{r}<\infty$, the Owen's $T$ function (<ref>) can be written as

\begin{align}
&=\tfrac{\arctan r}{2\pi}-\tfrac{r}{2\pi\,(1+r^2)}\sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!}\, P(k+1,q)\left(\tfrac{r^2}{1+r^2}\right)^k\label{Owen6a}\\
&=\tfrac{r}{2\pi\,(1+r^2)}\sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!}\, Q(k+1,q)\left(\tfrac{r^2}{1+r^2}\right)^k\label{Owen6b}\\
&=U(h,r)-\tfrac{\arctan r}{2\pi}+\tfrac{r}{2\pi\,(1+r^2)} \sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!}\,P(k+1,q) \left(\tfrac{1}{1+r^2}\right)^k\label{Owen6c}\\
&=U(h,r)-\tfrac{r}{2\pi\,(1+r^2)}\sum_{k=0}^\infty\tfrac{(2k)!!}{(2k+1)!!}\, Q(k+1,q)\left(\tfrac{1}{1+r^2}\right)^k,\label{Owen6d}
\end{align}

where $q=\frac{1}{2}(1+r^2)h^2$, $U(h,r)=\frac{1}{2}(\Phi(h)+\Phi(rh))-\Phi(h)\,\Phi(rh)-\beta$, and $\beta=0$ if $r\ge 0$ and $\beta=\frac{1}{2}$ if $r<0$.

Using (<ref>) and (<ref>), the equations (<ref>) and (<ref>) imply (<ref>) and (<ref>) respectively, which, together with (<ref>), imply (<ref>) and (<ref>) respectively.

Note that when truncating the series (<ref>) and (<ref>), the truncation error in absolute value can be estimated by (<ref>) and (<ref>) respectively. The same upper bounds with $\overline{r}=\frac{1}{r}$ instead of $r$ also apply to (<ref>) and (<ref>).

Because of (<ref>) and $0\le Q(k+1,q)\le 1$, the Owen's $T$ function can be seen as a modified Euler's arctangent series (<ref>) divided by $2\pi$, which is its majorant series.

§ RECURSION AND ASYMPTOTIC SPEED OF CONVERGENCE

Let $\tilde{S}=\frac{r\,S}{2\pi\,(1+r^2)}$ and $\tilde{e}_k=B_k$, where $S$ and $B_k$ relate to Theorem <ref>. Using (<ref>) and (<ref>), the series $\tilde{S}$ can be calculated by Recursion <ref> (tildes are omitted). Analogously, let $\hat{S}=\frac{r\,S}{2\pi\,(1+r^2)}$ and $\hat{e}_k=B_k$, where $S$ and $B_k$ relate to Theorem <ref>. The series $\hat{S}$ can be calculated by Recursion <ref> (hats are omitted). In practice, both recursions can be easily combined, since only four steps, marked with $\triangleright$, are different.

In both cases, given $\epsilon>0$, the recursion is terminated when the condition $e_{k+1}<\epsilon$ is met and $S_{k+1}$ is taken as its result. However, if $\epsilon>0$ is too small, or even intentionally negative, the calculation proceeds up to the actual accuracy threshold – the recursion is terminated when the condition $S_{k+1}\le S_k$ is met. This criterion was used in the calculations presented in the next section.

In both recursions, all variables are positive and are bounded upward. If $q\ne 0$, the sequence $\{d_k\}$ strictly increases and, due to (<ref>), converges to $1$. The sequence $\{S_k\}$ is strictly increasing if $r\ne 0$. If $r=0$, the recursion terminates with $S_1=0$. Note that (<ref>) implies $\abs{T(h,r)}\le\frac{\arctan\,\abs{r}}{2\pi}$. If step $7$ in Recursion <ref> is replaced by $S_0\leftarrow a_0\,(1-d_0)-\frac{\arctan\,\abs{r}}{2\pi}$, the sequence $\{S_k\}$ changes its sign but keeps the direction. However, the modified Recursion <ref> computes $\hat{S}=\frac{r\,S}{2\pi\,(1+r^2)}-\frac{\arctan r}{2\pi}$.

[1]
$g\gets p$ $(1-e^{-q})$ $a_0\gets\frac{\abs{r}}{2\pi\,(1+r^2)}$
$b_0\gets e^{-q}$
$d_0\gets b$
$S_0\gets a\,d_0$ $()$ $e_0\gets\frac{\abs{r}\,p}{3\pi}$ $k\gets 0$
$d_{k+1}\gets d_k+b_{k+1}$
$S_{k+1}\gets S_k+a_{k+1}\,d_{k+1}$ $()$ $e_{k+1}\gets\frac{(2k+4)\,g\,e_k}{2k+5}$
$k\gets k+1$
$e_{k+1}<\epsilon\textbf{ or }S_{k+1}\le S_k$

[1]
$b_0\gets e^{-q}$
$d_0\gets b$
$S_0\gets a\,(1-d_0)$
$k\gets 0$
$d_{k+1}\gets d_k+b_{k+1}$
$S_{k+1}\gets S_k+a_{k+1}\,(1-d_{k+1})$
$k\gets k+1$
$e_{k+1}<\epsilon\textbf{ or }S_{k+1}\le S_k$

Let $\tilde{S}_n$ and $\hat{S}_n$ be the results of Recursion <ref> and modified Recursion <ref> respectively, calculated with the same $h$ and $r$ if $\abs{r}\le 1$, or with $\overline{h}$ and $\overline{r}$ if $\abs{r}>1$. Let $\tilde{\Phi}_\varrho^2(h,0)$ be calculated by (<ref>) or (<ref>), where unknown remainders are ignored. Analogously, let $\hat{\Phi}_\varrho^2(h,0)$ be calculated by (<ref>) or (<ref>), where also $\frac{\arctan r}{2\pi}$ is ignored. If $\abs{r}\le 1$, then $\tilde{\Phi}_\varrho^2(h,0)\le\Phi_\varrho^2(h,0)\le\hat{\Phi}_\varrho^2(h,0)$ if $r\ge 0$ and $\tilde{\Phi}_\varrho^2(h,0)\ge\Phi_\varrho^2(h,0)\ge\hat{\Phi}_\varrho^2(h,0)$ if $r<0$. Even if $\abs{r}>1$, one of $\tilde{\Phi}_\varrho^2(h,0)$ and $\hat{\Phi}_\varrho^2(h,0)$ underestimates and the other overestimates $\Phi_\varrho^2(h,0)$. Consequently, the upper bounds (<ref>) and (<ref>) are not needed if $\tilde{S}$ and modified $\hat{S}$ are calculated in the same recursion which terminates when $\abs{\hat{S}_{k+1}-\tilde{S}_{k+1}}$ become small enough. Upper bounds are also not needed if in Recursions <ref> and <ref> step $17$ is replaced so that $S_{k+1}\le S_k$ is the only termination condition. In this case, steps $3$, $8$ and $15$ are unnecessary.

Given $q$, which one of the series (<ref>) and (<ref>) converges faster? Since $P(k+1,x)$ is a cdf of the gamma distribution with the median $\nu\in(k+\frac{2}{3},k+1)$, it follows $P(k+1,q)\le Q(k+1,q)$ if $q\le\nu<k+1$. The answer depends on the number of considered terms, which depends on the prescribed accuracy. On this basis, it seems that series (<ref>) is more suitable than (<ref>) for “small” $q$ and vice versa for “large” $q$. Tables <ref> and <ref> in the next section show how this is reflected in practice.

Gautschi's inequality $x^{1-s}<\frac{\Gamma(x+1)}{\Gamma(x+s)}<(x+1)^{1-s}$ <cit.> with $s=\frac{1}{2}$ implies $\frac{\sqrt{k\,\pi}}{2k+1}< \frac{(2k)!!}{(2k+1)!!}<\frac{\sqrt{(k+1)\,\pi}}{2k+1}$. Combining it with (<ref>), sending $k\to\infty$ and neglecting non-essential constants imply that the series (<ref>) asymptotically converges as $\sum_k\frac{q^{k+1}e^{-q}}{(k+1)!}\,\frac{\varrho^{2k}}{\sqrt{k+1}}$. Assuming calculation with $\overline{h}$, $\overline{r}$ and (<ref>) instead of $h$, $r$ and (<ref>) respectively if $\abs{r}>1$, hence $\varrho^2>\frac{1}{2}$, and recalling that $q$ is invariant on the transformation of $rh\mapsto\overline{h}$ and $\frac{1}{r}\mapsto\overline{r}$, we get an asymptotic convergence as $\sum_k\frac{\overline{q}^{k+1} e^{-\overline{q}}}{(k+1)!}\,\frac{\overline{\varrho}^{2k}}{\sqrt{k+1}}= \sum_k\frac{q^{k+1}e^{-q}}{(k+1)!}\,\frac{(1-\varrho^2)^k}{\sqrt{k+1}}$. The series (<ref>) or its corresponding transformed version, unwritten and based on (<ref>), asymptotically converges as $\sum_k\frac{q^{k+1}e^{-q}}{(k+1)!}\,\frac{(\min\{\varrho^2,1-\varrho^2\})^k} {\sqrt{k+1}}=\sum_k\frac{q^{k+1}e^{-q}}{(k+1)!}\, \frac{\min\{r^{2k},1\}}{\sqrt{k+1}\;(1+r^2)^k}$, where $q=\frac{1}{2}(1+r^2)h^2=\frac{h^2}{2\,(1-\varrho^2)}$. On the other hand, since (<ref>) implies $\lim_{k\to\infty}Q(k+1,x)=1$, the series (<ref>) or its corresponding transformed version, unwritten and based on (<ref>), asymptotically converges as $\sum_k\frac{(\min\{\varrho^2,1-\varrho^2\})^k}{\sqrt{k+1}}= \sum_k\frac{\min\{r^{2k},1\}}{\sqrt{k+1}\;(1+r^2)^k}$, hence in the worst case as $\sum_k\frac{1}{2^k\,\sqrt{k+1}}$.

Based on both asymptotic convergences, comparing (<ref>) to (<ref>) would be unfair because the calculation of $\arctan\,\abs{r}$ should also be taken into account for (<ref>). Since it is implicitly embedded in Recursion <ref>, it is calculated on the fly at the cost of worse speed of convergence. Each of the recursions is executable with the basic four arithmetic operations, two calls of elementary functions (if the calculation of $r$ from $\varrho$ is also considered), and $\arctan\abs{r}$ in the modified Recursion <ref>. However, the computational cost of $\Phi(h)$ (and $\Phi(rh)$ if needed) should also be taken in account, but it is the same in (<ref>) and (<ref>).

§ NUMERIC ISSUES AND TESTING RESULTS

A common cause of numerical instability, where two absolutely large values give an absolutely small result when added or subtracted, cannot occur in Recursions <ref> and <ref>. The calculation by them is stable, but something may happen that needs to be pointed out. The side effect of transforming $rh\mapsto\overline{h}$ and $\frac{1}{r}\mapsto\overline{r}$ is that $\abs{\overline{h}}$ can become unusually large. When computing $\Phi(h)$, according to <cit.>, “[...] the region of interest may not be the statistician's $0 < x < 5$, say, but rather $10 < x < 14$ [...]”. Assume the same needs when computing $\Phi_\varrho^2(h,0)$, and using transformed parameters if $\varrho^2>\frac{1}{2}$. If $h=14$ and $\varrho=\frac{10}{\sqrt{1+10^2}}\approx 0.995$, then $r\approx 10$, $\overline{r}\approx 0.1$ and $\overline{h}\approx 140$, implying that the region of interest could be $\abs{h}<140$. Even the usual $\abs{h}<5$ transformed to $\abs{\overline{h}}<50$ is well beyond the established limits.

The “problem” of too large $\abs{h}$ is mentioned because we have to be aware of it if we are going to perform high-precision computation. In any case, the recursion must be complemented by a suitable computation of $\Phi(h)$ and $\arctan r$ if needed. As an example that not all are like that, the computation by the well known series $\Phi(x)=\frac{1}{2}+\frac{1}{\sqrt{2\pi}}\sum_{k=0}^\infty \frac{(-1)^k\,x^{2k+1}}{2^k k!\,(2k+1)}$ <cit.> and the standard double precision arithmetic is numerically unstable even for moderately large $\abs{h}$, e.g. using the programming language R [R Core Team, 2023], for $h=10$ we get a completely wrong result, approximately $-1683$. A similar problem can lead to a significant decrease in accuracy if $T(h,r)$ is calculated with (<ref>), as will be seen in Subsection <ref>.

In double precision arithmetic, the minimal positive non-zero number is $\eta\approx 2.23\cdot 10^{-308}$. Assuming $b_0=e^{-q}\ge\eta$, it follows $q_{max}=-\log(\eta)\approx 708.40$ and $h_{max}(r)=\sqrt{-\,\frac{2\log(\eta)}{1+r^2}}$. If $\abs{r}\le 1$, then $h_{max}(r)$ is between $h_{max}(1)\approx 26.62$ and $h_{max}(0)\approx 37.64$. If $e^{-q}<\eta$, then all elements of the sequence $\{b_k\}$ are zero. Consequently, all elements of $\{d_k\}$ are one and all elements of $\{S_k\}$ in Recursion <ref> are zero.

From the integral in (<ref>), it follows $\abs{T(h,r)}\le\varphi(h)\,\frac{\Phi(\abs{r}\abs{h})-\frac{1}{2}}{\abs{h}}\le \frac{\varphi(h)}{2\,\abs{h}}$ if $h\ne 0$, where the upper bound is a decreasing function of $\abs{h}$. If $e^{-q}<\eta$ and $\abs{r}\le 1$, then $\abs{h}>h_{max}(1)$ and $\abs{T(h,r)}\le\frac{\varphi(h_{max}(1))}{2\,h_{max}(1)}\approx 1.12\cdot 10^{-156}$. Analogously to the case when $r=0$, the recursion terminates with $S_1=0$, which is practically a correct result for $T(h,r)$, as well as $\Phi_\varrho^2(h,0)=\frac{1}{2}$ if $h>0$, and $\Phi_\varrho^2(h,0)=0$ if $h<0$.

The Recursions <ref> and <ref> with adjustment of their results according to (<ref>)–(<ref>) are implemented in the package Phi2rho [Komelj, 2023]. It contains the functions OwenT(), hereafter Phi2rho::OwenT(), and Phi2xy(), which compute $T(h,r)$ and $\Phi_\varrho^2(x,y)$ respectively. They also enable computation with the tetrachoric and Vasicek's series, original or accelerated. The latter means that $h$ and $r$ are transformed if $r>1$ and $r<1$ for the tetrachoric and Vasicek's series respectively, and the equation (<ref>) is used, so (<ref>) instead of (<ref>) in the case of tetrachoric series. For the sake of comparison, all series were considered during testing in Subsection <ref>. In Subsection <ref>, the original tetrachoric and Vasicek's series were excluded because they converge too slowly for unfavorable $\varrho$.

During testing, the values of $T(h,r)$ and $\Phi_\varrho^2(x,y)$ were computed for a selected and randomly generated set of parameters. Reference values were computed using the accelerated tetrachoric series and quadruple ($128$-bit) precision. Absolute differences between tested and reference values are declared as absolute errors. The maximum absolute error is considered a measure of accuracy. Since the mean values and some quantiles are also important, they are shown in the tables in Subsection <ref>. An integral part of the testing was also a comparison of the results obtained with the Phi2rho package and competing packages available on CRAN.

§.§ Numeric issues and results when computing $T(h,r)$

Numerical calculations of $T(h,r)$ were performed on $39,999$ grid points $(h,r)$, where $h$ ranges from $-10$ to $10$ in increments of $0.1$, $r=\frac{\varrho}{\sqrt{1-\varrho^2}}$ and $\varrho$ ranges from $-0.99$ to $0.99$ in increments of $0.01$. Since for $\varrho=\pm 0.99$ we get $r\approx\pm 7$, the actual range of $h$ for which $\Phi(h)$ must be computed is approximately from $-70$ to $70$, so well beyond the $h_{max}(0)\approx 37.64$.

Note atanExt = no in tables means that Recursion <ref> was executed, hence (<ref>) was used if $\abs{r}\le 1$ and (<ref>) if $\abs{r}>1$, and atanExt = yes means that Recursion <ref> was executed, hence (<ref>) was used if $\abs{r}\le 1$ and (<ref>) if $\abs{r}>1$. Both recursions are well suited for vectors, so each of the parameters $h$ and $r$ can be a scalar or a vector. If they are both vectors, they must be of the same length, if one of $h$ and $r$ is a vector and the other is a scalar, the latter is replicated to the length of the former. Since the number of iterations is determined by the “worst” component, we would expect execution to be faster if function calls are done in a loop for each component of the vector separately, but usually the run time is significantly increased due to the overhead of the loop, more initialization, etc.

For the consistency test, $T(h,r)$ was calculated on the test grid with $h$ and $r$ both scalars, one a scalar and the other a vector. For all three calculations with Recursion <ref>, the results matched perfectly, as well as for all three with Recursion <ref>. The maximum absolute difference between calculations with Recursions <ref> and <ref> was approximately $9.49\cdot 10^{-17}$.

When testing accuracy, the tetrachoric and Vasicek's series, original and accelerated, were also tested, using $T(h,r)=\Phi_\varrho^2(h,0)-\frac{1}{2}\Phi(h)$. By transforming the parameters, if appropriate, both comparison methods gained a lot, as can be concluded from Table <ref>. The absolute errors, the maximum values of which are collected in this table, were calculated against the reference values computed with the accelerated tetrachoric series and $128$-bit precision.

Double precision computation of $T(h,r)=\Phi_\varrho^2(h,0)-\frac{1}{2}\Phi(h)$ on the test grid

1l   2cTetrachoric   1c    2cVasicek   1c   2cNovel series(lr)2-3 (lr)5-6 (lr)8-9
1l   2caccelerated   1c    2caccelerated   1c   2catanExt1l   1cno   1cyes   1c   1cno   1cyes   1c   1cno   1cyes$h$ is a vector                          Max abs. errorab   9.42e-16   1.13e-16      9.37e-14   1.68e-16      9.68e-17   1.43e-16  Average iter. 1ab   123.2   47.2      3,575.4   65.3      66.9   22.8$r$ is a vector                          Max abs. error ba   4.77e-15   1.13e-16      9.37e-14   1.68e-16      9.68e-17   1.43e-16  Average iter. 1 ba   1,601.3   65.2      188,834.8   66.6      75.5   36.4for both                          Average iter. 2 ba   84.2   29.8      3,355.8   31.2      38.5   18.6  Maximum iter. ba   2,447   113      199,790   141      125   50

Since the vector computations were performed, there are two average numbers of iterations. The first is empirical because the iterations were terminated according to the “worst” component of the vector, and the second is hypothetical, if the computations were performed for each component separately. Since all iterations ran to the accuracy limit and with an unbounded number of iterations, the original Vasicek's series needed the enormous number of iterations, not only for the worse case, but also in average, when calculating with scalar $h$ and vector $r$. However, it is not intended for cases with $\varrho^2<\frac{1}{2}$, and its error can be well estimated, which was ignored here.

Comparative computations with competing packages were also performed on the test grid, using $T(h,r)=\Phi_\varrho^2(h,0)-\frac{1}{2}\Phi(h)$. The pbinorm() function from the VGAM package, i.e. using (<ref>), had an average and maximum absolute error of $1.12\cdot 10^{-8}$ and $3.62\cdot 10^{-7}$ respectively. Similar maximum absolute errors were obtained by the function pbvnorm() from the pbv package [Robitzsch, 2020], which is based on [Drezner and Wesolowsky, 1990], and pmvnorm() from the mvtnorm package, but only with the parameter algorithm = Miwa.

§.§ Numeric issues and results when computing Phi2(x,y)

In all $\Phi_\varrho^2(x,y)$ computational methods based on (<ref>) and (<ref>), or (<ref>) and (<ref>), exceptional cases with $\abs{\varrho}=1$, $x\ne 0$ and $x-y\,\sgn(\varrho)=0$ are not problematic because (<ref>) can be used, however, nearby cases are their weak point. Recalling that $q_x=q_y=\frac{x^2-2\varrho xy+y^2}{2\,(1-\varrho^2)}=q$ and using (<ref>) and (<ref>), it follows $\frac{\partial}{\partial\varrho}\, \Phi_\varrho^2(x,y)=\frac{e^{-q}}{2\pi\,\sqrt{1-\varrho^2}}$ and $\lim_{\abs{\varrho}\to 1}\frac{\partial}{\partial\varrho}\,\Phi_\varrho^2(x,y)=0$, except if $x-y\,\sgn(\varrho)=0$. If $\abs{\varrho}$ is close to $1$ and $x-y\,\sgn(\varrho)\simeq 0$, then $x^2-2\varrho xy+y^2\simeq 0$ and $q$ may be small enough that $\frac{\partial}{\partial\varrho}\,\Phi_\varrho^2(x,y)$ is moderate or even large. As a result, even a small rounding error in the calculation of key variables that depend on $\varrho$ can cause a much larger error in the result. The calculation of $r_x$ and $r_y$ with (<ref>) is particularly problematic. In a similar situation, <cit.> recommends an equivalent but less sensitive calculation by $r_x=\frac{x-y}{x\,\sqrt{1-\varrho^2}}- \sqrt{\frac{1-\varrho}{1+\varrho}}$ if $\varrho\approx 1$ and $r_x=\frac{x+y}{x\,\sqrt{1-\varrho^2}}-\sqrt{\frac{1+\varrho}{1-\varrho}}$ if $\varrho\approx -1$, and analogously for $r_y$. In our case, his proposal improves accuracy, but does not solve the problem sufficiently.

The problem could be solved by a hybrid computation along the lines of some other methods, e.g. the series from <cit.> could be used if $\abs{\varrho}\approx 1$ and $x-y\,\sgn(\varrho)\simeq 0$. Using the initial idea to develop these series, i.e. axis rotation in a way that removes the cross product term $-2\varrho st$ in (<ref>), by setting $s=\frac{u+v}{\sqrt{2}}$ and $t=\frac{u-v}{\sqrt{2}}$ the equation (<ref>) transforms to

\begin{equation*}
\begin{split}
\Phi_\varrho^2(x,y)
&=\tfrac{1}{2\pi\,\sqrt{1-\varrho^2}}\int_{-\infty}^\frac{x-y}{\sqrt{2}} e^{-\frac{u^2}{2\,(1-\varrho)}}\left(\int_{-\infty}^{u+\sqrt{2}\,y} e^{-\frac{v^2}{2\,(1+\varrho)}}\, \mathrm{d}v\right)\mathrm{d}u\\
&+\tfrac{1}{2\pi\,\sqrt{1-\varrho^2}}\int_\frac{x-y}{\sqrt{2}}^\infty e^{-\frac{u^2}{2\,(1-\varrho)}}\left(\int_{-\infty}^{-u+\sqrt{2}\,x} e^{-\frac{v^2}{2\,(1+\varrho)}}\, \mathrm{d}v\right)\mathrm{d}u
\end{split}
\end{equation*}

and, using $\frac{u}{\sqrt{1-\varrho}}=\tilde{u}$ and $\frac{v}{\sqrt{1+\varrho}}=\tilde{v}$, to

\begin{equation*}
\begin{split}
\Phi_\varrho^2(x,y)&=\int_{-\infty}^\frac{x-y}{\sqrt{2\, (1-\varrho)}}\varphi(\tilde{u})\,\Phi \left(\tilde{u}\sqrt{\tfrac{1-\varrho}{1+\varrho}}+ y\sqrt{\tfrac{2}{1+\varrho}}\,\right) \mathrm{d}\tilde{u}\\
&+\int_\frac{x-y}{\sqrt{2\,(1-\varrho)}}^\infty\varphi(\tilde{u})\, \Phi\left(-\tilde{u}\sqrt{\tfrac{1-\varrho}{1+\varrho}}+ x\sqrt{\tfrac{2}{1+\varrho}}\,\right) \mathrm{d}\tilde{u}.
\end{split}
\end{equation*}

By transforming $\tilde{u}\mapsto-\tilde{u}$ in the second integral and using <cit.>, it follows $\Phi_\varrho^2(x,y)=\Phi_{\varrho^*}^2\Big(\frac{x-y}{\sqrt{2\,(1-\varrho)}},y\Big)+ \Phi_{\varrho^*}^2\Big(-\frac{x-y}{\sqrt{2\,(1-\varrho)}},x\Big)$, where $\varrho^*=-\,\sqrt{\frac{1-\varrho}{2}}$. This equation is suitable for $\varrho\approx 1$. From it and $\Phi_\varrho^2(x,y)=\Phi(x)-\Phi_{-\varrho}^2(x,-y)$ we obtain the equation

\begin{equation*}
\Phi_\varrho^2(x,y)=\Phi(x)-\Phi_{\hat{\varrho}}^2\left(\tfrac{x+y} {\sqrt{2\,(1+\varrho)}},-y\right)-\Phi_{\hat{\varrho}}^2 \left(-\tfrac{x+y}{\sqrt{2\,(1+\varrho)}},x\right),
\end{equation*}

where $\hat{\varrho}=-\,\sqrt{\frac{1+\varrho}{2}}$, that is suitable for $\varrho\approx-1$. Both are combined in

\begin{equation}\label{Phirxy6}
\Phi_\varrho^2(x,y)=\tfrac{1-\sgn(\varrho)}{2}\,\Phi(x)+\sgn(\varrho)\left(
\Phi_{\tilde{\varrho}}^2(z,\tilde{y})+\Phi_{\tilde{\varrho}}^2(-z,x)\right),
\end{equation}

where $\tilde{\varrho}=-\,\sqrt{\frac{1-\abs{\varrho}}{2}}$, $\tilde{y}=y\,\sgn(\varrho)$ and $z=\frac{x-\tilde{y}}{\sqrt{2\,(1-\abs{\varrho})}}$. Noting that

\begin{equation*}
\tilde{q}=\tfrac{\tilde{y}^2-2\tilde{\varrho}\tilde{y}z+z^2} {2\,(1-\tilde{\varrho}^2)}= \tfrac{x^2+2\tilde{\varrho}xz+z^2}{2\,(1-\tilde{\varrho}^2)}=\tfrac{x^2-2\varrho xy+y^2}{2\,(1-\varrho^2)}=q,
\end{equation*}

hence $\varphi_{\tilde{\varrho}}^2(z,\tilde{y})=\varphi_{\tilde{\varrho}}^2(-z,x)= \frac{e^{-q}}{\pi\,\sqrt{2\,(1+\abs{\varrho}})}$, we stop here because one potentially critical case is replaced with two non-critical ones.

If we continued to simplify the right side of (<ref>) by (<ref>), we would get (<ref>). This actually happens when calculating with the functions used in tests, but it turns out that the results are more accurate than if we start with (<ref>) if $\abs{\varrho}\approx 1$ and $x-y\,\sgn(\varrho)\simeq 0$. Of course, the equation (<ref>) could be used as a starting point instead of (<ref>) in all cases, but this is not rational because it requires $4$ calculations of the $T$ function instead of $2$, two of which cancel each other out, at least in theory, and as a result, in non-critical cases the error may be larger than when starting with (<ref>). When testing, (<ref>) was used as a starting point instead of (<ref>) if $\varphi_\varrho^2(x,y)>1$. However, since $\min\left\{\abs{\varrho},\sqrt{\frac{1-\abs{\varrho}}{2}}\right\}\le\frac{1}{2}$, using (<ref>) may also make sense in some other circumstances, e.g. when computing directly with (<ref>).

In this case, testing was performed on the randomly generated set of parameters. After setting the random generator seed to $123$ to allow replication and independent verification, one million $x$ values, one million $y$ values, and one million $\varrho$ values were drawn from a uniform distribution on the intervals $(-10,10)$, $(-10,10)$, and $(-1,1)$ respectively.

For each triplet $(x,y,\varrho)$, auxiliary values
$q=\frac{x^2-2\varrho xy+y^2}{2\,(1-\varrho^2)}$ and $\varphi_\varrho^2(x,y)=\frac{e^{-q}}{2\pi\,\sqrt{1-\varrho^2}}$ were calculated. For the new series, the tested $\Phi_\varrho^2(x,y)$ values were computed using (<ref>) as a branch point and (<ref>) for both branches if $\varphi_\varrho^2(x,y)>1$, and using (<ref>) directly otherwise. The same applies for the tetrachoric and Vasicek's series, only (<ref>) must be used instead of (<ref>). However, the described procedure is built in the Phi2xy() function. The tested $\Phi_\varrho^2(x,y)$ values needed for the preparation of Table <ref> were computed with:

row 1: Phi2xy( x, y, rho, fun = "tOwenT" );

row 2: Phi2xy( x, y, rho, fun = "vOwenT" );

row 3: Phi2xy( x, y, rho, fun = "mOwenT", opt = FALSE );

row 4: Phi2xy( x, y, rho, fun = "mOwenT", opt = TRUE ).

The absolute errors, which statistics are collected in Table <ref>, were calculated against reference values computed with the accelerated tetrachoric series as described, but with $128$-bit precision. From a user's point of view, only the parameters $x$, $y$ and $\varrho$ needed to be converted into $128$-bit “mpfr-numbers” before calling Phi2xy( x, y, rho, fun = "tOwenT" ), hence  and analogously for $y$ and $\varrho$. In this case, the computation took significantly more time, since all intermediate and final results were $128$-bit “mpfr-numbers”.

$\Phi_\varrho^2(x,y)$ computation – absolute error comparison

1lMethod   1cMedian   1cMean   1c3rd Qu.   1c99%   1cMax.Tetrachoric series (accelerated)   2,86e-17   3,68e-17   5,37e-17   1,51e-16   2,77e-16Vasicek's series (accelerated)   2,57e-17   3,78e-17   5,90e-17   1,65e-16   3,38e-16Novel series (atanExt = no)   2,86e-17   3,68e-17   5,38e-17   1,51e-16   2,77e-16Novel series (atanExt = yes)   3,36e-17   4,20e-17   6,11e-17   1,56e-16   3,45e-16

In this test, $\varrho_{min}\approx-0.99999991$, $\varrho_{max}\approx 0.99999775$, there are $5$ cases with $\varphi_\varrho^2(x,y)>1$ and $\max\varphi_\varrho^2(x,y)\approx 1.38$.

In order to better check the accuracy if $\abs{\varrho}$ is close to $1$, the computation was repeated with the same $x$ and $y$, and $\varrho^*=2\,\Phi(8\varrho)-1$. In this case, more than half of $\abs{\varrho^*}$ are greater than $0.9999$ and none is equal to $1$. The results are in Table <ref>.

$\Phi_{\varrho^*}^2(x,y)$ computation – absolute error comparison

1lMethod   1cMedian   1cMean   1c3rd Qu.   1c99%   1cMax.Tetrachoric series (accelerated)   2,92e-17   3,83e-17   5,55e-17   1,60e-16   2,79e-16Vasicek's series (accelerated)   1,72e-17   3,38e-17   5,38e-17   1,63e-16   3,01e-16Novel series (atanExt = no)   2,92e-17   3,83e-17   5,55e-17   1,60e-16   2,79e-16Novel series (atanExt = yes)   3,01e-17   3,91e-17   5,55e-17   1,56e-16   2,95e-16

In this test, there are $178$ cases with $\varphi_\varrho^2(x,y)>1$ and $\max\varphi_\varrho^2(x,y)\approx 376$. For both novel series in Table <ref>, the maximum absolute error is $1.54\cdot 10^{-14}$ if (<ref>) is not used, improving to $1.21\cdot 10^{-15}$ if [Meyer, 2013]'s advice is followed. Similarly applies to the other two series.

Due to the time-consuming $128$-bit precision computation, the presented tests are the only performed tests based on the $128$-bit precision reference values. However, additional tests were made by comparing the results of Phi2xy() with the results obtained by other functions from packages available on CRAN. For three functions, the accuracy was found to be worse than $2\cdot 10^{-7}$ as measured by the maximum absolute error on the test grid from Subsection <ref>. Regarding accuracy, only the function pmvnorm() with the parameter algorithm = TVPACK from the mvtnorm package proved to be equivalent to Phi2xy() and even slightly more accurate with the maximum absolute errors $2.58\cdot 10^{-16}$ and $2.19\cdot 10^{-16}$ on the test sets of triplets $(x,y,\varrho)$ and $(x,y,\varrho^*)$ respectively. The same function with algorithm = GenzBretz achieved $2.58\cdot 10^{-16}$ and $2.20\cdot 10^{-11}$.

The OwenQ::OwenT() function also proved to be equivalent to Phi2xy() as measured by the maximum absolute error on the test sets of triplets, but only when upgraded to compute $\Phi_\varrho^2(x,y)$ in the manner used in Phi2xy(). This function and Phi2rho::OwenT() are essentially wrappers for the internal functions tOwenT(), vOwenT() and mOwenT(), which as workhorses compute series.

In terms of reliability, stability and accuracy of the new methods, no problems or significant differences according to the parameters from the different data sets were detected, except those already described and related to $\abs{\varrho}\approx 1$ and $x-y\,\sgn(\varrho)\simeq 0$. To eliminate them, the equation (<ref>) was derived, which also proved to be successful when using the OwenQ::OwenT() function, upgraded to compute $\Phi_\varrho^2(x,y)$. Regarding parameters from different data sets, the only detected big difference is in the number of iterations, as can be seen from Figure <ref> in the next subsection.

Two million $\Phi_\varrho^2(x,y)$ computations on Windows 10 and 64-bit R 4.3.1 on the old Intel i7-6500U CPU @ 2.50 GHz and 8 GB RAM, using only one thread of four, lasted $26$, $31$ and $286$ seconds for Phi2xy(), upgraded OwenQ::OwenT() and pmvnorm() with algorithm = TVPACK respectively. The $128$-bit precision reference values computation took over $17$ hours, compared to $43$ seconds for the double precision one. Based on Tables <ref> and <ref>, we can conclude that the increased number of iterations has only an insignificant effect on the huge time extension factor.

§.§ High-precision computation

All functions in the Phi2rho package are ready to use the Rmpfr package, which enables using arbitrary precision numbers instead of double precision ones and provides all the high-precision functions needed. It interfaces R to the widely used C Library MPFR [Fousse et al., 2007]. Assuming that the Rmpfr package is loaded, the functions should only be called with the parameters, which are “mpfr-numbers” of the same precision. All $128$-bit precision benchmark values used in Subsections <ref> and <ref> were calculated using Rmpfr.

To get a sense of how the number of iterations depends on the required precision, the test from Subsection <ref> was partially repeated using the $128$-bit precision computation. The results are collected in Table <ref>. Note that $2^{-128}\approx 2.94\cdot 10^{-39}$ and that in this case the values in the table are not the maximum absolute errors because the comparison values are computed with the same precision.

$128$-bit precision computation of $T(h,r)=\Phi_\varrho^2(h,0)-\frac{1}{2}\Phi(h)$ on the test grid

1l   1cTetrachoric   1c   1cVasicek   1c   2cNovel series(lr)2-2 (lr)4-4 (lr)6-7
1l   1c(accelerated)   1c   1c(accelerated)   1c   1c(atanExt = no)   1c(atanExt = yes)$h$ is a vector                    Max abs. errorab   benchmark      6.61e-39      5.51e-39   3.88e-39  Average iter. 1ab   82.8      75.5      105.7   55.0$r$ is a vector                    Max abs. error ba   0.00e+00      6.61e-39      5.51e-39   3.88e-39  Average iter. 1 ba   132.9      120.9      149.2   72.5for both                    Average iter. 2 ba   63.2      57.7      73.4   40.3  Maximum iter. ba   191      144      199   119

From Tables <ref> and <ref> can be concluded that there is no significant difference in accuracy between the compared series if the accelerated series are considered, but one of the new series converges significantly faster than the others. It can also be concluded that the Vasicek's series with $\overline{\varrho}$ instead of $\varrho$ is a reasonable alternative to the tetrachoric series (<ref>) also for $\varrho^2<\frac{1}{2}$ and remains a reasonable alternative for $\varrho^2>\frac{1}{2}$ if compared to the transformed tetrachoric series (<ref>).

The number of iterations for the computation with the new series with double and $128$-bit precision can be compared in Figure <ref>. Only the first quadrant is presented because others are its mirror images.

$T(h,r)$ computation with the new series (atanExt = yes) with double (left) and $128$-bit (right) precision – number of iterations

We also provide examples with $r=1$ and $r=-1$ for which reference values can be computed alternatively. $A=\Phi_{\sqrt{2}/2}^2(2.1,0)$ and $B=\Phi_{-\sqrt{2}/2}^2(2.1,0)$ were computed by (<ref>) and (<ref>), and by the right side of (<ref>) by the pnorm() function from the package Rmpfr as a benchmark. In all cases, the latter was computed by double the precision expressed in bits and used for the former. Both reference values were used to compare how quickly the compared series converge when computing with more than $128$-bit precision. In Tables <ref> and <ref>, for each of them the absolute error and the number of iterations, which is the same for both computations, is presented. Since $\abs{r}=1$, there is no difference between the original and accelerated tetrachoric series, and the same applies to Vasicek's series.

A special case of high-precision computing – absolute error comparison

1c   1c   1c   3cTetrachoric   1c   3cVasicek(lr)4-6 (lr)8-10
1rBits   1rPrecision   1c    1c1A   1c1B    1r1n   1c   1c2A    1c2B   1r2n53   1.11e-16      6.98e-17   8.65e-17   46      6.98e-17   3.10e-17   4764   5.42e-20      6.89e-21   2.26e-20   57      6.89e-21   2.26e-20   57128   2.94e-39      5.47e-40   5.89e-41   119      5.47e-40   5.89e-41   120256   8.64e-78      6.29e-78   3.33e-78   246      6.29e-78   3.33e-78   246512   7.46e-155      3.65e-156   1.37e-155   500      4.09e-155   1.37e-155   5011024   5.56e-309      5.04e-309   2.41e-309   1011      2.25e-309   5.19e-309   1011

A special case of high-precision computing – absolute error comparison

1c   1c   1c   3cNovel series (atanExt = no)   1c   3cNovel series (atanExt = yes)(lr)4-6 (lr)8-10
1rBits   1rPrecision   1c    1c3A   1c3B    1r3n   1c   1c4A    1c4B   1r4n53   1.11e-16      6.98e-17   3.10e-17   54      6.98e-17   3.10e-17   2264   5.42e-20      3.40e-20   4.54e-21   65      6.89e-21   2.26e-20   25128   2.94e-39      5.47e-40   5.89e-41   128      9.23e-40   1.41e-39   41256   8.64e-78      6.29e-78   3.33e-78   256      6.29e-78   3.33e-78   69512   7.46e-155      4.09e-155   1.37e-155   511      4.09e-155   1.37e-155   1161024   5.56e-309      2.25e-309   5.19e-309   1023      5.04e-309   2.41e-309   199

Time-consuming high precision testing was less intensive than double precision one. Since there is no competing package on CRAN ready for high-precision computation of $\Phi_\varrho^2(x,y)$, it could not be supplemented by comparison computations with competing packages.

§ CONCLUSION

If the generally applicable Theorem <ref> is not considered, Theorems <ref> and <ref>, and Corollary <ref> are interesting mainly theoretically, because Theorem <ref>, together with (<ref>) and (<ref>), replaces and supplements them in practice. The series (<ref>)–(<ref>) enable fast and numerically stable computation, which is often more important than the speed of convergence. From Tables <ref> and <ref> can be concluded that the computation with (<ref>) or (<ref>) needs significantly fewer iterations than the computation with the accelerated tetrachoric and Vasicek's series, and from Tables <ref> and <ref> can be assumed that the advantage increases with increasing precision. For $53$-bit precision (double precision), it required a half, and for $1024$-bit precision only a fifth, of the iterations demanded by the competing series. However, the cost of calculating $\arctan r$ is not taken into account. In connection with this, let us just mention two more questions.

Recalling that $e^{-x}I_n(x)=P(n+1,x)$, the Taylor series of the arctangent function and (<ref>) imply

\begin{equation*}
T(h,r)=\tfrac{\arctan r}{2\pi}-\tfrac{1}{2\pi}\sum_{k=0}^\infty\tfrac{(-1)^k r^{2k+1}}{2k+1}\;P(k+1,\tfrac{1}{2}h^2)=\tfrac{1}{2\pi}\sum_{k=0}^\infty\tfrac{(-1)^k r^{2k+1}}{2k+1}\;Q(k+1,\tfrac{1}{2}h^2),
\end{equation*}

which <cit.> have already noticed. If the above series, (<ref>) and (<ref>) are viewed as functions of the variable $h$, they have a similar structure, but the coefficients, depending on $r$, belong to different arctangent series. Whether the similarity could be explained by a similar Euler transform as the one which transforms the Taylor arctangent series to the Euler's arctangent series seems interesting question, but was not deeply explored.

In Recursion <ref>, the external calculation of the arctangent function is avoided, assuming that the Euler's series (<ref>) is used for it, and the expectation that there would be no difference in the number of iterations for Recursions <ref> and <ref> if those for calculating the arctangent were also taken into account. However, the use of the external arctangent calculation allows the use of already known and possible future faster converging series. Such is the case with

\begin{equation}\label{Eq:arctan}
\arctan r=\tfrac{r}{\sqrt{1+r^2}}\sum_{k=0}^\infty\tfrac{(2k-1)!!}{(2k)!!\, (2k+1)}\left(\tfrac{r^2}{1+r^2}\right)^k\quad(r\in\mathbb{R}),
\end{equation}

which is based on the Taylor series of $\arcsin r$ and $\arctan r=\arcsin\frac{r}{\sqrt{1+r^2}}$. Indeed, the number of iterations is not significantly smaller if $\abs{r}\le 1$ and the accuracy obtained by double precision computation is sufficient, and even the square root must be calculated, but it could be significantly smaller in high-precision computation. The Euler's series (<ref>) is included in [Abramowitz and Stegun, 1972] and [Olver et al., 2010], but (<ref>) is not, even though it converges faster. It can be found in <cit.> as

\begin{equation*}
\arctan r=\tfrac{r}{\sqrt{1+r^2}}\sum_{k=0}^\infty\tfrac{(2k)!}{2^{2k}(k!)^2(2k+1)} \left(\tfrac{r^2}{1+r^2}\right)^k
\end{equation*}

with a reference to the 1922 source. Due to equations $(2k-1)!!=\frac{(2k)!}{2^k k!}$ and $(2k)!!=2^k k!$, the coefficients are the same as those in (<ref>). Whether we could also find a corresponding series for the Owen's $T$ function, which would converge faster than (<ref>), remains an open challenge.

[Abramowitz and Stegun, 1972]
Abramowitz, M., Stegun, I. A. (eds.) (1972).
Handbook of Mathematical Functions.
9th edition. Dover Publications, Inc., New York.

[Aǧca and Chance, 2003]
Aǧca, Ş., Chance, D. M. (2003).
Speed and Accuracy Comparison of Bivariate Normal
Distribution Approximations for Option Pricing.
The Journal of Computational Finance, 6(4), 61–96.

[Alzer, 1997]
Alzer, H. (1997).
On Some Inequalities for the Incomplete Gamma Function.
Mathematics of Computation, 66(218), 771–778.

[Alzer, 1998]
Alzer, H. (1998).
Inequalities for the Chi Square Distribution Function.
Journal of Mathematical Analysis and Applications,
223(1), 151–157.

[Divgi, 1979]
Divgi, D. R. (1979).
Calculation of Univariate and Bivariate Normal Probability
The Annals of Statistics, 7(4), 903–910.

[Donnelly, 1973]
Donnelly, T. G. (1973).
Algorithm 462: Bivariate Normal Distribution.
Communications of the ACM, 16(10), 638.

[Drezner and Wesolowsky, 1990]
Drezner, Z., Wesolowsky, G. O. (1990).
On the Computation of the Bivariate Normal Integral.
Journal of Statistical Computation and Simulation,
35(1-2), 101–107.

[Fayed and Atiya, 2014]
Fayed, H. A., Atiya, A. F. (2014a).
An Evaluation of the Integral of the Product of the Error
Function and the Normal Probability Density With Application to
the Bivariate Normal Integral.
Mathematics of Computation, 83(285), 235–250.

[Fayed and Atiya, 2014]
Fayed, H. A., Atiya, A. F. (2014b).
AnNovel Series Expansion for the Multivariate
Normal Probability Integrals Based on Fourier Series.
Mathematics of Computation, 83(289), 2385–2402.

[Fousse et al., 2007]
Fousse, L., Hanrot, G., Lefèvre, V., Pélissier, P., Zimmermann, P.
MPFR: A Multiple-precision Binary Floating-Point
Library With Correct Rounding.
ACM Transactions on Mathematical Software, 33(2),

[Gai, 2001]
Gai, J. (2001).
A Computational Study of the Bivariate Normal
Probability Function.
Master's thesis, Queen's University, Kingston.

[Gautschi, 1998]
Gautschi, W. (1998).
The Incomplete Gamma Functions Since Tricomi.
In Tricomi's Ideas and Contemporary Applied
Mathematics, Atti Convegni Lincei, vol. 147, pp. 203–237. Accad. Naz.
Lincei, Rome.

[Genz, 2004]
Genz, A. (2004).
Numerical Computation of Rectangular Bivariate and
Trivariate Normal and $t$ Probabilities.
Statistics and Computing, 14(3), 251–260.

[Genz and Bretz, 2009]
Genz, A., Bretz, F. (2009).
Computation of Multivariate Normal and $t$
Lecture Notes in Statistics. Springer-Verlag, Heidelberg.
ISBN 978-3-642-01688-2.

[Gradshteyn and Ryzhik, 2015]
Gradshteyn, I. S., Ryzhik, I. M. (2015).
Table of Integrals, Series, and Products.
8th edition. Academic Press, Amsterdam.

[Gupta, 1963]
Gupta, S. S. (1963).
Bibliography on the Multivariate Normal Integrals and
Related Topics.
The Annals of Mathematical Statistics, 34(3),

[Kendall, 1941]
Kendall, M. G. (1941).
Proof of Relations Connected With the Tetrachoric Series
and Its Generalization.
Biometrika, 32(2), 196–198.

[Komelj, 2023]
Komelj, J. (2023).
Phi2rho: Owen's $T$ Function and Bivariate Normal
R package version 1.0,

[Laurent, 2023]
Laurent, S. (2023).
OwenQ: Owen Q-Function.
R package version 1.0.7,

[Maechler, 2023]
Maechler, M. (2023).
Rmpfr: R MPFR – Multiple Precision
Floating-Point Reliable.
R package version 0.9-3,

[Marsaglia, 2004]
Marsaglia, G. (2004).
Evaluating the Normal Distribution.
Journal of Statistical Software, 11(4), 1–11.

[Meyer, 2013]
Meyer, C. (2013).
Recursive Numerical Evaluation of the Cumulative Bivariate
Normal Distribution.
Journal of Statistical Software, 52(10), 1–14.

[Mitrinović, 1970]
Mitrinović, D. S. (1970).
Analytic Inequalities.
Springer-Verlag, Berlin.

[Olver et al., 2010]
Olver, F. W. J., Lozier, D. W., Boisvert, R. F., Clark, C. W. (eds.) (2010).
NIST Handbook of Mathematical Functions.
National Institute of Standards and Technology & Cambridge
University Press, Cambridge, New York.

[Owen, 1956]
Owen, D. B. (1956).
Tables for Computing Bivariate Normal Probabilities.
The Annals of Mathematical Statistics, 27(4),

[Owen, 1980]
Owen, D. B. (1980).
A Table of Normal Integrals.
Communications in Statistics – Simulation and Computation,
9(4), 389–419.

[Patefield and Tandy, 2000]
Patefield, M., Tandy, D. (2000).
Fast and Accurate Calculation of Owen's T Function.
Journal of Statistical Software, 5(5), 1–25.

[Plackett, 1954]
Plackett, R. L. (1954).
A Reduction Formula for Normal Multivariate Integrals.
Biometrika, 41(3/4), 351–360.

[R Core Team, 2023]
R Core Team (2023).
R: A Language and Environment for
Statistical Computing.
R Foundation for Statistical Computing, Vienna, Austria.

[Robitzsch, 2020]
Robitzsch, A. (2020).
pbv: Probabilities for Bivariate Normal
R package version 0.4-22,

[Terza and Welland, 1991]
Terza, J. V., Welland, U. (1991).
A Comparison of Bivariate Normal Algorithms.
Journal of Statistical Computation and Simulation,
39(1-2), 115–127.

[Vasicek, 1998]
Vasicek, O. A. (1998).
A Series Expansion for the Bivariate Normal Integral.
The Journal of Computational Finance, 1(4), 5–10.

[Yee, 2023]
Yee, T. W. (2023).
VGAM: Vector Generalized Linear and Additive
R package version 1.1-8,