Coverage and Rate of Downlink Sequence Transmissions with
Reliability Guarantees

arXiv:1704.05296v2 [] 24 Jun 2017

Jihong Park and Petar Popovski
Abstract—Real-time distributed control is a promising application of 5G in which communication links should satisfy
certain reliability guarantees. In this letter, we derive closedform maximum average rate when a device (e.g. industrial
machine) downloads a sequence of n operational commands
through cellular connection, while guaranteeing a certain signalto-interference ratio (SIR) coverage for all n messages. The result
is based on novel closed-form n-successive SIR coverage bounds.
The proposed bounds provide simple approximations that are
increasingly accurate in the high reliability region.
Index Terms—Sequence transmissions, ultra-reliable communication, coverage, rate, stochastic geometry.

I. M OTIVATION AND C ONTRIBUTION

previous
commands

current
commands

interference

static machine

mobile machine

Fig. 1. A sequence of 3 command receptions for mobile and static UEs.

UE locations, the n-successive SIR coverage has been known
as (Theorem 2, [4]):

One of the focus in the emerging 5G wireless systems
will be communication with reliability guarantees. The most
challenging variant is ultra-reliable low-latency communication (URLLC), with extreme reliability guarantees. A representative use case for URLLC is industrial Internet of Things
(IoT) [1], in which a certain user equipment (UE), such as an
actuator, is controlled in real time by downloading operational
commands for high-precision cloud manufacturing [2]. A
loss of a command sequence may cause critical malfunction
or shutdown of the process. Therefore, the network should
guarantee wireless coverage with very high reliability, e.g.
≥ 99% in order to decode the n successive commands.
In this letter we rely on stochastic geometry to derive nsuccessive SIR (Signal-to-Interference Ratio) coverage pn (t)
for a target SIR threshold t and n ≥ 1, defined as pn (t) :=
Pr(SIR1 ≥ t, SIR2 ≥ t, · · · , SIRn ≥ t). We thereby derive
the closed-form maximum average rate Rn (η) that ensures η
reliability for n-successive message receptions, i.e. pn (t) ≥ η.
Regardless of whether the reliability is extreme as in URLLC
or more relaxed, the introduction of a reliability guarantee
poses a challenging problem. This precludes us from directly
exploiting previously known closed-form SIR coverage and
rate expressions [3]–[5]. The reason comes from the following
two issues: spatial correlation during n-successive receptions
and the trade-off between coverage and rate.
Regarding spatial correlation, if UE locations are independent and identically distributed (i.i.d.) during n message
receptions, then pn (t) = p1 (t)n straightforwardly.1 However,
a URLLC UE should receive n messages within a short
time span, during which all UEs and the interferers can be
considered static, see Fig. 1. For a cellular downlink with fixed

Although (2) is a simple convex optimizaition with an inequality condition, it cannot be analytically solved due to the
intractable hypergeometric function of pn (t) in (1).
A suitable way to avoid the difficulties brought by the hypergeometric function is to utilize more tractable SIR coverage

J. Park and P. Popovski are with Department of Electronic Systems, Aalborg
University, Denmark (email: {jihong, petarp}@es.aau.dk).
This work has been in part supported by the European Research Council
(ERC Consolidator Grant Nr. 648382 WILLOW) within the Horizon 2020
Program, and partly by Innovation Fund Denmark via the Virtuoso project.

1 For a typical UE’s SIR, i.i.d. mobility is almost surely identical to
any Markovian mobility such as a random walk and Lévy flight thanks to
displacement theorem [6] as utilized in [7] (see [8] for general
mobility).

2
2
2 A special case is p (t) = 1/ F
1 R
2 1 1, − α ; 1 − α ; −t in [4], which is
identical to 1/(1 + t2/α u>t−2/α du/[1 + uα/2 ]) presented in [3].

pn (t) = 1/2 F1 (n, −2/α; 1 − 2/α; −t)

(1)

2

