GOSPERS ALGORITHM 
A DECISION PROCEDURE FOR PARTIAL
HYPERGEOMETRIC SUMS
by
Richard B. Reeves
B.A., University of Colorado at Denver, 1990
A thesis submitted to the
" Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
1992
This thesis for the Master of Science
degree by
Richard B. Reeves
has been approved for the
Department of
Mathematics
Date
Reeves, Richard B. (M.S., Applied Mathematics)
Gospers Algorithm
A Decision Procedure for Hypergeometric Sums
Thesis directed by Associate Professor William Cherowitzo
Given an integer function fk, does there exist a function S such that
m
52 ft = S(m) S(0) or, equivalently fk = S(lc) S(k1) ?
jfc=i
An algorithm created by R.W. Gosper in 1978 gives a complete answer to
this question. In creating this algorithm Gosper has made a decisive step
in the mechanization of hypergeometric identity proofs.
The form and content of this abstract are approved. I recommend its
publication.
Signed
m
To my parents Jim and Millie
for their unwavering love, support and encouragement
and
To Terry
for her friendship, intellect, humor, and for
her love that compels me to see the beauty in life.
ACKNOWLEDGEMENTS
I would like to thank Bill Cherowitzo, Stan
Payne, and Zenas Hartvigson for sharing their
expertise. Above all, I would like to thank
Dr. Cherowitzo for his guidance, support, and
patience.
v
TABLE OF CONTENTS
Introduction .............................................. 1
Chapter 1: Preliminaries ................................. 2
1.1 Hypergeometric Functions....................... 2
1.2 Hypergeometric Terms .......................... 6
1.3 General Method for Finding Hypergeometric
Representation ................................ 8
Chapter 2: Gospers Algorithm ........................... 14
2.1 Partial Hypergeometric Sums................... 14
2.2 Algorithmic Method ........................... 15
Chapter 3: Mechanized Proofs............................. 24
Appendix
Hypergeometric Program...............................29
References ................................................55
vi
Introduction
Let fk be a complex values function of the integer k. Given such
a function fk, does there exist a function S (within an additive constant)
such that
E ft > so)
*=1
or, equivalently
/* = m s(k1) ?
Such a sum of the function f in finite calculus would be called an
"indefinite sum", in analogy with the indefinite integral of infinite calculus.
An algorithm created by R. W. Gosper in 1978 gives a complete answer
to this question. In other words Gospers algorithm will either return the
function S (with the desired properties), or it will guarantee that no such
function exists. It never returns the statements "not sure" or "unknown" to
this question. This paper examines Gospers Algorithm.
1
UL
Chapter 1: Preliminaries
1.1 Hypergeometric Functions
Gosper has discovered a way to decide whether a given function is
partially (indefinitely) summable with respect to a general class of
functions called hypergeometric terms. Therefore, in order to investigate
Gospers algorithm, it is first necessary to understand hypergeometric
series or functions.
Hypergeometric (HG) series were first investigated by Euler,
Gauss, Pfaff, and Riemann. The study of HG series has evolved into a
systematic means of handling a great variety of binomial coefficient
identities [5].
The following notation will be useful in the discussion of HG series.
Definition 1.1: Rising Factorial Powers
xm = x(x + l)...(x + ml) integer mz0 .
Definition 1.2: Falling Factorial Powers
x = x(xl)...(xm + l) integer mz0.
2
The general hypergeometric series is a power series in z (z e c),
and is defined as follows in terms of rising factorial powers:
Definition 1.3:
1 k
*am E teO *1 am zk
A 1 Z y bl A* k\
where, to avoid division by zero, none of the b parameters may be zero or
a negative integer. Using standard conventions, we define = 1 for all
x and 0 = 1. Therefore, the first term (k = 0) of any HG series is equal
to 1.
Many important functions occur as some special case of the general
hypergeometric series. This is one of the reasons hypergeometric
functions are so powerful. The following examples illustrate a small
portion of this power.
Example 1.1 (The simplest situation no parameters):
771 = 0, 71 = 0
= =
to k\
This "blank" notation does make sense, when m or n equals zero, because
the notation indicates the absence of parameters. However, this form has
proven to be perplexing to many. In order to avoid confusion, it has
3
Hi
become standard to add an extra one above and below whenever m = 0
or whenever n = 0. So,
F
z
fcO JtaO 1* kl
It is now important to note that, when using HG series in Gospers
algorithm, we are not concerned with convergence of the series. We are
simply using z as a formal symbol. Therefore, any identities derived will
be formally true. An example of this is the HG function F
k\zk.
This series doesnt converge for any nonzero z but it is still very useful in
solving problems [5].
Example 1.2 (specific example):
m = 1, n = 0, ax = 1 written as
m = 2, ax = a2 = 1, and n = 1, bl = 1
F
1
1 z
This is the geometric series. Hypergeometric series are so named
because the geometric series is included as a very specific case.
4
Example 1.3 (general example):
m = 1, n = 0,
F Y z
JhO
E
*s0
k\
hi 0
()(q+l)~(g+*l) *
it!
a+k1]k _
k r 
(1 zr
Example 1.4:
_ (a)* (~zf = y (fl)(fll)(flk+l) k
Â£o faO
E
k>0
This utilizes the binomial theorem. Note, a negative integer as an upper
parameter causes the infinite series to become finite since (a)1 = 0
whenever k > a > 0 and a is an integer.
Example 1.5 (The Gauss Hypergeometric so named because
many of its properties were first proved by Gauss in his doctoral
dissertation in 1812; although Euler and Pfaff had already discovered
many of these properties [5] ):
a,b 
= E^
albl zl
ki 0 Ck k\
5
Before everything was generalized to arbitrary m and n (around 1870), this
was considered the hypergeometric series.
One important special case of the Gauss Hypergeometric is:
k\k\ (z)*
S (*+l)!
to (* + l)
1+ + + ...
z2 z3 Z4
= z+ + ...
2 3 4
Note: z'1 In (1 + z) is a hypergeometric function but In (1 + z) is not,
since by convention every hypergeometric function has a value of 1, when
z = 0.
1.2 Hypergeometric Terms
It is essential to realize which series can be expressed as a
hypergeometric series. This determination is made by examining the ratio
of consecutive terms of any hypergeometric series. Let
\ k l
a\,am
I z
al "am Z'
= y' t, where t. =
to * *
As previously stated t0 = 1, and all other terms have the following ratio:
6
# Jfc+1 *+l h+1 al am blb] z*+1 k\
t. k k aX"am bf'bY1 (*+1)! zk
(aj +k)(am+k) z
(^+*3(&+*) k+1 '
Therefore, a hypergeometric series is a series Â£ tt in which the ratio of
ki o
every two consecutive terms is a rational function of the summation
variable k. This form exhibits the fact that the as are the negatives of the
roots of the numerator and the bs are negatives of the roots of the
denominator of the term ratio. In order to maintain a canonical form for
the term ratio, if (k+1) is not a factor in the denominator, we add it to
both the numerator and the denominator.
To restate, hypergeometric series are those series whose constant
term (first term) is 1 and whose term ratio is a rational function of the
summand variable. Given a rational function, we may always find a
hypergeometric series whose term ratio is this function.
Example 1.6:
= k2n+12 = (fc3)(fc4)(l/8)
h 8*M Qcnimm)
adding the missing (k+1) we get
7
= (*3)(*4)(l/8)(fc+l)
(k+if4)(ki/4)(k+l)
This is the term ratio of
F
3,4,i.n
*/4, i/4 g I *
1.3 General Method for Finding Hypergeometric Representation of a
Given Function S.
Step 1) Write s = Â£ tk, t0 0.
ki 0
Step 2) Calculate the term ratio.
Step 3) If term ratio is not a rational function of k, the series
cannot be expressed in hypergeometric form.
Step 4) Otherwise, find the negatives of the roots of the
numerator and denominator and express in the
following form:
(at+Â£) (aw+Â£) z
(bx+K) (bn+k) k+1
This gives the parameters .. ,am; bt....bn; and the
argument z such that
8
S <0F
IVA1
Example 1.7
Step 1)
(S = sin (z)):
sin(z) = z
sin(z) = Y,
kiO
z'1 sin(z) =
(1)* z(2**1}
(2fc+l)!
r (1)* z2*
fa (2&+1)!
first term is z.
first term is 1.
Step 2)
Step 4)
= (1)**1 z2**2 (2fc+l)i = (1) z2
h (1)* z2* (2fc+l)! (2*+3)(2fc+2) '
4 (k+3/2)(k+l)
sin(z) = z F
1
1,3/2
z
\
4
9
One of the strengths of hypergeometric series is the ability to
express and manipulate a great variety of binomial summations in a
systematic way.
Example 1.8:
A well known binomial coefficient identity is the parallel
summation law:
E
fr+V + ... +
, * > ll , V
n
(r+n+l\
integer n .
If we wished to express the parallel summation law in its hypergeometric
form we would proceed as follows:
Step 1)
E
'r+k' 'r+n+r
, n ,
integer n .
If we replace k with nk we get
E E E
nkin kiO kzO
r+nk'
nk ,
(r+nk) 1
ho (nk)l r!
where the first term
_ (r+w)l
, n j n\ r!
10
Step 2)
h*\ (r+nfc1)! (nfc)l r! _ (nk)
h (r+nk)\ (nÂ£l)! r! (r+nk)
Step 4)
(nk) _ (nk) (fc+1)
(r+nk) (r+nk) (k+l)
Ir+n) f n, 1
{ n J
r+n+lN
*
Solving for the hypergeometric function in this equation, we obtain
F
~n, 1
rn
\
1
/
r+n+1
r+1
i/
r+n
n
* 0 .
Note that since the lower parameters of a hypergeometric cannot
be zero or negative integers, this represents a degenerate case of the
parallel summation law because parallel summation is usually applied
when r and n are positive integers. This limitation is bypassed in this and
other situations where it arises by considering:
lim
e0
F
1 ,~n
nr+z
1
11
Example 1.9 (Vandermondes Convolution Identity):
E
r
5
, n ) '
A combinatorial description of this identity is the sum of all
possible ways to pick k men from a total of r men and the number of ways
to pick nk women from a total of s women equals the number of ways to
choose n people from a group consisting of r men and s women.
Let
_________rl s\__________
k\ (r k)! (n k) \(sn+k)\
We can make several general observations about the process of converting
to hypergeometric form [5].
1) Whenever tk contains a factor like (a+ifc)!, with a plus sign before
the k, we get = (a+jt+i) in the term ratio. Which
(a +k)\
contributes parameter (a+1) to the corresponding HG function, as
an upper parameter if (a+Jfc)! was in the numerator of tk, or as a
lower parameter if it was in the denominator of tk.
2) Whenever tk contains a factor like (ait)!, with a minus sign before
the k, we get (a~^~1)? = XlH in the term ratio. Which contributes
6 (ait)! (Jta)
12
parameter (a) to the corresponding HG function, as an upper
parameter if (aJfc)! was in the denominator tk and as a lower
parameter if it was in the numerator of tk. The above also negates
the argument of the HG function.
3) Whenever tk contains a factor like s! in the current example, that
is independent of k, it remains in ^ but disappears from the term
ratio.
Applying the observations to Vandermondes Convolution Identity, we
easily obtain
k\ (rk)! (nk)! (sn+k) l
r! s!
so,
(rk)(nk) _ (kr)(kn)
tk (&+l)(sn+&+l) (fc+l)(fcn+s+l)
This implies,
where
13
Chapter 2: Gospers Algorithm
2.1 Partial Hypergeometric Sums
We have seen how hypergeometric sums can be used to
systematically express a large number of binomial coefficient sums. It
would appear that partial sums may be needed in many instances
(combinatorial solutions, physical problems, and computational complexity
determination to name a few). When dealing with such sums it is often
desirable to find a closed form solution that works over a general range.
Gospers Algorithm decides whether a given summand (function) is
partially summable (within an additive constant) with respect to a general
class of functions called hypergeometric terms by examining
E /* = *() 5()
*=i
or, equivalently
ft = S(k) S(k1).
Given fk the algorithm finds those S(k) with the property
S(]A
^4 = a rational function of k .
S(k1)
(1)
2.2 Algorithmic Method [4]
14
Step 1
if sm
S(k1)
is a rational function of k then the term ratio
Sjk) j
fk = S(k) S(k1) = S(k1)
fk_x S(k1) Sik2) l S(k2)
is also a rational function of k. Therefore, we can express the term ratio
of (1) as
fk Pk
fk i Pk i rk
(2)
where pk, qk, and rk are polynomials in k subject the following condition:
(fc+a)  Qk and (fc+P)  rk
a p is not a nonnegative integer.
(3)
It is always possible to achieve this condition with a rational
function with the following steps:
a) Since the term ratio is a rational function of k, we can
express this ratio in hypergeometric form,
fk (ax+k) ~{am+k) z
= (V*) (V*)
b) Set:
15
and
Pti = 1.
qk = (k+aj  (k+aj z,
rk = (k+bi> (k+bn>
c) If (3) is violated then go to d.
d) If q and r have factors (k + a) and (k + /?) respectively that
violate (3), then divide them out of q and r and replace pk
by
pt(k+a)s = pt(I+a)(I+al) "(k+p+1) where N = ap.
Pk
The new p, q, and r still satisfy (2) smce  is replaced by
Pk1
Pk (fc+q)
Pk.i (Â£+P) *
e) Repeat step c until condition (3) is satisfied.
Step 2
Recall the goal is to find some hypergeometric term S(k) such that
(1) is satisfied. Assume
fk= S(k) S(k1)
and now write
S(fc) = Â£(*)/* (4)
Pk
16
where g(k) is a function to be found.
Combining (1) and (4) we get
e(k) = _ P* 1
" qM S(k)S(k1) qM J S(k1)
S(k)
This implies g(k) is a rational function of k whenever S(k) / S(kl) is.
By substituting (4) into (1) we get
A g) A g1) A,
Pt pt1
Multiplying through by pk/fk and using (2) (to solve for rk) we get
Pk = S(k) ~ rk g(k1) . (5)
If we can find a g(k) that satisfies this recurrence, weve found S(k) and
we are done.
Theorem: If S(k) / S(kl) is a rational function of k, then g(k) is a
polynomial.
Proof: We know that g(k) is a rational function when S(k) / S(kl)
is, so suppose by way of contradiction
g(k) = where v(k) 1
v(k)
17
and u(k) and v(k) have no factors in common. Let N be the
largest integer such that (k+/3) and (k+/? + N) both occur as
factors of v(k) for some (3 e C. The value of N is
nonnegative, since N = 0 will always satisfies this condition.
Substituting into equation (5) we get
u(k) u(k1)
p*= m ~
Now, equation (5) can be rewritten
Pk v(k) v(k1) = qkn u(k) v(fcl) rk u(k1) v(k) .
If we set k = /? +1 and k = 13N we get
r(p+l) u(P) v(p+l) = 0
and
q(VN+l) u(pA0 v(pIVl) = 0 .
Now u(/3) f 0 and u(/3N) f 0, because u and v have no
factors in common. Also, v(/3 + l) + 0 and v(/3Nl) + 0,
because v(k) would otherwise contain the factor (k+/3l) or
(k+/3 + N+1) which is contrary to N being maximal.
Therefore, r()3 +1) = 0 = q(/3N+1). But this contradicts
condition (3). Hence g(k) must be a polynomial [5].
18
Step 3
All that remains is finding, given pk, % and rk, if there exists a
polynomial g(k) that satisfies equation (5). This turns out to be
surprisingly easy, once we know the degree d of g(k).
Let
g(k) = adkd + aJ_1*d"1 + + a0 ad 0
(6)
Rewrite equation (5) in the following form:
2p(k) = DiffQlb x (g(k) + 8(k~ 1)) SumQ(k) x (g(k) g(k1))
where
DiffQ{k) = 9(^+1) r(k) and SumQ(k) = q(k+1) + r(k).
(7)
If g(k) has degree d then sum (g(k) + g(kl)) = 2 adkd + ... also has
degree d. Whereas the difference (g(k) g(kl)) = dadkd_1 + ... has
degree d1 (assume the zero polynomial has degree 1.) Write deg(p) for
the degree of any polynomial p(k).
Case 1: The deg(DiffQ) > deg(SumQ). This implies the degree
of the righthand side of (7) is deg(DiffQ) + d, so d =
deg(p) deg(DiffQ).
Case 2: The deg(DiffQ) < deg(SumQ). This implies the we can
write
19
DiffQik) = $kd'' +
and
SumQ(k) = ykd> + where y 0.
The righthand side of (7) now has the form
(2f>ad + ydud)kJ*d>~1 + This leaves two possibilities:
1) 2p + yd 0  J = deg(p) deg{SumQ) + 1,
or
2) 2p+yJ = 0 = Â£?> degip) deg(SumQ) + 1.
This second case needs to be examined only if 2)0/y is an
integer d greater than deg(p) deg(SumQ) + 1.
If d < 0, the algorithm stops with the message that no function S,
in equation (1), exists. If a nonnegative d is discovered, we can plug
equation (6) into equation (5). The polynomial g(x) will satisfy the
recurrence if and only if the as of equation (6) satisfy certain linear
equations, because each power of k must have the same coefficient on
both sides of (5).
Example 2.1 (This example was computed using the Derive
software package on an IBM compatible 486 PC):
20
Use Gospers algorithm to find a closed form identity for the following
summation:
m
E
IJ aj2 + bj + c
Jzl_________________
*+i
IJ aj2 + bj + d
y=i
Step 1
JJ aj2 + bj + c
f(k) = Â£1
J k* 1
n aJ2 + bi+ d
Step 2
Let
k1 k
n aJ2 + bj + c n +bJ + d
f(k) = 21____________________ J=1__________
f(k1) ** *2
n aj2 + bj + d JJ aj2 + bj + c
;=i ;'=i
a(k1)2 + fr(fcl) + c
a(fc+l)2 + fc(lfc+l) + d
q(k) = a(A;l)2 + Â£(fcl) + c
r(k) = a(k+l)2 + b(jc+1) + d
21
and suppose q and r have no factors in common. This implies pk = 1.
Therefore,
deg(qk+1 rk) = 1 < 2 = deg(qk+i + rk)
which requires Case 2, where we discover
2P/y = 2 > 1 = deg(p) deg(SumQ) + 1 .
So g will be a quadratic function of k with three undetermined
coefficients.
Let g(k) = g2k2 + g^k + g0
Step 3
Substituting into equation (5)
1 = k2 ( g2(2a+b+cd) + ^(a) )
+ k ( g2(b+2d) + gt(a + cd) + g0(2a) )
+ ( g^abd) + g^a + b+d) + (ab + cd) )
and solving for the coefficients of g we get the augmented matrix
f2a+b+cd a 0 O'
b+2d a+cd 2a 0
, abd a+b+d ab+cd 1,
22
which reduces to
1 0 0
0 1 0
______________2^_______________
(a 2 +2a(c +d) (b+c d)(bc +d)){c d)
_________2a(2a+b+cd)__________
(a 2 +2a(c +d) (b+c d)(bc +d))(c d)
q q ^ 2az+a(2b+3cd)+(cd)(b+cd)
((a 2 +2a(c +d)(b+cd)(bc +d))(c d) t
Having found the above coefficients for g we can now find S(k).
S(k) = q(k+1) g(k) M)
(ak2+bk+c)(2a2(k2+2k+l)+(a(2b(k+l +k+3)
_________________________________________________________________i1__________
*1
(a2+2a(c+d)(b+cd)(b~c+d))(cd)J^(aj2+bj +d)
j=i
Finally we discover the closed form of this partial sum,
*i
n }2 bi *c
e *
i=l
k+1
S(m) 5(0)
JJ aj2 + bj + d
/=i
ml
(am2+bm +c)(2a2(m2 +2m+1) +a(2b(m+1) +c(2m+3)d(2m+1)) +(cd)(b+c(aj2+bj+c)
________________________________________________________________M___________
m+1
(a2 +2a(c +d)(b+c d)(bc +d))(c 
>1
_ 2a2+a(2b+3cd)+(cd)(b+cd)
(a2+2a(c+d)(b+cd)(bc+d))(a+b+d)(cd)
23
Chapter 3: Mechanized Proofs
Combinatorial problems, such as the previous example, arise in
many areas of mathematics, engineering, and theoretical physics. Their
solution often entails the evaluation of a sum of products of binomial
coefficients [7]. A great many distinguished mathematicians have spent an
immense amount of time and ingenuity proving virtually thousands of
binomial identities. We now know that much of that time and talent was
wasted. Recently it has become widely known that "almost all" known
binomial coefficient identities are special cases of relatively few
hypergeometric identities [3]. Therefore, there are a comparatively small
number of essentially different binomial identities. The problem in
recognizing this has been twofold. First, the difficulty in recognizing
equivalent binomial expression due to the notation used. The notation
used to express binomial coefficients makes it very easy to disguise
equivalent expressions in a such a way that no relationship is readily
apparent. The second problem has been the reluctance of many to see the
advantage of using hypergeometric series to serve as a canonical form for
these identities (even though this was the method used by Euler and
Gauss [7]). Although this reluctance has diminished over the last decade
24
and the use of hypergeometric series has increased, the need for identities
has not decreased. This need for hypergeometric identities is being
largely supported by the use of computers.
The increased use of hypergeometric series has directly lead to an
increased use of computers in creating hypergeometric identities. Wilf
asserts, "that more and more hypergeometric identities are still being
conjectured and proved, so the need remains for mechanized proof' [3,78].
It would have been difficult (and in practice infeasible) to create
the identity in example 2.1, without the aid of a computer. For this
particular example, the software program Derive was used to find the
appropriate identity. Gosper credits the program MACSYMA and the
support of its developers at MIT for giving him ... the experiences
necessary to provoke the conjectures that led to this algorithm" [4, 42].
Doron Zeilberger has derived other identities (using algorithms largely
based on Gospers algorithm) using the MAPLE software package.
Finding identities and mathematical proofs using computers is becoming
more the norm [10]. Gospers algorithm is example of mechanized
theorem proving.
Mechanized theorem proving has been a goal for a long time. The
desire dates back to Leibniz in the late 17th and early 18th centuries[2].
25
In 1930 Herbrand proposed a system for mechanized (or automated)
theorem proving based on logic rules[6]. His method was slow and
untimely, since it came before the advent of the digital computer. With
the arrival of the computer came a renewed interest in automated
theorem proving. In the 1960s Herbrands method was implemented on
a computer by Gilmore[2]. Until 1965, "success was largely confined to
problems found in logic books rather than mathematics books." [6, 5] In
1965, J.A. Robinson developed a major breakthrough. She produced the
resolution inference system. This system was based on a single inference
rule and did not require generalized logic axioms (only those specific to
the problem). In essence the inference rule is:
an algorithm that, when successfully applied
to some given set of hypotheses or premisses,
yields a conclusion that follows inevitably and
logically from the premisses. [1,197]
In other words, the algorithm is a finite automata that stops. Gospers
algorithm "infers" that a closed form can be found if a certain polynomial
can be found. If this polynomial doesnt exist, Gospers algorithm infers
that no closed form can be found. Since 1965 many refinements of the
resolution system have been made.
Wilf and Zeilberger have shown that a very large class of identities
can be proved by computers [9]. They do not imply that computers need
26
to be "trusted blindly"; in fact, it is just the contrary. Although the proofs
are "discovered" by a computer, all proofs presented by Wilf and
Zeilberger produce a proof certificate that can be checked by hand if
desired [9]. Traditionalists might complain that identity proofs generated
by a computer can give no insight to the mathematician. Zeilberger
cautions against being too chauvinistic, "The team humancomputer is a
mighty one, and an openminded human can draw inspiration from all
sources, even from a machine" [10].
I found it very enlightening to write my own program (see
addendum) partially implementing Gospers algorithm. This program will
take as input the term ratio (already factored) of a specific hypergeometric
series and output a partial sum or the fact that one doesnt exist. The
purpose in writing this program was not to improve or advance
implementations of Gospers algorithm already in use. The purpose was
strictly one of education. Perhaps most edifying were the attempts
(admittedly small) at symbolic manipulation. I gained not only an insight
into mechanized proofs, but an appreciation of the software techniques use
in packages such as Derive.
In 1978 Gosper made a definitive step in the direction of
mechanized proofs of hypergeometric series. His algorithm is now fully
27
implemented in MACSYMA and partially implemented in Mathematica.
It will no doubt become more widely implemented when its usefulness
becomes more fully recognized [8].
28
Appendix A
{Rick Reeves}
{1992}
Program Hyper (Input, Output);
{This program is an attempt to implement R.W. Gospers algorithm perhaps
not in its complete form but perhaps in a way that gives the user
an appreciation of what the algorithm entails. }
{The algorithm is not shown here due to its length. Please see
Gospers original paper.}
USES
Gauss, { a gaussian elimination unit for n by n matrices }
Poly_l, { a polynomial unit add, subtract, multiply, show, etc.}
Crt; { standard crt unit }
TYPE
Parameter = array [1..MAXP] of real; { used for factored form of poly}
{}
Procedure Pause;
{ Holds user screen until key is press then clears screen.}
VAR
Ch : Char;
Begin
GOTOXY (27,25);
Write (Press any key to continue);
Ch := ReadKey;
CLRSCR;
end; {pause}
29
{}
Procedure Initialize (VAR UpperP, LowerP, P_array : Parameter);
{ This procedure initializes all the parameter arrays. }
VAR
i : integer;
Begin
for i : = 1 to MAXP do
Begin
UpperP[i] := 999;
LowerP[i] := 999;
P_array[i]:= 999;
end; {for}
P_array[l] := 9911; {flag for step one}
end; {initialize}
{}
Procedure GetHyper (VAR UpperPara, LowerPara : Parameter;
VAR UpNumber, LowNumber : Integer;
VAR Arguement : Real);
{ This procedure gets the factored hypergeometric term ratio from
the user. Note: this maybe the same as the hypergeometric
parameters except for perhaps the k +1 factor in the denominator.}
VAR
i : integer;
Begin
GOTOXY(10,3);
Write (OUTPUT, Please enter NUMBER of UPPER parameters: );
Readln (INPUT, UpNumber);
30
For i : = 1 to UpNumber do
Begin
GOTOXY (15,5 +1);
Write (OUTPUT, Please enter UPPER parameter number );
Readln (INPUT, UpperPara[i]);
end; {for}
GOTOXY(10,14);
Write (OUTPUT, Please enter NUMBER of LOWER parameters: );
Readln (INPUT, LowNumber);
For i : = 1 to LowNumber do
Begin
GOTOXY (15,15 + i);
Write (OUTPUT, Please enter LOWER parameter number ,i,: );
Readln (INPUT, LowerPara[i]);
end; {for}
GOTOXY (10,24);
Write (OUTPUT, Please enter value of arguement z: );
Readln (INPUT, Arguement);
end; {gethyper}
{}
Procedure StepOne (VAR Q_array, R array, P_array : Parameter;
Up, Low : integer);
{Get p(k), q(k), and r(k) in right form, ie, put term ratio in the
proper form. }
VAR
i, j, p, k, x : integer;
GOOD : boolean;
check : real;
Chk : integer;
31
Begin
p := 1;
{check to see if (alpha beta) is a positive number}
Repeat
GOOD := TRUE;
for i: = 1 to (Up) do
for j := 1 to (Low) do
begin
if (Q_array[i] < > 999.0) and (R_array[j] < > 999.0) then
begin
check := Q_array[i] R_array[j];
if ( check > 0 ) and (int(check) = check) then
begin
Chk : = trunc(check);
GOOD := FALSE;
x := 1;
for k : = p to (p + Chk 2 ) do
begin
P_array[k] : = Q_array[i] (x);
x := x + 1;
end; {for}
p : = p + chk 1;
Q_array[i] := 999;
R_array[j] := 999;
end; {if}
end; {if}
end; {for}
Until GOOD;
for i : = 1 to up do
Q_array[i] := Q_array[i] + 1;
end; {stepone}
32
{}
Procedure Display ( Q_ay, R_ay, P_ay : Parameter;
Up, Low : Integer;
Argue, Which : Real);
{ This procedure displays the polynomials in their factored form. }
VAR
i,j : integer;
Begin
GOTOXY (10,10);
if which = 1.0 then
Write (q(k) = )
else
Write (q(k+l) = );
For i : = 1 to Up do
if Q_ay[i] < > 998 then
if Q_ay[i] = 0 then
write 0 (k) )
else if Q_ay[i] < 0 then
write ( (k ,abs(Q_ay[i]):l:l,) )
else
write ( (k + ,Q_ay[i]:l:l,) );
Writeln ( ( ,Argue: 1:1, ) );
GOTOXY (10,12);
Write (r(k) = );
For i := 1 to (Low) do
if R_ay[i] < > 999 then
if R_ay[i] = 0 then
write ( (k) )
else if R_ay[i] < 0 then
write ( (k ,abs(R_ay[i]):l:l,) )
else
write ( (k + ,R_ay[i]:l:l,) );
GOTOXY (10,14);
i := 1;
33
If P_ay[l] < > 9911 then
begin
Write (p(k) = );
While P_ay[i] < > 999 do
begin
if P_ay[i] = 0 then
write ( (k) )
else if P_ay[i] < 0 then
write ( (k ,abs(P_ay[i]):l:l,) )
else
write ( (k + ,P_ay[i]:l:l,) );
i := i + 1;
if (WhereX >70) then
Writeln;
end; {while}
end {if}
else
writeln (?p(k) = 1);
writeln;
end; {display}
{}
Procedure Degree ( Q_ay, R_ay, P_ay : Parameter;
Up, Low : Integer;
Argue : Real;
VAR P,Q,R : Poly Type;
VAR DegreeG,DegR,
DegQ, DegP : Integer);
VAR
i,
DegDiff, DegSum, IntPivot : Integer;
PTerm, Qterm, Rterm, Temp, QSum, QDiff : Poly_Type;
First : Boolean;
Beta, Gamma, Pivot : Real;
{ This procedure finds the degree of the unknown function G if it
exists. If it doesnt exist, a message is outputted that indicates
that. This procedure also finds the degree of all pertinent
polynomials and puts them in standard form.
This procedure requires units: Polyl and Crt}
34
{ initialize all variables to zero.}
Begin
DegR := 0;
DegQ : = 0;
DegP : = 0;
DegSum : = 0;
DegDiff: = 0;
DegreeG := 0;
First := TRUE;
Get (Temp);
Temp[0] := Argue;
{ find degrees of p, q, and r}
for i := 1 to (Low+0) do
if (R_ay[i] < > 999) then
DegR : = DegR + 1; (degree of r}
if DegR = 0 then
R_ay[l] := 1;
for i : = 1 to (Up) do
if (Q_ay[i] < > 998) then
DegQ : = DegQ + 1; (degree of q}
if DegQ = 0 then
Q_ay[l] := 1;
for i : = 1 to MAXP do
if (P_ay[i] < > 999) and (P_ay[i] < > 9911) then (degree of p}
DegP := DegP + 1;
Get (R);
Get (Q);
Get (P); ( get "zeroed" coefficient arrays }
GeT (Qterm);
Get (Rterm);
Get (Pterm);
( find standard polynomial form for p, q, and r.}
35
For i : = 1 to Up do
if (Q_ay[i] < > 998) then
if First then
Begin
First := FALSE;
Qterm[0] := Q_ay[i];
Qterm[l] := 1;
Multiply (Qterm,temp, Q);
temp := Q;
end {if}
else
Begin
Qterm[0] := Q_ay[i];
Qterm[l] := 1;
Multiply (Qterm,temp, Q);
temp := Q;
end; {else}
Write (q(k+1) = );
Show(Q,DegQ);
First := TRUE;
Get (Temp);
Temp[0] := 1;
For i : = 1 to Low + 0 do
if (R_ay[i] < > 999) then
if First then
Begin
First := FALSE;
Rterm[0] := R_ay[i];
Rterm[l] := 1;
Multiply (Rterm,temp, R);
temp := R;
end {if}
else
Begin
Rterm[0] := R_ay[i];
Rterm[l] := 1;
Multiply (Rterm,temp, R);
temp := R;
end; {else}
36
{ find q}
{display q}
{ find r}
{display r}
Write (r(k) = );
Show(R,DegR);
First := TRUE;
Get (Temp);
Temp[0] := 1;
For i: = 1 to MAXP do
if (P_ay[i] < > 999) and (P_ay[i] < > 9911 ) then
if First then
Begin
First: = FALSE;
Pterm[0] := P_ay[i];
Pterm[l] := 1;
Multiply (Pterm,temp, P); { find p}
temp := P;
end {if}
else
Begin
PtermfO] := P_ay[i];
Ptermfl] := 1;
Multiply (Pterm,temp, P);
temp := P;
end; {else}
Write (p(k) = );
if DegP > 0 then
Show(P,DegP) {display p}
else
begin
Writeln (T);
P[0] := 1;
end;
Pause;
{ find term sums and differences }
Add (R,Q,QSum);
For i : = 0 to MAXP do
if QSum[i] < > 0 then
DegSum := i;
Write (Qsum(k) = );
Show (QSum, DegSum);
Subtract (Q,R,QDiff);
37
For i : = 0 to MAXP do
if QDiff[i] < > 0 then
DegDiff: = i;
Write (QDiff(k) = );
Show (QDiff, DegDiff);
{Will the degree of our mystery polynomial stand up and sign in please.}
Pivot := 0;
Beta := QDiffJDegDiff];
Gamma := QSum[DegSum];
if Gamma < > 0 then
Pivot: = 2 beta / gamma;
IntPivot := trunc(Pivot);
if DegDiff > = DegSum then
DegreeG : = DegP DegDiff
else if ( Pivot > (degP DegSum + 1) ) and ( Pivot = IntPivot) then
DegreeG := IntPivot
else
DegreeG : = (degP DegSum + 1);
Writeln (Degree of G = DegreeG);
DegreeG := 2;
end; {degree}
{ }
Function Factorial (n : integer): integer;
{ This is a recursive function that finds factorials.}
begin
if n = 0 then
Factorial := 1
else
Factorial := n Factorial(nl)
end; {factorial}
38
{
Function Combo (n,k :integer): integer;
}
{ This function finds combinations it uses function factorial.}
Var
temp : integer;
begin
if n < = k then
temp := 1
else
temp : = Factorial(n) Div (Factorial(k) Factorial(nk));
Combo := temp;
end; {combo}
{ }
Procedure BuildSK(Var Mat : Matrix_Type; Degree : Integer);
{ This procedure finds the term (in matrix form) g(kl) }
Var
i, j : integer;
Begin
For i := 1 to Degree+1 do
For j : = i to Degree +1 do
if odd (ji) then
Mat[i,j] := Combo(jl,ji)*(l)
else
Mat[i,j] := Combo(jl,ji) ;
For i := 1 to Degree+1 do
begin
For j := 1 to Degree+1 do
Write (Mat [i,j]:8:l);
writeln;
end;
Pause;
end; {buildSK}
39
{ }
Procedure BuildS (Var Mat : Matrix_Type; Degree : Integer);
{ This procedure finds the term (in matrix form) g(k) }
Var
i, j : integer;
Begin
For i := 1 to Degree+1 do
Mat [i,i] : = 1;
For i := 1 to Degree+1 do
begin
For j := 1 to Degree +1 do
Write (Mat [i,j]:8:l);
writeln;
end;
Pause;
end; {buildS}
{ }
Procedure BuildTerm (Mat : Matrix_Type; Poly : Poly_Type;
DegM,DegP : Integer; Var Term : Matrixjype);
{ This procedure is used to find the terms (in matrix form) q(k+l) g(k)
and r(k) g(kl). }
Var
i, j, k : integer;
Begin
For i : = 0 to DegP do
writeln (poly[i]:3:0);
For i : = 0 to DegP do
For j : = 1 to DegM + 1 do
For k : = 1 to DegM + 1 do
Term[i+j,k] := Term[i+j,k] + (Poly[i] Mat[j,k]);
40
For i := 1 to DegM+l + DegP do
Begin
For j := 1 to DegM+ 1 + DegP do
Write (Term [i,j]:8:l);
writeln;
end;
Pause;
end; {Term}
{ }
Procedure Subtract (Matl, Mat2 : Matrix_type; Var Answer : Matrb
{ This procedure is used to subtract terms that are in matrix form }
Var
i, j : integer;
begin
for i := 1 to MAXMATRIX do
for j := 1 to MAXMATRIX do
Answer[i,j] := Matl[i,j] Mat2[i,j];
end;
{OUTPUT}
Procedure OUTPUT(ANSWER : Answer_Type; B : Matrix_Type;
integer);
{ This procedure outputs to a file (Defaulted to Screen for now) the
Gaussian reduced matrix and the solution set or a message that
no unique solution exits.}
var
Outfile : text; {allows programmer to pick output device}
I, J : Integer; {loop and matrix counters}
begin
assign(Outfile, CON);
rewrite(Outfile);
_Type);
IERR,N
41
{ If there are not errors output the reduced matrix and the solution set.}
if IERR = 0 then
begin
writeln(Outfile);
writeln(Outfile);
writeln(Outfile,The Gaussian reduced matrix is as follows:);
writeln(Outfile);
for I: = 1 to N do
begin
for J : = 1 to (N +1) do
write(OUTFILE, B[I,J] :12:2);
writeln(OUTFILE);
end; {for}
writeln(Outfile);
writeln(Outfile,The solution set is as follows:);
writeln(Outfile);
for I : = 1 to N do
writeln(Outfile,G,I, = , ANSWER[I]:12:2);
writeln(Outfile);
end {if}
{ Else, output that there are no unique solutions to this linear system.}
else
begin
writeln(Outfile);
writeln(Outfile);
writeln (Outfile,There is no unique solution to this linear system.);
writeln(Outfile);
for I : = 1 to N do
begin
for J : = 1 to (N +1) do
write(OUTFILE, B[I,J] :8:1);
writeln(OUTFILE);
end; {for}
writeln(Outfile);
end; {else}
close(OUTFILE);
Pause;
end; {output}
42
{}
Procedure FindClosed (PolyG, P,Q,R : Poly_Type; DegG,DegR,DegQ,DegP
Integer);
{ This procedure is used to find the closed form of the summation.}
Var
Temp : Poly_Type;
begin
Get(Temp);
Multiply (Q,PolyG,Temp);
Multiply (Temp,R,Temp);
Show(temp, DegG + DegQ + DegR);
Pause;
end; {findclosed}
{ }
Procedure FindPoly (P,Q,R : Poly_Type; DegG,DegR,DegQ,DegP : Integer);
{ This procedure finds the coefficients of G (if they exist).}
Var
SK, SK_one : Matrix_Type;
Answer,
FirstTerm,
SecondTerm : Matrix_Type;
i,j, N, Error : Integer;
Coef : Answer_Type;
Max : Integer;
PolyG : PolyType;
Begin
error : = 0;
Initialize_Matrix (SK);
Initialize_Matrix (SK_one);
Initialize_Matrix (FirstTerm);
Initialize_Matrix (SecondTerm);
Initialize_Matrix (Answer);
43
Get (PolyG);
BuildS (SK,DegG);
BuildSK (SK one, DegG);
BuildTerm(SK, Q, DegG, DegQ, FirstTerm);
BuildTerm(SK_one R, DegG, DegR, SecondTerM);
Subtract (FirstTerm, SecondTerm, Answer);
if DegQ > DegR then { find # of rows}
N := DegQ
else
N := DegR;
N := N + DegG + 1;
if DegG+ 2 > N Then { find if more rows than cols}
Max : = DegG + 2
else
Max : = N;
For i : = 1 to MAXMATRIX do
Answer[i, DegG+2] := P[i1];
For i : = 1 to MAX do
begin
for j : = 1 to Max +1 do
write (Answer[i,j]:8:l); (show matrix}
writeln;
end; {for}
Pause; { find coefficients if possible}
N := DegG+1;
Gaussian_Elim (Answer,ERROR,N Coef);
Output ( Coef, Answer,ERROR, N);
For i := 1 to DegG+1 do
PolyG[il] := Coef[i];
write (g(k) := );
show(PolyG, DegG);
Pause;
FindClosed (PolyG, P,Q,R,DegG,DegR,DegQ,DegP);
end; { FindS }
44
{================== =DRIVER= = = = = = = = = = = = = = =}
VAR
Upper,
Lower,
P_array : Parameter;
UpNum, {number of upper parameters}
LowNum : Integer; {number of lower parameters}
Argue : Real;
DegG,DegR,
DegQ,DegP : Integer;
P, Q, R : Poly Type;
Begin
CLRSCR;
Initialize (Upper, Lower,P_array);
GetHyper (Upper, Lower, UpNum, LowNum, Argue);
Pause;
Display (Upper, Lower, P_array, UpNum, LowNum, Argue, 1);
Pause;
StepOne (Upper, Lower, P_array, UpNum, LowNum);
Display (Upper, Lower, P_array, UpNum, LowNum, Argue,2);
Degree (Upper, Lower, P_array, UpNum, LowNum,
Argue,P,Q,R,DegG,DegR,DegQ,DegP);
Pause;
if DegG > = 0 then
begin
FindPoly (P,Q,R,DegG,DegR,DegQ,DegP);
end {if}
else
begin
writeln ( not summable with respect to hypergeometric terms );
writeln ( ie, sum t(k) is not a hypergeometric term);
end; {else}
end. {hyper}
45
{Rick Reeves)
Unit Gauss;
{ This unit performs gaussian elimination on n by n matrices. }
I***********************************************************************
INTERFACE
Uses
CRT;
Const
MAXMATRIX = 12;
MAXANSWER = 20;
Type
Matrix_Type = Array [1..MAXMATRIX, 1..MAXMATRIX] of real;
Answer_Type = Array [1..MAXANSWER] of real;
Procedure GAUSSIAN_ELIM(var B : Matrix_Type; var IERR,N : integer;
var Answer : Answer_Type);
Procedure Initialize_Matrix(VAR Mat : Matrix_Type);
IMPLEMENTATION
VAR
ZERO MATRIX : Matrix Type;
Procedure GAUSSIAN ELIM (Var B : Matrix_Type; Var IERR,N : integer;
Var Answer : Answer_Type);
{ This procedure determines whether a linear system (in the form of
a augmented matrix) has a unique solution.
B Augmented matrix from the driver.
46
IERR Error code passed from the driver (See driver variables above).
N Number of equations in system.
P The smallest integer in the column that needs to have its
elements eliminated.
H,
K,
and L These are used to manipulate the various elements of the matrix.
TEMP This variable is used in the switching of rows if the smallest
P is not in the same row as I.
MJI This is the multiple used for creating zeros in the Gaussian
eliminations.
TOTAL This variable is used to determine the correct solutions during
backward substitution.
I,
J Both used as row and column indices.
ANSWER This array will hold the solution set if one exits. }
var
P,H,K,L,I,J : integer;
TEMP, MJI, TOTAL : real;
begin
IERR := 0;
I := 1;
{ If there are no errors, repeat the following for 1 through (Nl) }
While (IERR = 0) and (I < = (Nl) ) do
begin
{ Find the smallest integer P (no greater than N) for the particular
column search. }
P := I;
While (abs(B[P,I]) = 0.0) and (P < = N ) do
P := P + 1;
47
{ If no such P can be found then there is no unique solution. }
if P = (N +1) then
IERR := 1
{ If the value, P, found was not in the top row, I, then switch rows:
Ep Swithed with Ei }
else
begin
if ( P < > I) then
begin
for H : = 1 to (N + 1) do
begin
Temp := B[I,H];
B[I,H] := B[P,H];
B[P,H] := TEMP;
end; {for}
end; {if}
{ Perform the Gaussian Elimination on the elements below the chosen I,J
element. }
H := I + 1;
{ Step 1 find the multiples that will eliminate the other terms in
in the particular column. }
for J := (1+1) to N do
begin
MJI := B[J,I] / B[I,I];
48
{ Step 2 Perform (Ej Mji*Ei) and switch with Ej }
for K : = H to (N + 1) do
B[J,K] := B[J,K] MJI B[I,K];
B[J,I] := 0.0;
end; {for}
end; {else}
I := I + 1;
end; {while}
if IERR = 0 then
begin
{ If the Nth element of the last row equals zero there is no unique element.}
if (abs(B[N,N] ) = 0.0) then
IERR := 1
{ If not then perform backward substitution.}
else
begin
{ Step 1 Solve for the single variable of the last row first. }
ANSWER[N] := B[N,(N + 1)] / B[N,N];
{ Step 2 Using the above answer (and all previous answers), solve
the next variable. }
for K := 1 to (N 1) do
begin
I := (N 1) K +1;
H := I + 1;
TOTAL := B[I,(N + 1)];
for L : = H to N do
TOTAL := TOTAL B[I,L] ANSWER[L];
ANSWER [I] := TOTAL / B[I,I];
end; {for}
end; {else}
end; {if}
end; {Gaussian}
49
{}
Procedure Initialize_Matrix (VAR Mat: matrix_type);
I
begin
Mat := ZEROMATRIX;
end; {initialize}
{==================
VAR 1
I
i, j : integer; j
Begin {Initialization section}
For i: = 1 to MAXMATRIX do
For j:= 1 to MAXMATRIX do
ZERO MATRIX [i,j] := 0.0;
end. {Gauss unit}
I
I
50
UNIT Poly l;
. 
{This unit is a polynomial machine.
See below for data structure}
:^***********************#**********************************************j.
' ; i
INTERFACE I
i
i
CONST .
MAXP = 70; j
!
TYPE 1
Poly_Type =f= ARRAY [0..MAXP] of Real;
! i
i
i
. Procedure Get (VAR poly_out: Poly_Type);
1 Procedure Show (poly_in : Poly_Type; Deg : integer);
Procedure Add (first_in, second_in:Poly_Type; VAR answer_out:Poly_Type);
Procedure Subtract (first_in, second_in:Poly_Type; VAR answer_out:Poly_Type);
I
j
Procedure Multiply (first_in, second_in:Poly_Type;VAR answer_out:Poly_Type);
Function Evaluate (poly_in : Poly_Type; x_value: REAL): REAL;
i
**********************************************************************
IMPLEMENTATION
j i
USES 1
Crt;
i
i VAR !
ZERO POLY : Poly Type;
I
51
}
Procedure Get (VAR poly_out: Poly_Type);
Begin j
Poly Out := ZERO POLY;
end; {show} :
{:}
Procedure Show (poly in : Poly Type; Deg : integer);
; I
VAR
x,y,i : integer;
Begin
writeln; 1
writeln; '
y := WhereY; ,
For i: = 0 to Deg do
Begin
if WhereX > 70 then
begin
writeln;
writeln;
writeln;
writeln;
end;

write (poly_in[degi]:5:2);
if deg i < > 0 then
write(k);
x : = Whei;eX;
GOTOXY (x+1, y1);
if deg; i < > 0 then
Begin
write (degi);
GOTOXY (x+5, y);
write ( + );
end; {if}
end; {for} j
writeln;
52
writeln;
writeln;
end; {show}
Procedure Multiply ( First_in, Second_in : Poly_Type;
 VAR Answer_out : Poly_Type);
VAR
i,j : integer; j
, i
Begin
Answer out := ZEROPOLY;
For i := 0 to MAXP do
For j : = 0 to MAXP do
Answer_out[i+j] := Answer_out[i+j] +
! ( First_in[i] Second_in[j]);
end;
{
Procedure Add ( First_in, Second_in
VAR Answer_out
}
: Poly_Type;
: Poly Type);
VAR
i : integer;
Begin j
Answer out: = ZERO POLY;
For i := 0 to MAXP do
Answer_out[i] := First_in[i] + Second_in[i];
end:
{}
Procedure Subtract ( First_in, Second_in : Poly_Type;
j VAR Answer_out : Poly_Type);
VAR j
i : integer; 
Begin 
Answer_out :=! ZERO_POLY;
For i: = 0 to MAXP do
Answer_out[i] := First_in[i] Second_in[i];
end; j
{i}
Function Evaluate1 (poly_in : Poly_Type; x_value: REAL): REAL;
i
Begin 
j
end;
{=========i========================:
VAR
i : BYTE;
Begin {Initialization section)
For i: = 0 to MAXP do
ZERO POLY i[i] := 0.0;
end. 1
I
54
References
1. ! J. Boyle, E. Lusk, R. Overeek, and L. Wos, Automated
Reasoning. McGrawHill, (1992).
i
2. ; C.IL Chang and R. Lee, Symbolic Logic and Mechanical
Theorem Proving. Academic Press, (1973).
r
3. I. Gessel and D. Stanton, "Strange evaluations of
j, hypergeometric series," SIAM Journal of Mathematical
; Analysis, 13 (1982), 295308.
4. R.W. Gosper, "Decision procedures for indefinite
hypergeometric summation," Proceedings of the National
Academy of Science U.S.A., 75 (1978), 4042.
!
5. , R.ll. Graham, D.E. Knuth, and O. Patashnik, Concrete
; Mathematics. AddisonWesley Publishing Co., (1989).
6. D.W. Loveland, Automated Theorem Proving: A logical
Basis. NorthHolland Publishing, (1978)/
7. R. Roy, "Binomial identities and hypergeometric series,"
American Mathematical Monthly, 97 (1987), 3646.
8. : H.S. Wilf, Generatingfunctionologv. Academic Press, Inc.,
: (1990).
9. H.S. Wilf and D. Zeilberger, "Towards computerized proofs
of identities", Bulletin of the American Mathematical Society,
23 (July 1990), 7783.
10. !' D. jZeilberger, "A fast algorithm for proving terminating
' hypergeometric identities," Discrete Mathematics, 80 (1990),
t 207211.
55

PAGE 1
GOSPER'S ALGORITHM A DECISION PROCEDURE FOR PARTIAL HYPERGEOMETRIC SUMS by Richard B. Reeves B.A., University of Colorado at Denver, 1990 A thesis submitted to the Faculty of the Graduate School of the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Applied Mathematics 1992
PAGE 2
I I I I I I I I I I I I I This thesis for the Master of Science degree by Richard B. Reeves has been approved for the Department of Mathematics by William Cherowitzo Stanley E. Payne Sylvia Lu Date
PAGE 3
Reeves, Richard B. (M.S., Applied Mathematics) Gosper's Algorithm A Decision Procedure for Hypergeometric Sums Thesis directed by Associate Professor William Cherowitzo Given an integer function fk, does there exist a function S such that m L fk = S(m) S(O) or, equivalently fk = S(k) S(k1) ? kl An algorithm created by R.W. Gosper in 1978 gives a complete answer to this question. In creating this algorithm Gosper has made a decisive step in the mechanization of hypergeometric identity proofs. The form and content of this abstract are approved. I recommend its publication. Signed ________ William Cherowitzo iii
PAGE 4
To my parents Jim and Millie for their unwavering love, support and encouragement and To Terry for her friendship, intellect, humor, and for her love that compels me to see the beauty in life.
PAGE 5
ACKNOWLEDGEMENTS I would like to thank Bill Cherowitzo, Stan Payne, and Zenas Hartvigson for sharing their expertise. Above all, I would like to thank Dr. Cherowitzo for his guidance, support, and patience. v
PAGE 6
TABLE OF CONTENTS Introduction . . . . . . . . . . . . . . . . . . . . 1 Chapter 1: Preliminaries . . . . . . . . . . . . . . 2 1.1 1.2 Hypergeometric Functions .................... Hypergeometric Terms 1.3 General Method for Finding Hypergeometric 2 6 Representation . . . . . . . . . . . . . . 8 Chapter 2: Gosper's Algorithm . . . . . . . . . . . 14 2.1 Partial Hypergeometric Sums ................... 14 2.2 Algorithmic Method . . . . . . . . . . . . 15 Chapter 3: Mechanized Proofs ........................ 24 Appendix Hypergeometric Program . . . . . . . . . . . . . 29 References . . . . . . . . . . . . . . . . . . . . 55 Vl
PAGE 7
Introduction Let fk be a complex values function of the integer k. Given such a function fk, does there exist a function S (within an additive constant) such that m I; fk = S(m) S(O) koJ or, equivalently fk = S(k) S(k1) ? Such a sum of the function f in finite calculus would be called an "indefinite sum", in analogy with the indefinite integral of infinite calculus. An algorithm created by R. W. Gosper in 1978 gives a complete answer to this question. In other words Gosper's algorithm will either return the functionS (with the desired properties), or it will guarantee that no such function exists. It never returns the statements "not sure" or "unknown" to this question. This paper examines Gosper's Algorithm. 1
PAGE 8
Chapter 1: Preliminaries 1.1 Hypergeometric Functions Gosper has discovered a way to decide whether a given function is partially (indefinitely) summable with respect to a general class of functions called hypergeometric terms. Therefore, in order to investigate Gosper's algorithm, it is first necessary to understand hypergeometric series or functions. Hypergeometric (HG) series were first investigated by Euler, Gauss, Pfaff, and Riemann. The study of HG series has evolved into a systematic means of handling a great variety of binomial coefficient identities [5). The following notation will be useful in the discussion of HG series. Definition 1.1: Rising Factorial Powers xm x(x+l) ... (x+m1) integer m ;,: 0 Definition 1.2: Falling Factorial Powers x111 x(xl) ... (xm+l) integer m ;,: 0 2
PAGE 9
The general hypergeometric series is a power series in z (z c), and is defined as follows in terms of rising factorial powers: Definition 1.3: k k al ... am zk k! where, to avoid division by zero, none of the b parameters may be zero or a negative integer. Using standard conventions, we define x0 = 1 for all x and 0 = 1. Therefore, the first term (k = 0) of any HG series is equal to 1. Many important functions occur as some special case of the general hypergeometric series. This is one of the reasons hypergeometric functions are so powerful. The following examples illustrate a small portion of this power. Example 1.1 (The simplest situationno parameters): m 0, n 0 F( lz) hO k. This "blank" notation does make sense, when m or n equals zero, because the notation indicates the absence of parameters. However, this form has proven to be perplexing to many. In order to avoid confusion, it has 3
PAGE 10
become standard to add an extra one above and below whenever m = 0 or whenever n = 0. So, F ( lz) It is now important to note that, when using HG series in Gosper's algorithm, we are not concerned with convergence of the series. We are simply using z as a formal symbol. Therefore, any identities derived will be formally true. An example of this is the HG function F (1'i'1iz) = k!zk. This series doesn't converge for any nonzero z but it is still very useful in solving problems [5]. Example 1.2 (specific example): m = 1, n = 0, a1 1 written as m = 2, a1 = a2 = 1, and n = 1, b1 1 F (I lz) = F (1,1 lz) = I:zk = _1 1 ho 1z This is the geometric series. Hypergeometric series are so named because the geometric series is included as a very specific case. 4
PAGE 11
Example 1.3 (general example): m = 1, n = 0, F (ail lz) = L a 1 1zk = L (a)(a+l);(a+kl)zk hO k. hO k. Example 1.4: F ( 1z) = L (a+kl)zk = 1 ho k (1z)" = L (al' (d = L (a)(a1) (ak+l)zk hO k! hO k! = L (a)z k = (1 +z)a hO k This utilizes the binomial theorem. Note, a negative integer as an upper parameter causes the infinite series to become finite since ( ai 0 whenever k > a ;o: 0 and a is an integer. Example 1.5 (The Gauss Hypergeometric so named because many of its properties were first proved by Gauss in his doctoral dissertation in 1812; although Euler and Pfaff had already discovered many of these properties [ 5] ) : 5
PAGE 12
Before everything was generalized to arbitrary m and n (around 1870), this was considered the hypergeometric series. One important special case of the Gauss Hypergeometric is: ln(l +z)=z F(1 11z) = z L k!k! <d 2 hO (k+ 1)! k! = z L (zi hO (k+l) 4 Note: z'1 In (1 + z) is a hypergeometric function but In (I + z) is not, since by convention every hypergeometric function has a value of 1, when z = 0. 1.2 Hypergeometric Terms It is essential to realize which series can be expressed as a hypergeometric series. This determination is made by examining the ratio of consecutive terms of any hypergeometric series. Let where t = k As previously stated t0 = 1, and all other terms have the following ratio: 6
PAGE 13
k+l k+l al ... am k k a1 ... am (a1 +k) ... (am +k) z = __; __ __;;,;._ (bl +k) ... (b. +k) k+ 1 Therefore, a hypergeometric series is a series .E tk in which the ratio of hO every two consecutive terms is a rational function of the summation variable k. This form exhibits the fact that the a's are the negatives of the roots of the numerator and the b's are negatives of the roots of the denominator of the term ratio. In order to maintain a canonical form for the term ratio, if (k + 1) is not a factor in the denominator, we add it to both the numerator and the denominator. To restate, hypergeometric series are those series whose constant term (first term) is 1 and whose term ratio is a rational function of the summand variable. Given a rational function, we may always find a hypergeometric series whose term ratio is this function. Example 1.6: adding the missing (k + 1) we get 7 = _,__( k_3')_,__( k_4..:..;) (o...!.l/.'8) (k+i/4)(ki/4)
PAGE 14
(k3)(k4)(1/8)(k+ 1) (k+if4)(kif4)(k+ 1) This is the term ratio of (3,4,1 1) F i/4,i/4fg 1.3 General Method for Finding Hypergeometric Representation of a Given Function S. Step 1) Step 2) Step 3) Step 4) Write S = L t,, t0 0. hO Calculate the term ratio. If term ratio is not a rational function of k, the series cannot be expressed in hypergeometric form. Otherwise, find the negatives of the roots of the numerator and denominator and express in the following form: (a1 +k) (am +k) z (b1 +k) (bn +k) k+ 1 This gives the parameters a1 .. ,a01; b1 .. ,bn; and the argument z such that 8
PAGE 15
Example 1.7 Step 1) Step 2) Step 4) (S = sin (z)): z3 sin(z) = z 3! __ ( l)k z<2kl)' sin(z) L.. hO (2k+ 1)! z1 + . 7! first term is z. (l)kz2k z I sin (z) = L first term is 1. hO (2k+ 1)! tk+i = (1r1 z2k 2 (2k+1)! = tk ( 1)k z2k (2k+ 1)! ( 1) z2 (2k+3)(2k+2) tk 4 (k+3/2)(k+ 1) ( 1 z2 ) sm(z) = z F 1 ,312 I 4 9
PAGE 16
One of the strengths of hypergeometric series is the ability to express and manipulate a great variety of binomial summations in a systematic way. Example 1.8: A well known binomial coefficient identity is the paral.lel summation law: integer n If we wished to express the parallel summation law in its hypergeometric form we would proceed as follows: Step 1) integer n If we replace k with nk we get I: nk!:.n L L (r+nk) = L (r+nk)l k
PAGE 17
Step 2) Step 4) tk+l = (r+nk1)! (nk)! r! = (nk) tk (r+nk)! (nk1)! r! (r+nk) (nk) (r+nk) (nk) (k+ 1) (r+nk) (k+ 1) ( r+n) F (_ n,l I 1 ) = (r+n+I) n n Solving for the hypergeometric function in this equation, we obtain F ( n,1 I 1) = rn r+n+1 r+1 Note that since the lower parameters of a hypergeometric cannot be zero or negative integers, this represents a degenerate case of the parallel summation law because parallel summation is usually applied when r and n are positive integers. This limitation is bypassed in this and other situations where it arises by considering: lim F ( 1 n I 1) "0 nr+e 11
PAGE 18
Example 1.9 (Vandermonde's Convolution Identity): L ( ') ( s ) = ( r+s) ho k nk n A combinatorial description of this identity is the sum of all possible ways to pick k men from a total of r men and the number of ways to pick nk women from a total of s women equals the number of ways to choose n people from a group consisting of r men and s women. Let t ( ') ( s ) r! s! k k nk k!(rk)!(nk)!(sn+k)! We can make several general observations about the process of converting to hypergeometric form [5]. 1) Whenever tk contains a factor like (a+k)!, with a plus sign before the k, we get (a+k+l)! = (a+k+l) in the term ratio. Which (a +k)! contributes parameter (a+!) to the corresponding HG function, as an upper parameter if (a +k)! was in the numerator of tk, or as a lower parameter if it was in the denominator of tk. 2) Whenever tk contains a factor like (a k)!, with a minus sign before the k we get (akl)! = (!) in the term ratio. Which contributes (a k)! (ka) 12
PAGE 19
Chapter 2: Gasper's Algorithm 2.1 Partial Hypergeometric Sums We have seen how hypergeometric sums can be used to systematically express a large number of binomial coefficient sums. It would appear that partial sums may be needed in many instances (combinatorial solutions, physical problems, and computational complexity determination to name a few). When dealing with such sums it is often desirable to find a closed form solution that works over a general range. Gasper's Algorithm decides whether a given summand (function) is partially summable (within an additive constant) with respect to a general class of functions called hypergeometric terms by examining m L fk = S(m) S(O) kl or, equivalently fk = S(k) S(k1). (1) Given fk the algorithm finds those S(k) with the property S(k) = a rational function of k S(k1) 2.2 Algorithmic Method (4] 14
PAGE 20
Step 1 If S(k) is a rational function of k then the term ratio S(k1) S(k) S(k1) S(k1) S(k2) S(k) 1 = S(k1) 1 S(k2) S(k1) is also a rational function of k. Therefore, we can express the term ratio of (1) as (2) where pk, qk, and rk are polynomials ink subject the following condition: I qk and (k+ f3) I rk (3) f3 is not a nonnegative integer. It is always possible to achieve this condition with a rational function with the following steps: a) Since the term ratio is a rational function of k, we can express this ratio in hypergeometric form, (a +k)(a +k) z 1 m (b1 +k) (b n +k) b) Set: 15
PAGE 21
Step 2 Pk1 ; 1, c) If (3) is violated then go to J d) If q and r have factors (k + a) and (k + {3) respectively that violate (3), then divide them out of q and r and replace Pk by pk(k+uf; pk(k+a)(k+a1}(k+P+1) where N; ap. The new p, q, and r still satisfy (2) since Pt is replaced by Pk1 Pk (k+a) h1 (k+P) e) Repeat step condition (3) is satisfied. Recall the goal is to find some hypergeometric term S(k) such that (1) is satisfied. Assume fk; S(k) S(k1) and now write (4) 16
PAGE 22
where g(k) is a function to be found. Combining (1) and (4) we get g(k) = Pk S(k) qk+t S(k) S(k1) I l S(k1) S(k) This implies g(k) is a rational function of k whenever S(k) I S(k1) is. By substituting ( 4) into (1) we get Multiplying through by Pklfk and using (2) (to solve for r k) we get pk = qk+t g(k) rk g(k1) (5) If we can find a g(k) that satisfies this recurrence, we've found S(k) and we are done. Theorem: If S(k) I S(k1) is a rational function of k, then g(k) is a polynomial. Proof: We know that g(k) is a rational function when S(k) I S(k1) is, so suppose by way of contradiction g(k) = u(k) where v(k) 1 v(k) 17
PAGE 23
and u(k) and v(k) have no factors in common. Let N be the largest integer such that (k+ fl) and (k+ fl + N) both occur as factors of v(k) for some fl
PAGE 24
Step 3 All that remains is finding, given pk, qk, and rk, if there exists a polynomial g(k) that satisfies equation (5). This turns out to be surprisingly easy, once we know the degree d of g(k). Let (6) Rewrite equation (5) in the following form: 2p(k) = Dif.!Q(k) x (g(k) + g(k1)) SumQ(k) x (g(k) g(k1)) where Dif.!Q(k) = q(k+ 1) r(k) and SumQ(k) = q(k+ 1) + r(k). (7) If g(k) has degree d then sum (g(k) + g(k1)) = 2 o:dkct + ... also has degree d. Whereas the difference (g(k) g(k1)) = do:dkdJ + ... has degree d1 (assume the zero polynomial has degree 1.) Write deg(p) for the degree of any polynomial p(k). Case 1: Case 2: The deg(DiffQ) 2: deg(SumQ). This implies the degree of the righthand side of (7) is deg(DiffQ) + d, so d = deg(p) deg(DiffQ). The deg(DiffQ) < deg(SumQ). This implies the we can write 19
PAGE 25
Dif!Q(k) ; pkd'l + ... and SumQ(k) ; y kd' + ... where y 0. The righthand side of (7) now has the form This leaves two possibilities: 1) 2P + yd 0 d; deg(p) deg(SumQ) + 1, or 2) 2P + yd ; 0 = d > deg(p) deg(SumQ) + 1. This second case needs to be examined only if 2fJ h is an integer d greater than deg(p) deg(SumQ) + 1. If d < 0, the algorithm stops with the message that no function S, in equation (1), exists. If a nonnegative d is discovered, we can plug equation (6) into equation (5). The polynomial g(x) will satisfy the recurrence if and only if the a's of equation (6) satisfy certain linear equations, because each power of k must have the same coefficient on both sides of (5). Example 2.1 (This example was computed usmg the Derive software package on an IBM compatible 486 PC): 20
PAGE 26
Use Gosper's algorithm to find a closed form identity for the following summation: Step 1 f(k) f(k1) Step 2 Let k1 m rr aj2 + bj + c I: i"1 k+1 k"l rr aj2 + bj + d i"1 k1 rr aj2 + bj + c f(k) i"' k+1 rr a/+ bj + d j=1 k1 k rr aj2 + bj + c rr aj2 +bj + d j=1 j=! k+! k2 rr aj2 + bj + d rr aj2 + bj + c j=1 j=! a(k1i + b(k1) + c a(k+ 1)2 + b(k+ I) + d q(k) a(k1)2 + b(k1) + c r(k) a(k+ 1)2 + b(k+ 1) + d 21
PAGE 27
and suppose q and r have no factors in common. This implies pk = l. Therefore, which requires Case 2, where we discover 2{3/y ; 2 > 1 ; deg(p) deg(SumQ) + 1 So g will be a quadratic function of k with three undetermined coefficients. Step 3 Substituting into equation (5) 1 ; k2 ( gz(2a +b +c d) + g1 ( a) ) + k ( gz(b+2d) + g1(a+cd) + g0(2a)) + ( g2(abd) + g1(a+b+d) + (ab+cd)) and solving for the coefficients of g we get the augmented matrix [2a+b+cd a b+2d a+cd abd a+b+d OJ 2a 0 ab+cd 1 0 22
PAGE 28
which reduces to 1 0 0 0 1 0 0 0 1 (a2+2a(c+d) (b+cd)(bc+d))(cd) 2a(2a+b+cd) (a 2+2a(c+d) (b+cd)(bc+d))(cd) 2a 2 +a(2b+ 3c d)+( cd)(b +c d) ((a 2 + 2a( c +d) (b+cd)(bc +d))( c d) Having found the above coefficients for g we can now find S(k). S(k) q(k+ 1) g(k) j(k) k1 (ak2+bk+c)(2a2(k2+2k+ 1)+(a(2b(k+ 1 +k+3)d(2k+ 1))+(cd)(b+cd)) (aj2+bj+c) bl (a 2 +2a(c+d) (b+c d)(bc+d))(c d) n (aj2 +bj +d) j=J Finally we discover the closed form of this partial sum, k1 In II ajz + bj + c :E j=l ; S(m) S(O) k+l k=l II a/ + bj + d j=l ; m1 (am 2+bm+c)(2a2(m2+2m+ !) +a(2b(m+ l)+c(2m+3) d(2m+ l))+(cd)(b+cd)) n (aj2+bj+c) j=l m (a2+2a(c+d) (b+cd)(bc +d))(c d) n (aj2+bj+d) j=l 2a 2 +a(2b+ 3c d)+( c d)(b+c d) (a 2 +2a( c +d) (b +c d)( be +d))( a +b +d)( c d) 23
PAGE 29
Chapter 3: Mechanized Proofs Combinatorial problems, such as the previous example, arise in many areas of mathematics, engineering, and theoretical physics. Their solution often entails the evaluation of a sum of products of binomial coefficients (7]. A great many distinguished mathematicians have spent an immense amount of time and ingenuity proving virtually thousands of binomial identities. We now know that much of that time and talent was wasted. Recently it has become widely known that "almost all" known binomial coefficient identities are special cases of relatively few hypergeometric identities (3]. Therefore, there are a comparatively small number of essentially different binomial identities. The problem in recognizing this has been twofold. First, the difficulty in recognizing equivalent binomial expression due to the notation used. The notation used to express binomial coefficients makes it very easy to disguise equivalent expressions in a such a way that no relationship is readily apparent. The second problem has been the reluctance of many to see the advantage of using hypergeometric series to serve as a canonical form for these identities (even though this was the method used by Euler and Gauss [7]). Although this reluctance has diminished over the last decade 24
PAGE 30
and the use of hyper geometric series has increased, the need for identities has not decreased. This need for hypergeometric identities is being largely supported by the use of computers. The increased use of hypergeometric series has directly lead to an increased use of computers in creating hypergeometric identities. Wilf asserts, "that more and more hypergeometric identities are still being conjectured and proved, so the need remains for mechanized proof' [3,78]. It would have been difficult (and in practice infeasible) to create the identity in example 2.1, without the aid of a computer. For this particular example, the software program Derive was used to find the appropriate identity. Gosper credits the program MACSYMA and the support of its developers at MIT for giving him ... the experiences necessary to provoke the conjectures that led to this algorithm" [4, 42]. Doron Zeilberger has derived other identities (using algorithms largely based on Gosper's algorithm) using the MAPLE software package. Finding identities and mathematical proofs using computers is becoming more the norm [10]. Gosper's algorithm is example of mechanized theorem proving. Mechanized theorem proving has been a goal for a long time. The desire dates back to Leibniz in the late 17th and early 18th centuries[2]. 25
PAGE 31
In 1930 Herbrand proposed a system for mechanized (or automated) theorem proving based on logic rules[6]. His method was slow and untimely, since it came before the advent of the digital computer. With the arrival of the computer came a renewed interest in automated theorem proving. In the 1960's Herbrand's method was implemented on a computer by Gilmore[2]. Until 1965, "success was largely confined to problems found in logic books rather than mathematics books." [ 6, 5] In 1965, J.A. Robinson developed a major breakthrough. She produced the resolution inference system. This system was based on a single inference rule and did not require generalized logic axioms (only those specific to the problem). In essence the inference rule is: an algorithm that, when successfully applied to some given set of hypotheses or premisses, yields a conclusion that follows inevitably and logically from the premisses. [1,197] In other words, the algorithm is a finite automata that stops. Gosper's algorithm "infers" that a closed form can be found if a certain polynomial can be found. If this polynomial doesn't exist, Gosper's algorithm infers that no closed form can be found. Since 1965 many refinements of the resolution system have been made. Wilf and Zeilberger have shown that a very large class of identities can be proved by computers [9]. They do not imply that computers need 26
PAGE 32
to be "trusted blindly"; in fact, it is just the contrary. Although the proofs are "discovered" by a computer, all proofs presented by Wilf and Zeilberger produce a proof certificate that can be checked by hand if desired [9]. Traditionalists might complain that identity proofs generated by a computer can give no insight to the mathematician. Zeilberger cautions against being too chauvinistic, "The team humancomputer is a mighty one, and an openminded human can draw inspiration from all sources, even from a machine" [10]. I found it very enlightening to write my own program (see addendum) partially implementing Gasper's algorithm. This program will take as input the term ratio (already factored) of a specific hypergeometric series and output a partial sum or the fact that one doesn't exist. The purpose in writing this program was not to improve or advance implementations of Gasper's algorithm already in use. The purpose was strictly one of education. Perhaps most edifying were the attempts (admittedly small) at symbolic manipulation. I gained not only an insight into mechanized proofs, but an appreciation of the software techniques use in packages such as Derive. In 1978 Gosper made a definitive step in the direction of mechanized proofs of hypergeometric series. His algorithm is now fully 27
PAGE 33
implemented in MACSYMA and partially implemented in Mathematica. It will no doubt become more widely implemented when its usefulness becomes more fully recognized (8]. 28
PAGE 34
Appendix A {Rick Reeves} {1992} Program Hyper (Input, Output); {This program is an attempt to implement R.W. Gosper's algorithm perhaps not in its complete form but perhaps in a way that gives the user an appreciation of what the algorithm entails. } {The algorithm is not shown here due to its length. Please see Gosper's original paper.} USES Gauss, Poly_l, Crt; { a gaussian elimination unit for n by n matrices } { a polynomial unit add, subtract, multiply, show, etc.} { standard crt unit } TYPE Parameter = array [l..MAXP] of real; { used for factored form of poly} {} Procedure Pause; { Holds user screen until key is press then clears screen.} VAR Ch: Char; Begin GOTOXY(27,25); Write ('Press any key to continue'); Ch : = ReadKey; CLRSCR; end; {pause} 29
PAGE 35
{ } Procedure Initialize (V AR UpperP, LowerP, P _array : Parameter); { This procedure initializes all the parameter arrays. } VAR i : integer; Begin fori:= 1 to MAXP do Begin UpperP[i] : = 999; LowerP[i] : = 999; P array[i]: = 999; end; {for} P_array[l] := 9911; {flagfor step one} end; {initialize} { } Procedure GetHyper (V AR Upper Para, Lower Para : Parameter; V AR Up Number, Low Number : Integer; V AR Arguement : Real); { This procedure gets the factored hypergeometric term ratio from the user. Note: this maybe the same as the hypergeometric parameters except for perhaps the k + 1 factor in the denominator.} VAR integer; Begin GOTOXY(10,3); Write (OUTPUT, 'Please enter NUMBER of UPPER parameters:'); Readln (INPUT, UpNumber); 30
PAGE 36
For i : = 1 to Up Number do Begin GOTOXY(15,5 +I); Write (OUTPUT, 'Please enter UPPER parameter number ',i,': '); Readln (INPUT, UpperPara[i]); end; {for} GOTOXY(l0,14); Write (OUTPUT, 'Please enter NUMBER of LOWER parameters: '); Readln (INPUT, LowNumber); Fori := 1 to LowNumber do Begin GOTOXY(15,15+i); Write (OUTPUT, 'Please enter LOWER parameter number ',i,': '); Readln (INPUT, LowerPara[i]); end; {for} GOTOXY(l0,24); Write (OUTPUT, 'Please enter value of arguement z: '); Readln (INPUT, Arguement); end; {gethyper} {} Procedure StepOne (V AR Q array, R array, P array : Parameter; Up, Low : integer); {Get p(k), q(k), and r(k) in right form, ie, put term ratio in the proper form, } VAR i, j, p, k, X GOOD check Chk integer; : boolean; real; integer; 31
PAGE 37
Begin p := 1; {check to see if (alpha beta) is a positive number} Repeat GOOD : = TRUE; for i : = 1 to (Up) do for j : = 1 to (Low) do begin if (Q array[i] < > 999.0) and (R arrayUJ < > 999.0) then begincheck : = Q_ array[i] R array[j]; if ( check > 0 ) and (int( check) = check) then begin Chk : = trunc( check); GOOD:= FALSE; X:= 1; for k : = p to (p + Chk 2 ) do begin P array[k] : = Q_ array[i] (x); X : = X + 1; end; {for} p:=p+chk1; Q_ array[i] : = 999; R array[j] : = 999; end; {if} end; {if} end; {for} Until GOOD; for i : = 1 to up do Q_ array[i] : = Q_ array[i] + 1; end; {stepone} 32
PAGE 38
{} Procedure Display ( Q ay, R ay, P ay : Parameter; Up, Low : Integer; Argue, Which : Real); { This procedure displays the polynomials in their factored form. } VAR i,j : integer; Begin GOTOXY (10,10); if which = 1.0 then Write ('q(k) = ') else Write ('q(k+ 1) = '); Fori : = 1 to Up do if Q ay[ i] < > 998 then if Q ay[i] = 0 then wnte (' (k) ') else if Q ay[i] < 0 then write (' (k',abs(Q_ ay[i]):l:l,') ') else write(' (k + ',Q ay[i]:l:l,') '); Writeln (' ( ',Argue:l:i;) '); GOTOXY (10,12); Write ('r(k) = ;); Fori:= 1 to (Low) do if R_ ay[i] < > 999 then if R ay[i] = 0 then write (' (k) ') else if R ay[i] < 0 then write(' (k',abs(R_ay[i]):l:l,') ') else write (' (k + ',R ay[i]:l:l,') '); GOTOXY (10,14); i := 1; 33
PAGE 39
If P ay[l] < > 9911 then begin Write ('p(k) = '); While P ay[i] < > 999 do begin if P ay[i] = 0 then write (' (k) ') else if P ay[i] < 0 then write (' (k',abs(P ay[i]):l:l,') ') else write (' (k + ',P ay[i]:l:l,') '); i : = i + 1; if (WhereX >70) then Writeln; end; {while} end {if} else writeln ('p(k) = l'); writeln; end; {display} {} Procedure Degree ( Q ay, R ay, P ay : Parameter; Up, Low : Integer; Argue : Real; V AR P,Q,R : Poly_ Type; V AR DegreeG,DegR, DegQ, DegP : Integer); VAR I, DegDiff, DegSum, IntPivot : Integer; PTerm, Qterm, Rterm, Temp, QSum, QDiff Poly_ Type; First Boolean; Beta, Gamma, Pivot : Real; { This procedure finds the degree of the unknown function G if it exists. If it doesn't exist, a message is outputted that indicates that. This procedure also finds the degree of all pertinent polynomials and puts them in standard form. This procedure requires units: Polyl and Crt} 34
PAGE 40
Begin DegR := 0; DegQ := 0; DegP := 0; DegSum : = 0; { initialize all variables to zero.} DegDiff : = 0; DegreeG : = 0; First : = TRUE; Get (Temp); Temp[O] : = Argue; { find degrees of p, q, and r} for i : = 1 to (Low+ 0) do if (R ay[i] < > 999) then DegR := DegR + 1; if DegR = 0 then R_ay[1] := 1; for i : = 1 to (Up) do if (Q_ay[i] < > 998) then DegQ : = DegQ + 1; if DegQ = 0 then Q_ay[1] := 1; fori:= 1 to MAXP do {degree of r} {degree of q} if (P ay[i] < > 999) and (P _ay[i] < > 9911) then {degree of p} DegP : = DegP + 1; Get (R); Get (Q); Get (P); GeT (Qterm); Get (Rterm); Get (Pterm); { get "zeroed" coefficient arrays } { find standard polynomial form for p, q, and r.} 35
PAGE 41
For i : = 1 to Up do if (Q ay[i] < > 998) then if First then Begin First:= FALSE; Qterm[O] : = Q_ ay[i]; Qterm[1] :"' 1; Multiply (Qterm,temp, Q); temp:= Q; end {if} else Begin Qterm[O] : = Q_ ay[i]; Qterm[l] : = 1; Multiply (Qterm,temp, Q); temp:= Q; end; {else} Write ('q(k+ 1) = '); Show(Q,DegQ); First : = TRUE; Get (Temp); Temp[ OJ : = 1; For i : = 1 to Low + 0 do if (R ay[i] < > 999) then if First then Begin First : = FALSE; Rterm[O] : = R _ay[i]; Rterm[1] : = 1; Multiply (Rterm,temp, R); temp:= R; end {if} else Begin Rterm[O] : = R ay[i]; Rterm[1] : = 1;Multiply (Rterm,temp, R); temp:= R; end; {else} { find q} {display q} { find r} 36
PAGE 42
Write ('r(k) = '); {display r} Show(R,DegR); First : = TRUE; Get (Temp); Temp[O] : = 1; Fori:= 1 to MAXP do if (P ay[i] < > 999) and (P ay[i] < > 9911 ) then if First then Begin First : = FALSE; Pterm[OJ : = P ay[i]; Pterm[l] : = 1; Multiply (Pterm,temp, P); { find p} temp:= P; end {if} else Begin Pterm[O] : = P ay[i]; Pterm[1] : = 1; Multiply (Pterm,temp, P); temp:= P; end; {else} Write ('p(k) = '); if DegP > 0 then Show(P,DegP) {display p} else begin Writeln (' 1'); P[O] : = 1; end; Pause; { find term sums and differences } Add (R,Q,QSum); Fori:= 0 to MAXP do if QSum[i] < > 0 then DegSum := i; Write ('Qsum(k) = '); Show (QSum, DegSum); Subtract (Q,R,QDiff); 37
PAGE 43
Fori:= 0 to MAXP do if QDiff[i] < > 0 then DegDiff : = i; Write ('QDiff(k) = '); Show (QDiff, DegDiff); {Will the degree of our mystery polynomial stand up and sign in please.} Pivot : = 0; Beta : = QDiff[DegDiff]; Gamma : = QSum[DegSum]; if Gamma < > 0 then Pivot : = 2 beta / gamma; IntPivot : = trunc(Pivot); if DegDiff > = DegSum then DegreeG : = DegP DegDiff else if ( Pivot > ( degP DegSum + 1) ) and ( Pivot = IntPivot ) then DegreeG : = IntPivot else DegreeG := (degPDegSum + 1); Writeln ('Degree of G = ',DegreeG); DegreeG : = 2; end; {degree} {} Function Factorial (n : integer): integer; { This is a recursive function that finds factorials.} begin if n = 0 then Factorial : = 1 else Factorial : = n Factorial(n1) end; {factorial} 38
PAGE 44
{ } Function Combo (n,k :integer): integer; { This function finds combinations it uses function factorial.} Var temp : integer; begin if n < = k then temp:= 1 else temp : = Factorial(n) Div (Factorial(k) Factorial(nk)); Combo : = temp; end; {combo} {} Procedure BuildSK(Var Mat : Matrix_ Type; Degree : Integer); { This procedure finds the term (in matrix form) g(k1) } Var i, j integer; Begin Fori : = 1 to Degree+ 1 do For j :=ito Degree+1 do if odd (ji) then Mat[i,j] := Combo(j1,ji)*(1) else Mat[i,j] : = Combo(j1,ji) ; Fori : = 1 to Degree+ 1 do begin For j := 1 to Degree+ 1 do Write (Mat [i,j]:8:1); writeln; end; Pause; end; {buildSK} 39
PAGE 45
{ } Procedure BuildS (Var Mat : Matrix_ Type; Degree : Integer); { This procedure finds the term (in matrix form) g(k) } Var 1, J integer; Begin Fori : = 1 to Degree+ 1 do Mat [i,i] : = 1; Fori:= 1 to Degree+1 do begin For j : = 1 to Degree+ 1 do Write (Mat [i,j]:8: 1 ); writeln; end; Pause; end; {buildS} { } Procedure BuildTerm (Mat : Matrix_ Type; Poly : Poly_ Type; DegM,DegP :Integer; Var Term: Matrix _type); { This procedure is used to find the terms (in matrix form) q(k+ 1) g(k) and r(k) g(k1). } Var i, j, k integer; Begin Fori:= 0 to DegP do writeln (poly[i]:3:0); Fori:= 0 to DegP do For j : = 1 to DegM + 1 do For k : = 1 to DegM + 1 do Term[i + j,k] : = Term[i + j,k] + (Poly[i] MatU,k]); 40
PAGE 46
For i : = 1 to DegM + 1 + DegP do Begin For j := 1 to DegM+l+DegP do Write (Term [i,j]:8:1); write In; end; Pause; end; {Term} {} Procedure Subtract (Mat1, Mat2 : Matrix_type; Var Answer: Matrix_Type); { This procedure is used to subtract terms that are in matrix form } Var i, j integer; begin for i : = 1 to MAXMA1RIX do for j : = 1 to MAXMATRIX do Answer[i,j] : = Mat1[i,j] Mat2[i,j]; end; { 0 UTPUT } Procedure OUTPUT( ANSWER : Answer_ Type; B : Matrix_ Type; IERR,N integer); { This procedure outputs to a file (Defaulted to Screen for now) the Gaussian reduced matrix and the solution set or a message that no unique solution exits.} var Outfile I, J begin : text; : Integer; {allows programmer to pick output device} {loop and matrix counters} assign(Outfile, 'CON'); rewrite(Outfile); 41
PAGE 47
{ If there are not errors output the reduced matrix and the solution set.} if IERR = 0 then begin writeln(Outfile); writeln(Outfile ); writeln(Outfile,'The Gaussian reduced matrix is as follows:'); writeln(Outfile ); for I : = 1 to N do begin for J : = 1 to (N + 1) do write(OUTFILE, B[I,J] :12:2); writeln(OUTFILE); end; {for} writeln(Outfile ); writeln(Outfile,'The solution set is as follows:'); writeln(Outfile ); for I : = 1 to N do writeln(Outfile,'G',I,' = ', ANSWER[I]:12:2); writeln(Outfile ); end {if} { Else, output that there are no unique solutions to this linear system.} else begin writeln(Outfile ); writeln(Outfile); writeln (Outfile,'There is no unique solution to this linear system.'); writeln(Outfile ); for I : = 1 to N do begin for J : = 1 to (N + 1) do write(OUTFILE, B[I,J] :8:1); writeln(OUTFILE); end; {for} writeln(Outfile ); end; {else} close(OUTFILE); Pause; end; {output} 42
PAGE 48
{} Procedure FindClosed (PolyG, P,Q,R : Poly Type; DegG,DegR,DegQ,DegP : Integer); { This procedure is used to find the closed of the summation.} Var Temp begin Get(Temp); Poly_Type; Multiply (Q,PolyG,Temp ); Multiply (Temp,R,Temp); Show(temp, DegG+DegQ+DegR); Pause; end; {findclosed} {} Procedure FindPoly (P,Q,R : Poly_ Type; DegG,DegR,DegQ,DegP : Integer); { This procedure finds the coefficients of G (if they exist).} Var SK, SK_one Answer, First Term, SecondTerm i,j, N, Error Coef Max PolyG Begin error : = 0; Matrix Type; Matrix_ Type; Integer; Answer_ Type; Integer; Poly_Type; Initialize Matrix (SK); Initialize Matrix (SK one); Initialize Matrix (FirstTerm); Initialize Matrix (SecondTerm); Initialize =Matrix (Answer); 43
PAGE 49
Get (PolyG); BuildS (SK,DegG); BuildSK (SK_one, DegG); BuildTerm(SK, Q, DegG, DegQ, FirstTerm); BuildTerm(SK _one R, DegG, DegR, SecondTerM); Subtract (FirstTerm, SecondTerm, Answer); if DegQ > DegR then N := DegQ else N := DegR; N : = N + DegG + 1; if DegG+2 > N Then Max:= DegG+2 else Max:= N; For i : = 1 to MAXMATRIX do Answer[i, DegG + 2] : = P[i1]; Fori:= 1toMAXdo begin for j : = 1 to Max + 1 do write (Answer[i,j]:S:l); writeln; end; {for} { find # of rows} { find if more rows than cols} {show matrix} Pause; { find coefficients if possible} N:=DegG+1; Gaussian_ Elim (Answer,ERROR,N Coef); Output ( Coef, Answer,ERROR, N); Fori:= 1 to DegG+1 do PolyG[i1] : = Coef[i]; write ('g(k) : = '); show(PolyG, DegG); Pause; FindClosed (PolyG, P,Q,R,DegG,DegR,DegQ,DegP); end; { FindS } 44
PAGE 50
{===================DRIVER===============} VAR Upper, Lower, P _array : Parameter; UpNum, LowNum : Integer; Argue : Real; DegG,DegR, DegQ,DegP : Integer; P, Q, R : Poly_ Type; Begin CLRSCR; {number of upper parameters} {number of lower parameters} Initialize (Upper, Lower,P array); GetHyper (Upper, Lower,UpNum, LowNum, Argue); Pause; Display (Upper, Lower, P_array, UpNum, LowNum, Argue,l); Pause; StepOne (Upper, Lower, P array, UpNum, LowNum); Display (Upper, Lower, P array, UpNum, LowNum, Argue,2); Degree (Upper, Lower, P array, UpNum, LowNum, Argue,P,Q,R,DegG,DegR,DegQ,DegP); Pause; if DegG > = 0 then begin FindPoly (P,Q,R,DegG,DegR,DegQ,DegP); end {if} else begin writeln (' not summable with respect to hypergeometric terms '); writeln (' ie, sum t(k) is not a hypergeometric term'); end; {else} end. {hyper} 45
PAGE 51
{Rick Reeves} Unit Gauss; { This unit performs gaussian elimination on n by n matrices. } {***********************************************************************} INTERFACE Uses CRT; Canst MAXMATRIX = 12; MAXANSWER = 20; Type Matrix_ Type = Array [l..MAXMATRIX, l..MAXMATRIX] of real; Answer_Type =Array [l..MAXANSWER] of real; Procedure GAUSSIAN_ ELIM(var B : Matrix_ Type; var IERR,N : integer; var Answer : Answer_ Type); Procedure Initialize_ Matrix(V AR Mat : Matrix Type); {**********************************************************************} IMPLEMENTATION VAR ZERO_ MATRIX : Matrix_ Type; Procedure GAUSSIAN ELIM (Var B : Matrix Type; Var IERR,N: integer; Var Answer: Answer Type); { This procedure determines whether a linear system ( in the form of a augmented matrix) has a unique solution. B Augmented matrix from the driver. 46
PAGE 52
IERR N Error code passed from the driver (See driver variables above). Number of equations in system. p H, K, The smallest integer in the column that needs to have its elements eliminated. and L These are used to manipulate the various elements of the matrix. TEMP This variable is used in the switching of rows if the smallest P is not in the same row as I. MJI This is the multiple used for creating zeros in the Gaussian eliminations. TOTAL This variable is used to determine the correct solutions during backward substitution. I, J Both used as row and column indices. ANSWER This array will hold the solution set if one exits. } var P,H,K,L,I,J : integer; TEMP, MJI, TOTAL :real; begin IERR := 0; I : = 1; {If there are no errors, repeat the following for 1 through (N1)} While (!ERR = 0) and (I < = (N1) ) do begin { Find the smallest integer P (no greater than N) for the particular column search. } p :=I; While (abs(B[P,IJ) = 0.0) and (P < = N ) do p: = p + 1; 47
PAGE 53
{ If no such P can be found then there is no unique solution. } if P = (N + 1) then IERR := 1 { If the value, P, found was not in the top row, I, then switch rows: Ep Swithed with Ei } else begin if ( P < > I) then begin for H : = 1 to (N + 1) do begin Temp : = B[I,H]; B[I,H] : = B[P,H]; B[P,H] : = TEMP; end; {for} end; {if} { Perform the Gaussian Elimination on the elements below the chosen I,J element. } H:=I+1; { Step 1 find the multiples that will eliminate the other terms in in the particular column. } for J : = (I+ 1) to N do begin MJI : = B(J,I] / B(I,I]; 48
PAGE 54
{ Step 2 Perform (Ej Mji*Ei) and switch with Ej } for K : = H to (N + 1) do B[J,K] : = B[J,K] MJI B[I,K]; B[J,I] : = 0.0; end; {for} end; {else} I:= I + 1; end; {while} if IERR = 0 then begin { If the Nth element of the last row equals zero there is no unique element.} if ( abs(B[N,N] ) = 0.0) then IERR := 1 { If not then perform backward substitution.} else begin { Step 1 Solve for the single variable of the last row first. } ANSWER[N] := B[N,(N + 1)] I B[N,N]; { Step 2 Using the above answer (and all previous answers), solve the next variable. } for K : = 1 to (N 1) do begin I:=(N1)K+1; H:=I+1; TOTAL:= B[I,(N + 1)]; for L : = H to N do TOTAL:= TOTALB[I,L] ANSWER[L]; ANSWER[!] : = TOTAL I B[I,I]; end; {for} end; {else} end; {if} end; {Gaussian} 49
PAGE 55
{ } Procedure Initialize Matrix (V AR Mat: matrix_ type); begin Mat:= ZERO_MATRIX; end; {initialize} {==========================================} VAR i, j integer; Begin {Initialization section} Fori:= 1 to MAXMATRIX do For j: = 1 to MAXMATRIX do ZERO_ MATRIX [i,j] : = 0.0; end. {Gauss unit} 50
PAGE 56
UNIT Poly _1; {This unit is a polynomial machine. See below for data structure} {**********************************************************************} INTERFACE CONST MAXP = 70; TYPE Poly_Type = ARRAY [O .. MAXP] of Real; Procedure Get (VAR poly out : Poly Type); Procedure Show (poly _in :Poly_ Type; Deg : integer); Procedure Add (first in, second in:Poly Type; V AR answer out:Poly Type); Procedure Subtract (first_in, second _in: Poly_ Type; V AR answer _out:Poly _Type); Procedure Multiply (first in, second in:Poly Type;V AR answer out:Poly Type); Function Evaluate (poly _in : Poly _1ype; Xvalue: REAL): REAL; {**********************************************************************} IMPLEMENTATION USES Crt; VAR ZERO POLY : Poly Type; 51
PAGE 57
{} Procedure Get (V AR poly out : Poly Type); Begin Poly_ Out:= ZERO_POLY; end; {show} { } Procedure Show (poly in : Poly Type; Deg : integer); VAR x,y,i : integer; Begin write In; write In; y := WhereY; For i : = 0 to Deg do Begin if WhereX > 70 then begin write In; write In; writeln; writeln; end; write (poly _in[degi]:5:2); if deg i < > 0 then write('k'); x := WhereX; GOTOXY (x+l, y1); if deg i < > 0 then Begin write ( degi); GOTOXY (x+5, y); write (' + '); end; {if} end; {for} writeln; 52
PAGE 58
writeln; writeln; end; {show} { } Procedure Multiply ( First in, Second in : Poly Type; VAR Answer_out :Poly_ Type); VAR i,j integer; Begin Answer_out := ZERO_POLY; Fori:= 0 to MAXP do For j := 0 to MAXP do Answer out[i + j] : = Answer out[i + j] + ( First_in[i] Second _in[j]); end; { } Procedure Add ( First_in, Second _in : Poly_ Type; VAR Answer_out :Poly_ Type); VAR integer; Begin Answer out:= ZERO POLY; Fori:; o to MAXP do Answer_ out[i] : = First_in[i] + Second in[i]; end; 53
PAGE 59
{} Procedure Subtract ( First in, Second in : Poly Type; VAR Answer_out : Poly_Type); VAR 1 integer; Begin Answer out:= ZERO POLY; For i : ,;o to MAXP do Answer out[i] : = First in[i] Second in[i]; end; {} Function Evaluate (poly_in: Poly_Type; x_value: REAL): REAL; Begin end; {============================= VAR i: BYTE; Begin {Initialization section} For i: = 0 to MAXP do ZERO _POLY [i] : = 0.0; end. 54 } 
PAGE 60
References 1. J. Boyle, E. Lusk, R. Overeek, and L. Was, Automated Reasoning, McGrawHill, (1992). 2. C.L. Chang and R. Lee, Symbolic Logic and Mechanical Theorem Proving, Academic Press, (1973). 3. I. Gessel and D. Stanton, "Strange evaluations of hypergeometric series," SIAM Journal of Mathematical Analysis, 13 (1982), 295308. 4. R.W. Gasper, "Decision procedures for indefinite hypergeometric summation," Proceedings of the National Academy of Science U.S.A., 75 (1978), 4042. 5. R.L. Graham, D.E. Knuth, and 0. Patashnik, Concrete Mathematics, AddisonWesley Publishing Co., (1989). 6. D.W. Loveland, Automated Theorem Proving: A logical Basis, NorthHolland Publishing, (1978)/ 7. R. Roy, "Binomial identities and hypergeometric series," American Mathematical Monthly, 97 (1987), 3646. 8. H.S. Wilf, Generatingfunctionology, Academic Press, Inc., (1990). 9. H.S. Wilf and D. Zeilberger, "Towards computerized proofs of identities", Bulletin of the American Mathematical Society, 23 (July 1990), 7783. 10. D. Zeilberger, "A fast algorithm for proving terminating hypergeometric identities," Discrete Mathematics, 80 (1990), 207211. 55
