%appswz.tex:
%%a Plain TeX file by Moa Apagodu and Doron Zeilberger (8 pages)

%begin macros
\def\L{{\cal L}}
\baselineskip=14pt
\parskip=10pt
\def\halmos{\hbox{\vrule height0.15cm width0.01cm\vbox{\hrule height
  0.01cm width0.2cm \vskip0.15cm \hrule height 0.01cm width0.2cm}\vrule
  height0.15cm width 0.01cm}}
\font\eightrm=cmr8 \font\sixrm=cmr6
\font\eighttt=cmtt8
\magnification=\magstephalf
\def\P{{\cal P}}
\def\Tilde{\char126\relax}
\def\Q{{\cal Q}}
\def\1{{\overline{1}}}
\def\2{{\overline{2}}}
\parindent=0pt
\overfullrule=0in
\def\frac#1#2{{#1 \over #2}}
%\headline={\rm  \ifodd\pageno  \RightHead  \else  \LeftHead  \fi}
%\def\RightHead{\centerline{
%Title
%}}
%\def\LeftHead{ \centerline{Doron Zeilberger}}
%end macros

\centerline
{
\bf
FIVE Applications of Wilf-Zeilberger Theory to
Enumeration and Probability
\footnote{$^1$}
{\eightrm  \raggedright
Oct. 18, 2006.
Accompanied by Maple packages {\eighttt AppsWZ} and
{\eighttt AppsWZmulti} downloadable from Zeilberger's website.
Sample input and output can be viewed in: \hfill\break
{\eighttt http://www.math.rutgers.edu/\Tilde zeilberg/mamarim/mamarimhtml/
appswz.html} .
The research of the second author was supported in part by the NSF.
}
}
\bigskip
\centerline{Moa APAGODU and Doron ZEILBERGER}

\bigskip
\qquad\qquad\qquad\qquad\qquad\qquad\qquad
{\it Dedicated to Sergei Abramov on becoming {\bf FIVE}-dozen
years old}


{\bf Explicit Formulas vs. Algorithms}

In the old days, when one had to find some sequence, $a(n)$,
there were two extremes. In the lucky case, one had an
{\bf explicit formula}. For example, the probability
of tossing a fair coin $2n$ times and getting exactly $n$ Heads,
equals $(2n)!/(2^{2n} n!^2)$.
Sometimes, cheatingly,
one considered as `explicit' expressions in terms
of sums (or multisums) or integrals (or multi-integrals).
The other extreme was to just have a {\it numerical} algorithm,
that for each (numeric!) input $n$, found the output.
In that case the algorithm was rated by its efficiency.

Another compromise was an {\bf asymptotic} formula, valid 
(approximately!) for large $n$.

But what's a formula?, it is a kind of algorithm. Of course,
it is more than that, theoretically, but from a practical
point of view it should be judged by the efficiency of the
implied algorithm.

{\bf The Holonomic Ansatz}

Let's look at the explicit formulas that are called
`closed-form', or more precisely {\it hypergeometric} sequences.
A sequence $a(n)$ is called hypergeometric if
the ratio $a(n+1)/a(n)$ is a {\it rational} function of $n$,
i.e. a quotient $P(n)/Q(n)$ where $P(n)$ and $Q(n)$ are polynomials.
For example for the above-mentioned probability of getting exactly
$n$ Heads when tossing a fair coin $2n$ times,
$p(n):=(2n)!/(2^{2n} n!^2)$, we have
$p(n+1)/p(n)=(2n+1)/(2(n+1))$, or, by cross-multiplying
$$
2(n+1)p(n+1)-(2n+1)p(n)=0 \quad.
$$
This is an example of a {\it first}-order {\bf linear recurrence
equation with polynomial coefficients}.
Once you have the trivial value $p(0)=1$ you can use it
to compile a table of $p(n)$ for  $n<N$, for any desired
$N$ in $O(N)$ operations.

The same is true for solutions of any linear recurrence
equation with polynomial coefficients, 
$$
\sum_{i=0}^{L} a_i(n) p(n+i) =0 \quad ,
$$
of order $L$. The only difference is that we need $L$
{\it initial conditions}, $p(0),p(1), \dots, p(L-1)$.
We also assume that $a_L(n)=0$ has no positive integer
roots.