at a randomly picked (i.e. a typical) UE. The denominator
2 F1 (a, b; c; z) is a Gauss hypergeometric function that can
only be numerically computed, defined as 2 F1 (a, b; c; z) :=
P∞ Γ(a+n)Γ(b+n)Γ(c) zn
n=0 Γ(a)Γ(b)Γ(c+n) n! where Γ(x) denotes the gamma
function. It is thus difficult to establish a direct relation
between n and pn (t) due to the analytical intractability of
the hypergeometric function.
Next, the average rate expression E[log(1+SIR)] in previous
works [3]–[5] bore no relationship to the reliability constraint
pn (t) ≥ η. The reason is because in those works it was
implicitly assumed that each BS knows the full channel state
information (CSI) of all the desired and interfering links. This
enables the BS to adaptively adjust its modulation and coding
scheme (MCS) such that t = SIR resulting in pn (SIR) = 1.
The full CSI unfortunately cannot be exchanged within a low
target latency, and is therefore not affordable for URLLC.
Without CSI, the SIR threshold is kept constant during n
transmissions, while satisfying pn (t) ≥ η. This induces the
trade-off between a typical UE’s coverage pn (t) and its quasiconcave average rate pn (t) log(1 + t) with respect to t. We
focus on the maximum average rate Rn (η) given as
Rn (η) := arg max pn (t) log(1 + t) s.t. pn (t) ≥ η.

(2)

t≥0

n=3, a=4

n=3, a=4
11

1

0.7
0.6

+

0.5

LBa , UBa Lemma 1
LBb , UBb Lemma 2
LBc , UBc Theorem 1
LB⇥
Corollary 1
LB+

0.99

Exact
Bound in [5]

0.4
0.3

0.97
0.6
0.96

0.93

0.91
0.9
0

0.5
0.5 0.8

SIR Coverage, pn (t)

0

1

1.51

1.5 2.5

3 2

3.5

2.5 4

4.5 3

0
0.9
00

5

LB
0.5
0.005

0.01

1

0.015

1.5 0.02

LBc

LB

+

b

2
0.025

LB
2.5
0.03

⇥

0.035 3

0.04

SIRSIR
Threshold,
t t
Threshold,

SIR
Threshold,
t t
SIR
Threshold,

0.7
0.6

2

Exact

0.4
0.94

1
0.1

UBc

0.95

0.2
0.92

0.2

UBa

0.8
0.98

Bound in [5]

SIR
pnn(t)
(t)
SIR Coverage,
Coverage, pp
(t)
n

0.8

+

(t)
SIR Coverage, ppn(t)
n

0.9

(a) for the range within [0, 1].

(b) for the range within [0.9, 1].

0.5

Fig. 2. 3-successive
SIR coverage bounds (Theorem 1, α = 4), derived via the bounds after Pfaff transformation (Lemma 1) and the asymptotic bounds
0.4
(Lemma 2). They are compared to the existing bounds 1/(1 + cn t2/α ) in [5] and cn −1 t−2/α (= UBa ) in [9], [10]. The exact curve follows from (1).
0.3

+

0.2
0.1
bounds, as done
in our recent work [11]. Nonetheless, the
0
0
1
1.5
3
4
bounds in [11]
are0.5limited
only
to2 the 2.5
case n
= 3.51. Other
SIR Threshold,
existing SIR coverage bounds [5], [9],
[10] aret too loose for
ultra-reliable scenarios, within the range [0.9, 1] for low SIR
threshold t, as Fig. 2 shows.
In this letter we propose novel closed-form bounds of nsuccessive SIR coverage pn (t) (see LBc and UBc in Theorem 1), and thereby derive the closed-form bounds of maximum rate Rn (η) guaranteeing pn (t) ≥ η (Proposition 1).
For high reliability η ≥ 0.9, the closed-form results show that
the impact of spatial UE correlation during n receptions is
negligibly small. In this case the proposed lower bound under
i.i.d. UE locations tightly approximates the exact rate (LB× in
Corollaries 1 and 2). For moderate reliability 0.3 < η < 0.9,
a slightly modified lower bound becomes accurate as an
approximation (LB+ in Corollaries 1 and 2).

