(* The HolonomicFunctions package can be downloaded from 
  http://www.risc.jku.at/research/combinat/software/HolonomicFunctions/ *)
<<RISC`HolonomicFunctions`;

HolonomicFunctions Package version 1.7.1 (09-Oct-2013)
written by Christoph Koutschan
Copyright 2007-2013, Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

(* Define the functions used below. *)
(* Determine the holonomic rank. *)
HolonomicRank[ops_]:=
With[{u=UnderTheStaircase[ops]},If[u===Infinity,Throw["Not holonomic"],Length[u]]];
(* Deduce coefficient list of an OrePolynomial w.r.t. the given support. *)
mycfl[op_,uts_]:=Total[Cases[op[[1]],{a_,LeadingExponent[#]}:>a]]&/@uts;
(* Convert a holonomic system into a Pfaffian system (linear system of differential equations). *)
PfaffianSystem[ann_,op_]:=
With[{uts=UnderTheStaircase[ann]},mycfl[OreReduce[op**#,Together[ann]],uts]&/@uts];
(* Support for handling expressions in dB. *)
dB[x_]:=10^(x/10);

SNR moment generating function (m.g.f.)

(* Write M_{n_1}(s) as defined in Eq.(55). Upon changing index n_1 into m,  M_m corresponds to H_m in Eq.(67). *)
Mm=Sum[Binomial[(Subscript[n, 1]),(Subscript[m, 1])]*(-1)^(Subscript[m, 1])/(1-s*(Subscript[Γ, 1]))^(N+(Subscript[n, 1])-(Subscript[m, 1])),{(Subscript[m, 1]),0,(Subscript[n, 1])}]/.(Subscript[n, 1])->m

(* Write the factor that appears in the definition of G_n(z) in Eq.(67), except H_m. Note that this factor is used several times below. *)
TraditionalForm[facGn=Binomial[n,m]*Pochhammer[N,m]/Pochhammer[(Subscript[N, R])+n-m,m]*(Subscript[c, 1])^m]

(* Deduce the annihilator of G_n(z) from Eq.(67), written for the SNR m.g.f. *)
(* An annihilator is a set of partial differential equations (PDEs)/recurrences, in operator notation. *)
(* Note how creative telescoping is employed to deal with this summation. *)
annGn1=CreativeTelescoping[facGn*Mm,S[m]-1,{S[n],Der[s],Der[z]}][[1]]

(* Derive the annihilator of the SNR m.g.f., by multiplying the remaining terms from Eq.(67), and summing w.r.t. n. Again creative telescoping is employed to deal with the summation. *)
annMGF=CreativeTelescoping[DFiniteTimes[annGn1,Annihilator[Exp[-z]*z^n/n!,{S[n],Der[s],Der[z]}]],S[n]-1][[1]]

(* Deduce the size of the Pfaffian systems for the m.g.f. *)
HolonomicRank[annMGF]

(* Deduce the vector of derivatives of m(s; z), which is shown in Section V.C., Step (1). *)
UnderTheStaircase[annMGF]

(* From the m.g.f. annihilator, display the matrices required in Eq.(72), as lists of elements. *)
PSMGFs = PfaffianSystem[annMGF,Der[s]]
PSMGFz = PfaffianSystem[annMGF,Der[z]]
(* From the m.g.f. annihilator, display the matrices required in Eq.(72). *)
(Subscript[Θ, s])==MatrixForm[PSMGFs]
(Subscript[Θ, z])==MatrixForm[PSMGFz]

SNR probability density function (p.d.f.)

(* Write p_{n_1}(t) as defined in Eq.(58). Then, change index n_1 into m, so that p_{m}(t) corresponds to H_m from Eq.(67). The fact that Mathematica outputs p_m(t) as a closed form in terms of HypergeometricU for us can be disregarded. *)
pm=Sum[Binomial[(Subscript[n, 1]),(Subscript[m, 1])]*(-1)^(Subscript[m, 1])*t^(N+(Subscript[n, 1])-(Subscript[m, 1])-1)*Exp[-t/(Subscript[Γ, 1])]/(N+(Subscript[n, 1])-(Subscript[m, 1])-1)!/(Subscript[Γ, 1])^(N+(Subscript[n, 1])-(Subscript[m, 1])),{(Subscript[m, 1]),0,(Subscript[n, 1])}]/.(Subscript[n, 1])->m

(* Deduce the annihilator of G_n(z) from Eq.(67), written for the SNR p.d.f. *)
annGn2=CreativeTelescoping[facGn*pm,S[m]-1,{S[n],Der[t],Der[z]}][[1]]

(* The annihilator for the SNR p.d.f., using Eq.(67) as definition. *)
annPDF=CreativeTelescoping[DFiniteTimes[annGn2,Annihilator[Exp[-z]*z^n/n!,{S[n],Der[t],Der[z]}]],S[n]-1][[1]]

(* Alternatively, the annihilator for the SNR p.d.f., obtained from the m.g.f. by Laplace transform. *)
(* I.e., multiply the expression with Exp[-s*t] and integrate w.r.t. s. *)
ops={Der[s],Der[t],Der[z]};
annPDF1=CreativeTelescoping[DFiniteTimes[ToOrePolynomial[Append[annMGF,Der[t]],OreAlgebra@@ops],Annihilator[Exp[-s*t],ops]],Der[s]][[1]]

(* Compare the two results obtained for the p.d.f., they agree. *)
GBEqual[annPDF,annPDF1]

HolonomicRank[annPDF]

(* The vector of derivatives of p(t;z) needed for HGM, see Section V.C., Step (2). *)
UnderTheStaircase[annPDF]

(* Display the matrices required in the Pfaffian systems from Eq.(73), as lists of elements. *)
PSPDFt = PfaffianSystem[annPDF,Der[t]]
PSPDFz = PfaffianSystem[annPDF,Der[z]]
(* Display the matrices in matrix form. *)
(Subscript[Ξ, t])==MatrixForm[PSPDFt]
(Subscript[Ξ, z])==MatrixForm[PSPDFz]

Outage Probability

(* Note that our Eq.(63) involves the incomplete gamma function γ(k,x) defined as in Footnote 7 on page 7. On the other hand, Mathematica's Gamma[k,x] function is different as revealed by the following: *)
Integrate[t^(k-1)*Exp[-t],{t,0,x},Assumptions->k>0]

(* Write the summand in P_{o,n_1} defined in Eq.(60), by using Eq.(64) and by accounting for the Mathematica definition of Gamma[k,x]. *)
TraditionalForm[Pon1Smnd=Binomial[(Subscript[n, 1]),(Subscript[m, 1])]*(-1)^(Subscript[m, 1])*(1-Gamma[N+(Subscript[n, 1])-(Subscript[m, 1]),τ/(Subscript[Γ, 1])]/Gamma[N+(Subscript[n, 1])-(Subscript[m, 1])])]

(* Derive the annihilator of P_{o,n_1} w.r.t. n_1. Note that the result is a difference operator, which corresponds to a difference equation w.r.t. n_1 satisfied by P_{o,n_1}. *)
Annihilator[Sum[Pon1Smnd,{(Subscript[m, 1]),0,(Subscript[n, 1])}],S[(Subscript[n, 1])]]

(* Note that the above command internally uses creative telescoping to derive the annihilator. *)
(* Next, we show an alternative approach that explicitly employs creative telescoping; Finally, upon replacing n_1 with m, for Eq.(67), we obtain the annihilator of H_m (corresponding to P_{o,m}) from Eq.(67). Note that this is again a difference operator. *)
annPom=CreativeTelescoping[Pon1Smnd,S[(Subscript[m, 1])]-1,S[(Subscript[n, 1])]][[1]]/.(Subscript[n, 1])->m

(* Define the operators used below. *)
ops={S[m],S[n],Der[z]};
(* First, derive the annihilator for the summand of G_n in Eq.(67), written for the outage probability. *)
annGnSmnd=DFiniteTimes[ToOrePolynomial[Join[annPom,{S[n]-1,Der[z]}],OreAlgebra@@ops],Annihilator[facGn,ops]];
(* Then, do summation w.r.t. m, using a special command with a special option. *)
annGn3=FindCreativeTelescoping[annGnSmnd,S[m]-1,Denominator->(n+(Subscript[N, R]))*Pochhammer[n-m,5]][[1]]

(* Derive the annihilator for the outage probability, using Eq.(67) as definition. *)
annOP=CreativeTelescoping[DFiniteTimes[annGn3,Annihilator[Exp[-z]*z^n/n!,{S[n],Der[z]}]],S[n]-1][[1]]

(* Alternatively, compute the annihilator of the outage probability directly, *)
(* by integrating the SNR p.d.f. from 0 to τ, as in Eq.(59). *)
(* Because of the non-natural upper bound we have to consider inhomogeneous parts. *)
ct=CreativeTelescoping[annPDF,Der[t]];
annOP1=OreGroebnerBasis[Flatten[MapThread[Function[{p,q},(#**p)&/@DFiniteSubstitute[DFiniteOreAction[annPDF,q],{t->τ},Algebra->OreAlgebra[Der[z]]]],ct]]]

(* Test that the two different methods to compute this annihilator yield the same result. *)
GBEqual[annOP,annOP1]

HolonomicRank[annOP]

(* Reveal that annihilator of P_o(z) consists of a single ordinary differential equation of order 5 w.r.t. z. *)
Support[annOP]

(* The vector of derivatives needed in HGM. *)
UnderTheStaircase[annOP]

(* Display the companion matrix Φ of the corresponding Pfaffian system from Eq.(74) as list of elements. *)
PSOP = PfaffianSystem[annOP1,Der[z]]
(* Display the companion matrix Φ of the corresponding Pfaffian system from Eq.(74). Unfortunately, the matrix is fully viewable only under Mathematica (.nb), or under the Wolfram CDF player (.cdf). *)
MatrixForm[PSOP]

(* The following allows us to easily copy-paste the results into Matlab. *)
(*InputForm[PSOP]*)

Ergodic Capacity

(* Derive the annihilator for the ergodic capacity C(z), using the definition in Eq.(61). *)
annC=FindCreativeTelescoping[DFiniteTimes[Annihilator[Log[2,t+1],{Der[t],Der[z]}],annPDF],Der[t]][[1]]

(* Reveal that C(z) satisfies an ordinary differential equation of order 7. *)
Support[annC]

(* The vector of derivatives needed for HGM. *)
UnderTheStaircase[annC]

(* Display the companion matrix Ψ of the corresponding Pfaffian system, Eq.(74), as a list of elements. *)
PSC = PfaffianSystem[annC,Der[z]]
(* Display the companion matrix Ψ of the corresponding Pfaffian system, Eq.(74). Unfortunately, the matrix is fully viewable only under Mathematica (.nb), or under the Wolfram CDF player (.cdf). *)
MatrixForm[PSC]

(* The following allows us to easily copy-paste the results into Matlab. *)
(*InputForm[PSC]*)

HGM computation

(* Evaluate the series expression for the outage probability. *)
OutProb[Γb_,Γ1_,K_,der_,NT_,NR_,c1_,τ_,trunc_,prec_]:=
With[{
NN=NR-NT+1,
(* If necessary, differentiate w.r.t. z, then plug in its value. *)
zpart=D[Exp[-z]*z^n,{z,der}]/.z->K*NR*(NT-1)},
N[
(* The part from G_n(z) in Eq.(67). *)Sum[Sum[Binomial[n,m]*Pochhammer[NN,m]/Pochhammer[NR+n-m,m]*c1^m*
(* The part P_{o,m} from Eq.(60). *)
Sum[Binomial[m,m1]*(-1)^m1*(1-Gamma[NN+m-m1,τ/Γ1]/Gamma[NN+m-m1]),{m1,0,m}],
(* The remaining parts from Eq.(67). *)
{m,0,n}]*zpart/n!,{n,0,trunc}],prec]];

(* Here we produce the results shown in Fig. 1 using HGM. *)
RTKinv11=1.035478792935417;
c1=0.289132807530209;
With[{NT=4,NR=6,tau=dB[82/10],K0=dB[-7]},
Timing[data=Table[
z0=K0*NR*(NT-1);
z1=dB[K]*NR*(NT-1);
Gamma1=2*dB[gb]/(dB[K]+1)/RTKinv11;
inits=Table[Derivative[d][Po][z0]==OutProb[dB[gb],Gamma1,K0,d,NT,NR,c1,tau,20,20],{d,0,4}];
ode=ApplyOreOperator[First[annOP],Po[z]]/.{(Subscript[N, R])->NR,N->NR-NT+1,(Subscript[c, 1])->c1,τ->tau,(Subscript[Γ, 1])->Gamma1};
{gb,Po[z1]}/.First[NDSolve[Join[{ode==0},inits],Po,{z,z0,z1},AccuracyGoal->20]]
,{K,0,21,7},{gb,15,25,5/2}];]]

(* The plot shown in Fig. 1. *)
ListLogPlot[data,Joined->True,Ticks->{{15,17.5,20,22.5,25},Automatic}]

(* Here are the numeric values. *)
TableForm[Map[Last,data,{2}]]

(* Compute the values of P_o for K = 0dB using the truncated series. *)
Timing[test1=Table[OutProb[dB[gb],2*dB[gb]/(dB[0]+1)/RTKinv11,dB[0],0,4,6,c1,dB[82/10],50,10],{gb,15,25,5/2}];]

(* Compare with HGM. *)
test1/(Last/@data[[1]])

(* Compute the values of P_o for K = 7dB using the truncated series. *)
Timing[test2=Table[OutProb[dB[gb],2*dB[gb]/(dB[7]+1)/RTKinv11,dB[7],0,4,6,c1,dB[82/10],140,10],{gb,15,25,5/2}];]

(* Compare with HGM. *)
test2/(Last/@data[[2]])

(* Compute the values of P_o for K = 14 dB using the truncated series. *)
(* Note that the duration of this computation is substantial. *)
Timing[test3=Table[OutProb[dB[gb],2*dB[gb]/(dB[14]+1)/RTKinv11,dB[14],0,4,6,c1,dB[82/10],560,10],{gb,15,25,5/2}];]

(* Compare with HGM. *)
test3/(Last/@data[[3]])