Such sequences were dubbed {\it P-recursive} by Richard
Stanley in his seminal paper [St], but we prefer the
name {\it holonomic}. 
Zeilberger[Z] famously showed that many
sequences that arise in combinatorics, probability,
and elsewhere are holonomic, and this was made
into a full-fledged algorithmic theory by
Wilf and Zeilberger [WilZ]. 

In our humble opinion, a holonomic representation
of a sequence is to be considered {\it explicit}, since
it is almost as good as a {\it closed-form} (i.e. hypergeometric).
In the previous literature on WZ theory,
there were few scattered examples of potential applications,
but the focus was on the theory and the algorithms, not 
on specific applications. 

{\bf Why this Paper?}

The purpose of this paper is to fill this gap.
We only list five such applications, but the reader
can doubtless find many others. It is hoped that
our implementation of these five applications will
aid the reader to implement other ones that he or she might
be interested in.

{\bf The Maple packages AppsWZ and AppsWZmulti}

This article is accompanied by two Maple packages.
{\tt AppsWZ} that does applications of the single-variable
case ([AlZ]), and {\tt AppsWZmulti} that does applications of the 
multi-variable case ([WilZ],[ApZ]).

{\bf Asymptotics}

Another nice feature of being a solution of a linear recurrence
equation with polynomial coefficients is that using
the Birkhoff-Trijinski method (see [WimZ] for a lucid exposition),
one can deduce the {\it asymptotics} to any desired order.
This algorithm has been implemented by us in Maple, and is part of
both packages.


{\bf First Application: Rolling a Die}

If instead of tossing a coin $n$ times, you roll
a $k$-faced die, marked with positive or negative numbers,
and you win the amount shown on the landed face
(or lose, if it is a negative number). What is the probability
that after $n$ rolls, you break even? More generally,
how likely are you to win exactly $d$ dollars? 

If the $i^{th}$ face ($i=1 \dots k$) shows the amount
$m_i$, and lands with probability $p_i$, let the
{\it probability generating function} of the die be defined by 
$$
P(x):=\sum_{i=1}^{k} p_i x^{m_i} \quad.
$$
It is {\it very} well known (and {\it very very} easy to see), that
the probability, let's call it $a_d(n)$, of winding up with
$d$ dollars after $n$ rolls is the constant term
of
$$
{{P(x)^n} \over {x^{d}}} \quad .
$$
But that's exactly grist for the Almkvist-Zeilberger
mill! In the package {\bf AppsWZ} this is 
accompanied by the commands {\tt RecProbVisit} and  
{\tt RecProbVisitE}.

For example, let $a(n)$ be the probability of rolling
a fair (standard) die $2n$ times and having the total score
being {\it exactly} the expected value $7n$. Then
$a(n)$ satisfies the following third-order linear recurrence
equation with polynomial coefficients
$$
-4(5+2n)(2n+3)(2n+1)(7n+19)(5n+11)(7n+20)(7n+13)(n+2)(n+1)
a(n)
$$
$$
+4(7n+20)(5+2n)(2n+3)(n+2)(25480n^5+223496n^4+755066n^3+1223233n^2+
946889n+279936)a(n+1)-
$$
$$
(5n+6)(5+2n)(6+7n)(499359n^6+6777015n^5+38079431n^4
$$
$$
+113390385n^3+188723986n^2+166469280n+60800544)a(n+2)+
$$
$$
30(5n+14)(5n+13)(
5n+12)(7n+12)(5n+11)(5n+6)(7n+13)(6+7n)(n+3)a(n+3)=0 \quad ,
$$
and the Birkhoff-Trijinksi method implies that the asymptotics
is:
$$
(.197833497170804)n^{-1/2}
(1-(111/1400)/n-(12037/5488000)/n^2+(1367631/1097600000)/n^3 + \dots) \quad .
$$
Readers can produce their own output for the scenario of their choice.

{\bf Second Application: How many ways to have r people 
chip in to pay a bill of n cents}