II. S YSTEM M ODEL
Consider a downlink cellular network in a two-dimensional
plane. The network comprises uniformly located: BSs with
density λb and UEs with density λu . This leads to two
independent homogeneous Poisson point processes Φ(λb ) and
Φ(λu ), representing BS and UE locations. Each UE associates
with its nearest BS, and requests n messages. All messages
are of unit size. UE locations are fixed during the n-successive
message receptions, and randomly and independently changed
for the next n-successive message receptions.
Each BS serves a single UE with unit transmission power
and frequency bandwidth per time slot. We assume that
λu  λb such that any BS has at least a single associated
UE to be served. Transmitted signals experience path loss
attenuation with the exponent α > 2 and Rayleigh fading. A
typical UE is located at the origin, capturing any locations
without SIR changes according to Slyvnyak’s theorem [6].
The i-th message, 1 P
≤ i ≤ n, at the typical UE is received
with SIRi := g i r−α / rj >r gji rj −α , where rj := |xj | denotes
the j-th nearest interfering BS distance for its location xj ∈
Φ(λb ). Fading power gji at the i-th reception is exponentially
distributed with unit mean. We drop the subscript j for the

desired signal link. When i changes from 1 to n, the fading
5
powers
g i and gji vary with each i, while the desired BS and
interferer distances r and rj are fixed. In this case, pn (t) has
been known as (1) according to Theorem 2 in [4].

4.5

III. n-S UCCESSIVE SIR C OVERAGE B OUNDS
In this section we propose closed-form bounds of nsuccessive SIR coverage pn (t). The bounds are not only much
tighter for URLLC but also as simple as the existing bounds.
The desired bounds are derived by utilizing a transformed
function of pn (t) and pn (t)’s asypmtotic curves as follows.
We define the lower and the upped bounds of pn (t) through
UBi
i
pLB
n (t) ≤ pn (t) ≤ pn (t) for i ∈ {a, b, c, ×, +}. We then
apply a hypergeometric function transform [12] to (1), yielding
the following lemma.
Lemma 1. (Pfaff Transformed Bounds) At a static typical UE,
pn (t) has the following lower and upper bounds.
2

−1
a
pLB
(1 + t)− α
n (t) = cn
a
pUB
n (t)

(3)

2
−α

= (1 + t)
(4)

where cn := Γ 1 − α Γ n + α2 /Γ(n)
Proof: By using Pfaff’s hypergeometric
function transform

2
[12], 2 F1 n, − α2 ; 1− α2 , −t = (1 + t) α G(α, n,
 t) where
2

2
2
t
G(α, n, t) := 2 F1 1 − α + n , − α ; 1 − α ; 1+t that has
the range [1, cn ]. Applying this to (1) leads to the result. 

2

The term cn ≥ 1 is a monotone increasing/decreasing
function of n and
 α respectively. Its minimum value is thus