In George P\'olya's classic `picture writing' paper[Po], he considers
the problem of figuring out how many ways can {\it one} person
pay a bill of $n$ cents using any number of coins.
If the denominations are $\{d_1, \dots , d_k\}$
($\{1,5,10,25,50,100\}$ in the US), then the required number is
the coefficient of $x^n$ in the {\it generating function}
$$
\prod_{i=1}^{k} {{1} \over {1-x^{d_i}} } \quad .
$$
Calling this number $a(n)$, this entails, trivially,
a {\it linear recurrence equation} with {\bf constant} coefficients
of {\it order} $\sum_{i=1}^{k} d_i$.
But by allowing {\bf polynomial} coefficients, one can get,
thanks to Almkvist-Zeilberger, a recurrence of order
$ \leq lcm (d_1, \dots, d_k)$.

More generally, of $a_r(n)$ is the number of ways of breaking
$n$ cents with (up to) $r$ people participating (or
equivalently, one person with $r$ pockets in his or her pants or
dress), the generating function is
$$
\prod_{i=1}^{k} {{1} \over {(1-x^{d_i})^r} } \quad ,
$$
and applying Almkvist-Zeilberger to
$$
{{1} \over {x^{n+1}}} \prod_{i=1}^{k} {{1} \over {(1-x^{d_i})^r} } \quad ,
$$
produces a linear recurrence, still with polynomial coefficients
(but now these are polynomials in both $n$ and $r$) of the
above order. On the other hand if you stick to {\it constant} coefficients
then the order would be $r(\sum_{i=1}^{k} d_i)$, and would only make
sense for (small!) {\it numeric} $r$, while within the
{\it holonomic ansatz}, one can have {\it symbolic} $r$ without any
increase in the order.
The procedure that takes care of this problem
in the Maple package {\tt AppsWZ} is
{\tt TeamEffortMoneyChanging}. For example,
entering:

{\tt TeamEffortMoneyChanging($\{1$,5,10,25$\}$,n,N,r);}

would tell you that if $a_r(n)$ is the number of ways of breaking
$n$ cents with (up to) $r$ people chipping-in, using
only pennies, nickels, dimes, and quarters,
and denoting by $N$ the shift operator in $n$:
($Nx(n):=x(n+1)$), it turns out that $a_r(n)$ is annihilated
by the following $30^{th}$-order
linear recurrence operator with polynomial coefficients:
$$
(-41r-n)+(-r)N+(-r)N^2+(-r)N^3+(-r)N^4+
$$
$$
(-47r-n-5)N^5+(-2r)N^6+(-2r)N^7+(-2r)N^8+(-2r)N^9+(-22r)N^{10}
$$
$$
+(-2r)N^{11}+(-2r)N^{12}+(-2r)N^{13}+(-2r)N^{14}
$$
$$
+(-22r)N^{15}+(-2r)N^{16}+
(-2r)N^{17}+(-2r)N^{18}+(-2r)N^{19}+(-22r)N^{20}+
$$
$$
(-2r)N^{21}+(-2r)N^{22}+(-2r)N^{23}+(-2r)N^{24}
+(-6r+n+25)N^{25}
$$
$$
+(-r)N^{26}+(-r)N^{27}+
(-r)N^{28}+(-r)N^{29}+(n+30)N^{30} \quad .
$$

{\bf Third Application: Hidden Markov Models}

In contemporary bioinformatics (see for example
the fascinating and lucid expository article by Lior
Pachter and Bernd Sturmfels in the wonderful book
that they edited([PaStu]), one has to {\it estimate} probabilities
from output. The usual approach is {\it maximum likelihood},
that entails solving {\it huge} systems of polynomial equations
that are handled via Buchberger's Gr\"obner Basis Algorithm.
But, following Laplace, we can also use a {\it Bayesian} approach
and get quotients of integrals. 
These integrals can be handled via the multi-Almkvist-Zeilberger
([ApZ]) algorithm. For any {\it specific} instance
one may want to use {\it numerical integration}. If one wants
to compile a table anticipating all possible outputs (of course
up to a certain length), then the recurrences supplied
by Almkvist-Zeilberger are much more efficient.

Let's illustrate it with the simple case of (only) two {\it dice},
that enables us to stay within the single-variable
Almkvist-Zeilberger ([AlZ]).

In the proverbial {\it dishonest casino}, suppose that
there are two kinds of (identically looking) dice, with the
same possible outcomes (faces), let's call them
$\{1, \dots, m \}$, each with a {\it known} probability
distribution $(p_1, \dots , p_m)$ and
$(q_1, \dots , q_m)$. What is unknown is the probability,
$x$, of using the first die (and hence $1-x$ for the other die).
Judging from the sequence of outcomes (or rather by the relative
frequencies of the landed faces),
one has to {\it estimate} $x$.
If the (hidden) probability of using the first die was
$x$, (whatever it is), then the probability of
the output distribution being $a_1, \dots, a_m$ is
$(a_1 + \dots + a_m)!/(a_1! \cdots a_m!)$ times
$$
L(x):=\prod_{i=1}^{m} (xp_i+(1-x)q_i)^{a_i} \quad.
$$
The {\it maximum-likelihood} estimate is to {\it maximize}
$L(x)$ by solving $L'(x)=0$ (with more dice one gets
partial derivatives and systems of equations in several variables).
But, following the more democratic approach of Laplace,
that considers {\it all scenarios} (and that 
famously tells you that if so far you succeeded $m$ times and failed
$n$ times, then your estimated probability of success in your next try
is not $m/(m+n)$ but rather $(m+1)/(m+n+2)$),
we would have not the root of $L'(x)=0$ but rather
$$
{{\int_{0}^{1} x\prod_{i=1}^{m} (xp_i+(1-x)q_i)^{a_i} \, dx}
\over
{\int_{0}^{1} \prod_{i=1}^{m} (xp_i+(1-x)q_i)^{a_i} \, dx}
} \quad.
$$
Calling the top quantity $T(a_1, \dots, a_m)$ and the
bottom quantity $B(a_1, \dots, a_m)$, we see that
Almkvist-Zeilberger can manufacture {\it pure} linear recurrences
(with polynomial coefficients) in {\it each} of the 
discrete variables $a_1, \dots , a_m$ (in Nature
$m=4$, or $m=20$). The multi-dice analog of this
(which would employ the multi-variable Almkvist-Zeilberger
, done in [ApZ]) may be of interest to bioinformaticians.
So far we have only implemented the interface in the
two-dice case.

The procedure that takes care of this problem in
{\tt AppsWZ} is {\tt ExpProbFirstDie(Lp1,Lp2,a0)},
where {\tt Lp1}, and {\tt LP2} are the probability distributions
and {\tt a0} is the list of respective outcome. The
novelty here is that both top and bottom integrals
are not computed directly but via the recurrences outputted
by Almkvist-Zeilberger. With the {\it option remember} it should
be much more efficient if one wants to pre-compute a table
of estimated probability for each $a0$ of size less than some
(lagre) pre-assigned value.

For example, if there are two coins, one fair
and one loaded with probability of a Head being 
$1/3$, and the outcome was $10$ and $10$ (which
ML would say that the fair coin was used all the time),

{\tt evalf(ExpProbFirstDie([1/2,1/2],[1/3,2/3],[10,10])):}

would estimate that the fair coin was used only
\%58.76 of the time.

{\bf Fourth Application: Lattice Paths Counting}

We all know that the number of ways of walking from
the origin $(0,0)$ to the point $(m,n)$ in the
square lattice, with unit northbound and unit eastbound steps
$\{ (0,1), (1,0) \}$ is
$$
F(m,n)={{(m+n)!} \over {m!n!}} \quad .
$$
This immediately implies (and is equivalent to) the fact
that
$$
{{F(m+1,n)} \over {F(m,n)}}=
{{m+n+1} \over {m+1}} \quad, \quad
{{F(m,n+1)} \over {F(m,n)}}=
{{m+n+1} \over {n+1}} \quad.
$$
and cross-multiplying yields
$$
(m+1)F(m+1,n)-(m+n+1)F(m,n)=0 \quad , \quad  
(n+1)F(m,n+1)-(m+n+1)F(m,n)=0 \quad .
$$
In other words the discrete function $F(m,n)$ satisfies {\it pure}
linear recurrences equations with {\it polynomial} coefficients,
that happen, in this simple case, to be {\it first-order}.
Recall that a recurrence is {\it pure} if 
only one of the variables changes at a time.
For example, $F(m,n)$ trivially satisfies the ``mixed''
({\it partial}) recurrence $F(m,n)=F(m-1,n)+F(m,n-1)$.

An amazing consequence of Wilf-Zeilberger theory ([WilZ]) is that
this is still true for an {\it arbitrary} set of (positive) steps,
and in {\it arbitrary} dimension. Of course the pure recurrences
are no longer (usually) first-order, but as above, this
is a minor computational disadvantage.

Indeed, if we are walking in the $d$-dimensional (hyper)cubic lattice,
starting at the origin, and with a set of steps $S$
(all with non-negative coordinates, excluding the step {\bf 0}
[staying in place]),
the generating function is trivially seen to be
$$
\sum_{m} F(m_1, \dots , m_d) 
{x_1}^{m_1} \cdots {x_d}^{m_d}
$$
$$
={{1} \over 
{\left ( 1- {{\sum} \atop  {(s_1, \dots, s_d) \in S} }
{x_1}^{s_1} \cdots {x_d}^{s_d} \right )
}} \quad .
$$
So it follows that our discrete function of interest,
$F(m_1, \dots , m_d)$ equals the formal residue of
$$
{
{x_1^{-m_1-1} \cdots x_d^{-m_d-1} } \over 
{ \left (1- {{\sum} \atop  {(s_1, \dots, s_d) \in S} }
{x_1}^{s_1} \cdots {x_d}^{s_d} \right )
}} \quad .
$$
If you are only interested in getting to points
on the diagonal, then $f(n):=F(n,n, \dots, n)$ is given
by the formal residue of
$$
{{x_1^{-n-1} \cdots x_d^{-n-1} } \over 
{ \left (1- {{\sum} \atop {(s_1, \dots, s_d) \in S} }
{x_1}^{s_1} \cdots {x_d}^{s_d} \right )
}} \quad ,
$$
and once again thanks to Wilf-Zeilberger theory, satisfies
a linear recurrence equation with polynomial coefficients.

This is implemented in the second Maple package
accompanying this paper, {\tt AppsWZmulti} that is
``powered'' by the Maple package {\tt MultiAlmkvistZeilberger}
that accompanied [ApZ]. The relevant procedures
(in {\tt AppsWZmulti}) are {\tt LatticePaths} and
{\tt LatticePathsDiagonal} and for the verbose versions
{\tt LatticePathsStory} and {\tt LatticePathsDiagonalStory}.

For example, if you type
{\tt LatticePaths($\{$[0,1],[1,0],[1,1]$\}$,m,M);} you would
get that the following two operators annihilate 
$F(m_1,m_2)$, the number of ways of going from $(0,0)$ to
$(m_1,m_2)$ using the steps $\{ (0,1),(1,0),(1,1)\}$
(where $M_1,M_2$ are the shift operators in the $m_1, m_2$
variables respectively)
$$
[-(m_1+1)/(2+m_1)-((1+2m_2)/(2+m_1))M_1+M_1^2, 
-(m_2+1)/(2+m_2)-((1+2m_1)/(2+m_2))M_2+M_2^2] \quad,
$$
which in everyday parlance means that
$F(m_1,m_2)$ satisfies
$$
(m_1+2)F(m_1+2,m_2)-(2m_2+1)F(m_1+1,m_2)-(m_1+1)F(m_1,m_2) =0 \quad,
$$
$$
(m_2+2)F(m_1,m_2+2)-(2m_1+1)F(m_1,m_2+1)-(m_2+1)F(m_1,m_2) =0 \quad.
$$

If you type 
{\tt LatticePathsDiagonal($\{$[0,1],[1,0],[1,1]$\}$,n,N);} you would
get that $f(n)$, the number of ways of getting from $(0,0)$ to
$(n,n)$ using the same set of steps is:
$$
(n+2)f(n+2)-3(2n+3)f(n+1)+(n+1)f(n)=0 \quad ,
$$
subject to the initial conditions  $f(0)=1$,$f(1)=3$.
Thanks to Birkhoff-Trijinski, its asymptotics is
$$
C(3+2^{3/2})^n n^{-1/2}(1+
((3/32) 2^{1/2}-1/4)/n+(113/1024-(9/128) 2^{1/2})/n^2+
((1545/32768) 2^{1/2}-245/4096)/n^3) \quad,
$$
for some constant $C$.

A more interesting example are the (old-time) {\bf basketball numbers},
which is the number of ways a basketball game that
ended with the score $n:n$ can proceed. Recall that
in the old days (before 1961), an atom of basketball-scoring
could be only of one or two points.
Equivalently, this number is the number of ways of walking, in
the square lattice, from $(0,0)$ to $(n,n)$ using the
atomic steps $\{(1,0),(2,0),(0,1),(0,2)\}$. Entering this
into {\tt LatticePathsDiagonalStory} yields that,
calling this number $F(n)$, it satisfies the third-order
linear recurrence:
$$
(16/5)(2n+3)(11n+26)(1+n)/((n+3)(2+n)(11n+15))F(n)
$$
$$
-(4/5)(121n^3+649n^2+1135n+646)/((n+3)(2+n)(11n+15))F(1+n)
$$
$$
-(2/5)(
176n^2+680n+605)/((11n+15)(n+3))F(2+n)+F(n+3) = 0 \quad ,
$$
subject to the initial conditions:
$$
F(0) = 1, F(1) = 2, F(2) = 14 \quad .
$$

For the record, the first few terms are:


[1, 2, 14, 84, 556, 3736, 25612, 177688, 1244398, 8777612, 62271384, 
443847648, 3175924636, 22799963576, 164142004184, 
1184574592592, 8567000931404, 62073936511496, 450518481039956, 
3274628801768744, 23833760489660324].

The asymptotics is:
$$
(.37305616)(4+2 (3^{1/2}))^n n^{-1/2}(1+(67/1452)3^{1/2}-
(119/484) )/n
$$
$$
+((6253/117128)-(7163/234256) 3^{1/2})/n^2+(-(32645/15460896)
3^{1/2}+(129625/10307264))/n^3) \quad,
$$
or in floating-point:
$$
(.37305616) (7.464101616)^n n^{-1/2}\cdot
$$
$$
(1.-.1659453140/n+.42398086 \cdot 10^{-3} /n^2
+.8918933381 \cdot 10^{-2} /n^3) \quad .
$$

\vfill\eject

{\bf Fifth Application: Random Walk in Higher Dimensions}

This is the multivariate analog of the First Application.
The relevant procedures in {\tt AppsWZmulti} are
{\tt RandomWalkRecurrence} and {\tt RandomWalkRecurrenceE}.
By now readers should be able to generate their own examples.
A few sample input and output files are given in the
webpage of this article.

{\bf References}

[AlZ] G. Almkvist and D. Zeilberger,
{\it The method of differentiating under the
integral sign}, J. Symbolic Computation {\bf 10}, 571-591 (1990).

[ApZ]  M. Apagodu and D. Zeilberger
{\it Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the
Sharpening of Wilf-Zeilberger Theory },
Adv. Appl. Math. {\bf 37} (2006), 139-152.

[PaStu] L. Pachter and B. Sturmfels, {\it Statistics},
in: ``{\it Algebraic Statistics for Computational Biology}''
[L. Pachter and B. Sturmfels, eds.], Cambridge University Press,
2005, 3-33.


[Po] G. P\'olya , {\it On picture writing},
Amer. Math. Monthly {\bf 63} (1956), 689-697.

[St] R. Stanley, {\it Differentiably finite power series},
Europ. J. Combinatorics {\bf 1} (1980), 175-188.

[WilZ] H.S. Wilf and D. Zeilberger,
{\it An algorithmic proof theory for hypergeometric
(ordinary and ``q") multisum/integral identities}, Invent. Math. 
{\bf 108}(1992), 575-633 .

[WimZ] J. Wimp and D. Zeilberger,
{\it Resurrecting the asymptotics of linear
recurrences}, J. Math. Anal. Appl. {\bf 111}(1985), 162-177 .
  

[Z] D. Zeilberger, {\it A Holonomic systems approach to special functions
identities}, J. of Computational and Applied Math. {\bf 32} (1990),
321-368 .



-----

Moa Apagodu, Department of Mathematics, Virginia Commonwealth University,
Richmond, VA 23284-2014,  USA.
mapagodu  at vcu dot edu.

Doron Zeilberger,
Department of Mathematics, Rutgers University (New Brunswick),
Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway,
NJ 08854-8019, USA. zeilberg at math dot rutgers dot edu .


\end