2π α→∞
csc
= 1.
c1 = 2π
α
α
It is remarkable that both bounds identically share the
2
simple term (1 + t)− α only with the constant multiplication
difference. Nevertheless, they have the following limitations.
−1
a
First, pLB
] of which the
n (t) has the limited range [0, cn
maximum is much lower than our target minimum reliability
a
0.9 under 2 < α ≤ 6 for all n ≥ 1. Second, pUB
n (t) is
independent of n, and thus becomes less accurate for larger n.
For these reasons, we aim at further improving the bounds
in Lemma 1 such that the bounds behave similar to the original

pn (t) with respect to n and t under the entire range [0, 1] while
converging to the asymptotic curve of pn (t) for t → 0. To
this end, we first derive the asymptotic bounds as follows.
Lemma 2. (Asymptotic Bounds) At a static typical UE, nsuccessive SIR coverage pn (t) follows the asymptotic curves:
2n
t
α−2
2
−1 − α
b
pUB
t .
n (t) = cn
b
pLB
n (t) = 1 −

(5)
(6)

Proof: Applying Taylor expansion to (1) with the hypergeometric function definition, its first derivative satisfies:
limt→0 pn (t)0 = −(2nt)/(α − 2) and limt→∞ pn (t)0 =
2
t− α /cn . By the hypergeometric function definition, pn (0) = 1
and limt→∞ pn (t) = 0. The asymptotic curves thus become
UBb
b
pLB
respectively as t → 0 and ∞.
n (t) and pn
To prove the asymptotic curves become the bounds of
pn (t), we further identify the behavior of pn (t) by directly considering its first derivative given as: pn (t)0 =
−2pn (t) [1 − (1 + t)−n pn (t)] /(αt). This shows pn (t)0 < 0
for all t < ∞. Since pn (t) is continuous by the hypergeometric function definition, its functional and derivative
values at the extreme t’s indicate pn (t) is a function that is
(strictly) decresing as well as (strictly) convex for all t (strict
b
conditions hold for t < ∞). The asymptotic curve pLB
n (t)
LBb
shares pn (0) = pn (0) = 1 while having limited domain,
b
thus becoming the lower bound. The curve pUB
n (t) satisfies
UBb
limt→∞ pn (t) = limt→∞ pn (t) = 0 while having infinite
maximum range as t → 0, thus becoming the upper bound. 
b
Note that pUB
n (t) in Lemma 2 generalizes its special case
that has been derived in [9], [10]. It is also noted
a
that the previous Lemma 1’s pLB
n (t) converges to the derived
UBb
asymptotic curve pn (t) as t → ∞.
b
Next, we combine Lemma 1 and pLB
n (t) in Lemma 2,
resulting in the following Theorem.
b
pUB
1 (t)

Theorem 1. (Coverage Bounds) At a static typical UE, the
n-successive SIR coverage pn (t) is bounded by:
−

2

c
α−2
pLB
n (t) = (1 + nt)

− α2
nα
UBc
pn (t) = 1 +
t
.
α−2

(7)
(8)
2

Proof: We choose the common term (1 + t)− α in Lemma 1
as the baseline, which is also a strictly decreasing and convex
function for t < ∞. We modify this, and aim at its convergence to the asymptotic curve 1 − 2nt/(α − 2) in Lemma 2 as
t → 0. For this purpose, we consider its generalized function
2B
(1 + At)− α for constants A, B > 0, and consider its Taylor
expanded form 1 − 2ABt/α for an infinitesimal t, which
should satisfy AB = nα/ (α − 2) for the desired convergence.
We select {A, B} = {n, α/ (α − 2)} and {nα/ (α − 2) , 1}
c
c
respectively yielding pLB
and pUB
n
n .
c
To verify the lower bound, showing pn (t) ≥ pLB
n (t) only
at a single point t suffices thanks to their strictly decreasing
and convex-shaped behaviors (see the proof of Lemma 2). In

this respect, consider a sufficiently large constant M > 0 such
b
that pUB
pn (M ). In this case,
n (M ) in Lemma 2 approximates
α
2
c
α ≥ c 2 −1 /n, which can
pn (M ) ≥ pLB
(M
)
becomes
M
n
n
always be satisfied by choosing a sufficiently large M .
Similarly, we prove the upper bound by showing
UBb
c
pUB
n (M ) ≥ pn (M ) ≈ p2n (M ). This is identical to proving

α
≥ 1, which is a stricly increasD(n, α) := cn α−2
nα
ing function of n ≥ 1 and a quasi-concave function of
α > 2 by the definition of cn . Therefore, it suffices to
2
verify D(1, α) = 1 − α2 α c1 ≥ 1 as α goes to either ∞
α→∞
or 2 as follows. Firstly, recall c1 = 1 discussed after
Lemma 1, thereby yielding limα→∞ D(1, α) = 1. Secondly,


 (b)
(a)
limα→2 D(1, α) = limα→2 1 − α2 Γ 1 − α2 Γ 1 + α2 =
limα→2 Γ 2 − α2 Γ 1 + α2 = 1, where (a) and (b) follow
2

from limα→2 1 − α2 α / 1 − α2 = 1 and Γ(x+1) = xΓ(x).
The results indicate that the minimum of D(n, α) is unity, and
consequently D(n, α) ≥ 1 for all α > 2 and n ≥ 1.

UBc
c
have the full range [0, 1],
As Fig. 2 illustrates, pLB
n and pn
LBb
and are even tighter than pn within the range [0.9, 1].
Along the same lines, we additionally consider different
choices of {A, B} leading to the following lower bounds.

Corollary 1. At a static typical UE, pn (t)’s lower bounds are:
2n
− α−2

×
pLB
n (t) = (1 + t)
+
pLB
n (t)

= 1 + cn

α
2

(9)

− 2
t α.

(10)
α

Proof: Applying {A, B} = {1, nα/(α − 2)} and {cn 2 , 1} to
LB
LB
the proof of Theorem 1 results in pn × (t) and pn + (t).
LB

The first bound pn × (t) is tight under high reliabilty, i.e. low
LB
LB
t as Fig. 2 shows. In addition, it satisfies pn × (t) = p1 × (t)n .
This equality indicates that static (LHS) and i.i.d. (RHS) UE
locations during n receptions identically affect SIR almost
LB
surely. Knowing that pn × (t) ≤ pn (t), we conclude spatial
UE correlation during n receptions increases successive SIR
coverage, compared to the coverage under i.i.d. mobility. Such
a spatial correlation gain diminishes under high reliability.
LB
The second bound pn + (t) accuracy increases with t as
Fig. 2 shows, which may complement Theorem 1 for less
reliable applications. This is achieved by utilizing the opposite
b
asymptotic curve pUB
n (t) in Lemma 2 for high t accuracy,
b
instead of pLB
(t)
exploited
in Theorem 1.
n
IV. M AXIMUM AVERAGE R ATE B OUNDS UNDER
n-S UCCESSIVE SIR C OVERAGE C ONSTRAINTS
In this section we derive the closed-form Rn (η), a typical
UE’s maximum average rate ensuring pn (t) ≥ η. For brevity,
rate is measured in nats/sec (1 nat = 1/ log 2 bits).
Recall the optimization (2). We define t̂ as the t value when
the reliability constraint’s equality holds, i.e. pn (t̂) = η. We
consider t̂ is the optimal SIR
 threshold maximizing average
rate, i.e. Rn (η) = pn t̂ log(1 + t̂). To show this, we
additionally define t̃ as the t maximizing Rn (η) without
the reliable constraint, i.e. t̃ := arg maxt≥0 pn (t) log(1 + t).

0.6

0.4

0.1
0.7

0.09
0.6
0.08

n=1
+

+

0.5

LBc , UBc Proposition 1
LB⇥
Corollary 2
LB+
LBb
Exact

Maximum Average Rate,R
Rnn(⌘)
( ) (nats/sec)

n=1
+

Maximum Average Rate,R
Rnn(⌘)
( ) (nats/sec)

0.7

+

0.5
0.07
0.06
0.4

LBc , UBc Proposition 1
LB⇥
Corollary 2
LB+
LBb
Exact

0.05

0.3

Baseline: n = 3, ↵ = 4
n = 3, ↵ = 4

0.3
0.04

0.2
0.1
0
0.3

↵ = 2.5
0.4

0.5

0.6

0.7

0.8

0.9

1

Reliability,

(a) for the reliability η within [0.3, 1].

0.03
0.2

Baseline: n = 3, ↵ = 4
n = 3, ↵ = 4

0.02
0.1
0.01

↵ = 2.5

00
0.3
0.9

0.91 0.4

0.92

0.5
0.93

0.6
0.94

0.95

0.70.96

0.8
0.97

0.98 0.9 0.99

11

Reliability,

(b) for the reliability η within [0.9, 1].

Fig. 3. The bounds of maximum average rate Rn (η) ensuring n-successive coveraege with probability η, under different n and the path loss exponent α > 2.

If t̂ ≤ t̃, then the quasi-concave objective function strictly
increases, having the maximum at t = t̂. This is true for
sufficiently high reliability, e.g. when η ≥ 0.29 for all n ≥ 1
and 2 < α ≤ 4, which is much less than our target η ≥ 0.9.
For this reason, Rn (η) derivation is identical to deriving t̂.
For URLLC, it can be analytically derived via Theorem 1. The
original Rn (η) is then bounded by their corresponding results
UBi
i
RLB
n (η) ≤ Rn (η) ≤ Rn (η) given as follows.
Proposition 1. (Rate Bounds) At a static typical UE, maximum
average rate Rn (η) satisfying pn (t) ≥ η is bounded by:
!
α
η −( 2 −1) − 1
LBc
Rn (η) = η log 1 +
(11)
n




α
2 η− 2 − 1
UBc
Rn (η) = η log 1 + 1 −
.
(12)
α
n
Proof: The optimal SIR threshold t̂ is derived by taking the
LBc
c
and pUB
inverse pn (t) into η, i.e. t̂ = p−1
n
n (η). Applying pn
in Theorem 1 completes the proof.

Exploiting Corollary 1, we also derive the following bounds.
Corollary 2. At a static typical UE, Rn (η)’s lower bounds
are:
 1 α

− n ( 2 −1)
×
(η)
=
η
log
η
RLB
(13)
n




α
2
−α
+
RLB
η− 2 − 1 .
(14)
n (η) = η log 1 + cn
As reliability increases, Fig. 3 shows that the exact Rn (η)
curve from (1) asymptotically converges to RnLBc (η), RnUBc (η),
UB
and Rn × (η), as well as RnLBb (η) = α−2
2n (1 − η) similarly
b
derived by applying pLB
(t)
in
Lemma
2.
Especially, the conn
UB×
vergence to Rn (η) implies the spatial UE correlation impact
becomes negligible in URLLC as stated after Corollary 1.
LB
For moderate reliability 0.3 < η < 0.9, Rn + (η) becomes
b
accurate. Note that RLB
n (η) can also be an approximation if
n is large; otherwise, a huge gap exists as shown by the curve
in Fig. 3 when n = 1 (yellow circle).
Additionally, it is noted that the proposed closed-form
bounds are also applicable for conventional average rate
calculation under the full CSI assumption as in [3]–[5]. For

when n = 1, average rate
+ SIR)] =
Rexample,
R ∞ is E[log(1
∞
x
Pr
(log(1
+
SIR)
>
x)
dx
=
p
(e
−
1)
dx.
Apply1
0
0
UBa
c
ing the simplest bounds pLB
and
p
to
this
yields
the
1
1
average rate bounds (α − 2)/2 and α/2 respectively. Their
linear combination with the equal weight factor leads to the
rate approximation (α − 1)/2. For α = 4, this corresponds
to 1.5 nats/sec. The result is fit well with the exact calcu simplifying the√known expression

Rlation 1.49 nats/sec
√ x while
dx/
1
+
e
−
1
π/2 − arctan(1/ ex − 1) in [3].
x>0
V. C ONCLUSION
We proposed closed-form n-successive SIR coverage
bounds, and thereby derived the closed-form maximum average rate guaranteeing such a successive coverage requirement.
R EFERENCES
[1] Ericsson Business Review, Manufacturing Reengineered: Robots, 5G
and the Industrial IoT, no. 4, 2015.
[2] 5G-PPP, 5G and the Factories of the Future, White Paper, Oct. 2015.
[3] J. G. Andrews, F. Baccelli, and R. K. Ganti, “A Tractable Approach
to Coverage and Rate in Cellular Networks,” IEEE Trans. Commun.,
vol. 59, pp. 3122–3134, Nov. 2011.
[4] X. Zhang and M. Haenggi, “A Stochastic Geometry Analysis of InterCell Interference Coordination and Intra-Cell Diversity,” IEEE Trans.
Wireless Commun., vol. 13, pp. 6655–6669, Dec. 2014.
[5] J. Park, S.-L. Kim, and J. Zander, “Tractable Resource Management with
Uplink Decoupled Millimeter-Wave Overlay in Ultra-Dense Cellular
Networks,” IEEE Trans. Wireless Commun., vol. 15, pp. 4362–4379,
Jun. 2016.
[6] M. Haenggi, Stochastic Geometry for Wireless Networks. Cambridge
Univ. Press, 2013.
[7] J. Park, S. Jung, S.-L. Kim, M. Bennis, and M. Debbah, “User-Centric
Mobility Management in Ultra-Dense Cellular Networks under SpatioTemporal Dynamics,” Proc. IEEE GLOBECOM, Washington, DC, USA,
Dec. 2016.
[8] S. Krishnan and H. S. Dhillon, “Spatio-Temporal Interference Correlation and Joint Coverage in Cellular Networks,” [Online]. Available:
https://arxiv.org/abs/1606.05332.
[9] B. Blaszczyszyn, M. K. Karray, and H. P. Keeler, “Using Poisson
Processes to Model Lattice Cellular Networks,” Proc. IEEE INFOCOM,
Turin, Italy, Apr. 2013.
[10] H. S. Dhillon, R. K. Ganti, F. Baccelli, and J. G. Andrews, “Modeling
and Analysis of K-Tier Downlink Heterogeneous Cellular Networks,”
IEEE J. Sel. Areas Commun., vol. 30, pp. 550–560, Apr. 2012.
[11] J. Park, D. Kim, P. Popovski, and S.-L. Kim, “Revisiting Frequency
Reuse towards Supporting Ultra-Reliable Ubiquitous-Rate Communication,” Proc. IEEE WiOpt Wksp. SpaSWiN, Paris, France, May 2017.
[12] G. E. Andrews, R. Askey, and R. Roy, Special Functions. Cambridge
Univ. Press, 1999.

